Sub-millennial climate variability from high-resolution water isotopes in the EPICA Dome C ice core

The EPICA Dome C (EDC) ice core provides the longest continuous climatic record, covering the last 800 000 years (800 kyr). A unique opportunity to investigate decadal to millennial variability during past glacial and interglacial periods is provided by the high-resolution water isotopic record (δ¹⁸O and δD) available for the EDC ice core. We present here a continuous compilation of the EDC water isotopic record at a sample resolution of 11 cm, which consists of 27 000 δ¹⁸O measurements and 7920 δD measurements (covering, respectively, 94 % and 27 % of the whole EDC record), including published and new measurements (2900 for both δ¹⁸O and δD) for the last 800 kyr. Here, we demonstrate that repeated water isotope measurements of the same EDC samples from different depth intervals obtained using different analytical methods are comparable within analytical uncertainty. We thus combine all available EDC water isotope measurements to generate a high-resolution (11 cm) dataset for the past 800 kyr. A frequency decomposition of the most complete δ¹⁸O record and a simple assessment of the possible influence of diffusion on the measured profile shows that the variability at the multi-decadal to multi-centennial timescale is higher during glacial than during interglacial periods and higher during early interglacial isotopic maxima than during the Holocene. This analysis shows as well that during interglacial periods characterized by a temperature optimum at the beginning, the multi-centennial variability is strongest over this temperature optimum.

1 Introduction

Water stable isotopes (oxygen, δ¹⁸O, and hydrogen, δD) in ice cores are valuable proxy records that can be used to reconstruct past temperatures in polar regions. The water isotopic composition from an ice core sample is classically measured with the delta notation (δ), which expresses the variations of the isotopic ratio of heavy to light isotopes in the water molecule (i.e. ¹⁸O $/$ ¹⁶O and D $/$ H for δ¹⁸O and δD). During air mass transportation, distillation of moisture from the low-latitude regions of evaporation to the polar regions leads to a preferential loss of heavy isotopes (H${}_{\mathrm{2}}^{\mathrm{18}}$O and HD¹⁶O vs. H${}_{\mathrm{2}}^{\mathrm{16}}$O) during successive precipitation events and hence to a decrease of δ¹⁸O and δD toward cold regions. Despite known limitations due to temporal changes in intermittency of precipitation (Casado et al., 2020), vapor origin and transport (Helsen et al., 2006) and sea ice extent (Noone, 2004), and changes in condensation vs. surface temperatures (Buizert et al., 2021) or deposition and post-deposition effects (Casado et al., 2018), the spatial relationship between surface temperature and surface snow δD and δ¹⁸O has long been used to establish an isotopic palaeothermometer to infer past temperature variations, at least qualitatively (Jouzel et al., 2013).

Today, the oldest continuous isotopic record from ice cores is that retrieved through the European Project for Ice Coring Antarctica (EPICA) Dome C ice core (EDC), which covers the last 800 000 years (800 kyr) (Jouzel et al., 2007). The first analyses of water isotopic composition (δD) over the EDC ice core were performed at ∼ 4 m resolution and unveiled the δD variations over eight glacial–interglacial cycles (EPICA community members, 2004). Several years later, measurements of δD in bag samples (continuous 55 cm pieces of the EDC ice core) evidenced millennial-scale variability over the glacial periods in Antarctica (Jouzel et al., 2007; Stenni et al., 2010). In order to explore potential changes in high-frequency variability between interglacial periods, Pol et al. (2010, 2011, 2014) used 11 cm resolution δD measurements over interglacial periods during Marine Isotopic Stages (MIS) 5, 11 and 19, i.e. the periods between 112 and 134 ka (i.e. kyr before present, where “present” is taken to be 1950 AD), between 392 and 427 ka, and between 747 and 800 ka, respectively. Landais et al. (2015) focused on 11 cm resolution δ¹⁸O over the last glacial period back to 60 ka.

Pol et al. (2014) used the high-resolution water isotopic signals over the MIS 5 and 11 interglacial periods to estimate the relative variations of decadal to centennial climate variability during these interglacial periods with respect to the Holocene (Pol et al., 2011, 2014). Over the last glacial period, the high-resolution δ¹⁸O record showed an enhanced amplitude of the multi-decadal to centennial variability during the warm phases of the Antarctic Isotopic Maxima, AIM (Landais et al., 2015). These AIM events are key climatic features of the last glacial period: they are counterparts to the abrupt temperature increases in the Northern Hemisphere that were first identified in the Greenland ice cores (Dansgaard, 1985; Blunier and Brook, 2001; EPICA community members, 2004).

