The Antarctic Ice Core Chronology 2023 (AICC2023) chronological framework and associated timescale for the European Project for Ice Coring in Antarctica (EPICA) Dome C ice core

The EPICA (European Project for Ice Coring in Antarctica) Dome C (EDC) ice core drilling in East Antarctica reaches a depth of 3260 m. The reference EDC chronology, the AICC2012 (Antarctic Ice Core Chronology 2012), provides an age vs. depth relationship covering the last 800 kyr (thousands of years), with an absolute uncertainty rising up to 8000 years at the bottom of the ice core. The origins of this relatively large uncertainty are twofold: (1) the δ¹⁸O_atm, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and total air content (TAC) records are poorly resolved and show large gaps over the last 800 kyr, and (2) large uncertainties are associated with their orbital targets. Here, we present new highly resolved δ¹⁸O_atm, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and δ¹⁵N measurements for the EDC ice core covering the last five glacial–interglacial transitions; a new low-resolution TAC record over the period 440–800 ka BP (ka: 1000 years before 1950); and novel absolute ⁸¹Kr ages. We have compiled chronological and glaciological information including novel orbital age markers from new data on the EDC ice core as well as accurate firn modeling estimates in a Bayesian dating tool to construct the new AICC2023 chronology. For the first time, three orbital tools are used simultaneously. Hence, it is possible to observe that they are consistent with each other and with the other age markers over most of the last 800 kyr (70 %). This, in turn, gives us confidence in the new AICC2023 chronology. The average uncertainty in the ice chronology is reduced from 1700 to 900 years in AICC2023 over the last 800 kyr (1σ). The new timescale diverges from AICC2012 and suggests age shifts reaching 3800 years towards older ages over marine isotope stages (MISs) 5, 11 and 19. But the coherency between the new AICC2023 timescale and independent chronologies of other archives (Italian Lacustrine succession from Sulmona Basin, Dome Fuji ice core and northern Alpine speleothems) is improved by 1000 to 2000 years over these time intervals.

1 Introduction

1.1 Building age scales for deep polar ice cores

1.1.1 Motivation

Deep polar ice cores are unique archives of past climate, and their investigation is valuable to study mechanisms governing the Earth's climate variations. Precise chronologies are key to identifying the successions and lengths of climatic events, along with exploring phase relationships between the external forcing (changes in the Earth's orbit) and the diverse climatic responses (variations in temperature and atmospheric greenhouse gas concentrations). To date ice cores, we need to construct two separate chronologies: one for the ice and one for the younger air trapped in bubbles. Due to the thinning of horizontal ice layers as depth increases, a wide time span of paleoclimatic information is stored within the deepest part of the ice sheet. Therefore, many of the ice core community's ongoing efforts focus on improving deep-ice-core timescales for ice and gas phases, as well as extending them further back in time (Crotti et al., 2021; Oyabu et al., 2022). Ice cores drilled at sites characterized by a high accumulation rate of snow at the surface (10 to 30 cm a⁻¹) can be dated by counting ice layers deposited year after year (Svensson et al., 2008; Sigl et al., 2016). In contrast, East Antarctica sites are associated with very low accumulation rates (1 to 5 cm a⁻¹), which prevent annual layers from being identified and counted. As a consequence, chronologies of ice cores at low-accumulation sites are commonly established using ice flow and accumulation models (Nye, 1959; Schwander et al., 2001), then tied up with chronological and glaciological constraints (Veres et al., 2013; Bazin et al., 2013; Parrenin et al., 2017).

1.1.2 Glaciological modeling

Glaciological modeling has been historically used to date Greenlandic and Antarctic ice cores. A unidimensional ice flow model was first applied to the Camp Century ice core (Dansgaard and Johnsen, 1969) and later to other ice cores such as the ones drilled at EPICA (European Project for Ice Coring in Antarctica) Dome C (EDC) and Dome Fuji (EPICA Community Members, 2004; Parrenin et al., 2007a). First, water isotope (δD or δ¹⁸O) measurements provide estimates of past evolution of the accumulation rate of snow and temperature at the surface. Then, an ice flow model (Parrenin et al., 2004) takes as inputs past accumulation together with a vertical velocity–depth profile through the ice sheet to determine the thinning of annual snow/ice layers in time and therefore the ice timescale. This approach is very sensitive to some poorly known parameters including boundary conditions such as bedrock topography, geothermal properties or subglacial sliding. For this reason, the glaciological modeling approach is complemented with chronological constraints (gas or ice age known at certain depth levels).

1.1.3 Chronological constraints derived from measurements

Chronological constraints obtained either by measurement of radionuclides or by synchronization to a curve of reference are established for both ice and gas timescales. For building long chronologies, some time constraints can be obtained from the ¹⁰Be series measured in ice. The ¹⁰Be cosmogenic nuclide is produced at different rates depending on the solar activity, and its arrival on Earth is modulated by the strength of the Earth's magnetic field (Yiou et al., 1997; Raisbeck et al., 2007; Heaton et al., 2021). Some links hence exist between ¹⁰Be flux and precisely dated magnetic events such as the Laschamp excursion, an abrupt decline in the geomagnetic field magnitude occurring at about 41 ka BP and visible as a positive excursion in the ¹⁰Be flux records in ice cores (Lascu et al., 2016; Raisbeck et al., 2017). ⁴⁰Ar measurements in the gas phase of Antarctic ice cores also provide dating constraints for old ice, especially for non-continuous stratigraphic sequences (Yan et al., 2019). ⁴⁰Ar is produced in solid earth by the radioactive decay of ⁴⁰K, leading to an increasing concentration of ⁴⁰Ar in the atmosphere at a rate of 0.066±0.006 ‰ Myr⁻¹ (Bender et al., 2008). Recently, the possibility of measuring ⁸¹Kr in ice samples of a few kilograms gave a new absolute dating tool for ice cores (Jiang et al., 2020). ⁸¹Kr is a radioactive isotope that is suitable for dating ice cores in the range from 0.03 to 1.3 Ma BP (Ma: million years before 1950), making it perfectly adapted for Antarctic ice core dating (Buizert et al., 2014; Crotti et al., 2021).

To further constrain the oldest ice core chronologies, the so-called “orbital dating” tools are also used. These tools consist of aligning some tracers measured in ice cores to the Earth orbital series, called targets, whose fluctuations in time are accurately calculated from the known variations in orbital parameters (Berger, 1978; Laskar et al., 2011). The synchronization of the orbital tracer with its target provides ice or gas age constraints. So far, three orbital dating tools have been developed: δ¹⁸O of O₂ (δ¹⁸O_atm), $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and total air content (TAC). The δ¹⁸O_atm was typically aligned with the precession parameter (or with the 21 June insolation at 65^∘ N) delayed by 5000 years because such a lag between δ¹⁸O_atm and its orbital target was observed during the last deglaciation (Shackleton, 2000; Dreyfus et al., 2007). However, variations in the phasing between δ¹⁸O_atm and precession have been suspected (Jouzel et al., 2002) and identified since (Bazin et al., 2016). In particular, millennial-scale events (as Heinrich-like events) occurring during deglaciations have been shown to delay the response of δ¹⁸O_atm to orbital forcing (Extier et al., 2018a). Because there was a significant unpredictability in the lag between δ¹⁸O_atm and its orbital target, a large uncertainty in the δ¹⁸O_atm-based tie points (up to 6000 years) was assigned in the construction of the AICC2012 (Antarctic Ice Core Chronology 2012; Bazin et al., 2013). To improve the accuracy of the gas timescale, Extier et al. (2018a) rather aligned the variations in δ¹⁸O_atm to the δ¹⁸O_calcite recorded in absolute-dated East Asian speleothems between 640 and 100 ka BP. Indeed, the two records show similar orbital (related to the 21 July insolation at 65^∘ N) and millennial variabilities, which may correspond to southward shifts in the Intertropical Convergence Zone (ITCZ) position, themselves linked to Heinrich-like events as supported by the modeling study of Reutenauer et al. (2015).

In parallel, Bender (2002) observed that the elemental ratio $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ of air trapped in the Vostok ice core appears to vary in phase with the 21 December insolation at 78^∘ S (Vostok latitude) between 400 and 160 ka BP. Subsequent observations led Bender (2002) to assert that local summer solstice insolation affects near-surface snow metamorphism and that this imprint is preserved as snow densifies in the firn and, later on, affects the ratio $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measured in air bubbles formed at the lock-in zone. Wiggle matching between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and local summer solstice insolation has been used to construct orbital timescales for the Dome Fuji, Vostok and EDC ice cores reaching back 360, 400 and 800 ka BP, respectively, with a chronological uncertainty for each $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ tie point estimated between 250 and 4000 years (Kawamura et al., 2007; Suwa and Bender, 2008; Bazin et al., 2013; Oyabu et al., 2022). Finally, Raynaud et al. (2007) found very similar spectral properties between the TAC record of EDC and the integrated summer insolation (ISI) at 75^∘ S obtained by a summation over a year of all daily local insolation above a certain threshold over the last 440 kyr. As for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, these similarities may be explained by the insolation imprint in near-surface snow well preserved down to the lock-in zone, where it could affect the air content in deep ice, although the physical mechanisms involved during the snow and firn densification for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC are likely different (Lipenkov et al., 2011). Lately, Bazin et al. (2013) made use of TAC to constrain Vostok and EDC ice core chronologies back to 430 ka BP with an uncertainty for each TAC tie point varying between 3000 and 7000 years. Although these three orbital tools complement each other (TAC and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ inferred ages agree within less than 1000 years between 390 and 160 ka BP for the Vostok ice core; Lipenkov et al., 2011), they hardly ever have been employed together. Plus, they are often associated with large uncertainties (reaching 7000 years), which lie in the choice of the appropriate orbital target, in its alignment with ice core records that can be ambiguous during periods of low eccentricity in the Earth's orbit (leading to low-amplitude insolation variations) and in the poor quality of the signals measured in the deepest section of the cores.

To connect ice and gas timescales, the estimation of the lock-in depth (LID), indicating the lowest depth where the air is trapped in enclosed bubbles, and diffusivity becomes effectively zero (Buizert et al., 2013), is used to calculate the ice–gas age difference. Measurements of δ¹⁵N from N₂ yield a first estimate of this depth, and the LID can also be calculated with firn densification modeling (Goujon et al., 2003; Bréant et al., 2017).

