the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Novel approach to estimate the water isotope diffusion length in deep ice cores with an application to Marine Isotope Stage 19 in the Dome C ice core

### Andrew M. Dolman

### Torben Kunz

### Vasileios Gkinis

### Thomas Laepple

Accurate estimates of water isotope diffusion lengths are crucial when reconstructing and interpreting water isotope records from ice cores. This is especially true in the deepest, oldest sections of deep ice cores, where thermally enhanced diffusive processes have acted over millennia on extremely thinned ice. Previous statistical estimation methods, used with great success in shallower, younger ice cores, falter when applied to these deep sections, as they fail to account for the statistics of the climate on millennial timescales. Here, we present a new method to estimate the diffusion length from water isotope data and apply it to the Marine Isotope Stage 19 (MIS 19) interglacial at the bottom of the EPICA Dome C (EDC, Dome Concordia) ice core. In contrast to the conventional estimator, our method uses other interglacial periods taken from further up in the ice core to estimate the structure of the variability before diffusion. Through use of a Bayesian framework, we are able to constrain our fit while propagating the uncertainty in our assumptions. We estimate a diffusion length of 31±5 cm for the MIS 19 period, which is significantly smaller than previously estimated (40–60 cm). Similar results were obtained for each interglacial used to represent the undiffused climate signal, demonstrating the robustness of our estimate. Our result suggests better preservation of the climate signal at the bottom of EDC and likely other deep ice cores, offering greater potentially recoverable temporal resolution and improved reconstructions through deconvolution.

- Article
(1724 KB) - Full-text XML
- BibTeX
- EndNote

Large ice sheets from the polar regions offer unique insights into the climate up to hundreds of thousands of years ago. The drilling of deep ice cores in Greenland and Antarctica enables measurements of water isotopic ratios (*δ*^{18}O, *δ*D, *δ*^{17}O) impacted by fractionation effects upon evaporation and condensation. These ratios have been shown to relate to atmospheric temperatures at the time of deposition as snowfall (Dansgaard, 1964) and therefore provide a valuable proxy of past-climate conditions. However, water isotope records are not perfectly preserved, partially due to the molecules dispersing over time, smoothing the profile by attenuating high-frequency variability. This displacement is known as diffusion and occurs both in the firn, due to snow–vapour exchange in the pores (Johnsen, 1977; Whillans and Grootes, 1985), and in ice through processes such as ice diffusivity (Ramseier, 1967) and liquid water veins (Nye, 1998). Diffusion in ice is a much slower process than in firn but can act over much longer time periods, with the oldest, deepest sections of deep ice cores being most affected. Since ice diffusion increases with temperature, the warming of deep ice due to geothermal heat from the bedrock further accelerates the process. Additionally, the effect is exacerbated on the temporal scale by extreme layer thinning from ice flow. Collectively, these conditions can result in the attenuation of variability up to millennial timescales in deep ice cores (Pol et al., 2010). Water isotope diffusion can be characterised by the diffusion length, defined as the average displacement of water molecules along the vertical axis relative to their initial position within the ice sheet. In addition to informing which frequencies climate variability is preserved, knowledge of the diffusion length can enable techniques used to recover some of the lost information, such as deconvolution (Johnsen, 1977). These reconstructions are extremely sensitive to the diffusion length, so obtaining accurate values is crucial for interpreting the isotopic data. Furthermore, the temperature dependence of the diffusion process reflected on the magnitude of the diffusion length constitutes the latter as a good candidate for a firn temperature proxy (Simonsen et al., 2011; Gkinis et al., 2014; van Der Wel et al., 2015).

By comparing the power spectrum of a diffused water isotope record with that of the undiffused isotopic climate signal, it is possible to estimate the diffusion length. Current estimation methods assume that the isotopic variability before diffusion is constant across all frequencies, i.e. white noise (Johnsen et al., 2000; Gkinis et al., 2014; Jones et al., 2017; Kahle et al., 2018; Holme et al., 2018). This assumption is justified by the strong noise generated by precipitation intermittency and the stratigraphic noise, which dominates the signal up to decadal or multi-centennial timescales depending on the accumulation conditions of the site (Johnsen et al., 2000; Casado et al., 2020). However, on the longer timescales observed at the bottom of deep ice cores, the isotopic profile is not accurately represented by white noise. This is evident from other ice cores, which demonstrate strong variability over millennia, such as Dansgaard–Oeschger events from Greenland cores (NGRIP members, 2004; NEEM community members, 2013) or Antarctic Isotope Maxima events visible in Antarctic cores (Petit et al., 1999; EPICA community members, 2004). To accurately represent the statistical properties of such millennial variability, the estimation method requires a modified approach.

In this study we estimate the diffusion length for the Marine Isotope Stage 19 (MIS 19) section of the EPICA Dome C (EDC, Dome Concordia) ice core (763–795 ka). Previous estimations derived from water isotope data suggest a diffusion length between 40 and 60 cm for the MIS 19 interglacial (Pol et al., 2010). In contrast, analytical diffusion models using the physical properties of the ice core predict values between 16 and 22 cm (Pol et al., 2010; Grisart et al., 2022), which poses the question whether the model physics is wrong and/or the statistical estimator is biased when applied to deep ice. We address the latter point, removing the assumption that the climate signal can be approximated by white noise. Instead, the initial climate is inferred from similar periods further up in the core where the water isotopes have not undergone significant ice diffusion. We use a Bayesian methodology that propagates uncertainty in this reference climate spectrum while also constraining parameters to physically realistic ranges. This paper explains the new approach and discusses possible future applications and improvements. The new diffusion length estimate is compared with previous estimations and will serve as an indicator for how significant deep-ice diffusion will be in future deep ice cores such as the Beyond EPICA – Oldest Ice Core (BE-OIC).