High-resolution water isotopic measurements over the EDC ice core are hence key to documenting the temporal patterns of climatic variability over the past 800 kyr. Unfortunately, the analytical load required to obtain the full 800 kyr record at 11 cm resolution is enormous, and would represent 35 000 measurements. Even though several individual studies have been published, a complete synthesis of EDC high-resolution δD and δ¹⁸O records over the last 800 kyr is still missing. This is an important limitation on the documentation of past changes in sub-orbital climatic variability in Antarctica and a comparison of the climatic variability features between glacial and interglacial periods or between different interglacial (glacial) periods. As an example, the interglacial periods before the Mid-Brunhes Transition (MBT, 430 ka) are cooler than the five most recent interglacial periods, but there is limited evidence available to document climate variability during interglacial periods before and after the MBT at high resolution (Barth et al., 2018; Past Interglacials Working Group of PAGES, 2016). A first challenge is thus to provide homogeneous high-resolution isotopic records.

A second challenge when attempting to characterize the past high-frequency climate variability in Antarctica is non-temperature-related variability in the water isotopes from the depositional process (Fischer et al., 1985; Laepple et al., 2018) and alteration and smoothing effects of post-deposition processes. Indeed, post-deposition processes (Casado et al., 2018, 2020; Steen-Larsen et al., 2014) and firn and ice diffusion (Gkinis et al., 2011, 2021a) strongly limit the interpretation of water isotopic variability in terms of climatic variability. In the case of old ice, the impact of diffusion, which increases with depth and age, can reach the multi-centennial timescale and affect the climate variability recorded in δ¹⁸O and δD. Pol et al. (2010) showed that the 11 cm resolution δD record of MIS 19 (3147–3190 m deep in the EDC ice core) did not provide more information than the 55 cm resolution record due to the large impact of diffusion at this depth of the core. This effect is particularly important to quantify for the 1.5 Ma ice core to be drilled in East Antarctica. Indeed, documenting the evolution of diffusion length with depth is key to anticipating the kind of information on climate variability that can be retrieved from the deepest part of this future ice core.

Here we address the two aforementioned challenges (high-resolution records and the influence of diffusion) by presenting a compilation of new high-resolution measurements of δ¹⁸O and δD in the EDC ice core. The first section presents the analytical methods used to perform high-resolution measurements of δD and δ¹⁸O in the different sections of the EDC ice core in recent decades, as well as methods for spectral analyses and the calculation of isotopic diffusion along the EDC ice core. The second section describes how the different measurements performed over the past 20 years in different institutes with different analytical methods can be compiled together into a single record. The third section uses the high-resolution measurements to investigate changes in sub-orbital climatic variability across the last 800 kyr and how diffusion can affect some of the observed features.

2 Materials and methods

2.1 The EPICA ice core

The Franco-Italian Concordia Station is located 3233 m above sea level on the continental plateau of Antarctica (75^∘06^′12^′′ S 123^∘21^′30^′′ E). The mean annual surface temperature is −54.5 ^∘C and the snow accumulation rate is ca. 25 mm water equivalent per year (EPICA community members, 2004; Le Meur et al., 2018).

The EDC ice core was drilled at Concordia Station on Dome C, where the ice flow was small (EPICA community members, 2004). The EDC drilling project started in 1996 and was completed in 2004. In 1999, a second ice core (EDC2) was drilled from the surface because the drill for EDC1 was stuck at depth of 788 m. Bedrock was reached in 2004 at a depth of 3190 m. From here onwards, we refer to EDC1 and EDC2 as “the EDC ice core”. By the time this second ice core had been retrieved, the full 788 m of EDC1 had been analysed. Later, EDC2 measurements started 19 m higher than the bottom end of the EDC1 ice core in order to have an overlap to reconnect the two cores without duplicating all the measurements in the common depth range.