For many years, each polar ice core was characterized by its singular timescale, which was not naturally consistent with other ice core timescales. To address this issue, other measurements provide relative dating constraints (stratigraphic links), improving the coherency between timescales of ice cores from both hemispheres. The synchronization of globally well-mixed atmospheric methane ice core records gives tie points with an accuracy of a few decades to several centuries (60–1500 years) (Lemieux-Dudon et al., 2010; Epifanio et al., 2020). Climate-independent events, such as large volcanic eruptions, can be observed in ice cores from Greenland and Antarctica via singular patterns of the distribution of sulfate. Identification of these deposits permits the precise synchronization of several ice cores (within 5 to 150 years) (Svensson et al., 2020).

1.1.4 Bayesian dating tools

In order to integrate stratigraphic matching, independent synchronization and absolute dating constraints as well as glaciological modeling to produce coherent ice core chronologies, researchers developed Bayesian dating tools such as Datice (Lemieux-Dudon et al., 2010), IceChrono1 (Parrenin et al., 2015) and Paleochrono (Parrenin et al., 2021). These tools use an inverse method combining all chronological information to provide a coherent age scale for several ice cores. These probabilistic tools adjust prior estimates of ice and gas chronologies built with a glaciological model (background scenario) so that they respect chronological constraints.

Here we focus on the chronology of the EDC deep ice core. The EPICA project provided two cores in East Antarctica, including one at Dome C (EDC). The second (and final) drilling attempt at Dome C gave the 3260 m long EDC99 core, whose drilling was willingly stopped at 15 m above bedrock due to the expected presence of meltwater. EDC furnishes the oldest continuous ice core record so far, covering the last 800 kyr (EPICA Community Members, 2004; Jouzel et al., 2007).

1.2 The AICC2012 chronology

Bazin et al. (2013) and Veres et al. (2013) used the probabilistic dating tool Datice to establish the coherent AICC2012 chronology back to 800 ka BP for five ice cores, including EDC, Vostok, EPICA Dronning Maud Land ice core (EDML), North GReenland Ice core Project (NGRIP) and Talos Dome Ice core (TALDICE). To determine the EDC age scale, they used various orbital dating constraints including 39 tie points attached to a 6000-year uncertainty derived from δ¹⁸O_atm tuning to a 5000-year delayed precession between 800 and 300 ka BP, 20 tie points associated with a 4000-year uncertainty from $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ alignment to local summer solstice insolation between 800 and 300 ka BP, and 14 tie points linked to an uncertainty between 3000 and 7000 years using TAC synchronized to integrated summer insolation between 430 and 0 ka BP. However, due to the lack of data for the orbital dating approach, the AICC2012 1σ uncertainty is of 1700 years on average, reaching 8000 years at the bottom of the core. The origins of AICC2012 uncertainty can be divided in the following points: (i) some inherent dissimilarities between δ¹⁸O_atm, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC series and their curve-shaped orbital target; (ii) discontinuity and poor quality of the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC records; (iii) uncertainty in the phasing between δ¹⁸O_atm and precession; (iv) poor constraint on the LID scenario due to a disagreement between δ¹⁵N data and firn modeling estimates (Bréant et al., 2017).

1.3 The new AICC2023 chronology

It is now possible to address each source of uncertainty thanks to recent advances. (i) Since AICC2012, the δ¹⁸O_atm and TAC records have been extended, now covering the last 800 kyr (Extier et al., 2018b). In addition, new highly resolved δ¹⁸O_atm and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements are available over several glacial terminations (TII, III, IV, V and VI) (Grisart, 2023). (ii) Extier et al. (2018a) recently suggested a δ¹⁸O_atm-based timescale using δ¹⁸O_calcite of East Asian speleothems as an alternative tuning target to precession. This choice reduces the chronological uncertainty between 640 and 100 ka BP. (iii) Finally, new highly resolved δ¹⁵N data covering the terminations II to VI are available (Grisart, 2023). In parallel, firn densification models have been progressively improved, and the model described in Bréant et al. (2017) can be employed to estimate LID evolution in the past when δ¹⁵N data are still missing.

In this work, we implement new absolute age constraints spanning the last 800 kyr derived from ⁸¹Kr measured in air trapped in the EDC ice core as well as new orbital age constraints obtained by synchronizing up-to-date EDC records with their orbital target. We combine these data with recent volcanic matching and methane record synchronization, which provide additional stratigraphic links, relating EDC to other ice cores over the past 122 kyr (Baumgartner et al., 2014; Svensson et al., 2020). Finally, we propose the new chronology AICC2023 with reduced chronological uncertainties. AICC2023 is recommended as the new official age scale for the EPICA ice cores by the EPICA Scientific Steering Committee (Wolff, 2023).

2 Methods

2.1 Dating strategy

The Paleochrono Python software is a probabilistic dating tool similar to Datice and Icechrono1 with improved mathematical, numerical and programming capacities (Parrenin et al., 2021). The dating strategy of Paleochrono relies on the Bayesian inference of three glaciological functions forming the input background scenario: accumulation rate (A), thinning of annual ice layers (τ) and lock-in depth (LID). The three variables evolve along the ice core depth z and are used to estimate the ice (ψ) and gas (χ) age profiles as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathit{\psi }\left(z\right)=\underset{\mathrm{0}}{\overset{z}{\int }}{\displaystyle \frac{D\left(z^{\prime }\right)}{\mathit{\tau }\left(z^{\prime }\right)A\left(z^{\prime }\right)}}\mathrm{d}{z}^{\prime }\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\chi }\left(z\right)=\mathit{\psi }(z-\mathrm{\Delta }\mathrm{depth}(z\left)\right)\text{(3)}& {\displaystyle }& {\displaystyle }\underset{z-\mathrm{\Delta }\mathrm{depth}\left(z\right)}{\overset{z}{\int }}{\displaystyle \frac{D\left(z^{\prime }\right)}{\mathit{\tau }\left(z^{\prime }\right)}}\mathrm{d}{z}^{\prime }=\mathrm{LID}\left(z\right)\times {\left.{\displaystyle \frac{D}{\mathit{\tau }}}\right|}_{\mathrm{firn}}^{\mathrm{0}},\end{array}$$

where D is the relative density of the snow/ice and ${\left.\frac{D}{\mathit{\tau }}\right|}_{\mathrm{firn}}^{\mathrm{0}}$the average value of $\frac{D}{\mathit{\tau }}$ in the firn when the air particle was at the lock-in depth (this parameter is usually ∼0.7; Parrenin et al., 2012). The age scales are further constrained to respect chronological constraints identified from observations. To specify the credibility of the background scenario for the age scales and the chronological constraints, the glaciological functions (accumulation, thinning and LID) and the chronological information can be mathematically expressed as probability densities, which are presumed to be Gaussian and independent (i.e., decorrelated between them). Thus, the inference is based on the least square optimization method (implying all probability densities are Gaussian). It is numerically solved using the trust region algorithm (assuming that the model is roughly linear around the solution), and the Jacobian of the model is evaluated analytically for an improved computation time. As a result, the best adjustment between the background scenario and chronological observations is found, providing the most probable scenario as a posterior evaluation of the three glaciological functions and hence chronologies for ice and air. For each ice core, the input files for Paleochrono are the following: (i) the background values of the three glaciological functions with depth; (ii) gas and ice stratigraphic links; (iii) gas- and ice-dated horizons, which are tie points derived for one core from absolute and synchronization dating methods; (iv) gas and ice intervals of known durations; and (v) depth difference estimates between the same event recorded in the gas and ice matrix (Δdepth). Specific relative or absolute uncertainties are attached to each of these parameters in each input file.

In this study, we added numerous gas- and ice-dated horizons for EDC as well as an updated background scenario for the LID. Then, to construct a new chronology for the EDC ice core that is consistent with the timescales of the Vostok, TALDICE, EDML and NGRIP ice cores, we followed the same strategy as for the construction of AICC2012. Glaciological background parameters and dating constraints for the Vostok, TALDICE, EDML, NGRIP and EDC drillings are compiled in one run of Paleochrono to obtain AICC2023. The Vostok, TALDICE, EDML and NGRIP background parameters and dating constraints are extracted from Bazin et al. (2013), except for (i) new Vostok gas age constraints determined from the alignment of δ¹⁸O_atm and East Asian δ¹⁸O_calcite records as for EDC (see Fig. S10), (ii) new TALDICE background parameters and age constraints from Crotti et al. (2021), and (iii) corrected LID background scenarios for the Vostok and EDML sites (see Fig. S11). In order to prevent any confusion with reference ice core timescales, the new AICC2023 chronology for NGRIP is compelled to respect exactly the layer-counted GICC05 timescale through absolute tie points placed at 1 m intervals over the last 60 kyr (Andersen et al., 2006). For this reason, we did not use the methodology described by Lemieux-Dudon et al. (2015), who implemented layer counting as a constraint on the duration of events in the dating tool, inducing a slight shift (maximum 410 years) in the AICC2012 timescale. The resulting Paleochrono experiment provides the new official chronology AICC2023 for the EDC ice core. The contingent timescales obtained for the four other sites are not the subject of this study but are also provided (see “Data availability” section). We acknowledge the exclusion of the WAIS (West Antarctic Ice Sheet) Divide ice core (WDC) from the construction of the AICC2023 age scale as for the AICC2012 age scale. Over the last 60 kyr, though, we recommend the use of timescales tied to the WAIS Divide 2014 age model (WD2014; Buizert et al., 2015; Sigl et al., 2016). A correspondence between AICC2012, AICC2023 and WD2014 age models based on the volcanic synchronization of WDC and EDC using sulfate data (Buizert et al., 2018) is provided over the 0–58 ka BP period (that is to say for the section above the depth of 915 m for the EDC ice core; see “Data availability” section).

2.2 Analytical method

2.2.1 δ¹⁸O_atm, δO${}_{\mathrm{2}}/$N₂ and δ¹⁵N

The measurements of the isotopic and elemental compositions of O₂ and N₂ were performed by Grisart (2023) at the Laboratory for Sciences of Climate and Environment (LSCE) following the method described by Bréant et al. (2019) and Extier et al. (2018a). The air trapped in the EDC ice core is extracted using the semi-automatic line, which eliminates CO₂ and H₂O. Samples of 30 to 40 g are prepared in a cold environment (−20 ^∘C), their exterior layer (3–5 mm) is removed so that there is no exchange with atmospheric air, and each sample is cut in two replicates. Each day, three ice samples (and replicates) are placed in six flasks, and the atmospheric air is evacuated from the flasks. Samples are then melted and left at ambient temperature for approximately 1 h 30 min in order to extract the air trapped in ice samples. The extracted air is then cryogenically trapped within a dedicated manifold immersed in liquid helium (Bazin et al., 2016). Along the way to the cryogenic trap, the air goes through cold traps to remove CO₂ and H₂O. Two additional samples containing exterior modern air are processed through the same line every day for calibration and for monitoring the analytical set-up. Lastly, the δ¹⁵N, δ¹⁸O of O₂ and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ of each sample are measured by a dual-inlet Delta V plus (Thermo Electron Corporation) mass spectrometer.