## 2.1 Water isotope data

We use discrete *δ*^{18}O data from the EDC ice core, measured at the Niels Bohr Institute, University of Copenhagen, with a water–CO_{2} equilibration mass spectrometry system (Finnigan MAT 251) (Grisart et al., 2022). The data have a resolution of 11 cm and a reported accuracy of 0.07 ‰ and were dated using the Antarctic ice-core chronology (AICC2012) (Veres et al., 2013). For the diffusion length estimate, we define MIS 19 from a depth of 3147 m (748 ka) to 3190 m (802 ka), which marks the beginning of the MIS 18 glacial period and the termination of the MIS 20 glacial period respectively (Pol et al., 2010). Depth ranges for the more recent interglacial periods of MIS 1 (the Holocene), MIS 5, and MIS 9 are given in Table 1, along with depth ranges where data are missing. All section selections are a trade-off between using a long section allowing for a more precise statistical estimation and using a shorter section only containing warm periods.

## 2.2 Method overview

In summary, the diffusion length in deep ice is estimated using a modification of the existing statistical model by representing the power spectrum of the water isotope record before diffusion as a power law. An appropriate power law is estimated from selected water isotope data from shallower sections of the same ice core, where diffusion will not have had time to affect the relevant frequencies, spanning time periods of a similar climate state and length to the deep-ice section.

## 2.3 Spectral analysis

The diffusion length and the properties of the isotopic variability are estimated in the frequency domain. For this, we use the raw periodogram on linearly detrended data with a split cosine bell taper of 10 % as this estimator is unbiased and we are only using it for parametric fits.

The isotope time series are irregular in time due to ice flow thinning, which is incompatible with classical power spectra analysis. Therefore, when working with data on the time domain, the respective records were linearly interpolated to equidistant time points. Although this introduced power loss at high frequencies in the power spectra (Schulz and Stattegger, 1997; Hébert et al., 2021), our analysis did not include these frequencies and so remains largely unbiased. Gaps were only present in the temporal data used for spectral analysis and did not exceed 100 data points (at most corresponding to 863 years). The effect on the spectra was small and did not significantly impact the final result (see Appendix A).

For comparison, we also calculated diffusion lengths using the conventional white-noise model for a running window over the entire EDC ice-core record. Gaps in the data were resolved through linear interpolation, and large or frequent linear interpolated gaps can significantly affect the power spectra. Therefore, the estimate was only performed over windows in which less than 25 % of the *δ*^{18}O values were missing to minimise estimation biases, explaining the gaps seen in Fig. 5b.

## 2.4 Spectral representation of diffusion length

The diffusion of water molecules averages out variability, with rapidly increasing effectiveness as frequency increases. While the lowest frequencies remain mostly unchanged, information at higher frequencies can be greatly affected. This relationship allows for the diffusion length of a water isotope time series to be estimated through analysis of its spectral properties combined with a representative model.

Mathematically, the effect of diffusion on a time series can be represented by a convolution with a Gaussian filter (Johnsen, 1977; Johnsen et al., 2000):

where *z* is the depth and *σ* is the diffusion length in the same units as *z*. By applying the convolution theorem, we can represent the effect of diffusion in the frequency domain with the transfer function:

where *k*=2*π**f*_{z} and *f*_{z} is frequency in the depth (units of m^{−1}). From this the power spectral density (PSD) is

where *P*_{0} is the PSD of the isotope profile before diffusion and *P* is the PSD of the signal after diffusion. In real isotope records the measurement process will add some noise *ϵ* to the signal:

Therefore it is possible to estimate the diffusion length of a water isotope record by taking its power spectrum and fitting Eq. (4) as a model of the effect of diffusion.

The white-noise climate assumption of conventional methods relates to the *P*_{0} term, prescribing it as frequency independent. A constant *P*_{0} means the model contributes all of the power difference between frequencies to diffusion. Also, when water isotope measurements are made of discrete samples, the measurement noise *ϵ* is often considered to be white noise, as each measurement should be independent of the previous measurement. This does not hold true in a continuous-flow analysis system, where the noise can have positive autocorrelation due to mixing within the water lines (Gkinis et al., 2014), which would also influence the shape of the power spectrum. All the data used in this report were discretely measured, so the measurement noise was considered to be white.

Discrete sampling will also introduce a block-averaging effect, further smoothing the data and biasing our diffusion length estimate if not taken into consideration. The magnitude of this smoothing depends on the sampling resolution and can be computed using the method described in Gkinis et al. (2014). For the 11 cm resolution data used here, we compute the corresponding rectangular filter with a width of 11 cm to be equivalent to a diffusion length of ∼ 3.3 cm. Such small-scale smoothing is unlikely to have a significant effect on the heavily diffused deep ice, so the resulting bias is considered negligible.

## 2.5 Fitting methods

There are two common approaches for fitting a spectrum with Eq. (4), a linear regression method applied only to the lower frequencies (Jones et al., 2017) and a non-linear method applied over the entire spectrum (Gkinis et al., 2014; Holme et al., 2018; Kahle et al., 2018).

