On the relationship between δO₂∕N₂ variability and ice sheet surface conditions in Antarctica

While the processes controlling pore closure are broadly understood, the physical mechanisms driving the associated elemental fractionation remains ambiguous. Previous studies have shown that the pore closure process leads to a depletion in small-sized molecules (e.g. H₂, O₂, Ar, Ne, He) in ice core bubbles relative to larger-sized molecules like N₂. This size-dependent fractionation, identified using ice core $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ records, exhibits a clear anti-correlation with local summer solstice insolation, making $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ a valuable ice core dating tool. Mechanisms controlling this relationship are attributed to the physical properties of deep firn. In this study, we compile $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ records from 15 polar ice cores and show a new additional link between $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ and local surface temperature and/or accumulation rate. Using the Crocus snowpack model, we perform sensitivity tests to identify the response of near-surface snow properties to changes in insolation intensity, accumulation rate, and air temperature. These tests support a mechanism linked to firn grain size, such that the larger the grain size for a given density, the stronger the pore closure fractionation and, hence, the lower the $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ values archived in the ice. Based on both snowpack model outputs and data compilation, our findings suggest that local accumulation rate and temperature should be considered when interpreting $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ as a local insolation proxy.

1 Introduction

Ice cores store crucial information for our understanding of past climate variability and atmospheric composition. Interpreting ice core gas records first requires an understanding of the evolution of snow into ice via the firnification processes. Firn is the name given to the layer of consolidated snow which makes up the top 50–120 m of ice sheets. Atmospheric air moves through porous networks within the firn until a critical depth (known as the lock-in depth) where vertical diffusion effectively stops, and pores gradually become sealed off from the atmosphere. The lock-in depth and the depth at which all pores are closed (close-off depth) are largely determined by local accumulation rate, temperature, and possibly the degree of density layering (Schwander et al., 1997; Martinerie et al., 1994; Mitchell et al., 2015). The region between the lock-in depth and close-off depth, where gases are progressively occluded into bubbles, is known as the lock-in zone (LIZ).

Gas records from ice cores provide a vital dating tool, especially at low-accumulation sites where other methods are unsuitable. One such tool is orbital dating, which enables certain gas records to be tuned to insolation curves directly calculated from astronomical variables (Laskar et al., 2004). For example, ice core δ¹⁸O_atm records (δ¹⁸O of atmospheric O₂) are strongly correlated with precession (mid-June 65° N insolation) (Bender et al., 1994). δ¹⁸O_atm provides a direct atmospheric signal which may ultimately be used to align different ice core records (Extier et al., 2018). Two alternative proxies for orbital dating are (1) total air content, which is anti-correlated with integrated summer insolation, i.e. the annual sum of daily insolation above a certain threshold (Raynaud et al., 2007; Eicher et al., 2016), and (2) $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$, which is anti-correlated with local summer solstice insolation intensity (e.g. Bender, 2002; Kawamura et al., 2007; Suwa and Bender, 2008b; Landais et al., 2012; Bouchet et al., 2023). Unlike δ¹⁸O_atm, total air content (TAC) and $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ reflect processes within the firn column, making the records site specific. The term $\mathit{\delta }({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}})$ (hereafter, simply $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$) describes the relative difference between the ratio of O₂ to N₂ molecules trapped within the ice and that of the standard atmosphere and is expressed in the delta notation commonly used for stable isotope ratios.

The use of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ for dating was first proposed by Bender (2002) after observations of an anti-correlation between local summer solstice insolation (hereafter SSI) and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ from the Vostok ice core. A similar relationship was then observed at other sites such as Dome Fuji (Kawamura et al., 2007) and EPICA Dome C (Landais et al., 2012) in Antarctica, as well as GISP2 (Suwa and Bender, 2008b) in Greenland. Over orbital timescales, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is in anti-phase with local SSI when drawn on the ice-age chronology, indicating that the firn properties controlling the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ fractionation are set near the surface. It is understood that increased temperature gradients under high SSI conditions enhance near-surface snow metamorphism, thus increasing near-surface grain size which persists throughout the firnification process down to the close-off depth (Bender, 2002; Severinghaus and Battle, 2006; Suwa and Bender, 2008a; Fujita et al., 2009).

Parallel firn air studies of the open porosity revealed an enrichment in O₂ and other small molecules (Ar, Ne, and He) in the LIZ, providing further evidence of size-dependent fractionation during pore closure (Battle et al., 1996; Huber et al., 2006; Severinghaus and Battle, 2006). While the physical mechanisms controlling the amount of fractionation are not fully understood, it is believed that smaller molecules (<3.6 Å diameter) escape during pore closure via molecular diffusion through the ice lattice (or permeation), driven by pressure gradients between recently closed pores and neighbouring open pores (Ikeda-Fukazawa et al., 2004; Huber et al., 2006; Severinghaus and Battle, 2006). This process is facilitated by the pore network's capacity to export the fugitive gases back to the atmosphere, which is required for the observed depletion of O₂ in bubbles (Fujita et al., 2009).

In addition, Fujita et al. (2009) proposed that enhanced layering in the deep firn under high SSI conditions would also facilitate $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ depletion. They argue that large-grained, relatively low-density “summer” layers in the deep firn close off deeper and take longer to do so than neighbouring “winter” layers (smaller-grained, relatively high density). As a result, summer layers remain permeable for longer, allowing the O₂-enriched air in the open porosity to be exported to the atmosphere, thereby reducing bulk ice $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ under high SSI conditions (Fujita et al., 2009). While the proposed mechanisms are posited to explain the SSI imprint on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, they are also influenced by local climate conditions such as temperature and accumulation rate. Indeed, there is a substantial amount of evidence linking local climate conditions with both firn physical properties (Casado et al., 2021; Inoue et al., 2024) and deep firn layering (Hörhold et al., 2011).

There is a growing body of evidence for a local climatic imprint on ice core $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records. Firstly, spectral analysis has revealed climate-related 100 kyr cyclicity at EPICA Dome C corresponding to glacial cycles (Bazin et al., 2016). However, this 100 kyr cyclicity is not apparent in the Dome Fuji $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record, which Kawamura et al. (2007) attribute to the idea that temperature and accumulation effects cancel each other out. Secondly, millennial-scale variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records from the GISP2 core appeared in-phase with local temperature fluctuations driven by Dansgaard–Oeschger events (Suwa and Bender, 2008b). In contrast, Kobashi et al. (2015) evidenced an anti-phase effect of accumulation rate on $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ records at GISP2 on decadal timescales over the last 6000 years. Like O₂, Ar is a smaller molecule than N₂, invoking the same driving mechanisms for $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variations, albeit the $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ anomaly tends to be half as large as the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ anomaly (Bender et al., 1995; Buizert et al., 2023). Kobashi et al. (2015) proposed a direct effect of accumulation rate or temperature on the $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ($\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$) variations through the firn depth. Higher accumulation rates (or lower temperatures) increase firn weight, leading to greater overloading pressures in microbubbles, which results in the preferential expulsion of Ar (O₂) in the LIZ. Alternatively, Severinghaus and Battle (2006) proposed that higher accumulation rates result in a more rapid burial of bubbles, allowing less time for gases to escape during pore closure.

Figure 1 Locations of each ice core site initially included in our study (Matsuoka et al., 2018; Moon et al., 2023).

Using a combination of data compilation and snowpack modelling, we aim to develop our understanding of the formation of the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records by first determining the role of local climate parameters (accumulation rate and temperature), on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability and subsequently identifying potential mechanisms related to snow physical properties using snowpack sensitivity tests. We use a compilation of datasets from 15 ice cores from both Antarctica and Greenland to identify spatial and temporal patterns in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ depending on local surface conditions. The impacts of SSI and local climate on snow properties are then assessed using the SURFEX-ISBA-Crocus detailed snowpack model (Vionnet et al., 2012). We aim to constrain the influence of near-surface snow properties on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability, potentially contributing to a mechanistic explanation for elemental fractionation during pore closure.

2 Methods

2.1 Ice core sites