Classical corrections are applied on the measurements (pressure imbalance, chemical slopes, as per Landais et al., 2003). In addition, δ¹⁵N data are used to get the values of atmospheric δ¹⁸O of O₂ and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ after gravitational fractionation occurred in the firn so that δ¹⁸O_atm=δ¹⁸O of ${\mathrm{O}}_{\mathrm{2}}-\mathrm{2}\times {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}$ and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}\left(\mathrm{corr}\right)}=\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}\left(\mathrm{raw}\right)}-\mathrm{4}\times {\mathit{\delta }}^{\mathrm{15}}\mathrm{N}$ (Landais et al., 2003; Bazin et al., 2016; Extier et al., 2018a). Note that our samples were stored at −50 ^∘C after drilling to minimize the gas loss effect. As a consequence, no correction for gas loss was applied (see Sect. S1 in the Supplement), and if gas loss may explain a slight scattering in the data, the peak positions are not affected.

Existing and new EDC data are compiled in Table 1. The dataset's resulting pooled standard deviations for the new measurements are of 0.006 ‰, 0.03 ‰ and 0.4 ‰ for δ¹⁵N, δ¹⁸O_atm and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, respectively.

Table 1 Information on isotopic and elemental compositions measured in air trapped in the EDC ice core. ^* Details on storage and measurement conditions of δO${}_{\mathrm{2}}/$N₂ are available in Sect. S1 in the Supplement.

2.2.2 Total air content

The TAC record was measured in the entire EDC ice core at the IGE (Institute of Environmental Geosciences) following the barometric method first described by Lipenkov et al. (1995). The TAC record measured in the younger part of the core (400–0 ka BP) has been published in Raynaud et al. (2007) (Table 2). TAC estimates need to be corrected for the cut-bubble effect. After correction, the uncertainty in TAC values is of about 1 %, and the analysis replicability is better than 1 %.

Table 2 Information on TAC measurements in the EDC ice core.

2.2.3 ⁸¹Kr extraction and analysis

The analytical method is the same as described by Crotti et al. (2021). Three 6 kg ice samples are taken from the bottom part of EDC, and a slight shaving (1 mm) of the external layer is performed before processing. The air extraction is performed through a manual extraction line following the protocol described in Tian et al. (2019). The ice sample is placed in a 40 L stainless-steel chamber. The atmospheric air is pumped while the chamber is kept at −20 ^∘C. The air is then slowly extracted, passing through a water trap, and compressed in a stainless-steel cylinder. The three cylinders are sent to the University of Science and Technology of China (USTC, Hefei, China) for krypton extraction and analysis. Krypton extraction is performed after the methodology of Dong et al. (2019), who set up an automated system for dual separation of argon and krypton, composed of a titanium getter module followed by a gas-chromatography separator module. The extracted krypton is analyzed by the Atom Trap Trace Analysis (ATTA) instrument set up at the Laser Laboratory for Trace Analysis and Precision Measurement (LLTAPM, USTC, Hefei, China), giving the ⁸¹Kr abundance R₈₁ in the sample. R₈₁ is determined by the number of counted ⁸¹Kr atoms in the sample as compared to the atmospheric reference. The anthropogenic ⁸⁵Kr is measured simultaneously with ⁸¹Kr to control any present-day air contamination. Here, the ⁸⁵Kr abundance measured in ice samples is inferior to the detection limit, so contamination has occurred.

From the ⁸¹Kr abundance, it is possible to estimate ⁸¹Kr radioactive decay and to calculate the ice sample's age. As a noble gas isotope, ⁸¹Kr is globally mixed in the atmosphere, and its decay cannot be affected by complex chemical reactions (Lu et al., 2014). The ⁸¹Kr half-life ($t_{\mathrm{1}/\mathrm{2}}$) is estimated to be $\simeq \mathrm{229}\pm \mathrm{11}$ kyr (Baglin, 2008). ⁸¹Kr age can be calculated as per the following equation:

$$\begin{array}{}\text{(4)}& \mathrm{age}=-{\displaystyle \frac{t_{\frac{\mathrm{1}}{\mathrm{2}}}}{\ln \left(\mathrm{2}\right)}}\times \ln \left({\mathrm{R}}_{\mathrm{81}}\right).\end{array}$$

The atmospheric abundance of ⁸¹Kr is not constant in the past, and its value is corrected using reconstruction of the geomagnetic field intensity (Zappala et al., 2020). The error in ⁸¹Kr age estimates is estimated from the statistical error in atom counting, from the uncertainty in ⁸¹Kr half-life (inducing a systematic age error) and from the size of the sample (larger sample resulting in a smaller uncertainty).

2.3 Firn model

Firn densification models have been progressively improved over the years (Herron and Langway, 1980; Alley, 1987; Arthern et al., 2010; Ligtenberg et al., 2011; Kuipers Munneke et al., 2015; Oraschewski and Grinsted, 2022). While these models generally explain the evolution of δ¹⁵N in time well through changes in the LID, they fail to reproduce values of δ¹⁵N in some regions, including coastal areas and cold and low-accumulation sites such as EDC (Capron et al., 2013). This disagreement can be explained by an inaccurate estimate of glacial temperature and accumulation rate at the surface (Buizert, 2021) and/or by the impossibility of tuning empirical firn models to sites with no present-day equivalent in terms of temperature and accumulation rate (Dreyfus et al., 2010; Capron et al., 2013). Recently, the firn model described in Bréant et al. (2017) was developed from the IGE model (Pimienta and Duval, 1987; Barnola et al., 1991; Arnaud et al., 2000; Goujon et al., 2003) by implementing a dependency of the firn densification rate on temperature and impurities. The temperature dependence is added to the classical formulation of the densification rate following an Arrhenius law with an activation energy Q as per $\exp (-Q/RT)$, with R the perfect gas constant and T the firn temperature. Rather than using a constant activation energy (Goujon et al., 2003), Bréant et al. (2017) stated that the value of the activation energy should be contingent on the firn temperature value as observed in material science, where the temperature dependency exhibits the predominance of one physical mechanism among others for a material compaction at a specific temperature. Through several sensitivity tests, Bréant et al. (2017) adjusted three values for activation energy in three different temperature ranges to reproduce best the δ¹⁵N evolution over the last deglaciation at East Antarctic sites. The firn model also considers that firn densification is facilitated by the dissolution of impurities within the snow (Freitag et al., 2013). If the impurity content in snow (i.e., concentration of calcium ions) is greater than a certain threshold, the densification rate dependence on impurities is traduced by a relationship between the new activation energy Q^′ and the concentration of calcium ions [Ca²⁺]: $Q^{\prime }=f_{\mathrm{1}}\times \left(\mathrm{1}-\mathit{\beta }\ln \left(\frac{\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]}{{\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]}_{\mathrm{threshold}}}\right)\right)\times Q$ (Freitag et al., 2013). Bréant et al. (2017) assumed the impurity effect equal for all physical mechanisms and tuned β and f₁ constants so that the modeled-δ¹⁵N data mismatch is minimized over the last glacial termination at cold East Antarctic sites.

As a consequence, and in addition to our new extensive δ¹⁵N dataset, we have chosen to use here the firn model approach of Bréant et al. (2017). In order to make a correct calculation of uncertainties linked to firn modeling at EDC, we ran two tests of the model with and without including the impurity concentration parameter (see Sect. S3.1).

The firn densification model takes as input scenarios of temperature and accumulation rate at the surface. It computes both the LID and the thermal gradient in the firn (ΔT) and then deduces the δ¹⁵N${}_{\mathrm{therm}}=\mathrm{\Omega }\cdot \mathrm{\Delta }T$ with Ω the thermal-fractionation coefficient (Grachev and Severinghaus, 2003). The final δ¹⁵N is calculated as δ¹⁵N = δ¹⁵N_therm + δ¹⁵N_grav and δ¹⁵N_grav ≃ LID ⋅ $\frac{g}{RT}$ (first-order approximation), with g the gravitational acceleration (9.8 m s⁻²), R the gas constant (8.314 J mol⁻¹ K⁻¹) and T the mean firn temperature (K).

Figure 1 Alignment of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and insolation between 800 and 100 ka BP. (a) Old raw EDC $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data between 800 and 100 ka BP (black circles for data of Extier et al., 2018c, and purple squares for data of Landais et al., 2012), outliers (gray crosses) and filtered signal (black and purple lines). New raw EDC $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data (blue triangles, this study) and filtered signals (blue line). The $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data are plotted on the AICC2012 ice timescale. Zooms between 270 and 100 ka BP and between 570 and 300 ka BP are shown in Fig. S2. (b) Extrema in the compiled filtered $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ dataset (plain blue line) are identified and matched to extrema in the (c) 21 December insolation at 75^∘ S plotted on a reversed y axis and on the age scale given by Laskar et al. (2004) (black line). The peaks are matched by vertical black bars. (d) The zero value in the time derivative of insolation (black line) and of the filtered $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ dataset (blue line) corresponds to extreme values in the signals. The determined tie points between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and insolation are depicted by markers on the horizontal line. Green circles are attached to a 3 kyr 1σ uncertainty (horizontal green error bars show 2σ in panel c), purple squares are associated with a 6 kyr 1σ uncertainty (horizontal purple error bars show 2σ in panel c), and red markers are associated with a 10 kyr 1σ uncertainty (horizontal red error bars show 2σ in panel c). Between 390 and 475 ka BP, all extrema are not tuned to the target due to the poor resemblance between the signal and insolation.

3 Age constraints and background scenarios

3.1 ⁸¹Kr age constraints

Three ice samples from the bottom part of EDC have been analyzed and provide three age estimates displayed in Table 3: 629, 788 and 887 ka BP, with statistical age uncertainties between 30 and 50 kyr and a 4.8 % systematic error due to the uncertainty in the half-life of ⁸¹Kr. The deepest sample suggests the presence of ice older than 800 ka BP below the 3200 m depth level, and further studies would be valuable in exploring whether the stratigraphy of the lowermost EDC section is continuous (Tison et al., 2015), although this is beyond the scope of this work.