The linear approach (Jones et al., 2017) focuses only on the lower frequencies, where the power of the climate signal is much greater than the measurement noise term (${P}_{\mathrm{0}}\left(k\right){e}^{-{k}^{\mathrm{2}}{\mathit{\sigma}}^{\mathrm{2}}}\gg \mathit{\u03f5}\left(k\right)$). In this frequency range, Eq. (4) can be approximated as Eq. (3), and by taking the natural log of both sides we get

Assuming ln *P*_{0} is independent of frequency (Johnsen et al., 2000), we can model *P* on a logarithmic scale as a function of *k* with a linear regression which will have a slope proportional to the square of the diffusion length. The cut-off frequency (that is, the upper frequency limit after which the *ϵ* term is no longer negligible) is usually manually defined, which makes this method difficult to generalise and automate. This method relies on the frequency independence of *P*_{0}, as any changes with frequency will create a non-linear relationship between *P* and *k*.

Another possible method of statistically estimating the diffusion length involves modelling ln *P* as a function of *f*_{z} and fitting all parameters of Eq. (4) including the noise (Gkinis et al., 2014). Unlike the linear approach, this fit can be applied over all available frequencies, removing the subjective step of choosing a cut-off frequency. Variants of this approach exist for situations where the measurement noise is not white, such as with data measured using continuous-flow analysis (Gkinis et al., 2014; Kahle et al., 2018).

While previous studies assumed *P*_{0} is equivalent to white noise, the latter fitting method, unlike the linear one, also allows for a frequency dependence of *P*_{0}, and so it was the preferred choice for our new approach going forwards.

## 2.6 Estimating diffusion length when *P*_{0} is not constant

To get a realistic alternative model for *P*_{0} suitable for longer timescales, we use water isotope data less affected by diffusion from the same site. Assuming the statistics of similar climate states over time are comparable, we estimate *P*_{0} empirically from a shallower section of the same ice core spanning a time period of a similar length and climate state. The ice diffusion in this shallow record is still negligible, as it is younger, colder, and less thinned than the deepest ice. Additionally, any significant firn diffusion is on timescales much shorter than we are analysing. For the MIS 19 case, we take water isotope data from interglacial periods further up in the EDC ice core.

To parameterise the variability across a large range of frequencies, we use a power law,

which has been shown to be a good approximation of the spectrum of climate variability (Pelletier, 1998; Huybers and Curry, 2006). Here, *α* and *β* are constants, and *f*_{t} is the frequency in the time domain. We work in time units only for the *P*_{0} fit because similar climate variability during different past-climate states results in comparable power spectra in the time domain but not in the depth domain, as the annual layer thickness will vary over time and ice depth.

For the case of the MIS 19 record from EDC, an appropriate time period to estimate the climate signal before long-term ice diffusion is another, more recent interglacial. Using the same EDC ice core, we selected sections of water isotope data from the MIS 1 (the Holocene), MIS 5, and MIS 9 interglacials, which were retrieved from depths shallow enough to remain unaffected by lower-frequency diffusion. Large data gaps are present over the MIS 7 and MIS 11 interglacials, so a reliable power spectral estimate could not be acquired. Interglacial records from deeper in the ice core (MIS 13/15/17) were not suitable for our analysis as diffusion has attenuated the frequencies over which we are inferring the spectrum of *P*_{0}. Using different interglacial periods to estimate *P*_{0} allows us to evaluate the sensitivity of our diffusion length estimate on our choice of the interglacial.

## 2.7 Bayesian fitting approach

We used a Bayesian approach for all our power spectrum fitting as it has a several advantages over classical methods. Most importantly, rather than having to either set parameters to specific fixed values or leave them free to assume any value, the prior distributions of a Bayesian model allow us to put physically realistic restrictions on the values of parameters while also allowing uncertainty in their true values to be propagated through into uncertainty in the final estimate of diffusion length. Additionally, the Bayesian method allows us to specify a gamma distribution which is the true distribution of the errors in power spectral density estimated by Fourier methods (Bloomfield, 2004). The Gaussian approximation assumed by classical fitting in log space produces a positive bias in the estimated power, which here would affect estimated diffusion lengths.

Using a reduced model without diffusion (Eq. 6 from earlier), we fit our undiffused climate spectra (*P*_{0}) for MIS 1, MIS 5, and MIS 9 with largely uninformative, truncated normal priors (see Appendix B) of *α* ∼ *N*(0.1 ‰^{2} kyr, 1 ‰^{2} kyr) and *β* ∼ *N*(1.5, 1). Here, *N*(*x*,*y*) refers to a normal distribution of mean “*x*” and standard deviation “*y*”. This gives a large range of possible values to find the best fit for the initial MIS 19 climate signal. The posterior means and standard deviations of *α* and *β* from these fits were then used to parameterise informative priors for subsequent fits of the full model including diffusion (Eq. 4) to MIS 19.