After drilling and core logging, the EDC ice core was cut into 55 cm long sections, and each section was further cut longitudinally on-site for several measurements (e.g. water isotopes, physical properties, ¹⁰Be, chemistry, and gas analysis). The archival piece (∼ one-quarter of the section) was stored in polystyrene boxes in the EPICA snow cave at Concordia Station at −50 ^∘C. Two types of contiguous samples were dedicated to the analyses of water isotopes in the EDC ice core. First, a 55 cm long stick with a 1 cm² cross section was melted and stored on-site in plastic bottles for the low-resolution measurements. The second was a 55 cm long stick with a 2 cm² cross section that was cut into 11 cm long samples. Each sample was placed in a plastic sheath cut to obtain a plastic bag of the right dimensions, and then the bag was thermally sealed. The sample was stored at −20 ^∘C for a few months prior to being melted and transferred into plastic bottles that were kept at −20 ^∘C.

2.2 Measurement techniques and coherency of the dataset

Several analytical techniques have been used to measure δ¹⁸O and δD in the EDC1 and EDC2 ice cores (Tables 1 and 2). Initial analytical techniques included a uranium reduction method for δD (Vaughn et al., 1988) and a CO₂–H₂O equilibrium method for δ¹⁸O (Meyer et al., 2000), while the most recent method used to determine δ¹⁸O and δD on the EDC2 ice core is cavity ring-down spectroscopy (CRDS) (Kerstel and Gianfrani, 2008; Busch and Busch, 1999). When uranium reduction was replaced by CRDS measurements at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE), an extensive series of comparisons were performed which showed that there was excellent agreement between the two methods within the uncertainty ranges of the instruments. The analytical precisions for the methods are comparable, with 2σ values ranging between 1 and 1.4 ‰ for δD and between 0.1 and 0.4 ‰ for δ¹⁸O (Table 2). Figure 1 displays the full high-resolution (11 cm) datasets for δ¹⁸O and δD of water over the EDC ice core.

Table 1 Summary of available δD measurements of the EDC ice core and associated analytical methods. 2σ values come from the instrumental measurement uncertainty as provided in the original studies.

Table 2 Summary of available δ¹⁸O measurements of the EDC ice core and associated analytical methods. 2σ values come from instrumental measurement uncertainty as provided in the original studies.

Figure 1 EDC ice core and other palaeoclimate records as well as variations in Milankovitch cycles over the past 800 kyr. (a) Precession (pink) and obliquity (blue) from Laskar et al. (2004); (b) composite EDC and Vostok CO₂ record over the last 800 kyr (Lüthi et al., 2008; Bereiter et al., 2015); (c) the EDC δD record at 11 cm (orange) and 55 cm (black) resolution; and (d) the EDC δ¹⁸O record at 11 cm (light blue) and 55 cm (dark blue) resolution. All ice core records are presented on the AICC2012 scale (Bazin et al., 2013; Veres et al., 2013). Grey rectangles indicate the positions of interglacial periods.

2.3 Multi-resolution analysis (MRA)

A discrete wavelet analysis was used to identify the contributions to the overall isotopic variability from signals of different periodicities (i.e. corresponding to decadal to multi-millennial signal variability). We produced a multi-resolution analysis (MRA) using R software with the waveslim wavelet package (Whitcher, 2020) containing the MRA function with a Daubechies orthonormal wavelet filter. MRA is a mathematical analysis tool which decomposes a signal at different resolution levels. An important feature of MRA is its ability to capture temporally localized changes at its nearest neighbour. A low (high) resolution level corresponds to a coarse (detailed/high-frequency) component of the original signal. Each MRA level can thus be used to interpret the temporal variability within a frequency range. Adding all MRA levels exactly reproduces the original undecomposed signal. As the wavelet analysis needed to be applied at time intervals with a uniform sample resolution, we divided the EDC isotopic record on the AICC2012 age scale (Bazin et al., 2013) into six intervals. These ranged from the youngest interval, between 0 and 56 ka (where the longest time span covered by 11 cm is 10 years), to the oldest interval (at the bottom of the core), between 651 and 800 ka (where the longest time span covered by 11 cm is 320 years) on the AICC2012 age scale (Bazin et al., 2013) (Table 3). Over each interval, we performed an interpolation with a uniform resolution corresponding to the longest time span covered by 11 cm of ice (i.e. interpolation at 10 years between 0 and 56 ka; at 20 years between 56 and 144 ka; see the details for all periods given in Table 3).