Table 3 Ice sample details and radio krypton dating results. Reported errors are 1σ errors. Upper limits have a 90 % confidence level. The average ⁸⁵Kr activity in the Northern Hemisphere is about 75 dpm cm⁻³ (decays per minute per cubic centimeter) in 2017. The measured ⁸⁵Kr concentrations are inferior to the detection limit, verifying that no relevant contamination with modern air has occurred. In addition to the statistical error in the ⁸¹Kr age from atom counting, a systematic error due to the uncertainty in the half-life of ⁸¹Kr is considered. This error would shift the calculated ⁸¹Kr ages up or down for all ice samples. ^a dpm cm⁻³: decays per minute per cubic centimeter standard temperature and pressure (STP) of krypton (conversion: 100 dpm cm⁻³ corresponds to ⁸⁵Kr/Kr = 3.03×10¹¹). ^b pMKr: percent modern krypton.

3.2 Determination of orbital age constraints using new data

3.2.1  $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$

In this work, new highly resolved $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data on the EDC ice core are presented over terminations II, III, IV, V and VI (Fig. 1). As these novel $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements have been performed on ice samples stored at −50 ^∘C, there is little storage effect, and they can directly be merged with the 800 kyr long record of Extier et al. (2018c) (Table 1). The new dataset improves the resolution of the long EDC record, reaching sub-millennial-scale accuracy over marine isotope stages (MISs) 5, 7 and 9 and in particular over MIS 11 and MIS 13, periods of sparsity in the ancient record (Extier et al., 2018c). Although the two datasets agree well over recent periods (last 350 kyr), they show some discrepancies during older periods (between 550 and 375 ka BP; see Fig. 1). Such dissimilarities are observed over MIS 11 (between 424 and 374 ka BP), where the sampling resolution of the previous dataset is particularly low (2500 years). In addition, MIS 11 is a period characterized by a low eccentricity in the Earth orbit, inducing subdued variations in insolation, causing $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ changes of smaller magnitude and leading to a lower signal-to-noise ratio. Data by Landais et al. (2012) (shown by purple squares in Figs. 1 and S2) are consistent with the highly resolved data presented here, supporting the relevance of the new dataset over this period. Over termination VI (from 550 to 510 ka BP), the old dataset continuously increases, while the novel dataset shows a brief maximum at around 525 ka BP followed by a minimum at around 520 ka BP. These newly revealed variations seem in phase with insolation variations, suggesting that the new dataset shows improved agreement with insolation. Still, highly resolved measurements are needed in the lowermost part of the ice core, where noise is significantly altering the temporal signal.

Following a data processing treatment consistent with the method described in Kawamura et al. (2007), the compiled dataset is linearly interpolated every 100 years and then smoothed using a finite-duration impulse response (FIR) filter with a KaiserBessel20 window (cut-off from 16.7 to 10.0 kyr period, 559 coefficients for the 800 kyr long record) designed with the software Igor Pro in order to reject periods inferior to 10 000 years and erase the noise present in the data. Note that using a low-pass (rejecting periods below 15 kyr) or a band-pass filter (keeping periods between 100 and 15 kyr periods, used by Bazin et al., 2013) does not alter the peak positions in the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ curve (see Fig. S2). The noise is particularly significant for highly resolved $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data, and without preliminary filtering, it becomes ambiguous to identify the exact peak position (which needs to be subjectively placed on a 1000- to 2000-year interval; see Sect. S2.1).

The filter is then applied to the local summer solstice insolation curve to check that it does not induce the shift in extrema positions by more than 100 years. This condition is verified over the last 800 kyr, except for the peaks located at the endpoints of the record (respectively around 107 and 788 ka BP), which are then not used for tie point determination. Outliers in the raw $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ dataset are discarded if they show an anomaly greater than 3.2 ‰ when compared to the low-pass-filtered signal. Five outliers are rejected out of 294 points. The $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is interpolated and filtered again after removal of the outliers.

The orbital target chosen is the 21 December insolation at 75^∘ S, which is calculated every 100 years over the last 800 kyr (Laskar et al., 2004). The peak positions in the compiled filtered $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ signal and in the summer solstice insolation are detected via an automated method using the zero values of the time derivatives of the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and its orbital target. Each $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ maximum is matched to an insolation minimum and each $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ minimum to an insolation maximum. The data treatment and tie point identification method used here are consistent with the approach recently conducted by Oyabu et al. (2022) on a novel 207 kyr long $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record of the Dome Fuji (DF) ice core.

Some periods, such as between 450 and 350 ka BP (encompassing MIS 11) and older ages (before 600 ka BP), are characterized by a poor resemblance between the signal and the target. For instance, two or three peaks in the insolation curve only correspond, respectively, to one or two peaks in the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data. This could be explained by a low eccentricity-induced subdued variability in the insolation target and hence in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ signal over MIS 11 and by the poor resolution of the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements before 600 ka BP. In such cases, the uncertainty (1σ) associated with each tie point ranges from 6 to 10 kyr (precession half period), and some extrema in the target are not used to tune the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record (5 extrema over MIS 11 out of 63 over the last 800 kyr). Otherwise, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ seems to evolve in phase with the inverse summer solstice insolation variations, and the tie point uncertainty (1σ) is set at 3 kyr. A 3–4 kyr uncertainty was evaluated by Bazin et al. (2016) on the following arguments. They examined three $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records from the Vostok, DF and EDC ice cores over MIS 5 and detected some site-specific $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ high-frequency variability that could not be explained by a timescale issue. This observation, along with the presence of a 100 kyr periodicity in the EDC $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record and the difficulty of identifying $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ mid-slopes and maxima because of a scattering of the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ signal at the millennial scale, led them to recommend the use of a 3–4 kyr uncertainty. Because our higher-resolution $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data give the possibility of filtering the signal with more confidence and hence reduces the uncertainty in the identification of $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ tie points, we propose taking a 3 kyr uncertainty. The orbital tuning results in 58 new tie points over the last 800 kyr (displayed in Fig. 1 and compiled in Table S5), replacing the 20 tie points used to constrain AICC2012 between 800 and 300 ka BP that were derived from synchronizing mid-slopes of band-pass-filtered $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ with the insolation (Bazin et al., 2013).

Figure 2 Alignment of TAC and insolation between 800 and 0 ka BP. (a) Raw EDC TAC data (blue circles; Raynaud et al., 2007), outliers (gray crosses) and filtered signal (blue line) on the AICC2012 ice timescale. The raw data between 800 and 440 ka BP are not shown here and will be published in a separate study. (b) ISI375 at 75^∘ S on a reversed axis. The peaks and mid-slopes are matched by vertical bars. (c) Temporal derivative of insolation (black line) and TAC (blue line). Its zero value corresponds to extreme values in insolation and TAC. The determined tie points between TAC and insolation are depicted by markers on the horizontal line. Green circles are attached to a 3 kyr 1σ uncertainty (horizontal green error bars show 2σ in panel c), purple squares are associated with a 6 kyr 1σ uncertainty (horizontal purple error bars show 2σ in panel c), and red markers are associated with a 10 kyr 1σ uncertainty (horizontal red error bars show 2σ in panel c).

The uncertainty arising from the filter used and from the tie point identification method can be estimated by a comparison of the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ peak positions identified before and after filtering of the signal with two different methods (see Sect. S2.1). The resulting uncertainty is of 700 years on average (with a standard deviation of 250 years), reaching 2100 years around 230 ka BP.

The new highly resolved data presented here enable a better description of the signal variability and a reduction in the uncertainty associated with orbital tie points.

3.2.2 Total air content

The TAC record is extended over the last 800 kyr with a mean sampling resolution of 2000 years (Fig. 2). The raw data between 800 and 440 ka BP are not shown here and will be published in a separate study. The TAC series shows a good resemblance with the integrated summer insolation (ISI, obtained by a summation over a year of all daily insolation at 75^∘ S above a chosen threshold). After comparing the EDC TAC record, within its frequency domain, with ISI curves obtained using different thresholds, the ISI curve calculated for a threshold of 375 W m⁻² (ISI375) exhibits the finest spectral agreement with the EDC TAC record over the past 800 kyr. The coherency between the TAC record and ISI is deficient over MIS 11 (between 430 and 370 ka BP) and in the deepest part of the core (prior to 700 ka BP), where the signal-to-noise ratio is low.

Following a data processing treatment consistent with the method described by Lipenkov et al. (2011), the 800 kyr long TAC dataset is interpolated every 100 years and then filtered with a band-pass filter rejecting periods below 15 000 and above 46 000 years (Igor Pro FIR filter with a KaiserBessel20 window, cut-off from 15 to 14 kyr and from 46 to 47 kyr, 559 coefficients). Outliers in the raw TAC dataset are discarded if they show an anomaly greater than 1.0 mL kg⁻¹ (standard deviation of TAC record) when compared to the band-pass-filtered signal. A total of 45 outliers are rejected out of 399 data points (among which 16 outliers are identified between 100 and 0 ka BP). The TAC is interpolated and filtered again after removal of the outliers.

Tie points are mostly determined by matching variation extrema of TAC and integrated summer insolation at 75^∘ S (Fig. 2). Indeed, in the case of a non-linear relationship between TAC and insolation, extrema are better indicators of TAC response to insolation forcing. Moreover, filtering the dataset induces a bias in the mid-slope position. The method employed to determine extrema position is the same as for $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ insolation tie points. Only one of the tie points is identified by matching mid-slopes (i.e., derivative extremum) at 362 ka BP rather than minima at 375 ka BP due to the flatness of the insolation minimum, which precludes identification of an accurate tie point. Not all extrema are tuned to the target due to the poor resemblance between the signal and insolation, and 42 unambiguous tie points were kept out of 64 detected by the automated method. The tie point uncertainty finds its origin in the age errors associated with the filtering (∼700 years), tie point identification and outlier rejection (∼900 years). The 1σ uncertainty is evaluated to be 3 kyr when there is good agreement: (i) between the signal and its target, meaning that one peak in ISI375 is reflected by a singular peak in the TAC record, and (ii) between the tie points identified by the automated method and manually (age shift < 1300 years, average value) (see green circles, Fig. 2). A 6 kyr uncertainty (1σ) is attached to the tie points if the latter condition is not respected (age shift > 1300 years) (see purple squares, Fig. 2), and a 10 kyr uncertainty (1σ) (precession half period) is ascribed to the tie points if the ISI375 variations are not reflected by the TAC record, meaning that one peak in ISI375 could be associated with two peaks in the TAC record, or if the signal-to-noise ratio of the TAC record is too large (see red markers, Fig. 2). The choices of filter and orbital target have no significant impact on the chronological uncertainty; a further detailed study is thus beyond the scope of this work.