For the power of the noise we used a somewhat informative prior, as the uncertainty of the *δ*^{18}O data is well defined from the lab-based measurements, while for the key parameter of interest, the diffusion length, we used a weak, uninformative prior. We set a lower limit of 0 for all parameters in all fits to prevent negative, non-physical values for *α*, *σ*, and the noise *ϵ* and to constrain *β* to power-law red-noise-like behaviour. For the error in spectral estimation we used a gamma distribution, *γ*(*ϕ*, $\mathit{\varphi}/\widehat{P}\left({f}_{z}\right))$, where *ϕ* represents the scale parameter and $\mathit{\varphi}/\widehat{P}\left({f}_{z}\right)$ represents the shape parameter. For this distribution, the magnitude of the error is proportional to the estimated power. We fit an additional final model with a prior for *P*_{0} defined by taking the mean and standard deviation of the *α* and *β* estimates from the individual interglacials (“Mean” row in Table 2).

The models were defined in the Stan language (Carpenter et al., 2017) and fit using the No-U-Turn sampler (Hoffman and Gelman, 2014) with the R package CmdStanR from Gabry and Češnovar (2021). For each fit, we sampled 2000 values from 4 independent chains, with the first 1000 iterations discarded as a warmup. Inspection of Rhat values and trace plots of the posterior model parameters indicated that chains were well mixed for all parameters and all models converged.

## 2.8 Conventionally estimated EDC diffusion lengths

In order to determine the significance of this new method, we first estimate the diffusion length using the conventional model, where the climate variability *P*_{0} is considered constant across all frequencies, which is the second method in Sect. 2.5. However, we still use our Bayesian approach to circumvent the bias of the least-squares estimator. We use the same *σ* and noise priors as before, as well as a weak *α* prior of *α* ∼ *N*(0.5 ‰^{2} m, 0.1 ‰^{2} m). We apply this to a running window of 500 data points across the full high-resolution EDC record.

## 3.1 New *P*_{0} estimate

The full power spectra of the selected interglacial time periods used to constrain *P*_{0} proved difficult for a direct power-law fit. Different timescales can have different power laws (Pelletier, 1998), and this is evident in the power spectra, which seem to level off at frequencies above 3 kyr^{−1} (Fig. 2). Considering also that the lowest few frequencies in the spectra are biased due to tapering, we did not apply the *P*_{0} fit over the whole spectra. We heuristically selected a frequency range between 0.25 kyr${}^{-\mathrm{1}}<{f}_{t}<\mathrm{2.5}$ kyr^{−1} which can be reasonably modelled with a power law for each interglacial. We also devised a technical approach to compute this range, taking care to include all frequencies relevant for the MIS 19 diffusion length estimate, as explained in Appendix C. This computed range closely matched our manually chosen range, so we proceeded with the latter for simplicity. This frequency range also excluded the significantly diffused higher frequencies, with the highest included frequency having lost only 0.2 % of its initial power due to firn diffusion, which makes a negligible difference in the power-law fit.

Using the Bayesian sampling method, the power-law fit is applied to each of the three more recent interglacials within these frequency bounds, with best fits using the mean *α* and *β* values per interglacial shown in Fig. 2. The corresponding average parameter values are shown in Table 2, with errors representing the standard deviation of the Bayesian estimates. We also use the mean of both parameters across the three interglacial fits to represent the average interglacial. After an age-to-depth conversion using the MIS 19 record, this gave Gaussian prior distributions of *α* ∼ *N*(0.0176 ‰^{2} m, 0.0025 ‰^{2} m) and *β* ∼ *N*(1.19, 0.28), with the strength taken as the standard deviation of the interglacial parameter means, accounting for differences between the statistics of the similar climate states.

## 3.2 Improved diffusion length estimate for MIS 19

While prior distributions of *α* and *β* could be derived from data, we still need to choose priors for the measurement noise and the diffusion length itself for our MIS 19 diffusion length estimate. A Gaussian prior distribution of *N*(0.07 ‰, 0.02 ‰) is sufficient for the *ϵ* term as the uncertainty of the *δ*^{18}O data is well defined from the lab-based measurements. The diffusion length *σ* was given a very weak Gaussian prior distribution of *N*(0.4 m, 0.4 m), allowing the model almost complete flexibility to find the best fit. Using the mean *α* and *β* priors, a diffusion length of 31±5 cm (95 % confidence) was estimated. A full table of the results for each interglacial is shown in Table 3.

## 3.3 Comparison of the conventional and new approach

Estimating diffusion length using an incorrectly parameterised best-fit model can significantly bias the result. To demonstrate this, we use the conventional method in Sect. 2.8 and apply it to the theoretical spectrum of a diffused power law with measurement noise, with *α*=1, *β*=1.5, *σ*=0.3, and *ϵ*=0.07. The computed best fit greatly misrepresents the lower frequencies (Fig. 4) and estimates a diffusion length of 0.415, an overestimation of almost 40 %.

Aware of this bias, we applied the conventional method to the MIS 19 spectrum to evaluate the difference between the two estimates, using the Bayesian approach with a very weak *P*_{0} prior of *N*(1 ‰^{2} m, 10 ‰^{2} m) and the same *σ* and noise priors as before. Directly comparing the two fits, our new method demonstrates a clear qualitative improvement (Fig. 5a). The conventional approach estimated a diffusion length of 55±6 cm, over 77 % larger than the new method.

To get an understanding about the evolution of diffusion length at Dome C, we applied the conventional method over a running window across the full Dome C record (Fig. 5b). Gaps in the diffusion length profile arise due to missing high-resolution *δ*^{18}O values, as we excluded diffusion length estimates where more than 25 % of the data in the window were missing. Adding our new result (pink) gives a new impression of the scale of diffusion in the deepest section of Dome C.