Table 3 Time resolutions (in years) of the different MRA decompositions for specific intervals (0–56, 56–144, 144–305, 305–420, 420–651 and 651–800 ka). Letters a, b, c, d, e and f represent segments that have the same time resolution and can be combined.

The resolution of the MRA was chosen to increase by a factor of 2 between neighbouring intervals i and i+1, i being a number between 1 and 5. As a consequence, the second MRA decomposition of the interval i has the same resolution as the first MRA of the interval i+1 (Table 3). We then concatenated the MRAs with the same temporal resolution, leading to nine successive composites (named a, b, c, d, e, f, g, h and i in Table 3). The longest (composites f, g, h and i) corresponded to the variability of the signal at 320 years resolution and covered the whole 800 ka, and the shortest (composite a) corresponded to the variability of the signal at 10 years resolution and covered only the last 56 kyr.

2.4 Effect of isotopic diffusion

The effect of isotopic diffusion with depth is convolved using a function G(z) of associated diffusion length σ_z (Gkinis, 2011; Laepple et al., 2018; Gkinis et al., 2021a):

$$\begin{array}{}\text{(1)}& G\left(z\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_z\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left({\displaystyle \frac{-z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_z^{\mathrm{2}}}}\right),\end{array}$$

where z is the depth along the ice core and σ_z is the diffusion length.

We quantify the amplitude decay of the signal between the initial amplitude A₀ and the measured amplitude A at a certain depth as described in Johnsen et al. (2000) and Gkinis et al. (2021a) for a given period λ with the following equation:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{A}{A_{\mathrm{0}}}}=\exp \left(-\mathrm{2}{\left({\displaystyle \frac{\mathit{\pi }\times {\mathit{\sigma }}_z}{\mathit{\lambda }}}\right)}^{\mathrm{2}}\right).\end{array}$$

For our purpose, the diffusion length along the EDC ice core is calculated by considering the firn diffusion (i.e. due to water vapour diffusion in the open porosity) and the ice diffusion (i.e. due to water molecule diffusion in the ice matrix).

We used two different estimates for the firn diffusion length, σ_firn, along the EDC ice core. In a first approach, we assumed a constant σ_firn all along the EDC ice core and took the value of 0.07 m estimated by Johnsen et al. (2000) for the EDC. In a second, refined approach, we considered a changing σ_firn between interglacial and glacial periods, as described in Gkinis et al. (2021a). This causes σ_firn to vary from 0.075 m in interglacial periods to 0.065 m in glacial periods. The ice thinning, S, also affects the visible effect of the firn diffusion length along the ice core, so the thinned firn diffusion length should be ${\mathit{\sigma }}_{\mathrm{thinned}\mathrm{\_}\mathrm{firn}}=S\times {\mathit{\sigma }}_{\mathrm{firn}}$. In this study, for consistency, we used the thinning function for the EDC ice core corresponding to the AICC2012 chronology (Bazin et al., 2013).

The ice diffusion depends on the thinning and the temperature. The following formulation permits the calculation of the diffusion length associated with ice diffusion, σ_ice, as a function of the age (and depth) of the ice (Gkinis et al., 2011):

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{ice}}^{\mathrm{2}}\left(\mathit{\tau }\right)=S{\left(\mathit{\tau }\right)}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathit{\tau }}{\int }}\mathrm{2}D\left(t\right)S{\left(t\right)}^{-\mathrm{2}}\mathrm{d}t,\end{array}$$

where S is the thinning of the ice layers at the considered age τ. In order to estimate the ice diffusion coefficient D(t), we used the classical formulation of Ramseier (1967):

$$\begin{array}{}\text{(4)}& D=D_{\mathrm{0}}\times \exp \left({\displaystyle \frac{-Q}{RT}}\right),\end{array}$$

with D₀ = 9.13 cm² s⁻¹ and Q = 59.820 kJ mol⁻¹. At −50 ^∘C, D is equal 8.866 × 10⁻¹⁴ cm² s⁻¹; at −10 ^∘C, D is equal to 1.1993 × 10⁻¹¹ cm² s⁻¹. T represents the ice temperature from the borehole.

The total calculated diffusion length that was expected to be measured in the ice core was estimated using the diffusion length associated with the firn diffusion and the diffusion length associated with ice diffusion in a quadratic addition, so that

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_z=\sqrt{\left({\mathit{\sigma }}_{\mathrm{ice}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{thinned}\mathrm{\_}\mathrm{firn}}^{\mathrm{2}}\right)}.\end{array}$$