The orbital tuning results in 42 new tie points over the last 800 kyr (displayed in Fig. 2 and compiled in Table S5). They replace the 14 tie points used to constrain the EDC ice timescale in AICC2012 between 425 and 0 ka BP that were derived by directly matching mid-slope variations in unfiltered TAC and ISI target and attached to an uncertainty varying between 2.9 and 7.2 kyr.

3.2.3 δ¹⁸O_atm

In this work, new highly resolved δ¹⁸O_atm data on the EDC ice core are presented over terminations II, III, IV, V and VI (Fig. 3). The available δ¹⁸O_atm data can be sorted out in two groups: new δ¹⁸O_atm data (Grisart, 2023) at high temporal resolution (between 333 and 1375 years; see Table 1) and old measurements compiled by Extier et al. (2018b), characterized by a lower sampling resolution (between 1000 and 1500 years; see Table 1), except between 374 and 153 ka BP (resolution between 160 and 700 years; see Table 1). The new dataset improves the resolution of the long EDC record over MISs 5, 7 and 9 and in particular over MISs 11 and 13, periods of sparsity in the ancient record (Extier et al., 2018b). Although the two datasets agree well globally over the last 800 kyr, the new highly resolved dataset refines the signal between 255.5 and 243 ka BP where a lot of noise is present in the record of Extier et al. (2018b) (see inset in Fig. 3). This noise may be explained by the fact that highly resolved (mean sampling resolution of 381 years) measurements were performed on ice samples stored at −20 ^∘C in the compilation by Extier et al. (2018b), while the new measurements are performed on ice stored at −50 ^∘C. Therefore, we chose to remove the noisy dataset of Extier et al. (2018b) between 255.5 and 243 ka BP before combining the novel dataset with the remaining 800 kyr long record of Extier et al. (2018b).

Figure 3 Evolution of EDC δ¹⁸O_atm record between 800 and 100 ka BP. (a) Old raw EDC δ¹⁸O_atm data (black circles; Extier et al., 2018b) and new raw EDC δ¹⁸O_atm data (blue triangles; Grisart, 2023) on the AICC2012 gas timescale. (b) Compilation of the two datasets after removal of old measurements between 255.5 and 243 ka BP. (c) δ¹⁸O_calcite composite record from speleothems from Sambao, Dongge, Hulu (red line) and Yongxing (brown line) caves (Cheng et al., 2016; Zhao et al., 2019) on U–Th age scales. (d) Climatic precession from Laskar et al. (2004) delayed by 5000 years. Inset is a zoom between 290 and 190 ka BP. Vertical gray rectangles highlight the improved agreement between new data of Grisart (2023) (blue triangles) and δ¹⁸O_calcite (red line) compared to old data (gray circles) and δ¹⁸O_calcite.

Following the dating approach proposed by Extier et al. (2018a), δ¹⁸O_atm and δ¹⁸O_calcite are aligned using mid-slopes of their variations over the last 640 kyr. To do so, the compiled EDC δ¹⁸O_atm record and the Chinese δ¹⁸O_calcite signal are linearly interpolated every 100 years and smoothed (25-point Savitzky–Golay) and extrema in their temporal derivative are aligned. It should be specified that synchronizing δ¹⁸O_atm and East Asian δ¹⁸O_calcite is not always obvious due to the long residence time of oxygen in the atmosphere (1–2 kyr), which may not be compatible with abrupt δ¹⁸O_calcite variations over glacial inceptions and terminations. In particular, the slow increase in the δ¹⁸O_atm record from 370 to 340 ka BP does not resemble the evolution of δ¹⁸O_calcite, which is first moderate then abrupt over the same period (Fig. 4, red area). For this reason, we chose not to use the two tie points identified by Extier et al. (2018a) at 351 and 370.6 ka BP. The new highly resolved data enable the identification of five new tie points and shifting of five tie points that have been determined beforehand by Extier et al. (2018a) (Fig. 4). Between 248 and 244 ka BP, the new δ¹⁸O_atm measurements do not coincide with the δ¹⁸O_calcite variations, and we decided to remove the tie point identified by Extier et al. (2018a) at 245.4 ka BP (Fig. 4, red area). Between 480 and 447 ka BP, the δ¹⁸O_atm variations are characterized by a low resolution (1.1 kyr) and a weak amplitude, which prevents unambiguous matching of δ¹⁸O_atm and δ¹⁸O_calcite. The four tie points identified by Extier et al. (2018a) at 447.3, 449.9, 455.9 and 462.8 ka BP are thus rejected (Fig. 4, red area). The remaining 39 tie points defined by Extier et al. (2018a) are preserved and used here to constrain the EDC gas age. Their uncertainty (1σ) varies between 1.1 and 7.4 kyr.

Figure 4 Alignment of EDC δ¹⁸O_atm and Chinese δ¹⁸O_calcite records over time periods where new tie points are defined. (a) New and old EDC δ¹⁸O_atm datasets on the AICC2012 gas age scale. (b) Compiled EDC δ¹⁸O_atm. (c) Chinese δ¹⁸O_calcite on U–Th age scale (Cheng et al., 2016). (d) Temporal derivatives of compiled EDC δ¹⁸O_atm (blue curve) and of the old δ¹⁸O_atm dataset (black curve). (e) Temporal derivative of Chinese δ¹⁸O_calcite (red curve). Extrema in temporal derivatives are matched. Tie points represented by vertical black bars are determined by Extier et al. (2018a), and those by vertical blue bars are determined by this study. Both are used in the AICC2023 chronology. Dashed vertical bars show tie points identified by Extier et al. (2018a) that are not used in AICC2023; 2σ uncertainties attached to the tie points are shown by the horizontal error bars in panel c). Vertical red areas frame periods of lacking resemblance between δ¹⁸O_atm and δ¹⁸O_calcite variations.

Figure 5 Alignment of EDC δ¹⁸O_atm and climatic precession between 810 and 590 ka BP. (a) Compiled EDC δ¹⁸O_atm on the AICC2012 gas timescale. (b) Precession delayed by 5000 years (dashed gray line) and not delayed (black line) (Laskar et al., 2004). (c) Temporal derivative of precession (black line), delayed precession (dotted gray line) and the compiled δ¹⁸O_atm record (plain purple line). (d) Ice-rafted debris at ODP983 site (North Atlantic Ocean, southwest of Iceland) by Barker et al. (2019, 2021). The gray rectangles indicate periods where IRD counts are greater than the 10 counts per gram threshold shown by the dotted horizontal blue line. Vertical gray bars illustrate new tie points between EDC δ¹⁸O_atm and delayed precession mid-slopes (i.e., derivative extrema) when IRD counts are greater than the threshold. Vertical black bars illustrate new tie points between EDC δ¹⁸O_atm and precession mid-slopes (i.e., derivative extrema) when no Heinrich-like events are shown by the IRD record. The 12 kyr 2σ uncertainty attached to the tie points is shown by the horizontal error bars in panel (b).

Between 810 and 590 ka BP, the δ¹⁸O_atm–δ¹⁸O_calcite dating uncertainty becomes larger than 6 kyr, and no East Asian speleothem δ¹⁸O_calcite records are available before 640 ka BP. Over this time interval, we updated the following approach of Bazin et al. (2013): EDC δ¹⁸O_atm and 5 kyr delayed climatic precession are synchronized using mid-slopes of their variations. However, from the findings of Extier et al. (2018a), δ¹⁸O_atm should rather be aligned to precession without delay when no Heinrich-like events occur. Indeed, δ¹⁸O_atm is sensitive to both orbital- and millennial-scale variations in the low-latitude water cycle (Landais et al., 2010; Capron et al., 2012), and Heinrich-like events occurring during deglaciations delay the response of δ¹⁸O_atm to orbital forcing through southward ITCZ shifts (Extier et al., 2018a). We thus chose to align δ¹⁸O_atm to precession when no ice-rafted debris (IRD) peak is visible in the studied period in the ODP983 record (Barker et al., 2019, 2021) and keep a 5 a delay when IRD peaks are identified. This results in shifting 12 tie points of Bazin et al. (2013) by 5000 years towards older ages (Fig. 6). The eight remaining tie points of Bazin et al. (2013) that coincide with peaks in the IRD record are kept (Fig. 6). To confirm the validity of our approach, we tested three methodologies to align δ¹⁸O_atm and precession over well-dated periods when δ¹⁸O_atm–δ¹⁸O_calcite matching was done (see Sect. S2.2.2). These tests support our approach, but in order to account for potential errors associated with this tuning method (Oyabu et al., 2022), a 6 kyr uncertainty (1σ) is attributed to the δ¹⁸O_atm-derived tie points over the period between 810 and 590 ka BP.

A total of 69 new δ¹⁸O_atm tie points are determined over the last 810 kyr (displayed in Figs. 4 and 5 and compiled in Table S5). They replace the 39 tie points used to constrain the EDC gas timescale in AICC2012 between 800 and 363 ka BP (Bazin et al., 2013). The age constraints are attached to an uncertainty varying between 1.1 and 7.4 kyr, which is the sum of the uncertainties in the speleothem ²³⁰Th dating, the δ¹⁸O_atm response to orbital forcing (1 kyr) and the δ¹⁸O_atm–δ¹⁸O_calcite matching (0.5 kyr). The same alignment method is applied between Vostok δ¹⁸O_atm (Petit et al., 1999) and Chinese δ¹⁸O_calcite, and 36 new tie points are determined (see Sect. S4.1.2), replacing the 35 tie points used to constrain the Vostok gas timescale in AICC2012.