## 4.1 Overcoming the estimation bias of the conventional method in deep ice

Applying the conventional method to a running window over EDC, the estimated diffusion lengths rapidly increase from ∼ 20 cm at a depth of ∼ 2800 m to nearly 60 cm at the bottom of the ice core, a depth of almost 3200 m (Fig. 5b). The average 500-data-point window in this depth range covers ∼ 50 kyr, a timescale over which variability cannot be considered uniform across all frequencies. Furthermore, this interval in the record spans four glacial–interglacial cycles, with some windows containing the transition between these states, contributing to more low-frequency variability. This is inconsistent with the common assumption of approximating the original, undiffused record as white noise, strongly biasing conventional estimates of deep-ice diffusion. Adding frequency dependence to *P*_{0} in the estimator model through prescribing power-law scaling helps to eliminate such biases.

## 4.2 Isotopic variability before diffusion and interglacial assumption

The downside of the proposed method is that it needs information about the isotope variability before diffusion (*P*_{0}). In our estimate of this undiffused isotopic signal, we assumed the spectral properties of more recent interglacial periods are comparable to that of the MIS 19 interglacial. It could be argued this is a logical leap, as the interglacials vary in duration and amplitude. However, our diffusion length estimate was insensitive to the interglacial chosen (Fig. 3b, Table 3). For our final result estimated using the mean of the *P*_{0} interglacial fits, the strength of the priors of *α* and *β* are taken as the standard deviation across the individual interglacial fits, corresponding to the uncertainty of our average interglacial. Regardless of whether these precautions represent the full uncertainty in our assumption, it is certainly a more accurate representation of isotopic variability over multi-millennial timescales than white noise. It is usually assumed that isotopes in ice cores are faithful recorders of multi-centennial and longer temperature variability. Therefore, *P*_{0} is dominated by the climate variability at these timescales. This would also allow for the use of other proxy information, for example from marine sediment cores (Shakun et al., 2015), in future studies to test how appropriate it is to use *P*_{0} from more recent interglacials to approximate the variability in MIS 19.

## 4.3 Advantages of the Bayesian approach

Our implementation of a Bayesian fitting approach enables our prior knowledge of the physically realistic values of the parameters to be incorporated into the fit (Gelman et al., 2013). Using a fixed power-law slope to represent our estimation of the past-climate statistical properties would be restrictive, and its accuracy would influence our diffusion length estimate. On the other hand, allowing complete freedom of the *P*_{0} spectrum fit could produce a biased diffusion length estimate as the shape of *P*_{0} may account for some of the loss of variability due to diffusion. The Bayesian method offers a middle ground, where the fit parameters are suggested through priors and the strength of the suggestion depends on the confidence of our knowledge of the parameter. Given our inference for *P*_{0} is an assumption, this is ideal as it allows for some flexibility during the fitting process to account for differences between the MIS 19 interglacial and the more recent interglacials, without overcompensating for diffusive attenuation effects. In other words, the priors enabled our knowledge and confidence in the parameter values to be incorporated into the fitting process, allowing for the uncertainty of the parameters to propagate through to the fit. A further advantage of this method was its capability to correctly treat the residuals of the model as a gamma distribution. Standard non-linear fitting methods assume a Gaussian distribution of the residuals, but this is not the case for power spectral density estimators.

## 4.4 Implications of the new diffusion length estimate for MIS 19

Our estimate of 31±5 cm for the diffusion length over the MIS 19 interglacial period, relative to earlier estimates of 40–60 cm (Pol et al., 2010), has important implications for the interpretation of this part of the EDC ice core. This is because the remaining power of certain frequencies after diffusion is extremely sensitive to the diffusion length. In the EDC record, millennial variability would be reduced to around 1 % of the record before diffusion for our new estimate, compared to $\mathrm{1}/\mathrm{100}$ 000 for the conventionally estimated 55±6 cm. The former offers the opportunity to reconstruct millennial variability using deconvolution, whereas this is not possible for the 55 cm case given the magnitude of measurement noise. The significant decrease in our diffusion length estimate vs. conventional estimates can be explained by the introduction of a “red-noise” behaviour from our *P*_{0} fit (Fig. 4). Red-noise contains more variability in low frequencies than high frequencies and so may account for some of the power difference across frequencies which would otherwise be attributed entirely to diffusion, as in the conventional case. The new estimate is much closer to the diffusion lengths at the bottom of the EDC ice core modelled from physical first principles, which average around 20 cm (Pol et al., 2010; Grisart et al., 2022), reducing the discrepancy between the empirical and analytical results. This new result offers a much more optimistic outlook for recovering millennial-scale climate features from future deep-ice-core projects such as the BE-OIC.

## 4.5 Future outlook

Possible improvements to the precision of our estimate could be made with a higher sampling resolution or data spanning longer time periods. Increasing the resolution would provide some opportunity to reduce the uncertainty, but the effect will be limited as the high frequencies do not contribute to the diffusion length estimate, with the main improvement arising from the effective measurement noise reduction at a given frequency. It would also lessen the smoothing effect of block averaging due to the rectangular sampling scheme, although in the future this could be directly incorporated into the fit. Alternatively, increasing the time span of our selected water isotope record would improve the reliability of our *P*_{0} estimate. However, for MIS 19 the record could not be extended without including data from glacial periods and different climate states with different variability structures, which would complicate our *P*_{0} model.