The increase of the diffusion length with increasing depth in the ice core is shown in Fig. S1 in the Supplement. It is mainly due to the increase in temperature. Indeed, the borehole temperature evolves almost linearly from −53.5 to −2.6 ^∘C along the 3255 m ice core (Buizert et al., 2021). The variation of the calculated diffusion length around 3000 m is explained by the variability of the thinning function (Parrenin et al., 2007).

3 Coherency of different analytical measurements

Different analytical instruments and techniques have been used to determine δ¹⁸O and δD in the EDC1 and EDC2 ice cores at different laboratories (Table 1). To determine the coherency of the different datasets, two different comparisons were performed; (1) a comparison of the isotopic values from the same samples measured by different analytical techniques; and (2) a comparison of the 55 cm sample resolution data with the 11 cm sample resolution data using a five-point average. We averaged the 11 cm resolution in a five-point window to compare it with the 55 cm resolution measurements in exactly the same window.

3.1 Comparison of isotopic data using different analytical techniques

The CRDS analysis in 2019–2020 measured previously analysed samples from 2004–2010 by uranium reduction for δD in MIS 5.5 (1670–1693 m) and by H₂O–CO₂ equilibration for δ¹⁸O (1670–1793 m) (Table 1). Figures 2 and 3 provide two examples of analyses performed for all overlapping intervals. Additional comparisons of isotopic data measured by different analytical techniques on the same samples are also presented in the Supplement (Figs. S5 and S6). The samples stayed refrozen between the different measurements and they were refrozen immediately after analysis. Tests were performed by storing low δ¹⁸O and δD internal standards for several years in the freezer. In some cases, but not systematically and not significantly compared to the analytical precision, small increases of δ¹⁸O and δD were obtained. In a comparison of the old and new records, we did not observe systematic increases of δ¹⁸O and δD for the samples analysed recently compared to the analyses performed 15 years ago, so unfortunately we cannot give a solid explanation for the small differences between the series of measurements.

Figure 2 (a) EDC δD measurements versus depth (m) over Termination 2: measurements completed in 2010 at LSCE (uranium reduction method; Pol et al., 2014) are shown in blue, and measurements completed in 2019 at LSCE (CRDS method) are shown in red. (b) Difference between the δD values measured in 2010 and 2019. (c) Probability density function for the difference between the first (uranium reduction) and the new (CRDS) δD measurements. A Gaussian curve (red) is fitted to the data. Another Gaussian curve (green) with the standard deviation equal to the classically displayed 1σ uncertainty of δD measurements with the CRDS method at LSCE (1σ = 0.7 ‰) is also displayed.

Figure 3 (a) EDC δ¹⁸O measurements versus depth (m) over Termination 6: measurements completed in 2010 at the University of Triestre with the CO₂ equilibration method are shown in blue, and measurements completed in 2019 at LSCE (CRDS method) are shown in red. (b) Difference between the δ¹⁸O values measured in 2010 and 2019. (c) Probability density function for the difference between the old (University of Triestre) and the new (LSCE) δ¹⁸O measurements. A Gaussian curve (red) is fitted to the data. Another Gaussian curve (green) with the standard deviation equal to the classically displayed 1σ uncertainty of δ¹⁸O measurements by CRDS at LSCE (1σ = 0.2 ‰) is also displayed.

3.2 δD comparison

The difference between analytical techniques (Fig. 2) seems to depend on the absolute value for δD (there is a negative difference for low δD values). This was confirmed by a statistical test of the correlation between the absolute value of δD and the δD difference between the two series of measurements, which led to a Pearson coefficient of 0.13 and a p-value of 0.003. Such an isotopically dependent feature may arise from a possible calibration effect, despite the fact that exchanges of internal laboratory water standards and regular intercalibrations between the laboratories measuring water isotopes of the EDC ice cores were performed.

Despite such a tendency for the δD differences between the two series, the absolute value of the difference remains small. We used a Welch t-test to show that the 2010 and 2019 time series have equal means at the 99.9 % confidence level (t = 3.5, N = 1000) with respect to the experimental uncertainties.