We compiled $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records from 18 ice cores from Antarctica and Greenland but use data from 15 of those sites for reasons explained in Sect. 2.2.3. Previously published data were measured on ice cores from the following: Dome Fuji (DF), EPICA Dome C (EDC), Greenland Ice Core Project 2 (GISP2), Law Dome (LD) DE08 and DSSW20k, North Greenland Ice core Project (NGRIP), Roosevelt Island Climate Evolution (RICE), Siple Dome (SD), South Pole (SP), Talos Dome (TALDICE), Vostok, and the West Antarctic Ice Sheet Divide (WAIS) ice cores (references for all datasets are presented in Table S1 in the Supplement). We also present unpublished data from Berkner Island (BI), EPICA Dronning Maud Land (EDML), Fletcher Promontory (FP), GISP2, James Ross Island (JRI), North Greenland Eemian Ice Drilling (NEEM), Skytrain Ice Rise (SIR), and Talos Dome (TALDICE) ice cores (an overview of unpublished datasets is in Table S2). Table 1 provides an overview of the site characteristics.

Table 1 Overview of ice core site characteristics. The acronym BCTZ stands for brittle–clathrate transition zone, and COD stands for close-off depth.

¹ Mulvaney et al. (2007)${}^{\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d}}$; ² Oyabu et al. (2021)^a,d, Fujita et al. (1998)^b,c; ³ Parrenin et al. (2012), ^a EPICA community members (2004)^b,c, Landais et al. (2006)^d; ⁴ Ueltzhöffer et al. (2010)^a, EPICA community members (2006)${}^{\mathrm{b},\mathrm{c},\mathrm{d}}$; ⁵ Mulvaney et al. (2014)${}^{\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d}}$; ⁶ Gow et al. (1997)^a, Alley et al. (1993)^b, Alley and Koci (1988)^c, Cuffey and Paterson (2010)^d; ⁷ Mulvaney et al. (2014)^a,d, Capron et al. (2013)^b,c; ⁸ Etheridge et al. (1996)${}^{\mathrm{a},\mathrm{b},\mathrm{c},\mathrm{d}}$, Etheridge and Wookey (1988)^b,c, Etheridge et al. (1996)^d; ⁹ Rubino et al. (2019)^a, Morgan et al. (1997)^b,c, Buizert and Severinghaus (2016)^d; ¹⁰ Warming et al. (2013)^a, NEEM Community Members (2013)^b,c, Buizert et al. (2012)^d; ¹¹ NGRIP project members (2004)${}^{\mathrm{a},\mathrm{b},\mathrm{c}}$, Martinerie et al. (2009)^d; ¹² Lee et al. (2020)^a,, Bertler et al. (2017)${}^{\mathrm{b},\mathrm{c},\mathrm{d}}$; ¹³ Gow and Meese (2007)^a, Severinghaus et al. (2001)${}^{\mathrm{b},\mathrm{c},\mathrm{d}}$; ¹⁴ Mulvaney et al. (2021)${}^{\mathrm{a},\mathrm{b},\mathrm{c}}$, Hoffmann et al. (2022)^d; ¹⁵ Lazzara et al. (2012)^a, Severinghaus et al. (2001)${}^{\mathrm{b},\mathrm{c},\mathrm{d}}$; ¹⁶ Neff (2014)^a, Stenni et al. (2002)^b,c; ¹⁷ Uchida et al. (1994)^a, Petit et al. (1999)^b,c, Cuffey and Paterson (2010)^d; ¹⁸ Fitzpatrick et al. (2014)^a, Fegyveresi et al. (2011)^b,c, Battle et al. (2011)^d.
^*Reported bubble–clathrate transition depth corresponds to the maximum depth of RICE core (Lee et al., 2020).

2.2 Analytical techniques for previously unpublished data

The previously unpublished $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ datasets were measured at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE), with the addition of new GISP2 data measured at Scripps Institution of Oceanography (Scripps). At both LSCE and Scripps, gases are extracted from the ice using a melt–refreeze technique based on the method described by Sowers et al. (1989), with modifications to the LSCE method as described by Landais et al. (2003) and to the Scripps method as outlined in Petrenko et al. (2006). In short, ice samples are placed into glass flasks at −20 °C and the atmospheric air is evacuated. Samples are then left to slowly melt to release the trapped gases, before being refrozen. The extracted air samples are then passed through a CO₂ and water vapour trap, before being trapped in the stainless-steel dip-tubes submerged in liquid helium. An alternative semi-automated extraction technique was also used at LSCE which removes the need for refreezing the samples (Bazin et al., 2016; Bouchet et al., 2023).

2.2.1 New $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements

Samples from the following sites were measured at LSCE, France, on a 10-collector Thermo Delta V Plus isotope ratio mass spectrometer (IRMS), unless otherwise stated. This instrument allows for simultaneous measurements of elemental and isotopic ratios, including $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, δ¹⁵N-N₂, and δ¹⁸O-O₂. We note that δ¹⁵N-N₂ data from the BI, EDML, JRI, and TALDICE samples were previously published by Capron et al. (2013). In most cases, the values used in this study are the average of at least two replicate measurements and have an average analytical uncertainty of 0.5 ‰ for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$.

• Berkner Island. Measurements were performed on bubbly ice from the BI ice core every 55 cm (every bag) between 631 and 680 m, corresponding to 10.27 to 21.35 ka (Capron et al., 2013). Replicate samples were prepared at LSCE using the melt–refreeze method between March 2010 and March 2011.

• EDML. Nine samples were measured on bubbly ice from the EDML ice core over five depth levels between 328 and 473 m (327.8 m (4.51 ka), 354.2 m (4.95 ka), 381 m (5.43 ka), 467 m (7.04 ka), and 473 m (7.16 ka)) (Bazin et al., 2013). The samples were prepared using the melt–refreeze method at LSCE. Corresponding measurements of δ¹⁵N of N₂ were previously published in Capron et al. (2013).

• Fletcher Promontory. In January 2015, 39 depth levels were measured from the FP ice core, retrieved in 2012. Measurements were performed approximately every 3 m starting at 289 m down to 388 m. There is currently no published age scale for the FP ice core. All samples were prepared using the melt extraction method.

• James Ross Island. Between February and March 2011, measurements were performed at 16 depth levels on the JRI ice core. The depth resolution varied between 2–50 m starting at 52 m until 363 m corresponding to 0.03–14.3 ka (Mulvaney et al., 2012). Samples were prepared using the melt extraction method.

• NEEM. Clathrate ice from the NEEM ice core was measured between February and April 2011, a year after the core was retrieved. A total of 119 depth levels were sampled at varying resolutions over the following intervals: 55 cm intervals (every bag) between 1757–1773 m (38.127–39.735 ka b2k), 5.5 m intervals (every 10 bags) between 2205 and 2370 m (108.56–120.237 ka b2k), every 2 bags from 2375–2434 m (no published age scale available below these depths), and 5.5 m intervals (every 10 bags) between 2436 and 2519 m (Gkinis et al., 2021; Rasmussen et al., 2013). Samples were prepared using the melt–refreeze method.

• Skytrain Ice Rise. Measurements were performed on bubbly ice from the SIR ice core between March and April 2021. Samples were taken sporadically (1–15 m intervals) at 16 depth levels between 307 and 436 m depth (4.707–11.696 ka) (Mulvaney et al., 2023). Each sample was prepared at LSCE using the melt extraction method.

• TALDICE. Numerous measurements have been performed on bubbly and clathrate ice from TALDICE between 2008 and 2022 at LSCE. A total of 308 depth levels were measured at varying intervals starting at 155 m down to 1617 m. Published age scales reach 1548 m, giving an age range of 1.55–343 ka for TALDICE samples (Buiron et al., 2011; Crotti et al., 2021). All samples were prepared using the melt extraction technique. Some measurements between 1356–1620 m depth have been published previously and are available in Crotti et al. (2021).

Samples from GISP2 were measured at Scripps on a 3 kV Thermo Finnigan Delta V Plus dual-inlet IRMS. All measurements were measured in replicate.