While the focus of this new method was to improve the diffusion length estimate for deep ice, it should also be considered for application to shallower ice cores. Power-law behaviour is observable in the climate across multi-millennial to sub-annual timescales (Pelletier, 1998), which suggests younger ice cores and firn cores could also benefit from this modified approach. Specifically at core sites with higher accumulation rates and therefore less stratigraphic noise, such as the West Antarctic Ice Sheet (WAIS) Divide ice core (Jones et al., 2018), a power law *P*_{0} may offer a more accurate estimation than white noise.

Perhaps the main disadvantage of this new method compared to conventional estimation methods is the necessary inference of the undiffused isotopic profile from elsewhere in the ice core. It would be much more practical if it were possible to generalise the approach for an entire ice core, acquiring a full diffusion length profile with depth. This is potentially achievable using a different parameterisation for *P*_{0} that is valid across a broader frequency range. Using white noise (reflecting depositional noise) superposed on a piecewise power-law scaling (reflecting internal climate variability vs. multi-millennial orbitally forced variability) might be one option. Another idea is to test if non-diffused parameters such as the dust content from the same section of the ice core can be related to the climate variability and thus *P*_{0}.

We have described a new approach for water isotope diffusion length estimations in deep ice cores, resolving the biased assumption of a white-noise undiffused climate signal. Our method instead implements a power-law slope inferred from water isotope sections of a similar climate state in the shallower parts of the ice core, better representing the climate on millennial scales. Incorporating Bayesian statistics enables us to use priors, chosen based on our knowledge of the parameters and propagating our uncertainties into the fitting procedure. We applied our new method to the MIS 19 water isotope record from the bottom of the EDC ice core, estimating a diffusion length of 31 ± 5 cm, a 23 % reduction to the smallest previously estimated value of 40 cm (Pol et al., 2010). A smaller diffusion length offers a more optimistic outlook for the preservation of millennial climate signals in the oldest-ice projects such as the BE-OIC. Future work could be made to generalise the process for entire ice-core records, possibly through adjusting the Bayesian priors for different climate states or, alternatively, inferring *P*_{0} from non-diffused parameters.

The data for the *P*_{0} fits are gappy and unevenly spaced (time axis). To resolve this for finding the power spectra we create an evenly spaced age axis with the same number of data points and then linearly interpolate the *δ*^{18}O data onto this new axis.

To see the effect this might have on the resulting *P*_{0} fit, for each interglacial (MIS 1/5/9) we simulate a very high-resolution time series with the same parameters as the interglacial. We then bin the data to the unevenly spaced axis of the interglacial record and remove data where the interglacial data have gaps. We compare the spectra and fits with evenly binned data with no gaps, representing the ideal, unbiased spectra. The results for each interglacial are shown below.

We found the linear interpolation has a negligible effect in the frequency range the fit is applied over, as it only introduces error in much higher frequencies. For MIS 1 and MIS 9 there is not much data missing, so the interpolated fit matches almost perfectly with the “true” fit. This is not the case for MIS 5, as enough data are missing that it impacts the accuracy of the fit. However, the resulting fit still falls within the 90 % confidence interval, and given it only contributes partially to a suggestive prior, the MIS 19 fit is not strongly affected (less than 1 cm bias).

Priors for *α* and *β* in the *P*_{0} fits are uninformative (Fig. B1), as are the priors for the diffusion length and noise in the MIS 19 diffusion length estimate fit (bottom plots in Fig. B2). The more informative *α* and *β* priors in the MIS 19 diffusion length estimate (top plots in Fig. B2) are defined from the *P*_{0} fits.

The *P*_{0} power-law fit need only be applied over the frequencies most sensitive to small changes to the diffusion length, as lower frequencies are unaffected by diffusion and higher-frequency variability is attenuated to such an extent that measurement noise dominates. Therefore, an estimate is only necessary over the frequencies in between these two extremes. To find the relevant frequency range, we rearrange Eq. (3) for frequency:

where $\overline{\mathit{\lambda}}$ represents the mean annual layer thickness and is required to convert the diffusion length into the time domain. Here, the term $P\left({f}_{t}\right)/{P}_{\mathrm{0}}\left({f}_{t}\right)$ represents the fraction of power remaining after diffusion at the frequency *f*_{t}. Therefore, we can calculate a frequency range for the fit by defining an upper and lower limit of power loss, assuming an approximate diffusion length. Given we do not know the diffusion length, we take a conservative guess to prevent excluding important frequencies. For the MIS 19 case, we approximate a minimum diffusion length of 15 cm and a maximum diffusion length of 60 cm and use these values to estimate the upper and lower frequency limit respectively. The lower limit is defined as the frequency where the remaining power first drops below $\mathrm{1}/e$, while the upper limit is defined as the frequency after which the remaining power is below 10 % of the measurement noise level. For the latter definition, the power of the measurement noise can be reliably estimated from the flat tail of the diffused power spectrum, while the assumed *P*_{0} is estimated from the highest frequencies in the power spectra (giving a conservative, larger frequency range as this is a larger *P*_{0} than would be expected for the higher frequencies).