Finally, the distribution of the differences between the 2010 and the 2019 δD measurements is not Gaussian and is not centred around zero (Fig. 2c). Still, this distribution is narrow and is encompassed within a Gaussian distribution where 2σ = 1.4 ‰ is associated with the classical analytical uncertainty of the δD measurements. Note that the analytical uncertainty associated with this CRDS measurement series was evaluated from the analysis of the difference between the same samples (1000 samples, which represent 10 % of the whole series) measured twice, 1 to 3 months apart. We thus conclude that the δD difference between the uranium reduction vs. CRDS datasets is smaller than the uncertainty associated with CRDS measurements, and thus that we can combine the different datasets if we consider a 2σ uncertainty of 1.4  ‰ in the final δD data.

3.3 δ¹⁸O comparison

No significant statistical differences are observed between the δ¹⁸O measurements performed using the CO₂ equilibrium and CRDS methods. The standard deviation of the difference series between 2010 and 2019 δ¹⁸O measurements (2σ = 0.2 ‰) is smaller than the classical analytical uncertainty of the δ¹⁸O measurements by CRDS (2σ = 0.4 ‰) (Fig. 3). A statistical test of the correlation between the absolute value of δ¹⁸O and the δ¹⁸O difference between the two series of measurements led to a Pearson coefficient of 0.0049 and a p-value of 0.9. In addition, we did a Welch t-test of the equality between two averages in δ¹⁸O to check whether the difference between the 2010 and 2019 δ¹⁸O data is significantly different within the experimental margin of error. Upon doing so, the result shows that the two series have equal means at the 70 % confidence level (t = 0.557, N = 1000).

Second, we compared the low-resolution (55 cm) and the high-resolution (11 cm) δ¹⁸O series after calculating the average δ¹⁸O of five 11 cm samples that overlap with the same sample depth as the 55 cm samples (Figs. S2 to S4). The difference between the two time series is 0.008 ± 0.001 ‰ (Fig. S4). A comparison between low- and high-resolution series for δD has already been performed by Pol et al. (2011, 2014). In the 2011 paper, the coherency between the 55 and 11 cm samples was studied by calculating the average signal over five 11 cm data. They observed that the signal from the 55 cm samples is similar to the average signal but has a lower statistical accuracy (1σ = 0.5 ‰) than the average signal (1σ = 0.23 ‰).

The two comparisons performed above suggest there is no statistically significant difference in δ¹⁸O and δD between the datasets compiled here (Fig. 1).

4 Results and discussion

4.1 Recorded multi-decadal to multi-millenial isotopic variability over the last 800 kyr

The compiled high-resolution EDC water isotope record is presented in Fig. 1. We applied the MRA decomposition to each of the six selected intervals (see the “Methods” section) and we present the decadal to multi-millennial variability across the last 800 kyr (Fig. 4). We calculated the running standard deviation (1σ) in a 3 kyr window and we use this value as an estimate of the level of variability. For the first MRA composite at 10 years resolution (Fig. 4a), we observe stronger isotopic variability during the Holocene than during the Last Glacial Maximum (average 1σ of 0.46 ‰ and 0.24 ‰, respectively). The 20-year variability (Fig. 4b) inferred from the second composite shows a globally uniform pattern over the last 150 kyr. The 80-year variability (Fig. 4d) is smaller during the interglacial periods (1σ = 0.18 ‰) than during the glacial periods (1σ = 0.30 ‰) over the last 400 kyr. The 160-year to 640-year variability (Fig. 4e to g) also shows a small decrease of variability over interglacial periods and decreasing variability for the oldest ice core sections. For the lower-frequency variabilities (composite at resolutions of 1280 and 2560 years, Fig. 4h to i), the amplitude of the variability envelope increases during glacial inceptions and glacial periods, with a notable strong 2,560-year variability at the onset of MIS 9 (1σ = 1.13 ‰, compared to an average of 1σ = 0.20 ‰ over the whole series). The large centennial to multi-centennial water isotope variability in glacial periods is linked tothe successive Antarctic Isotopic Maxima (AIM) during glacial periods (EPICA community members, 2004; Jouzel et al., 2007). Finally, the decreasing amplitude of the signal variability toward older ages is probably the result of the diffusion of water isotopes in open-porosity firn and ice crystals. While we can disentangle the effects of diffusion and climate-driven isotopic variability for low-frequency signals and at greater depths, the respective influences of diffusion and climate are less obviously identified at shallower depths and for high frequencies.