• GISP2. Two datasets were measured on GISP2 ice. The first set was measured in 2009 at 266 depth intervals from 318.67–1201.22 m. The second set of measurements was performed on samples from 643 depths between 1740 and 2400 m (13–50 ka) over several measurement campaigns between 2017–2020 The δ¹⁵N of N₂ data from these samples were previously reported (Martin et al., 2023). In both cases, gases were extracting using the melt–refreeze technique, and all measurements were referenced to La Jolla pier air.

2.2.2 Corrections

Chemical slope and pressure imbalance corrections are applied to the measurements during data processing (Landais et al., 2003). In addition, all data are corrected for gravitational fractionation in the firn using δ¹⁵N of N₂ from the same samples.

$$\begin{array}{}\text{(1)}& \mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\text{grav. corr.}=\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}-\mathrm{4}\cdot {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}\end{array}$$

Gas loss effects during ice core storage are well documented to modify $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, causing significant depletion in O₂ in clathrate ice stored above −50 °C (Ikeda-Fukazawa et al., 2005; Kawamura et al., 2007; Landais et al., 2012). Ikeda-Fukazawa et al. (2005) proposed a correction for storage gas loss effects as a function of temperature and time. However, given the incomplete storage history for all ice cores, we do not attempt to correct for storage gas loss but rather define rejection criteria outlined in Sect. 2.2.3.

2.2.3 Data rejection criteria

Ice core storage histories need to be considered before interpreting the data to account for post-coring gas loss effects which disturb the signal (Sect. 2.2.2). Successive $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements from TALDICE and GISP2 clathrate ice samples stored above −50 °C show strong depletion of O₂ through time (Sect. S1 in the Supplement), which is consistent with observations from EDC (Bouchet et al., 2023). We systematically reject measurements from clathrate ice stored at −20 °C for over 3 years or at −36 °C for more than 4 years. Bubbly ice stored at these same temperatures appears to be mostly unaffected by gas loss (Sect. S1.2), with the exception of Vostok (Bender, 2002).

Several measurements were performed on ice within the bubble–clathrate transition zone (BCTZ), where high hydrostatic pressure in bubbles cause entrapped gases to form clathrate hydrates (Schaefer et al., 2011). At these depths, elemental fractionation occurs due to some gas species being preferentially incorporated into the clathrate structures (Ikeda-Fukazawa et al., 2001), thus making the interpretation of gas measurements unreliable (Bender, 2002). Measurements from the BCTZ may either have increased mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (usually in excess of 0 ‰) or strong data scattering, expressed as a high standard deviation (Oyabu et al., 2021). To avoid adding biases to our analysis, measurements from BCTZ are removed (Sect. S1.1). Additional scattering in elemental ratios, characterised by a standard deviation of 6.2 ‰ compared to 1.8 ‰ in pure bubble and clathrate ice, is observed below the BCTZ in the WAIS record between 1300–1500 m (Shackleton, 2019). Similar effects have been documented in the EDC and TALDICE ice cores (Lüthi et al., 2010), as well as the Dome Fuji ice core (Oyabu et al., 2021). Data influenced by this scattering effect were also removed from our analysis, followed by the removal of outliers (2σ) from the cleaned datasets.

We note that these criteria result in the removal of all data from NGRIP, as well as sections of data from other sites (Sect. S1 for overview). We also exclude FP and SIR due to limited availability of accumulation rate and temperature records over the depth/age ranges covered by the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ date. The remaining 15 datasets are presented in Table 3 and were used to analyse the drivers of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability.

2.3 Modelling near-surface snow properties

The second component of our study addresses the modelled response of snow physical properties to perturbations in SSI, accumulation rate, and temperature with the aim of identifying which properties may be influencing elemental fractionation during pore closure. We use the SURFEX-ISBA-Crocus detailed snowpack model (Crocus hereafter) to simulate snowpack evolution (Vionnet et al., 2012). Crocus simulates changes in snow physical properties induced by surface metamorphism and the evolution of these properties with depth. The model is forced by ERA5 reanalysis data (Hersbach et al., 2020), and the snowpack is initialised with measurements of snow density, effective optical radius of snow grains, and snow temperature. Optical radius is defined as the radius which snow grains would have for their surface area-to-volume ratio if they were spherical (Domine et al., 2006). Optical radius is thus directly linked to specific surface area (SSA), defined as the surface area of snow at the ice–air interface per unit mass (units m² kg⁻¹) (Legagneux et al., 2002), via the following equation:

$$\begin{array}{}\text{(2)}& \text{SSA}={\displaystyle \frac{\mathrm{3}}{r_{\mathrm{opt}}\cdot {\mathit{\rho }}_{\mathrm{ice}}}},\end{array}$$

where r_opt is the optical radius and where ρ_ice is the density of ice (Gallet et al., 2014). We use this model to assess changes in snow physical properties near the surface which are invoked to explain $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability. Dome C is used as the test site, given the abundance of snowpack observations as well as high-resolution $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data.

2.3.1 Crocus model description

Crocus is a one-dimensional model which simulates the evolution of snow properties with time and depth on a layer-by-layer basis, i.e. in a Lagrangian framework (Vionnet et al., 2012). A detailed description of the model can be found in Vionnet et al. (2012). Briefly, the initial number of layers is defined by the user, with the thickness of each layer allowed to change along the simulation (layer thickness ranging from 2 mm at the surface to metres thick). The maximum number of layers available in the model was increased from 50 (Libois et al., 2014) to 80 (this study) to maximise the resolution with depth, owing to the higher number of thin layers forming at Dome C than at Alpine sites. Once the simulated snowpack consists of 80 layers, the aggregation scheme merges internal neighbouring layers with similar properties, allowing a new surface layer to form. The key physical processes incorporated into Crocus for dry snow conditions are accumulation of snowfall, snow metamorphism, compaction of snow by the wind, compaction due to the weight of the overlaying layers, absorption of solar radiation, heat diffusion, and surface energy budget.

For our study, two fundamental user-defined model components are the snow metamorphism and radiative transfer schemes. We use the semi-empirical model from Flanner and Zender (2006) (F06) for the metamorphism scheme which describes the evolution of optical radius with time. F06 was found to be the most appropriate formulation for Dome C conditions (Carmagnola et al., 2014; Libois et al., 2014). To successfully reproduce the snow temperature profile (vital for realistically simulating snow metamorphism), the Two-streAm Radiative TransfEr in Snow (TARTES) model is used to account for vertical distribution of absorbed solar radiation in the snowpack (Libois et al., 2013). TARTES also considers the effect of impurities on snow temperature via albedo. For Dome C, we include black carbon content, which is set to 3 ng g⁻¹ (Warren et al., 2006; Libois et al., 2015). Here we assess the simulated snow density, snow temperature, and snow SSA from the Crocus model outputs.

2.3.2 Dome C-specific Crocus configuration

Crocus was initially developed for Alpine or sub-polar regions with seasonal snowpacks. Libois et al. (2014) modified multiple components of the Crocus model to improve its suitability to high-latitude sites with low accumulation rates – specifically for Dome C. The modifications are extensively described in Libois et al. (2014) and were implemented into the current version of Crocus for this study. The changes are as follows:

1. Fresh snow properties. The parameterisation of fresh snow density is based on temperature and wind speed, which results in an unrealistically low density for Dome C (50 kg m⁻³). Fresh snow density is fixed to a minimum of 170 kg m⁻³, the lowest fifth percentile from Dome C observations (Libois et al., 2014). Similarly, fresh snow SSA is set to 100 m² kg⁻¹ instead of the 65 m² kg⁻¹ used in the standard version of the model (Grenfell et al., 1994; Libois et al., 2014).

2. Wind-induced compaction. At low-accumulation sites, snow remains at the surface for prolonged periods of time. The long exposure time to surface winds facilitates compaction and, hence, increases density. The maximum surface snow density is increased from 350 kg m⁻³ to 450 kg m⁻³ to account for this effect (Albert et al., 2004; Libois et al., 2014).