The measurement noise was estimated from the MIS 19 power spectrum, which tails off over the highest frequencies at a PSD of $\mathrm{5.4}\times {\mathrm{10}}^{-\mathrm{4}}$ ‰^{2} m, corresponding to a measurement noise of ±0.07 ‰. Then, given $\overline{\mathit{\lambda}}=\mathrm{8.4}\times {\mathrm{10}}^{-\mathrm{4}}$ m yr^{−1} or 0.84 m kyr^{−1} and assuming an upper diffusion length estimate of 60 cm (Pol et al., 2010), a power drop to $\mathrm{1}/e$ occurs at a frequency of *f*_{t}=0.22 kyr^{−1}. Likewise, using the mean annual layer thickness and a lower diffusion length estimate of 15 cm from models (Pol et al., 2010), the power drops below 10 % of the estimated measurement noise at *f*_{t}=2.45 kyr^{−1}.

The data used in this study were previously published and can be found at https://doi.org/10.1594/PANGAEA.939445 (Gkinis et al., 2021). Code used in this study is available on Zenodo at https://doi.org/10.5281/zenodo.12773092 (Shaw, 2024).

TL and FS designed the study. TK and TL contributed to the spectral analysis and signal processing. AMD contributed to the Bayesian estimation. VG provided insights into diffusion and diffusion length estimation. FS performed the analysis and wrote the manuscript. All authors contributed to the interpretation and to the preparation of the final manuscript.

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

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

This article is part of the special issue “Oldest Ice: finding and interpreting climate proxies in ice older than 700 000 years (TC/CP/ESSD inter-journal SI)”. It is not associated with a conference.

This publication was generated in the frame of the DEEPICE project and received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 955750 and the ERC SPACE no. 716092. It contributes to Beyond EPICA, which received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 815384 (Oldest Ice Core). Vasileios Gkinis was supported by the Villum Foundation (“The whisper of ancient air bubbles in polar ice”, grant no. 00028061, and “Unraveling paleo-climate knots with lasers”, grant no. 00022995). Andrew M. Dolman was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project no. 468685498, and SPP 2299, project no. 441832482). This is Beyond EPICA publication number 37.

This research has been supported by the European Research Council through the European Union Horizon 2020 research and innovation programme Marie Skłodowska-Curie Actions (grant no. 955750), the European Research Council through the European Union Horizon 2020 research and innovation programme (grant no. 716092), the Villum Foundation (grant nos. 00028061 and 00022995), Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project no. 468685498, and SPP 2299, project no. 441832482), and European Union’s Horizon 2020 research and innovation programme under grant agreement no. 815384 (Oldest Ice Core).

The article processing charges for this open-access publication were covered by the Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung.

This paper was edited by Carlos Martin and reviewed by Christian Holme and one anonymous referee.

Bloomfield, P.: Fourier analysis of time series: an introduction, John Wiley & Sons, ISBN 0471889482, 2004. a

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M. A., Guo, J., Li, P., and Riddell, A.: Stan: A probabilistic programming language, J. Stat. Softw., 76, 1–32, https://doi.org/10.18637/jss.v076.i01, 2017. a

Casado, M., Münch, T., and Laepple, T.: Climatic information archived in ice cores: impact of intermittency and diffusion on the recorded isotopic signal in Antarctica, Clim. Past, 16, 1581–1598, https://doi.org/10.5194/cp-16-1581-2020, 2020. a

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468, 1964. a

EPICA community members: Eight glacial cycles from an Antarctic ice core, Nature, 429, 623–628, 2004. a

Gabry, J. and Češnovar, R.: cmdstanr: R Interface to “CmdStan”, https://mc-stan.org/cmdstanr, https://discourse.mc-stan.org (last access: 10 July 2024), 2021. a

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Aki, V., and Rubin, D. B.: Bayesian data analysis, Chapman and Hall/CRC, 3rd Edn., ISBN 9780429258411, 2013. a

Gkinis, V., Simonsen, S. B., Buchardt, S. L., White, J., and Vinther, B. M.: Water isotope diffusion rates from the NorthGRIP ice core for the last 16,000 years–Glaciological and paleoclimatic implications, Earth Planet. Sci. Lett., 405, 132–141, 2014. a, b, c, d, e, f, g

Gkinis, V., Dahl-Jensen, D., Steffensen, J. P., Vinther, B. M., Landais, A., Jouzel, J., Masson-Delmotte, V., Cattani, O., Minster, B., Grisart, A., Hörhold, M., Stenni, B., and Selmo, E.: Oxygen-18 isotope ratios from the EPICA Dome C ice core at 11 cm resolution, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.939445, 2021. a

Grisart, A., Casado, M., Gkinis, V., Vinther, B., Naveau, P., Vrac, M., Laepple, T., Minster, B., Prié, F., Stenni, B., Fourré, E., Steen-Larsen, H. C., Jouzel, J., Werner, M., Pol, K., Masson-Delmotte, V., Hoerhold, M., Popp, T., and Landais, A.: Sub-millennial climate variability from high-resolution water isotopes in the EPICA Dome C ice core, Clim. Past, 18, 2289–2301, https://doi.org/10.5194/cp-18-2289-2022, 2022. a, b, c, d

Hébert, R., Rehfeld, K., and Laepple, T.: Comparing estimation techniques for temporal scaling in palaeoclimate time series, Nonlin. Processes Geophys., 28, 311–328, https://doi.org/10.5194/npg-28-311-2021, 2021. a

Hoffman, M. D. and Gelman, A.: The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo, J. Mach. Learn. Res., 15, 1593–1623, 2014. a