Figure 4 Contributions to the original δ¹⁸O signal (red) of the MRA composites with resolutions of 10 (a), 20 (b), 40 (c), 80 (d), 160 (e), 320 (f), 640 (g), 1280 (h) and 2560 years (i). Marine Isotope Stage intervals are marked as grey bars. The black envelope presents the running standard deviation (1σ) in 3 kyr windows.

4.2 Effect of isotopic diffusion on the recorded signal variability

We evaluate the effects of diffusion on the isotopic signal recorded in the ice core by computing the decrease of Holocene variability from Eq. (2). The calculated $A/A_{\mathrm{0}}$ signal amplitude is hence scaled for each MRA composite to the mean amplitude of the variability of the MRA composite signal between 2 and 8 ka for each resolution (Fig. 6).

Figure 5 High-resolution record (light blue) and comparison of its variability (3 kyr standard deviation, black) to the variability (3 kyr standard deviation) of the diffused Holocene signal for the different periods (10, 20, 40, 80, 160 and 320 years for a to f). The diffused Holocene signal was calculated using two σ_firn estimates: a constant σ_firn of 7 cm (dark blue) and a variable σ_firn equal to 6.5 cm in glacial periods and 7.5 cm in interglacial periods (pink).

Figure 6 Contributions to the original δ¹⁸O signal (red) of the MRA composites with resolutions of 20 (b), 40 (c), 80 (d), 160 (e), 320 (f) and 640 (g) years for MIS 5 (left panel) and MIS 9 (right panel)​​​​​​​. The black envelope presents the running standard deviation (1σ) in 3 kyr windows. The red rectangles indicate periods with enhanced variability and the blue rectangles indicate periods with reduced variability.

As explained in the section “Methods”, we used the ice diffusion coefficient from Ramseier (1967) with two different estimates for σ_firn. The different estimates of σ_firn did not have a significant effect on the calculated amplitude of the variability (Fig. 5).

Diffusion has the expected effect of decreasing the amplitude of the variability of the isotopic signal for older and deeper ice core sections (Fig. 4). In the 10-year series (Fig. 4a), diffusion dampens the amplitude of the recorded variability in the last glacial period by half compared to that in the Holocene. The calculated amplitude of the variability due to diffusion is actually much smaller than the recorded one, which suggests that either the 10-year isotopic variability during the last glacial period is larger than the 10-year variability during the Holocene or that measurement noise dominates the 10-year variability.

For the deepest sections of the ice core, i.e. sections older than 600 ka, the diffusion model overestimates the damping of centennial and multi-centennial variability compared to what is retrieved from the ice core isotopic composition. This discrepancy calls for future reassessment of the isotopic diffusivity in the deepest sections of the EDC ice core.

4.3 Climate variability at different time intervals over the last 800 kyr

Combining our high-resolution water isotopic records with frequency analysis and the impact of diffusion, we can suggest some patterns for the decadal to millennial climate variability over the last 800 kyr.

First, at the decadal scale, our findings show larger variability during the last glacial period compared to the Holocene. The analysis of Jones et al. (2017) using the water isotopic record in the WAIS Divide record in West Antarctica supports this higher variability at the decadal scale during the last glacial maximum. At the high accumulation site of WAIS, diffusion has a minimal effect on the variability with a 4- to 15-year periodicity, and the higher water isotopic variability observed during this period is interpreted as an increase in the strength of the teleconnections between the tropical Pacific and West Antarctica (Jones et al., 2018). Jones et al. (2018) invoke the expansion of the Northern Hemisphere ice sheets during the LGM, leading to a shift in the location of tropical convection, to explain these characteristics. The same pattern is observed for the 20-year periodicity (Fig. 5b), i.e. the variability of the water isotopic composition after diffusion is smaller than the measured one during the last glacial period, while there is a good agreement between the diffused and measured signal over MIS 5e. For the 40- and 320-year periodicity (Fig. 5c to f), the variability of the last glacial period is also higher than the diffused Holocene variability. This is also the case for MIS 6 for the 80-year periodicity (Fig. 5d) and for MIS 8 and 10 for the 160- and 320-year periodicity (Fig. 5e to f). For these periods and frequency ranges, the impact of diffusion on the variability is limited and the isotopic signal in the ice core is preserved. The multi-centennial variability increase during glacial periods can be related to the presence of AIM.