3. Aggregation scheme. The formation of a new snow layer requires a minimum amount of snowfall. Due to the low accumulation rate at Dome C, the amount of snowfall needed to form a new layer was decreased from 0.03 mm h⁻¹ to 0.003 mm h⁻¹. In the instance when the snowpack has the maximum number of layers (80) at the time a new snow layer is formed, layers with similar properties will be aggregated, resulting in a smoothed signal. The aggregation scheme was disabled for the top six layers to resolve realistic near-surface snow temperature profiles, required to accurately simulate snow metamorphism.

Table 2 Overview of modifications made to forcing test parameters in Crocus snowpack sensitivity test scenarios. Perturbed forcing parameters are shown in bold.

2.3.3 Model initialisation

The snowpack was initialised with density and optical radius profiles measured in January 2010 at Dome C down to 20 m (Champollion et al., 2019), as well as snow temperature data from a probe installed at Dome C in 2012 with 5 cm resolution near the surface, coarsening with depth down to 12 m. ERA5 reanalysis data for Dome C was used to force the model at 3-hourly resolution over the period between 1 January 2000 and 1 December 2020 (Hersbach et al., 2020). The model requires atmospheric forcings for air temperature, accumulation rate, wind speed and direction, incoming shortwave and longwave radiation, and specific humidity. ERA5 gives a mean annual snowfall rate between 2000 and 2020 of 2.3 $\mathrm{cm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$, and as such, the snowfall rate was multiplied by 1.2 to match the observed mean annual accumulation rate of around 2.8 $\mathrm{cm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (Frezzotti et al., 2004; Libois et al., 2014). To ensure that at least the top 1 m consists of accumulated snow, a 100-year spin-up was used by running the forcing file 10 times between 2000 and 2010, followed by the period from 2000 to 2020. The outputs from 2010 to 2020 were then used for analysis.

2.3.4 Sensitivity tests

The sensitivity of snowpack properties to perturbations in surface forcings is tested by modifying one of three forcing parameters: incoming shortwave radiation, accumulation rate, or 2 m air temperature. The magnitude of the perturbation to each parameter corresponds to minimum and maximum values reconstructed over the last 800 kyr. We use shortwave radiation as a proxy of insolation and scale the values in proportion to the SSI values. A total of seven simulations are used to perform sensitivity analysis and are outlined in Table 2. The model configuration and initial snow profile were kept constant for each simulation; only the tested parameter in the atmospheric forcing file was modified as follows:

• Summer solstice insolation (SSI). Over the last 1000 years, the average SSI at 75.1° S was 544 W m⁻² (Laskar et al., 2004), compared to a minimum of 462 W m⁻² and a maximum of 601 W m⁻² over the past 800 kyr. To translate to forcing perturbations, the incoming shortwave radiation (SWR) is scaled by 85 % and 111 %, respectively, to reach the target values (462 and 601 W m⁻²). No additional modifications are applied to the annual distribution of SWR.

• Annual mean accumulation rate (A). Present-day accumulation rate at Dome C is set to 2.8 $\mathrm{cm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ following Libois et al. (2014). Hereafter, accumulation rates are expressed as water equivalent centimetres per year (in the form cm yr⁻¹). ERA5 snowfall was scaled by 36 % to reach the target accumulation rate of 1.0 cm yr⁻¹, representing the 800 kyr minimum, and 146 % to produce an accumulation rate of 4.1 cm yr⁻¹, which corresponds to the 800 kyr maximum (Bazin et al., 2013).

• Annual mean air temperature (T). Snowpack sensitivity to air temperature is tested by decreasing the 2 m air temperature by 10 °C for glacial conditions (Jouzel et al., 2007) and applying a 4 °C increase to represent the 800 kyr maximum temperatures. We note, however, that borehole temperature measurements and delta-age are more consistent with a 5 °C cooling (Buizert et al., 2021). Furthermore, these temperature modifications do not include changes in seasonal temperature variability but suffices for the purpose of identifying bulk changes in the snow properties. The average seasonal cycle is kept constant with an average amplitude of 35 °C.

It is important to highlight that, at polar sites, accumulation rate is dependent on temperature, and temperature is influenced by insolation, such that these parameters are not independent. However, we use the model to constrain the influences of each forcing parameter in an independent manner to understand the mechanisms, even if, in reality, these parameters are inter-dependent.

Figure 2 Scatterplots showing the relationship between SSI and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ at three Antarctic sites. Significant negative correlations (over 99 % confidence) are observed using relatively high-resolution data from (a) EDC (Bouchet et al., 2023), (b) Dome F (Kawamura et al., 2007; Oyabu et al., 2021), and (c) South Pole (Severinghaus, 2019). The slope (m, units ‰ W⁻¹ m²) and r values are presented for each site.

3 Results

3.1 Influence of SSI and local climate on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability in ice cores

Figure 2 shows $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ versus SSI for EDC, Dome F, and South Pole. These three sites are used owing to their temporal range and high-resolution $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements, largely uninfluenced by gas loss. The regression slopes vary between $-\mathrm{0.09}\pm \mathrm{0.006}$ ‰ W⁻¹ m² for South Pole to $-\mathrm{0.06}\pm \mathrm{0.005}$ ‰ W⁻¹ m² for EDC. The regression for Dome F falls within two standard deviations (2σ) of the regression for EDC, but the regression for South Pole falls outside this uncertainty. $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values from South Pole are higher than for EDC and Dome F for the same SSI, suggesting that additional factors are influencing the records, such as accumulation rate, which at South Pole is around three times that of both EDC and Dome F. In the following sections, we provide evidence for the influences of accumulation rate and air temperature on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, in addition to SSI, using both spatial (inter-site) and temporal (EDC ice core) variability of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$.

3.1.1 Inter-site comparison of mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$

In addition to the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ datasets, we compile SSI, accumulation rate, and temperature reconstructions for each site over the age/depth ranges which overlap in all records (Table 3). Drivers of mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ are then explored using regression analysis, as presented in Fig. 3. Due to limited data availability from the BI ice core, we use the average of present-day and Last Glacial Maximum temperature and accumulation rates from Capron et al. (2013). This is justified given the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data span approximately the Last Glacial Maximum to the start of the Holocene (Massam, 2018). Present-day values are used for accumulation rate and temperature at DE08 and DSSW20k due to the relatively young ages corresponding to the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data (∼0.09 ka).

Figure 3 Scatterplots showing the dependence of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ on (a) SSI, (b) accumulation rate (A), and (c) air temperature (T_air). Each point represents the mean values of each site over the depth interval in Table 3, with error bars showing the standard deviation (1σ). Sites from Greenland are identified by square markers. Linear regressions are shown in black, along with the associated slope (m, units ‰ W⁻¹ m²), correlation coefficient (r), and p value (p).

Table 3 Overview of data used in Fig. 3. Mean values are presented for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, accumulation rate (in water equivalent), temperature, and SSI at the current site latitude (Laskar et al., 2004). Where age scales are available, the SSI represents the mean values at the present-day site latitude. For sites without age scales, the average SSI over the last 1000 years is used. All data compiled for this study can be accessed at https://doi.org/10.5281/zenodo.12571888 (Harris Stuart, 2024).

¹ This study^a, Capron et al. (2013)^b,c
2 Kawamura et al. (2007)^a, Oyabu et al. (2021)^a, Kawamura et al. (2017)^b, Uemura et al. (2018)^c
3 Bouchet et al. (2023)^a, Bazin et al. (2013)^b, Jouzel et al. (2007)^c
4 This study^a, Bazin et al. (2013)^b, Stenni et al. (2010)^c
5 This study^a, Bender et al. (1994)^a, Smith (1998)^a, Suwa and Bender (2008b)^a, Cuffey and Clow (1999)^b, Clow (1999)^c
6 This study^a, Mulvaney et al. (2012)^b,c, Capron et al. (2013)^b,c
7 Buizert et al. (2020)^a, Etheridge and Wookey (1988)^b,c
8 Buizert et al. (2020)^a, Morgan et al. (1997)^b,c
9 This study^a, Rasmussen et al. (2013)^b
10 Lee et al. (2020)^a, Winstrup et al. (2019)^b
11 Severinghaus (2009)^a, Buizert (2021)^b,c
12 Severinghaus (2019)^a, Kahle et al. (2020)^b,c
13 This study^a, Crotti et al. (2021)^a, Bazin et al. (2013)^b,c
14 Bender (2002)^a, Suwa and Bender (2008a)^a, Bouchet et al. (2023)^b, Petit et al. (1999)^c
15 Severinghaus (2015)^a, Fudge et al. (2017)^b, White et al. (2019)^c.

Figure 3 indicates significant correlations between mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and all surface variables (SSI, accumulation rate, and temperature). Mean SSI appears to have the weakest influence on mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (Fig. 3a; $r=-\mathrm{0.46}$, p<0.1), but a slope of $-\mathrm{0.09}\pm \mathrm{0.05}$ ‰ W⁻¹ m² is statistically indistinguishable from the temporal slope (−0.07 ‰ W⁻¹ m²; Fig. 2). The weaker spatial correlation may be linked to the narrow range in SSI between sites compared to the temporal range of SSI at an individual site. Figure 3b and c show strong, significant correlations between mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and both the natural log of accumulation rate (r=0.86, p<0.0001) and temperature (r=0.79, p<0.005). The log-linear model in Fig. 3b indicates that a doubling of accumulation rate would result in a 1.4 ‰ increase in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ($\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}=\mathrm{2.6}\cdot \log \left(A\right)-\mathrm{12.6}$). However, it is important to note that temperature and the logarithm of accumulation rate are strongly correlated in Antarctica, such that the correlations seen in panels (b) and (c) of Fig. 3 are dependent on one another.

Deviations from the regression lines in Fig. 3 may have multiple origins, such as discrepancies in site latitude resulting from the ice flow speeds at different sites or simply uncertainties in the records. For BI and SD, we consider that the relatively low $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values may be the result of storage gas loss effects (Sects. 2.2.3 and S1), given that the measurements were carried out ∼7 and ∼8 years after coring for BI and SD, respectively. While clathrates are reportedly absent from BI (Mulvaney et al., 2007; Neff, 2014) and while SD data were removed below 600 m, i.e. the depth at which clathrate-like features first appear (Gow and Meese, 2007), the low values may link to gas loss during storage within increasingly brittle ice at depth in the core. Indeed, brittle ice is reported from 450 m at BI and from 400 m at SD (Neff, 2014).

3.1.2 Temporal variability of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ at Dome C

High-resolution data are required to investigate the temporal variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ as a function of accumulation rate and temperature. The δ-deuterium (δD) record from water isotope measurements is used as a qualitative proxy for accumulation rate and temperature, whereby higher δD values are generally associated with increased accumulation rate and temperature in ice cores from the East Antarctic plateau (Jouzel et al., 2007; Parrenin et al., 2007). The following analysis uses the longest period of relatively high-resolution $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements from the EDC ice core between 190–259 ka (1980–2350 m) on the AICC2023 ice-age chronology (Bouchet et al., 2023).

Figure 4 Evolution of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, SSI, and δD on the AICC2023 ice-age timescale from EDC (Bouchet et al., 2023). Panel (a) presents the 100-year interpolated $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ overlain with the low-pass-filtered curve (blue). The right y axis in panel (a) shows SSI at 75.1° S (grey). The predicted $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ as a function of SSI ($\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$–SSI predicted) is shown on the right axis of panel (b), and the residuals of the observed $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and predicted data are on the left axis. Panel (c) shows the 100-year interpolated δD overlain with the low-pass-filtered curve (green). The bottom two panels present the residuals of the low-pass-filtered curves for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (d) and δD (e). Grey bars highlight the pronounced millennial-scale variability.

Two approaches are used to extract the non-SSI signals in the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records. First, we interpolate SSI onto the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ages and take the deviations from the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$–SSI linear regression to isolate the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability not explained by SSI (hereafter, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$–SSI residual). The second approach directly investigates the millennial-scale variability by applying a low-pass filter to the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and δD records (interpolated onto a 100-year time step), using a 10 kyr cut-off to isolate the low-frequency signals associated with SSI. The filtered curves are then subtracted from the original curves to remove the orbital (SSI) signal, and a five-point moving average is applied to the residuals to reduce noise but retain millennial-scale variability in the records.

Figure 4 primarily shows the dominant SSI cyclicity in the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record, as has been documented previously (e.g. Landais et al., 2012). Superimposed onto this signal are millennial-scale peaks in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ which appear to coincide with peaks in δD, highlighted by grey bars. This high-frequency variability is more clearly identified from the filtered curves in Fig. 4d and e. The coherence between millennial-scale peaks in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and δD suggests a positive correlation with accumulation rate and temperature, which shares analogy with the positive spatial correlation between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and both accumulation rate and temperature (Fig. 3). We note a remnant 20 kyr variability in the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$–SSI residual, in phase with δD, which is present when using either the AICC2012 chronology (Bazin et al., 2013) or the AICC2023 chronology (Bouchet et al., 2023). Given that the new AICC2023 age scale (used here) has been orbitally tuned using these data (Bouchet et al., 2023), this signal is not expected to be the result of a phase difference between the EDC chronology and SSI.

3.2 Crocus model results

3.2.1 Crocus model evaluation for Dome C

The Crocus model outputs are first evaluated by comparing the reference simulation (“Reference” in Table 2) to observational data from Dome C. Simulated density and SSA profiles are compared to data from Libois et al. (2014), measured daily between 23 November 2012 and 16 January 2013 at two sites within 600 m of Concordia station, Dome C. Density was measured at 2.5 cm resolution down to 25 cm, while SSA was measured at 1 cm depth intervals down to 50 cm. Snow temperature at Dome C has been continuously measured since 2012 at 30 min intervals. Simulated snow density and SSA were interpolated to a fixed grid of 1 cm depth resolution and a 24 h time step between 1 January 2010 and 1 December 2020, while snow temperature was interpolated onto a 1 mm by 6-hourly grid over the same time period.

Density and SSA outputs presented in Fig. 5a and b are averaged values over the measurement period. For the most part, the observations fall within the range of the simulations. The simulated SSA profile is consistently within 1 standard deviation of the measurements below 10 cm; above this depth, simulated SSA is overestimated by up to 10 m² kg⁻¹. Small standard deviations associated with simulated density and SSA profiles suggest that variability is not well reproduced. As discussed in Libois et al. (2014), the standard version of Crocus is unable to reproduce density and SSA variability with depth due to its one-dimensional nature. For the purpose of our study, we consider the standard version sufficient to assess the overall sensitivity of snowpack properties to perturbations in forcings.

Figure 5 Comparison of observations and Crocus-simulated snowpack profiles. Density (a) and SSA (b) profiles represent the average from 23 November 2012 and 16 January 2013, covering the measurement period (Libois et al., 2014). Snow temperature profiles (c) both from Crocus (blue) and observations (grey) represent 23 November 2019 to 16 January 2020 averages. Faded grey lines represent individual profiles from observations, and the shaded bands show the standard deviations. Snow temperature distributions in panel (d) represent January and June temperatures between 2016 and 2020.

Snow temperature in Fig. 5c covers the period between 23 November 2019 and 16 January 2020. A 3 °C cold bias is apparent, but the mean falls within 1σ of the observations. Figure 5d presents distributions of the stacked January snow temperatures between 2016 and 2020, further highlighting this 3 °C cold bias during summer ($-\mathrm{37.1}\pm \mathrm{3.6}$ and $-\mathrm{34.4}\pm \mathrm{4.1}$ °C for Crocus outputs and observations, respectively). However, winter snow temperatures are well simulated with a 1 m mean of $-\mathrm{64.7}\pm \mathrm{4.5}$ and $-\mathrm{63.9}\pm \mathrm{5}$ °C from Crocus outputs and observations, respectively. The overestimation of simulated SSA in the top 10 cm may be linked to the cold bias in modelled snow temperature during summer via the reduced rate of near-surface snow metamorphism. Alternatively, high SSA may be the result of a recent precipitation event. During large precipitation events, the initial surface snow layer is buried rapidly, allowing less time for snow metamorphism (SSA decrease). This is particularly pronounced near the surface (top 5 cm) where the majority of metamorphism occurs in response to strong diurnal temperature gradients driven by solar radiation (Picard et al., 2012). This could result in sustained high SSA in the top 10 cm.

3.2.2 Simulated response of snowpack properties to surface forcing perturbations

Sensitivity tests were run using the Crocus model to assess the response of near-surface snowpack properties to perturbations in surface forcings. The following analysis uses optical radius as a measure of grain size, which is directly linked to SSA and the density of ice (Sect. 2.3.1), and focuses on the response of near-surface snow density and grain size to the six scenarios outlined in Table 2. We firstly assess the bulk changes in physical properties before looking at the variability with depth.

Figure 6 Comparison of density and grain size (r_opt) over the top 20 cm from Crocus sensitivity simulations. Distributions of (a) density and (b) grain size outputs from each test simulation are compared to the reference simulation (bin size is 5 kg m⁻³ and 0.005 mm). In panel (c), bars represent the percentage change in mean density and mean grain size for perturbations in SSI (orange), accumulation rate (green), and temperature (purple), with the decreased scenarios represented by the faded colour and the increased scenarios by the bold colour.

3.2.3 Bulk snowpack sensitivity

Numerous studies have suggested that modifications in near-surface density and grain size are key parameters influencing elemental fractionation during pore closure (e.g. Bender, 2002; Fujita et al., 2009). Figure 6 shows results from the sensitivity analysis, comparing differences in density and grain size between simulations outlined in Table 2. Overall, grain size is more sensitive to changes in surface forcing than density in the top 20 cm. Perturbations in SSI appear to have a minimal influence on density (<1 % change compared to the reference simulation): a 15 % decrease in SSI (SSI min) corresponds to an 8 % decrease in grain size, and an 11 % increase in SSI (SSI max) causes a 12 % increase in grain size. The magnitudes of change in both density and grain size are much larger under accumulation rate and temperature perturbations than under SSI perturbations.

Accumulation rate and temperature have opposite effects on density and grain size: a 4 °C increase in temperature (T max; 7 % increase) increases density and grain size by 8 % and 29 %, respectively, while increasing the accumulation rate to 4.1 cm yr⁻¹ (A max; 46 % increase) results in a 7 % decrease in density and a 10 % decrease in grain size. Mean density and grain size respond non-linearly to perturbations in all forcing parameters, as is shown by the magnitude of the increase in grain size from decreased accumulation rate being 3 times greater than the magnitude of decrease induced by an increase in accumulation rate (right panel of Fig. 6c, green bars). This is in line with the dependence of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ to the logarithm of accumulation rate documented in Fig. 3.

Figure 7 Variability in density and grain size with depth over the top 1 m. Each pair of simulations is represented by a colour: SSI in orange, accumulation rate in green, and temperature in purple. The faded line of each pair represents the 800 kyr minimum simulation, and the bold line represents the 800 kyr maximum simulation.

3.2.4 Depth variability sensitivity

An alternative (although possibly complementary) explanation for the mechanistic control of snow properties on elemental fractionation links to stratigraphic layering due to seasonality (Fujita et al., 2009). Here we explore the influence of layering by looking at the depth variability in density and grain size as a qualitative measure of stratification (Hörhold et al., 2011), assuming that higher variability indicates stronger layering. Variability is defined as the standard deviation of each depth interval (denoted σ) over the period between 1 January 2010 and 1 December 2020.

The σ values for each simulation are presented in Fig. 7 for the top 1 m of snowpack. In all runs, σ peaks near the surface and decreases with depth for both density and grain size (σ_ρ and σ_gs). Density variability for four out of the six test simulations is largely similar to the reference. High-accumulation rate (A max) and low-temperature (T min) are the exception, with an increase in σ_ρ of up to 5 kg m⁻³ compared to all other runs. The spread in σ_gs between simulations is broader than for σ_ρ over the top 50 cm, with T min, SSI min, A min, and T max all resulting in reduced variability compared to the reference run. The σ_gs values increase with an increase in SSI and accumulation rate throughout the top metre. At 100 cm, there is no significant difference in the σ_ρ and σ_gs for the different forcing scenarios, which is discussed in Sect. 4.4.

To summarise, outputs from the Crocus model indicate that insolation modifies mean grain size but has negligible effect on mean density in the top 20 cm. On the other hand, higher air temperature and lower accumulation rate result in increases in both mean density and mean grain size, although this change is not necessarily proportional. The depth variability of the two snowpack parameters is strongly influenced by increased accumulation rate, with decreased temperature also having a large effect on density variability, and increased SSI causing an increase in grain size variability.

4 Discussion

4.1 Evidence for non-SSI dependence of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$

The compilation of deep ice core $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records in Fig. 2 reinforces the widely documented anti-correlation between SSI and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (e.g. Oyabu et al., 2022; Bouchet et al., 2023). However, a comparison of EDC, Dome F, and South Pole (Fig. 2) reveals an additional influence of site conditions which cannot be explained by SSI alone. Indeed, our compilation of records from 15 sites shows that mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is strongly correlated with accumulation rate and temperature (Fig. 3). While similar spatial and temporal slopes characterise the SSI dependence of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, accumulation rate and temperature appear to be better predictors of site mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$. However, a key distinction between the temporal and spatial variations in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ lies in the range of variations in the three forcing parameters. For example, spatial variations in accumulation rate are much larger than the variations at a single site through time (2–110 cm yr⁻¹ compared to 1–4 cm yr⁻¹ in the EDC core). Conversely, the spatial range of SSI values is between 500–560 W m⁻², whereas the range at EDC, for instance, is between 470–600 W m⁻².

Figure 8 $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ as a function of accumulation rate and temperature. Contours show the predicted $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values calculated using the multiple linear regression of accumulation rate (A) and temperature (T) ($\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}=-\mathrm{0.03}\cdot T+\mathrm{2.96}\cdot \log \left(A\right)-\mathrm{14.6}$). Measured mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values are noted alongside each site. The accumulation rate vs. temperature regression (dashed, black line) is defined using a regression on temperature and accumulation rate for all sites in the data compilation ($A={\exp }^{\left(\right(\mathrm{0.069}\cdot T)+\mathrm{4.98})}$).

To further constrain the drivers of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability, we construct a multiple regression fit to parameterise mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ as a function of accumulation rate and temperature (Fig. 8). Combined, accumulation rate and temperature can explain up to 80 % of the total variance in mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (r²=0.81, adjusted r²=0.766). Including SSI to the multiple linear regression further improves the adjusted r² from 0.77 to 0.82. These results reiterate that $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is highly sensitive to accumulation rate, especially at low-accumulation sites. Although not explicitly explored, the site dependence of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ has been noted in previous studies. For example, Bazin et al. (2016) identified an offset in absolute $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values in the Vostok core compared to Dome F and EDC, similar to our observation in Fig. 2. By including records from numerous ice cores, we are able to definitively show that absolute values are systematically linked to site accumulation rate and temperature.

4.2 Millennial-scale variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records

Further investigation into drivers of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability in the EDC core reveals a millennial-scale climate signal in the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record. The timing of these anomalies in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ broadly corresponds to those in the δD record from EDC between 200–260 ka (Fig. 4), suggesting that millennial-scale variations can arise from changes in accumulation rate and/or temperature. The anomalies observed in the EDC $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record are not unique to this period but vary between time periods and sites (Sect. S2). Indeed, a similar effect was previously observed in the GISP2 core, where $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ was found to be positively correlated with local temperature over millennial timescales (i.e. Dansgaard–Oeschger (D–O) events) (Suwa and Bender, 2008b).

Previous studies suggest an overlap between the drivers of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability and those of TAC variability (Fujita et al., 2009; Lipenkov et al., 2011). Ice core $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC records exhibit slight differences in their spectral signals, whereby $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is dominated by precession (hence SSI pacing) and TAC is dominated by obliquity (hence the integrated summer insolation pacing). However, variability in both records is linked to local insolation via the modulation of near-surface snow properties, which ultimately influence pore closure processes (Lipenkov et al., 2011). Fujita et al. (2009) hypothesised that permeation, the process driving size-dependent fractionation of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, can explain about half the variation in TAC, with the rest being driven by effusion (not believed to substantially influence $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$). We therefore draw on TAC studies to aid our interpretation of millennial-scale variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$.

Several mechanisms have been proposed to explain millennial-scale variability in TAC and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, with apparent differences in the dominant mechanisms between hemispheres; TAC is anti-correlated with accumulation rate in the NGRIP core (Eicher et al., 2016) but positively correlated in the South Pole ice core (Epifanio et al., 2023). Eicher et al. (2016) attributed this anti-correlation between TAC and accumulation rate during D–O events at NGRIP, Greenland, to transient effects in the firn column. During the initial stage of a D–O event, rapid increases in accumulation rate increase overburden pressure, reducing the pore volume at close off and consequently decreasing TAC (Eicher et al., 2016). Similarly, Kobashi et al. (2015) suggested that rapid changes in overburden pressure could explain the negative correlation between $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ (and thus $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$) and accumulation rate over multi-decadal timescales at GISP2 via its influence on the permeation rate of gases in the deep firn. In contrast, a positive correlation between TAC and accumulation rate at South Pole is attributed to increased snow burial rates, leading to smaller grains and therefore increased pore volume at close-off (Epifanio et al., 2023). Such a grain size mechanism is understood to drive the integrated summer insolation signal in TAC records (Raynaud et al., 2007) as well as the SSI signal in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records (e.g. Bender, 2002).

Transient effects in the firn column in response to rapid climatic changes (linked to overburden pressure) are expected to be largely absent from Antarctic sites due to reduced magnitude in climate variability compared to Greenland. However, an accumulation-dependent grain size mechanism may explain the positive correlation between δD and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ at EDC. Epifanio et al. (2023) proposed that the contradictory behaviour of TAC between NGRIP and South Pole may be explained by varying responses of the firn to changes in accumulation rate, depending on the site's local climate conditions. They suggest that a grain size mechanism dominates TAC modulation at low-accumulation sites, while transient effects from rapid climatic changes are more important at warm, high-accumulation sites. We expect that the mechanisms driving $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ are also modulated by accumulation rate. Indeed, observations of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records from South Pole appear to support the in-phase coherence between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and δD in the EDC core, but only when δD is sufficiently low (Fig. S4 in the Supplement). A shifted, anti-phase relationship is apparent when δD is higher, which suggests that different mechanisms may be dominant under high-accumulation conditions.

4.3 Towards a mechanistic understanding of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability at low-accumulation sites

Previous studies have demonstrated the link between climate and snow metamorphism, showing that snow metamorphism is enhanced during summers with very low accumulation rates (Picard et al., 2012; Casado et al., 2021). However, the connection to pore closure processes has received less attention. One possible reason for this gap is the added complexity which would be required to implement surface snow metamorphism into firn models, which necessitates prior information such as surface grain size. The following sections propose some ideas to bridge this gap by examining the influence of local climate parameters (accumulation rate and temperature) and SSI on near-surface snow metamorphism and how these might modulate elemental fractionation during pore closure.

4.3.1 Surface snow properties

Variations in grain properties are thought to drive the SSI signal in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records (e.g. Bender, 2002; Fujita et al., 2009). Indeed, results from sensitivity tests in Fig. 6 show an increase in grain size under increased SSI, which is attributed to both stronger near-surface temperature gradients (Appendix B in Vionnet et al., 2012) and higher snow temperatures during summer (Fig. S4). Negligible change in simulated density under perturbations in SSI appears to support the hypothesis of a near-surface grain size mechanism (e.g. Bender, 2002), although this may be linked to the absence of vapour transport in the Crocus model. Increased snow metamorphism rate can also explain the increased grain size with increased temperature (Legagneux et al., 2003), while an increase in density results from increased compaction rates (Vionnet et al., 2012). The decrease in grain size with increased accumulation rate is the result of shorter residence time of a snow layer in the upper centimetres to metres of the snowpack, where temperature gradients are strongest (Picard et al., 2012).

It is widely accepted that permeation is the process by which small molecules escape during pore closure (e.g. Ikeda-Fukazawa et al., 2004; Huber et al., 2006; Severinghaus and Battle, 2006). Gregory et al. (2014) found that permeability is increased in high-density, large-grained firn due to a less complex pore structure. Indeed, Calonne et al. (2022) showed that grain size has a strong influence on permeability, such that for a given density permeability is increased with grain size. Our results show that high SSI and low accumulation rate are associated with larger grain size and lower $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$. This observation supports the hypothesis that increased SSI (decreased accumulation rate) ultimately leads to lower $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ due to the increased permeability associated with increased grain size (e.g. Severinghaus and Battle, 2006; Fujita et al., 2009), and suggests a similar mechanism may partially explain the accumulation rate dependence.

Identifying the mechanisms driving the accumulation rate dependence of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is complicated by the fact that, at polar sites, accumulation rate tends to covary with temperature (Fig. 8). Thus, an increase in grain size with increased near-surface residence time (decreased accumulation rate) will be countered by reduced snow metamorphism rates with a decrease in temperature (Fig. 6). These counteracting effects were invoked by Kawamura et al. (2007) to explain the absence of a 100 kyr periodicity in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records at Dome Fuji. However, a stronger correlation between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and accumulation rate than temperature in Fig. 3 suggests a dominant role of accumulation. Therefore, we propose that a grain size mechanism controlling pore structure at close-off, like that of SSI, may explain the positive correlation between accumulation rate and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ in the EDC core (Fig. 3b).

Previous studies propose that firn column characteristics, such as delta-age, will also modulate $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ via changes in accumulation rate and temperature (Severinghaus and Battle, 2006; Bazin et al., 2016). One mechanism links to the time for pores to close off. At low-accumulation sites, pores will take longer to fully close; thus, they experience more elemental fractionation (Severinghaus and Battle, 2006). Additional hypotheses point to LIZ thickness and pore space geometry (or tortuosity) (e.g. Fujita et al., 2009). Low-accumulation sites tend to have thinner LIZs (e.g. Landais et al., 2006; Witrant et al., 2012) and less tortuous pore structure in deep firn associated with large, rounded grains (Gregory et al., 2014). Such characteristics are associated with increased gas diffusivity that would facilitate the removal of O₂-enriched gas back to the atmosphere. In contrast, at high-accumulation sites, the associated decrease in gas diffusivity in the LIZ would lead to a build up of a stagnant air enriched in O₂, resulting in the trapping of gas with relatively high $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$. Alternatively, Hutterli et al. (2009) proposed that varying near-surface temperature gradients under different SSI intensities modulates anisotropy of the snow, leading to vertically elongated pores (as was confirmed by Leinss et al., 2020) under high SSI (low-accumulation) conditions. They argue that elongated pores facilitate vertical diffusivity in the LIZ, leading to greater fractionation of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (Hutterli et al., 2009), an alternative but complementary mechanism to LIZ thickness. A simplified and non-exhaustive illustration of these various interactions can be found in Fig. S6.

The aforementioned mechanisms relating to firn column characteristics indicate a positive correlation between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and accumulation rate. This is consistent with Fig. 3b but cannot be directly supported by our sensitivity results. Based on the sensitivity tests in Fig. 6, a link to surface density is not clear. We expect that the role of density in modulating $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ derives from the grain size dependence of densification rates invoked by Freitag et al. (2004), providing an additional link between near-surface grain size and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, as is outlined in the following section.

4.3.2 Depth-dependent variability as a proxy for layering

Density stratification in deep firn has also been invoked to modulate $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (Fujita et al., 2009). Results from the Crocus model are tentatively used to infer the sensitivity of near-surface density and grain size variability to perturbations in input forcing parameters. Figure 7 shows that grain size variability is increased with an increase in both SSI and accumulation rate. For SSI, this is expected given the enhanced snow metamorphism during summer due to higher snow temperatures, while winter conditions are largely unchanged (Fig. S5). Decreased grain size variability in both the increased and decreased temperature simulations is likely linked to the way in which temperature forcing is perturbed. By applying a constant increase in air temperature, the strength of metamorphism is increased during winter due to higher temperatures (Legagneux et al., 2003; Flanner and Zender, 2006) but largely unchanged during summer due to the dominant influence of insolation on summer snow temperature (Fig. S5), resulting in homogeneity within the snowpack.

Linking variability near the surface with variability in deep firn is not trivial. Indeed, a study by Hörhold et al. (2011) used a compilation of firn cores from numerous polar sites to show that density variability in deep firn is positively correlated with local accumulation rate and temperature but anti-correlated with near-surface density variability. However, this anti-correlation was not observed by Inoue et al. (2024) when comparing firn cores in the Dome Fuji area. Additional consideration, both from data and modelling approaches, is required for density inversions in the firn column – whereby relatively low-density layers in the upper part of the firn become relatively high-density layers below the density inversion depth, due to preferential deformation in the upper firn (Freitag et al., 2004; Fujita et al., 2009). Freitag et al. (2004) proposed that the inversion is largely driven by grain size, such that initially low-density layers associated with large grain size become relatively high-density layers with large grain size below the density inversion. This has since been supported by Gregory et al. (2014), who further show that near-surface grain size determines the density at which pores close off. They show that large-grained, high-density layers are more permeable and, hence, close off deeper than fine-grained, low-density layers. In addition, Hörhold et al. (2012) suggested that stratification at depth is influenced by impurity content, with impurity-rich layers being more susceptible to densification. Given the array of factors proposed to influence deep firn density variability, we focus on the role of grain size variability.

No inversion is observed for grain size variability. Deep firn layers retain evidence of near-surface grain size, which is determined by snow metamorphism and accumulation rate (Gregory et al., 2014), as mentioned previously. While our sensitivity tests are constrained to the near surface, we draw on previous studies to suggest potential links with deep firn properties and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ fractionation. Given the presence of a density inversion at Dome C (Hörhold et al., 2011), we may expect that relatively increased variability in near-surface grain size from our sensitivity tests would translate to increased grain size and density variability in the deep firn. As such, increased grain size variability in our SSI max simulation (Fig. 7) supports the conclusions of Fujita et al. (2009), that elemental fractionation is enhanced under high-SSI conditions due partially to layering.

The link to accumulation rate is less clear. As for SSI, we observe an increase in grain size variability with increased accumulation rate. However, according to the layering hypothesis, this should result in more $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ fractionation at close off, which contradicts our observations of a positive correlation between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and accumulation rate (Fig. 3b). Instead, our observations suggest that a decrease in mean grain size with increased accumulation rate is dominant over the layering effect. Future work would benefit from an inter-site comparison of deep firn properties alongside gas measurements in order to disentangle these effects.

Extracting concrete conclusions from our variability analysis and extrapolating these into the deep firn is inhibited by both the ascribed forcing perturbations and the aggregation scheme of the model, which is particularly sensitive to changes in accumulation rate. We thus conclude that the simulations performed with the Crocus model can support a mechanism linked to average grain size on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ but not a mechanism implying density or grain size variability. However, our conclusion does not rule out the effect of layering or grain size variability on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability but highlights a limitation in our study, as explained in the next section.

4.4 Limitations and perspectives

A number of limitations have been mentioned throughout our study and warrant some further discussion. As mentioned in Sect. 4.3.2, our study is unable to fully consider the sensitivity of layering to perturbations in forcing parameters. Firstly, near-surface density and grain size variability are not fully captured by the model which can be partially attributed to the absence of snow transport by the wind in the standard version of Crocus (Libois et al., 2014). While wind is known to strongly influence snow properties (e.g. Inoue et al., 2024), we are limited in our understanding of winds throughout the past 800 kyr; therefore, wind effects are not considered in our study. Secondly, the aggregation of layers with depth in Crocus make it difficult to focus on both the fine-layered near-surface snow and the propagation of stratified layers into the deep firn. Accurate assessment of the layering effect would require a new dedicated snow model preserving individual snow layers and properties over a wide depth range, from the surface to LIZ, and at high resolution.

While the single-parameter sensitivity tests presented here provide useful insights for understanding physical mechanisms, they do not account for complex compound effects associated with the covariance of accumulation rate and temperature at polar sites. Additional tests perturbing both accumulation rate and temperature simultaneously indicate that snowpack properties are very sensitive to the ascribed accumulation rate and temperature values, i.e. for glacial temperature reconstructions which are debated to have been overestimated by as much as 5 °C at Dome C (Buizert et al., 2021). Moreover, the additive effects from the single-parameter sensitivity tests for accumulation rate and temperature do not equal the effects of multi-parameter sensitivity simulations (i.e. perturbing both accumulation rate and temperature). This discrepancy highlights the complex interactions between quasi-covarying accumulation rate and temperature on snow properties. Indeed, a cancellation of the accumulation rate and temperature effects was invoked by Kawamura et al. (2007) to explain the absence of a 100 kyr periodicity in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ at Dome Fuji. Given the uncertainties in climate reconstructions and the limitations of the model, we refrained from using multi-parameter simulations for interpretation of empirical results. Future studies would benefit from obtaining $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements in open and closed porosity at sites with distinct characteristics to better constrain the physical mechanisms.

Finally, given that the high-resolution $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records tend to be from Antarctic sites (e.g. Fig. 2), our study is slightly biased towards low-accumulation conditions. Indeed, initial temporal analysis of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability at South Pole suggests that different mechanisms may be dominant at relatively high-accumulation sites, as has been previously suggested (Suwa and Bender, 2008b; Kobashi et al., 2015), but this was outside the scope of this study. Future studies focusing on the dominance of different mechanisms under different site conditions would therefore benefit from obtaining high-resolution measurements from sites with distinct characteristics. Furthermore, a comparison of deep firn properties alongside $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements from the open and closed porosity would aid a mechanistic understanding.

5 Conclusions

We present a compilation of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records measured on multiple ice cores from Antarctica and Greenland to develop a mechanistic understanding of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability. Analysis of both temporal (single-site) and spatial (multi-site) variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ presents new evidence of a dependence on local climate (accumulation rate and temperature), in addition to the well-documented SSI dependence. High-resolution measurements from the EDC ice core reveal millennial-scale variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ in phase with δD records when both parameters are plotted on the AICC2023 ice-age scale. The inter-site analysis showed an increase in mean $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ for sites with higher accumulation rate and temperature, which is analogous with the temporal analysis from EDC, showing $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ to increase together with δD.

We suggest that a mechanism relating to firn physical properties can partially explain both the influence of SSI and local climate on $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability. Sensitivity tests using the Crocus model show that grain size is very responsive to perturbations in SSI, accumulation rate, and air temperature, while density responds to all but SSI perturbations. Our findings support the hypothesis that a grain size mechanism is a dominant driver of elemental fractionation at low-accumulation sites, such that increased grain size for a given density facilitates O₂ expulsion via enhanced permeability. Furthermore, our results support the hypothesis that the presence, or lack thereof, of a local climatic signal in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ variability is due to the delicate balance between the counter-effects of accumulation rate and temperature on grain properties. However, the inter-site results suggest that low-accumulation, low-temperature sites experience stronger elemental fractionation, having a comparable effect to high insolation.

While our findings from the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data compilation are supported by the Crocus sensitivity tests, we acknowledge that additional mechanisms are at play. In particular, the influence of deep firn layering – itself linked to surface snow metamorphism – could not be tested fully in this study but is believed to play a major role in modulating $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$. Determining the relative influence of stratification, firn physical properties, and residence time in the lock-in zone using a combination of firn models and observations would be useful for future studies.