Holme, C., Gkinis, V., and Vinther, B. M.: Molecular diffusion of stable water isotopes in polar firn as a proxy for past temperatures, Geochim. Cosmochim. Ac., 225, 128–145, 2018. a, b

Huybers, P. and Curry, W.: Links between annual, Milankovitch and continuum temperature variability, Nature, 441, 329–332, 2006. a

Johnsen, S. J.: Stable isotope profiles compared with temperature 50 profiles in firn with historical temperature records, Isotopes and Impurities in Snow and Ice, 118, 210–219, 1977. a, b, c

Johnsen, S. J., Clausen, H. B., Cuffey, K. M., Hoffmann, G., Schwander, J., and Creyts, T.: Diffusion of stable isotopes in polar firn and ice: the isotope effect in firn diffusion, in: Physics of ice core records, Hokkaido University Press, 121–140, 2000. a, b, c, d, e

Jones, T., Cuffey, K., White, J., Steig, E., Buizert, C., Markle, B., McConnell, J., and Sigl, M.: Water isotope diffusion in the WAIS Divide ice core during the Holocene and last glacial, J. Geophys. Res.-Earth, 122, 290–309, 2017. a, b, c

Jones, T. R., Roberts, W. H., Steig, E. J., Cuffey, K., Markle, B., and White, J.: Southern Hemisphere climate variability forced by Northern Hemisphere ice-sheet topography, Nature, 554, 351–355, 2018. a

Kahle, E. C., Holme, C., Jones, T. R., Gkinis, V., and Steig, E. J.: A generalized approach to estimating diffusion length of stable water isotopes from ice-core data, J. Geophys. Res.-Earth, 123, 2377–2391, 2018. a, b, c

NEEM community members: Eemian interglacial reconstructed from a Greenland folded ice core, Nature, 493, 489–494, 2013. a

NGRIP members: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, 2004. a

Nye, J.: Diffusion of isotopes in the annual layers of ice sheets, J. Glaciol., 44, 467–468, 1998. a

Pelletier, J. D.: The power spectral density of atmospheric temperature from time scales of 10-2 to 106 yr, Earth Planet. Sc. Lett., 158, 157–164, 1998. a, b, c

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

Pol, K., Masson-Delmotte, V., Johnsen, S., Bigler, M., Cattani, O., Durand, G., Falourd, S., Jouzel, J., Minster, B., Parrenin, F., Ritz, C., Steen-Larsen, H. C., and Stenni, B.: New MIS 19 EPICA Dome C high resolution deuterium data: Hints for a problematic preservation of climate variability at sub-millennial scale in the “oldest ice”, Earth Planet. Sc. Lett., 298, 95–103, 2010. a, b, c, d, e, f, g, h, i, j

Ramseier, R. O.: Self-diffusion of tritium in natural and synthetic ice monocrystals, J. Appl. Phys., 38, 2553–2556, 1967. a

Schulz, M. and Stattegger, K.: SPECTRUM: Spectral analysis of unevenly spaced paleoclimatic time series, Comput. Geosci., 23, 929–945, 1997. a

Shakun, J. D., Lea, D. W., Lisiecki, L. E., and Raymo, M. E.: An 800-kyr record of global surface ocean *δ*^{18}O and implications for ice volume-temperature coupling, Earth Planet. Sc. Lett., 426, 58–68, 2015. a

Shaw, F.: deep_ice_diffusion, Zenodo [code], https://doi.org/10.5281/zenodo.12773092, 2024. a

Simonsen, S. B., Johnsen, S. J., Popp, T. J., Vinther, B. M., Gkinis, V., and Steen-Larsen, H. C.: Past surface temperatures at the NorthGRIP drill site from the difference in firn diffusion of water isotopes, Clim. Past, 7, 1327–1335, https://doi.org/10.5194/cp-7-1327-2011, 2011. a

van der Wel, G., Fischer, H., Oerter, H., Meyer, H., and Meijer, H. A. J.: Estimation and calibration of the water isotope differential diffusion length in ice core records, The Cryosphere, 9, 1601–1616, https://doi.org/10.5194/tc-9-1601-2015, 2015. a

Veres, D., Bazin, L., Landais, A., Toyé Mahamadou Kele, H., Lemieux-Dudon, B., Parrenin, F., Martinerie, P., Blayo, E., Blunier, T., Capron, E., Chappellaz, J., Rasmussen, S. O., Severi, M., Svensson, A., Vinther, B., and Wolff, E. W.: The Antarctic ice core chronology (AICC2012): an optimized multi-parameter and multi-site dating approach for the last 120 thousand years, Clim. Past, 9, 1733–1748, https://doi.org/10.5194/cp-9-1733-2013, 2013. a

Whillans, I. and Grootes, P.: Isotopic diffusion in cold snow and firn, J. Geophys. Res.-Atmos., 90, 3910–3918, 1985. a

- Abstract
- Introduction
- Data and methods
- Results
- Discussion
- Conclusions
- Appendix A: Interpolation effect
- Appendix B: Prior sensitivity
- Appendix C: Defining the frequency range of the fit
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Special issue statement
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Data and methods
- Results
- Discussion
- Conclusions
- Appendix A: Interpolation effect
- Appendix B: Prior sensitivity
- Appendix C: Defining the frequency range of the fit
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Special issue statement
- Acknowledgements
- Financial support
- Review statement
- References