Figure 6 EDC δ¹⁵N record and past LID evolution as a function of EDC depth. (a) New and highly resolved δ¹⁵N dataset (blue circles), dataset of Bréant et al. (2019) (purple circles), old dataset (black circles) and outliers (rejection criterion of 1σ) (gray crosses). (b) LID calculated as per LID ≃ δ¹⁵N${}_{\mathrm{grav}}\cdot \frac{RT}{g}$ for three cases: (1) δ¹⁵N_grav=δ¹⁵N with the δ¹⁵N record constructed by interpolation between data when no data are available (gray), (2) δ¹⁵N_grav=δ¹⁵N with the δ¹⁵N record constructed by normalization of the δD record when no data are available (black ), (3) δ¹⁵N${}_{\mathrm{grav}}={\mathit{\delta }}^{\mathrm{15}}\mathrm{N}-{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}_{\mathrm{therm}}$ with δ¹⁵N_therm estimated by the firn model (Bréant et al., 2017) and the δ¹⁵N record constructed by interpolation between data when no data are available (blue). (c) Modeled LID with impurity concentration (blue) and without impurity concentration (red). (d) Background scenarios of LID used to construct AICC2012 (black) and inputs in Paleochrono to obtain AICC2023 (blue). (e) Absolute difference between prior LID of AICC2012 and AICC2023. The gray line separates the top 5 % values from the rest. The gray rectangles cover areas when no δ¹⁵N data are available.

Finally, there was a redundancy in the dating of the bottom part of the EDC ice core in AICC2012, where both δ¹⁸O_atm orbital tie points and ¹⁰Be peaks corresponding to the Matuyama–Brunhes geomagnetic reversal event were used. Indeed, the two ¹⁰Be dating constraints at 780.3 and 798.3 ka BP were directly derived from the δ¹⁸O_atm orbital dating and not obtained independently (Dreyfus et al., 2008). We thus decided to remove the ¹⁰Be age constraints.

3.3 Background scenario of LID

In this work, new highly resolved δ¹⁵N data in the EDC ice core are presented over terminations II, III, IV, V and VI (Fig. 6a). The available δ¹⁵N data can be sorted out in two groups: δ¹⁵N measured by Grisart (2023) and Bréant et al. (2019) at high temporal resolution (between 333 and 1375 years; see Table 1) and the older measurements (Bazin et al., 2013) used to estimate LID in AICC2012, characterized by a lower sampling resolution (between 1400 and 2400 years; see Table 1). The measurements of Bazin et al. (2013) and Bréant et al. (2019) have been shifted down by 0.04 ‰ to account for calibration errors. The new dataset permits the extension of the record by around 1100 m and between 1700 and 2500 m and improvement in the resolution over terminations II to VI.

Outliers are discarded if they show an anomaly greater than 0.045 ‰ when compared to the smoothed record (Savitzky–Golay algorithm with 25 points). This results in the rejection of 25 data points out of 475 measurements for the new dataset (Fig. 6). The two δ¹⁵N datasets are merged, and the compiled record is interpolated every 100 years. Then, assuming that the firn is solely a diffusive zone (i.e., no convection layer at the top) at EDC during the last 800 kyr, in agreement with current observations (Landais et al., 2006), past LID is calculated as per the first-order estimate of the barometric equation:

$$\begin{array}{}\text{(5)}& \mathrm{LID}\simeq {{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}_{\mathrm{grav}}{\displaystyle \frac{RT}{g}},\end{array}$$

with T the temperature at EDC estimated from combined measurements of ice δ¹⁸O and δD after correction of the influence of the seawater δ¹⁸O (Landais et al., 2021).

Table 4 Method of determination of the background LID scenario according to the EDC depth range.

In the absence of a large thermal gradient within the firn (mostly present in Greenlandic ice cores during Dansgaard–Oeschger events), δ¹⁵N is mainly modulated by gravitational fractionation of N₂ molecules occurring from the surface down to the lock-in zone, and δ¹⁵N measured in bubbles hence approximately reflects the LID (Severinghaus et al., 1996; Landais et al., 2006) and δ¹⁵N_grav≃δ¹⁵N in Eq. (5) (gray and black lines in Fig. 6b). To account for a small temperature gradient in the firn in the Antarctic ice core, the thermal-fractionation term δ¹⁵N_therm can be estimated by the firn model (Bréant et al., 2017). Past LID is then calculated as per Eq. (5) with δ¹⁵N_grav=δ¹⁵N − δ¹⁵N_therm (blue curve in Fig. 6b). Thermal fractionation represents a maximum correction of 4.2 m to the LID at EDC.

When δ¹⁵N measurements are not available, Bazin et al. (2013) used a synthetic δ¹⁵N signal based on the correlation between δ¹⁵N and δD to estimate the LID background scenario at Dome C (black curve in Fig. 6b). Indeed, for different low-accumulation Antarctic sites, it has been observed that δ¹⁵N and δD are well correlated over the last termination on a coherent timescale (Dreyfus et al., 2010; Capron et al., 2013). Since then, Bréant et al. (2019) presented new high-resolution EDC measurements of δ¹⁵N extending the signal over termination III (around 2300 m, 250 ka BP). Their study unveiled the anatomy of this atypical deglaciation: the interplay between Heinrich-like events and the bipolar seesaw mechanism induced a strong warming of Antarctic temperature, resulting in divergent δ¹⁵N and δD records. Therefore, using δD to construct a synthetic δ¹⁵N scenario should be done carefully. For this reason, the firn densification model described in Bréant et al. (2017) is employed to estimate EDC LID evolution in the past when δ¹⁵N data are missing, rather than using the δ¹⁵N–δD relationship, as was done for AICC2012. After different sensitivity tests, we chose to keep the parameterization preferred by Bréant et al. (2017) (i.e., firn densification activation energy depending on the temperature and the impurity concentration) as it is believed to give the most probable evolution of LID over the last 800 kyr (see Sect. S3.1).

The final background LID scenario is calculated as a function of EDC depth (Table 4, Fig. 6d). It has been smoothed using a Savitzky–Golay algorithm (25 points) and then provided as an input file to Paleochrono.

The other necessary input files for Paleochrono, accumulation (A) and thinning (τ) background scenarios, are the same as in Bazin et al. (2013). A is estimated from water isotopes (Parrenin et al., 2007b) and τ from unidimensional ice flow modeling (Parrenin et al., 2007a).

Figure 7 EDC ice age difference between AICC2012 and different test chronologies obtained with Paleochrono over the last 800 kyr. The ice age difference is calculated as per test chronology − AICC2012. Two test chronologies are obtained by addition of either new Vostok δ¹⁸O_atm–δ¹⁸O_calcite age constraints (dotted green line) or stratigraphic and absolute TALDICE constraints between 470 and 129 ka BP from Crotti et al. (2021) (dotted red line). The other “test chronologies” are constructed by replacing AICC2012 constraints by either (1) ⁸¹Kr (green); (2) $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ (purple); (3) TAC (orange); (4) δ¹⁸O_atm (red); (5) CH₄ tie points with NGRIP, EDML, TALDICE and Vostok (black); or (6) volcanic matching points with EDML and NGRIP (gray). Vertical bars represent the corresponding age horizons. AICC2023 is obtained by implementing the new constraints all together (light-blue line). Vertical light-blue bars show new data collected by Grisart (2023) and presented in this work. The three largest inconsistencies between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, TAC and δ¹⁸O_atm chronologies are shown by red areas. Gray rectangles indicate interglacials from MIS 19 to MIS 1.

3.4 New stratigraphic links between the EDC and other ice cores

EDC can be linked to other ice cores via ice and gas stratigraphic links identified during abrupt climate changes recorded in Greenlandic and Antarctic ice cores. To establish AICC2012, Bazin et al. (2013) used 255 gas stratigraphic tie points coming from the matching of CH₄ (or δ¹⁵N when CH₄ is not available at NGRIP) or δ¹⁸O_atm variations between EDC, EDML, Vostok, NGRIP and TALDICE. Here we revise some of these tie points using the synchronization of CH₄ series of EDC, Vostok and TALDICE to up-to-date, highly resolved records from EDML and NGRIP ice cores over the last interglacial offset and the last glacial period (Baumgartner et al., 2014). From 122 to 10 ka BP, Baumgartner et al. (2014) identified 39 stratigraphic links between EDML and NGRIP by matching mid-points of the CH₄ abrupt changes with a precision of 300 to 700 years. When they also detected such rapid variations in lower-resolution CH₄ records of the TALDICE, Vostok and EDC ice cores, they extended the stratigraphic links to the five ice cores but assigned them a larger uncertainty (up to 1500 years). AICC2012 was further constrained by 534 ice stratigraphic links identified from volcanic matching and synchronization of cosmogenic isotopes between the five ice cores. Here we replace some of the stratigraphic links between NGRIP, EDML and EDC by highly resolved volcanic matching points (Svensson et al., 2020). The application of volcanic proxies and annual layer counting helped them identify large volcanic eruptions that left a specific signature in both Greenland and Antarctica. Such a signature is defined by sulfate patterns (indicating singular volcanic events separated by the same time interval in ice cores from both poles). Their study spotted 82 large bipolar volcanic eruptions over the second half of the last glacial period (from 60 to 12 ka BP), providing as many ice stratigraphic links synchronizing EDC with EDML and EDML with NGRIP within a small relative uncertainty (i.e., ranging from 1 to 50 years, 12 years on average). Between 43 and 40 ka BP, five cosmogenic tie points associated with the Laschamp geomagnetic excursion (Raisbeck et al., 2017) replace the volcanic matching over this period (Svensson et al., 2013), shifting the tie points by ∼30 years.

4 Discussion

4.1 New AICC2023 chronology

4.1.1 Impact of absolute age constraints

A large uncertainty is linked with ⁸¹Kr dating; therefore ⁸¹Kr age estimates do not significantly change the chronology (maximum 200 years) (Fig. 7). ⁸¹Kr age estimates are systematically older than the new timescale (by 25 to 36 kyr; see Fig. 8). This observation could also indicate an undervaluation of ⁸¹Kr half-life.

Figure 8 EDC ice age and uncertainty as a function of the depth. (a) EDC ice age (AICC2012 in black, AICC2023 in blue), (b) 1σ uncertainty (AICC2012 in black, AICC2023 in blue). Crosses and slashes represent new age constraints (ice stratigraphic links in black, gas stratigraphic links in gray, δ¹⁸O_atm in red, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ in blue, TAC in orange, ⁸¹Kr in green). Inset is a zoom in between 800 and 600 ka BP. Gray rectangles frame periods where the new AICC2023 uncertainty is larger than the AICC2012 uncertainty. See Fig. S13 for EDC gas age profile.

4.1.2 Consistency between orbital age constraints

To evaluate the consistency between the orbital age constraints, several “test chronologies” are produced. Each test chronology of the EDC ice core is obtained by running one multi-site (EDC, Vostok, EDML, TALDICE, NGRIP) experiment of Paleochrono. In each of these tests, we implemented one category of new age constraints presented in this work while keeping AICC2012 parameters for other categories. Several test chronologies are thus constructed: the ⁸¹Kr-, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$-, TAC-, δ¹⁸O_atm-, CH₄-matching- and volcanic-matching-based chronologies (Fig. 7). Two additional test chronologies are obtained by implementing and modifying age constraints either on Vostok or TALDICE with respect to the AICC2012 chronology, as explained in Sect. 2.1 (Fig. 7, dotted lines). The EDC ice age difference between each test chronology and the AICC2012 timescale is represented in Fig. 7 so that it is possible to read which type of dating tool suggests shifting the background chronology towards either older or younger ages.

Although the three orbital dating tools globally agree with each other over the last 800 kyr, meaning that they all tend to shift the background chronology towards either older or younger ages over a certain period of time, they are sometimes inconsistent (Fig. 7). The three largest inconsistencies involve age differences between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$-, TAC- and δ¹⁸O_atm-based chronologies reaching 4.15 to 8.3 kyr (Table 5). At 390 ka BP, a large 8.3 kyr discrepancy is observed between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$- and δ¹⁸O_atm-based chronologies. Over this period, the low-resolution $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record variations do not match its orbital target variations (two insolation minima against one $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ maximum; see Fig. 1). For this reason, the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ age constraints identified between 480 and 350 ka BP were attached to a 6 kyr uncertainty (quarter of a recession period; Fig. 1). In contrast, the δ¹⁸O_atm record agrees well with δ¹⁸O_calcite (Fig. 3), and the uncertainty attached to the δ¹⁸O_calcite inferred tie points over this interval is smaller. Hence, the new AICC2023 chronology suggests shifting AICC2012 towards older ages by 2.2 kyr, as per the δ¹⁸O_atm-based chronology (Fig. 7). Around 550 ka BP, the TAC- and δ¹⁸O_atm-based chronologies strongly diverge. This may be caused by the absence of TAC tie points due to the non-coincidence of TAC and ISI375 extrema (Fig. 3), while there is good agreement between δ¹⁸O_atm and δ¹⁸O_calcite records. Therefore, we decide to increase the uncertainty attached to the four TAC age constraints between 600 and 550 ka BP up to 6 kyr (Fig. 2), and AICC2023 follows the δ¹⁸O_atm-based chronology, inducing older ages than AICC2012. At 765 ka BP, the discordance between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$- (and TAC-) and δ¹⁸O_atm-based chronologies is likely due to the poor quality of the records from the lowermost part of the core. Over these oldest time periods, δ¹⁸O_atm, TAC and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ were tied up, respectively, with precession, integrated insolation and insolation with a large uncertainty (6 to 10 kyr). This leads to a final chronology AICC2023, suggesting a larger chronological uncertainty than AICC2012 as well as younger ages (as per TAC and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ chronologies) over MIS 18 and then older ages (as per δ¹⁸O_atm chronology) over MIS 19.

Table 5 Description of the inconsistencies between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$-, TAC- and δ¹⁸O_atm-based chronologies. The age shift suggested by each dating tool with respect to AICC2012 age is detailed. The age position of the disagreement is given as per AICC2023. We did not highlight inconsistencies between $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$- and TAC-based chronologies as they remain within their respective orbital uncertainty.

4.1.3 Final chronology and uncertainty

The new AICC2023 chronology suggests significant age shifts when compared to AICC2012 over old periods, including 3.8 and 5 kyr shifts towards older ages around 800 and 690 ka BP as well as a 2.1 kyr shift towards younger ages around 730 ka BP. The chronology is also modified over MIS 5 and MIS 11, where AICC2023 is about 2 kyr older than AICC2012. These 2 kyr shifts are induced by $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and δ¹⁸O_atm dating constraints and stratigraphic links over MIS 5 and by TAC and δ¹⁸O_atm constraints over MIS 11. When averaged over the past 800 kyr, the chronological uncertainty is reduced from 1.7 kyr for AICC2012 to 900 years here. Still, it remains significant (above 2 kyr) over MIS 11 and in the lowermost part of the core, between 800 and 650 ka BP. Specifically, between 800 and 670 ka BP, the uncertainty associated with the new AICC2023 timescale is sometimes larger than the AICC2012 uncertainty (Fig. 8). This is caused by a larger relative error attached to $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC age constraints as well as by the exclusion of the two redundant ¹⁰Be age constraints at 780.3 and 798.3 ka BP associated with the Matuyama–Brunhes geomagnetic reversal event.

The age difference between ice and gas timescales (Δage) is 3 kyr on average, reaching its largest values (∼4 kyr) during the cold eras of MISs 12, 8, 6 and 4 (at 440, 260, 145 and 70 ka BP, respectively; Fig. 9). A 4 kyr Δage is obtained at around 160 ka BP (Fig. 9), consistent with the use of new δ¹⁵N data of Bréant et al. (2019), leading to a background scenario of LID that is 13 m smaller than the prior LID scenario used in AICC2012 between the depths of 1900 and 2000 m (Fig. 6). Using the definition of an interglacial period implying an EDC δD value surpassing the threshold of −403 ‰ (EPICA Community Members, 2004), we identify 10 substages of interglacials (MISs 1, 5e, 7a, 7c, 7e, 9e, 11c, 15a, 15e and 19; Fig. 9). The average duration of these substages is reduced by 320 years with the new AICC2023 timescale in comparison with the AICC2012 chronology (Fig. 9). More specifically, MISs 5e to 15a are shorter, while only MIS 15e and MIS 19 are longer. The largest decreases in duration affect the Last Interglacial (MIS 5e) and MIS 11c, whose lengths are decreased from 16.3 to 15.1 kyr and from 31.1 to 30.1 kyr, respectively, in agreement with the durations of 14.8 and 29.7 kyr proposed by Extier et al. (2018a).

Figure 9 EDC gas and ice records on AICC2023 (blue) and AICC2012 (black) timescales over the last 800 kyr. (a) EDC CH₄ (Loulergue et al., 2008) on AICC2012 and (b) AICC2023 gas timescales. (c) Gas age difference AICC2023–AICC2012. Gray and blue envelopes are AICC2012 and AICC2023 chronological 1σ uncertainties, respectively. (d) EDC δD (Jouzel et al., 2007) on AICC2012 and (e) AICC2023 ice timescales. Gray and blue rectangles indicate interglacial periods defined when δD is greater than the threshold of −403 ‰ (horizontal lines) (EPICA Community Members, 2004). Interglacials are numbered from MISs 1 to 19 (Berger et al., 2016). (f) Ice age difference AICC2023–AICC2012. (g) Age difference between ice and gas AICC2023 timescales (Δage).

4.2 Comparison with other chronologies

4.2.1 MIS 5 (from 130 to 80 ka BP)

When Veres et al. (2013) presented the AICC2012 chronology over the last climatic cycle, they identified a disagreement with the Greenland timescale GICC05-modelext between 115 and 100 ka BP. The comparison between the Greenland δ¹⁸O_ice record and the δ¹⁸O_calcite from U–Th-dated Alpine speleothems showed a delay of up to 2.7 kyr during the Dansgaard–Oeschger (D–O) events 23, 24 and 25. Later, this disagreement between abrupt changes in δ¹⁸O_ice from NGRIP (Greenland surface temperature) and δ¹⁸O_calcite from the Alps was re-evaluated based on a different use of δ¹⁸O_atm in ice core chronology, and Extier et al. (2018a) presented better agreement between the two records with an older NGRIP timescale than AICC2012 by ∼2200 years for D–O 23 to 25.

In Fig. 11, the NGRIP δ¹⁸O_ice record is represented on the AICC2023 timescale and is compared to ancient and novel records of δ¹⁸O_calcite from speleothems that grew along the northern rim of the Alps (NALPS; Boch et al., 2011; Moseley et al., 2020). Thanks to new $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and δ¹⁸O_atm age constraints, the new AICC2023 chronology is also older than AICC2012 between 115 and 100 ka BP and leads to improved agreement between the records along with a reduction in the uncertainty. This amelioration is particularly visible over D–O warmings 23 and 24, where the difference between NALPS and NGRIP chronologies is reduced from ∼2000 years (AICC2012) to 430 and 325 years (AICC2023), respectively (Table 6).

Table 6 Timing of D–O warmings 23 and 24 and D–O 25 rebound event onset. The GICC05-modelext age uncertainty is undetermined.

The Greenland Interstadial (GI) 25 can be subdivided in three substages: GI-25a, GI-25-b and GI-25-c, where GI-25a is the earliest glacial so-called “rebound event” (Capron et al., 2010). This latter consists of a brief warm–wet excursion during the slow cooling trend of the longer GI-25 period before jumping back to a cool–dry climate. The GI-25a warm–wet interval corresponds to a temperature increase in Greenland and continental Europe and is hence identified by a positive excursion in NGRIP and NALPS δ¹⁸O records (D–O 25 rebound) (Boch et al., 2011; Capron et al., 2012). At lower latitudes, this rebound likely affected the rainfall amount variations, as exhibited by the abrupt decrease in the δ¹⁸O_calcite from a U–Th-dated Sardinian stalagmite from Bue Marino Cave (BMS1; Columbu et al., 2017). The 2 kyr shift in the new AICC2023 chronology towards older ages improves the coherency between NALPS, NGRIP and BMS1 timescales over the GI-25a onset (traceable in the δ¹⁸O series; Fig. 10). The age discrepancy is reduced from ∼3600 years (between AICC2012 and BMS1 timescale) to 1640 years (between AICC2023 and BMS1 timescale; Table 6).

Figure 10 Northern Alpine speleothems (NALPS) and Bue Marino stalagmite (BMS1) δ¹⁸O_calcite records and NGRIP δ¹⁸O_ice evolution between 114 and 100 ka BP. NGRIP δ¹⁸O_ice data by Andersen et al. (2004) on AICC2012 (gray) and AICC2023 (blue) chronologies. NALPS δ¹⁸O_calcite data by Moseley et al. (2020) (red) and Boch et al. (2011) (brown). BMS1 δ¹⁸O_calcite data by Columbu et al. (2017) (dark blue). Vertical bars indicate D–O 25 rebound, D–O 24 and D–O 23 warmings at the onset of the GI-25a warm–wet substage, GI-24, and GI-23. They correspond to abrupt increases in the NALPS δ¹⁸O_calcite and NGRIP δ¹⁸O_ice records and to a decrease in the BMS1 δ¹⁸O_calcite series (for the GI-25a onset). Dashed black bars and blue bars show increases in δ¹⁸O_ice, respectively, on AICC2012 and AICC2023 chronologies. Brown bars and dotted red bars show increases in NALPS δ¹⁸O_calcite datasets. The dotted blue bar indicates the decrease in BMS1 δ¹⁸O_calcite. GI–GS (Greenland stadial) boundaries and GI-25 subdivision are indicated on the new AICC2023 chronology by horizontal bars.

Figure 11 Evolution of EDC and DF records on AICC2023 and DF2021 chronologies between 180 and 100 ka BP. (a) δD records from DF (red; Uemura et al., 2018) and EDC (purple; Jouzel et al., 2007). (b) $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ records from DF (red triangles; Oyabu et al., 2022) and EDC (purple circles; this work). (c) δ¹⁸O_atm records from DF (red triangles; Kawamura et al., 2007) and EDC (purple circles; this work). DF and EDC records are represented on DF2021 and AICC2023 timescales. (d) IRD from ODP 983 (Barker et al., 2019, 2021). (e) Ice age difference between DF2021 and (i) AICC2023 (blue), (ii) Extier et al. (2018a) chronology (orange) and (iii) AICC2012 (black). The age difference is calculated as per EDC age − DF2021 age. DF2021 age is transferred onto the EDC ice core via the volcanic synchronization of Fujita et al. (2015). The gray rectangle indicates MIS 5e.

Between 128 and 103 ka BP, the comparison between the AICC2012 timescale and the novel Dome Fuji ice core DF2021 chronology indicates that AICC2012 is likely too young by up to 4 kyr. Here, thanks to new highly resolved $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data and to the alignment of δ¹⁸O_atm and δ¹⁸O_calcite records, we improve the consistency between AICC2023 and DF2021, now agreeing within 1.7 kyr over MIS 5e (Fig. 11). With the new chronologies, the records of δ¹⁸O_atm and $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ from the Dome Fuji and EDC ice cores show synchronous variations between 140 and 115 ka BP, although the $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements from EDC are more scattered than DF data due to the use of smaller samples (see Sect. S1 in the Supplement). However, δD records are still slightly discordant, and the EDC record lags DF by up to 1700 years over MIS 5e and at the onset of the Antarctic Isotope Maximum (AIM) 24 (Fig. 11), suggesting some remaining chronology problems (AIM 24 onset) or regional climatic differences (δD decrease over MIS 5e). Between 180 and 150 ka BP, AICC2012 shows better agreement with the DF2021 chronology than the new AICC2023 chronology, which suggests younger ages as per TAC and δ¹⁸O_atm dating constraints.

4.2.2 MIS 11 (from 425 to 375 ka BP)

Over the time interval from 430 to 360 ka BP, encompassing MIS 11, the new AICC2023 chronology predicts older ages than AICC2012 (by up to 2 kyr) with a diminished uncertainty (from 3.9 to 1.7 kyr). This shift towards older ages is induced by δ¹⁸O_atm–δ¹⁸O_calcite (Hulu, Sambao and Dongge caves) tie points at 377.3, 385.7 and 398.5 ka BP and by the TAC age constraint at 362.1 ka BP (Figs. 7 and 12). As a result, two major rises in the EDC atmospheric CO₂ and CH₄ concentration records (corresponding to carbon dioxide jumps CDJ+ 11a.3 and 11a.4, labeled as per Nehrbass-Ahles et al., 2020) occur at 385.6±1.4 and 389.8±1.5 ka BP (Fig. 12). These two rapid jumps in CO₂ and CH₄ are better aligned with two abrupt decreases in the highly resolved δ¹⁸O_calcite record of Zhao et al. (2019) from Yongxing cave (independently dated with ²³⁰Th at 386.4±3.1 and 390.0±3.0 ka BP) than when using the AICC2012 chronology (improvement by ∼800 years). Such millennial-scale synchronicity is expected between CH₄ and δ¹⁸O_calcite series from Chinese speleothems as they are both influenced by Asian monsoon area displacements (and associated methane emissions from wetlands) (Sánchez Goñi et al., 2008).

Figure 12 Evolution of climate tracers from the EDC ice core and Yongxing cave stalagmites between 405 and 370 ka BP. EDC records of (a) δD, (b) CH₄ (Nehrbass-Ahles et al., 2020), (c) CO₂ (Nehrbass-Ahles et al., 2020) and (d) δ¹⁸O_atm on AICC2012 (gray triangles) and AICC2023 (blue circles) chronologies. (e) Speleothem δ¹⁸O_calcite from Hulu, Dongge and Sambao cave, used to constrain AICC2023 (dashed red curve; Cheng et al., 2016), and from Yongxing cave, independently dated with ²³⁰Th dating and annual band counting (plain brown curve; Zhao et al., 2019). CDJ+ is labeled as per Nehrbass-Ahles et al. (2020). Vertical dashed black and blue bars show jumps in CO₂ on AICC2012 and AICC2023 chronologies, respectively; vertical red bars show corresponding decreases in δ¹⁸O_calcite. Black lines show the three tie points between δ¹⁸O_atm and δ¹⁸O_calcite (Cheng et al., 2016) used to constrain AICC2023.

4.2.3 MIS 19 (from 780 to 760 ka BP)

The Matuyama–Brunhes event (geomagnetic field reversal) is reflected by a globally synchronous event in the ¹⁰Be signal: an abrupt termination of the large ¹⁰Be peak following a long-term increasing trend recorded in both ice and sedimentary cores (Giaccio et al., 2023). The ⁴⁰Ar/³⁹Ar age-constrained chronology of a lacustrine succession from Sulmona basin (Giaccio et al., 2023) gives an age of 770.9±1.6 ka BP for the ¹⁰Be peak termination. The new AICC2023 chronology provides an estimate of 767.3±3 ka BP for the same ¹⁰Be peak termination, an age which is closer to the ⁴⁰Ar/³⁹Ar age evaluation than the AICCC2012 chronology estimate (766.2±3 ka BP; Fig. 13). The new AICC2023 chronology indeed indicates an increasingly older age than AICC2012 over MIS 19 (from 790 to 761 ka BP) due to the new δ¹⁸O_atm-based timescale (Fig. 7).

Figure 13 EDC ¹⁰Be record on AICC2012 and AICC2023 chronologies between 778 and 760 ka BP. Vertical gray and blue bars indicate the age of the abrupt EDC ¹⁰Be peak termination on AICC2012 (gray triangles) and AICC2023 (blue circles) chronologies, respectively. The horizontal gray and blue rectangles correspond to AICC2012 and AICC2023 2σ confidence intervals (±3 and ±2.7 ka, respectively). The vertical red bar and horizontal rectangle show the ¹⁰Be peak termination age and its 2σ confidence interval (770.9±1.6 ka BP; Giaccio et al., 2023).

We acknowledge that the Chiba composite section also provides a high-resolution ¹⁰Be record, as the Montalbano Jonico marine section (Simon et al., 2017), the Sulmona basin succession and the EDC ice core do. Although the ¹⁰Be flux records of Sulmona and EDC show a similar pattern and the same asymmetrical shape (i.e., slow increase followed by an abrupt ¹⁰Be peak termination), the sharp termination is less obvious in the Montalbano Jonico and Chiba records. In addition, the Chiba and Montalbano Jonico records are shallow marine deposits; hence expression of paleoclimatic proxies can be amplified and/or hampered by fluvial input (Nomade et al., 2019). Finally, substantial adjustments, up to 10.2±5.5 kyr (i.e., exceeding the related uncertainty), are required to fit the millennial-scale variability in the Chiba record within the Sulmona radioisotopic-based chronology. Giaccio et al. (2023) point out that, despite these relatively large temporal offsets for the Chiba record, the Sulmona-based age model is more linear and describes a simpler, and likely more realistic, history of sediment accumulation. Therefore, we rather use the Sulmona succession to compare with AICC2023.

5 Conclusions

In this study, we have established a new reference chronology for the EDC ice core, AICC2023, covering the last 800 kyr, which is consistent with the official GICC05 timescale over the last 60 kyr. A valuable update of the chronology construction has been the compilation of chronological and glaciological information, including new age markers from recent high-resolution measurements on the EDC ice core. As a result, the chronological uncertainty is reduced from 1.7 kyr in AICC2012 (standard deviation of 995 years) to 900 years on average in AICC2023 (standard deviation of 720 years). A total of 90 % of the new AICC2023 timescale is associated with an uncertainty lower than 2 kyr, compared to only 60 % in the AICC2012 chronology. First, the distinct orbital chronologies derived from $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, TAC and δ¹⁸O_atm are coherent within their respective uncertainties except over three periods, including MIS 11 and MIS 19. Second, new δ¹⁵N measurements along with new sensitivity tests with the firn densification model described by Bréant et al. (2017) and adapted for the EDC ice core provide the most plausible evolution of LID at EDC over the last 800 kyr.

The majority of the age disparities observed between the AICC2023 and AICC2012 chronologies are smaller than 500 years (median), hence minor given the average uncertainty in AICC2012 (1.7 kyr). Exceptions are significant age shifts reaching 3.4, 3.8 and 5 kyr towards older ages suggested over MIS 15, MIS 17 and MIS 19, respectively. However, most of these age discrepancies lead to an improved coherency between the new EDC timescale and independent absolute chronologies derived for other climate archives, especially over the following periods: MIS 5, MIS 11 and MIS 19.

We have identified time intervals where building the chronology is more complicated, such as termination VI (from 540 to 456 ka BP) and from 800 to 600 ka BP, corresponding to the lowermost section of the core, and we would like to draw attention to the requirement for new measurements over these periods. In particular, the links between the variability in $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC records and their orbital targets are not obvious over the 800–600 ka BP period (Fig. 1). This may be due to bad quality of the ice and/or diffusion of gases through the ice matrix (Bereiter et al., 2009). The imprecision of the signal may also be partially explained by the limited temporal resolution of the existing dataset in this deep section. To address these issues, highly resolved $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and TAC measurements are needed in the lowermost section of the EDC ice core. In addition, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ from ice samples over the period covering TVI should also be analyzed to investigate the mismatch between old and new datasets (Fig. 1).

A final important aspect would be to further extend the Paleochrono dating experiment by implementing other ice cores such as Dome Fuji, WAIS Divide and NEEM (North Greenland Eemian), for which a large amount of chronological and glaciological information is now available.