Our analysis shows that there is a clear enhanced isotopic variability during glacial periods at the multi-decadal to multi-centennial timescale in the EDC ice core, which could be attributed to climate variability (Fig. 4). This result is in agreement with the findings of Rehfeld et al. (2018), who used a worldwide data synthesis to show increased interannual to millennial climatic variability during the last glacial maximum with respect to the Holocene at all latitudes along with a factor of 2 increase in variance at the high latitudes of the Southern Hemisphere, a result in agreement with output of coupled model simulations (Rehfeld et al., 2018).

Second, while the effect of diffusion is important when we want to compare the variability in one interglacial period to that in another, it does not affect the evolution of the recorded variability during the course of an interglacial period.

A previous study focused on the warm phase of MIS 5 (115.5 to 132 ka), where wavelet analysis of the 11 cm resolution δD record showed three different isotopic phases with different levels of variability (Pol et al., 2014). The first phase, from 111 to 119 ka, has a low orbital forcing context, but the variability increases during the entry into glaciation, with centennial periodicities dominant. The second phase, from 119 to 123 ka, is a stable warm phase, warmer than the Holocene. The δD variability of the second phase is notably lower than those of the other phases, 3.7 ‰ compared to a 4.5 ‰ average. Finally, during the third phase, there is again a higher variability, with multi-centennial periodicities dominant between 123 and 133 ka.

When doing a similar analysis with our MRA decomposition, we find similar variability of the high-resolution signal (Fig. 6a), i.e. the maximum amplitude of the multi-decadal to multi-centennial variability of the signal is encountered over the optimum of MIS 5 (phase 3 in Pol et al., 2014; between 125 and 131 ka) and toward the end of this warm period (phase 1 of Pol et al., 2014). The minimum amplitude of the multi-decadal to multi-centennial variability of the signal is encountered between 119 and 123 ka (phase 2 in Pol et al., 2014), when the δ¹⁸O and δD signals are on a plateau.

Thus, during MIS 5, multi-decadal to multi-centennial variability of the water isotopic signal can be interpreted as climate variability at these multi-decadal to multi-centennial timescales. This can be compared to the variability over the interglacial period of MIS 9 (∼ between 325 and 338.5 ka), which is also characterized by a temperature optimum at its start. The amplitudes of the variability for the different MRA decompositions for the interglacial period of MIS 9 cannot be directly compared to those over MIS 5 because of the effects of diffusion and thinning (see Fig. 5). However, in Fig. 6b we observe the same pattern as for MIS5: higher amplitudes for the multi-decadal to multi-centennial variability are observed over the δ¹⁸O optimum (333–338 ka) and at the end of this warm period (321–326 ka), while the minimum amplitudes for the multi-decadal to multi-centennial variability are observed over the plateau of the interglacial period (326–332 ka). This result strengthens the conclusion of Pol et al. (2014) that the climate over the temperature optima of interglacial periods may also be more variable at the multi-decadal to multi-centennial timescale. A parallel can be drawn with the larger high-frequency water isotopic variability observed during the temperature optimum of the AIM of the last glacial period in the EDC ice core (Landais et al., 2015), since this temperature optimum at the beginning of the interglacial could also be the result of millennial-scale variability (Past Interglacials Working Group of PAGES, 2016).

5 Conclusion

Here, we compiled and presented a EDC ice core water isotopic record (δ¹⁸O and δD) using new and previously published 11 cm data spanning the last 800 kyr. This compilation and the comparison performed between different series of measurements showed that the differences between water isotopic data measured by different laboratories and techniques over the last 20 years for the same samples are not statistically significant and are within analytical uncertainty. As a result, all the available EDC water isotope data were combined to produce a continuous high-resolution dataset at mostly 11 cm sample resolution.

MRA decomposition of the water isotopic record at temporal resolutions varying between 10 and 2560 years shows that the variability during glacial periods at multi-decadal to multi-centennial timescales is higher than the variability during the Holocene and that the variability is enhanced over the early temperature optimum during MIS 5 and 9. These results are not influenced by diffusion in the open-porosity firn and in the ice matrix, but the interpretation of high-resolution δD and δ¹⁸O profiles needs to take this effect into account. Finally, our study calls for further analyses quantifying the diffusivity in EDC, which is essential for analyses of the Beyond EPICA – Oldest Ice core.