The impact of measurement precision on the resolvable resolution of ice core water isotope reconstructions

Stable water isotopes in ice cores serve as a valuable proxy for the climate of the past hundreds of thousands of years. Over time, water isotope diffusion causes significant attenuation of the isotopic signal, exacerbated in deep ice due to extreme layer thinning and increased temperatures from geothermal heat flux. This damping affects higher frequencies to a greater extent, erasing information on the shortest timescales. It is possible to restore some of the attenuated variability through deconvolution, a method which reverses the effect of diffusion. However, since the measured isotopic signal always contains noise from the measurement process, deconvolution inevitably amplifies this measurement noise along with the isotopic signal. Thus the effectiveness of deconvolution depends on the precision of the measurements, with noisier data limiting the ability to restore otherwise resolvable frequencies. Here, we quantify the upper frequency limit introduced by the magnitude of the measurement noise analytically for different climate states, and offer a numerical example using the Beyond EPICA Oldest Ice Core (BE-OIC). We also demonstrate the qualitative significance of measurement noise on simulated Antarctic isotopic profiles. The general resolution improvement for firn or upper ice records is on the order of 1.5 times for a 10-fold reduction in measurement noise. Similarly, throughout the BE-OIC, we find the deconvolution of δ¹⁸O records with measurement error of 0.1 ‰ contributes a 1.5 times increase in the maximum resolvable frequency, which rises to a factor of 2 improvement after reducing the measurement noise to 0.01 ‰. While progress is continuously being made towards improving precision of stable isotope measurements, further improvements using longer integration times should be considered when analysing limited and precious deep ice in order to obtain the most faithful climate reconstructions possible.

1 Introduction

Ice cores offer a valuable archive of climate proxies, including stable water isotopes (δ¹⁸O, δD, δ¹⁷O) which are known to relate to air temperature (Dansgaard, 1964). Thus, the drilling of deep ice cores enables the retrieval of long, continuous temperature proxy profiles hundreds of thousand of years in the past (EPICA community members, 2004; NEEM community members, 2013; NGRIP members, 2004; Petit et al., 1999). However, over time these isotopic records are subject to alteration, primarily due to the constant motion of the archived water molecules, a process known as diffusion, and quantified by the diffusion length (Johnsen, 1977; Johnsen et al., 2000). This has the undesirable effect of reducing high frequency variability, smoothing the signal. Diffusion occurs both in the porous firn through vapour-snow exchange (Johnsen, 1977; Whillans and Grootes, 1985), and in the solid ice matrix through ice diffusivity (Ramseier, 1967; Nye, 1998). While the former process ceases below the pore close-off depth, the slower ice diffusion is amplified in the warmer, deeper parts. This becomes increasingly problematic for deep ice cores, as due to the age, temperature and extreme thinning of the ice, the smoothing effect of the ice diffusion is strongest at the bottom of the ice sheet, where the oldest climate information is stored. Methods to retrieve this original signal are thus crucial for the use of water isotopes in deep ice cores for paleoclimate reconstructions.

A commonly employed technique for signal recovery is deconvolution, a method capable of reversing the effects of a linear impulse response on a desired signal. When applied to water isotope diffusion, deconvolution is often referred to as back-diffusion, and the impulse response takes the form of a Gaussian (Johnsen et al., 2000). If the variability of the signal before diffusion and the diffusion length are known, the attenuated frequencies can be amplified, restoring their original variability. However, noise introduced during instrumental measurements complicates this process. For strongly diffused ice, the variability of the measurement noise will mask the climate information at the highest, most-diffused frequencies, so the full restoration of these frequencies through deconvolution is not possible without greatly inflating the measurement noise. Therefore, a compromise must be found between recovering the climate signal and avoiding amplification of the measurement noise. Wiener deconvolution offers an optimal solution, by deriving the amplification factor for each frequency which minimises the mean square error between the original and deconvolved records (Wiener, 1949; Johnsen, 1977). Since the signal-to-noise ratio (hereafter SNR) of a given frequency is not improved, it is optimal to dampen the variability of frequencies where measurement noise dominates the signal. This establishes an effective limit on the highest resolvable frequency through Wiener deconvolution which depends on the SNR and the diffusion length. Consequently, reducing measurement noise enables the recovery of higher frequencies that were previously irretrievable. Figure 1 illustrates a simplified version of this concept through idealised power spectra. The variability of the diffused signal (black) decreases with increasing frequency, until it falls below the variability of the measurement noise (dashed lines). The corresponding frequencies at which SNR = 1 is marked by the vertical lines. It is clear to see that decreasing the measurement noise by some amount ΔP_n (orange to blue) increases the frequency at which SNR = 1 by Δf.

Figure 1 Conceptual diagram illustrating how reducing the measurement noise will enable the recovery of higher frequencies. Shown is an idealised power spectrum of a diffused signal (black), and measurement noise (dashed horizontal lines). Reducing the measurement noise (orange to blue) increases the frequency at which SNR = 1 (solid vertical lines).

We are particularly interested in the effect on the Beyond EPICA – Oldest Ice Core (hereafter BE-OIC) project, due to its aim in acquiring a climate record dating back 1.5 million years, unprecedented in ice core records. In these deepest sections, millennial variability will likely be attenuated beyond recovery, and even the 41 000 year (41 kyr) glacial-interglacial cycles before the Mid-Pleistocene Transition could be dampened. It is therefore crucial to be aware of the limitations introduced by our measurement process so attempts can be made to recover the sought-after temporal resolution.

In this paper, we determine the relationship between measurement precision and highest resolvable frequency. We begin by addressing the problem analytically, deriving general equations that describe the highest resolvable frequency and its dependencies across various scenarios. We then apply these equations to the Beyond EPICA – Oldest Ice Core project, estimating the measurement precision required to achieve specific temporal resolutions, particularly in the deepest and oldest sections of the core. We simulate surrogate BE-OIC records with differing levels of measurement noise, and apply Wiener deconvolution to visualise the qualitative improvement offered by very precise data. Finally, we discuss potential techniques for achieving such high precision through long-integration time measurements.

2 Data and Methods

2.1 Data

For our application to the BE-OIC, we use the age-depth model from Chung et al. (2023), based on a 1D numerical model, radargrams of the Dome C area and the EPICA Dome C (EDC) ice core age profile. This model is used for the transformation of frequencies between the time domain and the depth domain in the following analysis. We also use the discrete 11 cm resolution δ¹⁸O data from the EDC ice core (EPICA community members, 2004; Grisart et al., 2022) to estimate the isotopic variability before diffusion in the Dome C area.

2.2 Power spectrum describing the isotopic variability at Dome C

Over deep ice cores such as the BE-OIC we are working with vastly different frequencies at different depths, spanning decadal to multi-millennial timescales. Before diffusion, water isotope records from ice cores are commonly approximated as white-noise on sub-annual to centennial timescales (Johnsen et al., 2000; Gkinis et al., 2014; Jones et al., 2017) justified by precipitation intermittency and depositional processes obfuscating the temperature signal (Fisher et al., 1985; Johnsen et al., 1997; Münch and Laepple, 2018; Casado et al., 2020). However, these processes have a lesser effect on longer timescales, which tend to display greater low-frequency variability due to a relatively stronger climate signal. Instead, the power spectra of records on longer timescales often closely resemble a power-law (Pelletier, 1998; Huybers and Curry, 2006; Laepple and Huybers, 2014; Shaw et al., 2024). A simple model across all frequencies for the power spectrum before diffusion (hereafter P₀) of isotopic records from the Dome C region is therefore a combination of the models of both regimes,

$$\begin{array}{}\text{(1)}& P_{\mathrm{0}}\left(f\right)=\mathit{\alpha }{f}^{-\mathit{\beta }}+\mathit{ϵ},\end{array}$$

where αf^−β is the low frequency power-law component characterised by numerical constants α and β, and ϵ is a constant representing the high frequency white-noise component.

To best estimate P₀(f), we use δ¹⁸O data from the current oldest ice core to have been retrieved, the EDC ice core. We bin the data to equidistant spacing in time which is required to compute the power spectrum. Since we are interested in all timescales, we create multiple timeseries using different bin sizes. We then estimate the power spectrum of each of these timeseries using Thomson's multitaper method (Percival and Walden, 1993) and compute the mean across all individual spectra (see Appendix A). Using Eq. (1) we find a best fit in logarithmic space, to attribute proportional weighting to all frequencies (Fig. 2). Frequencies greater than 10 kyr⁻¹ were excluded due to the influence of the firn diffusion (see Appendix A), while frequencies less than 0.01 kyr⁻¹ were also excluded due to the known bias in the spectral estimation method (Thomson, 1982).

Figure 2 P₀(f) fit (red) over the grey shaded region of a composite power spectrum of the EDC water isotope record (black). By excluding high frequencies attenuated by firn diffusion, this fit serves as a good approximation of the initial isotopic variability of the BE-OIC.

3 Theory

3.1 Diffusion and Wiener deconvolution

The effect of diffusion on some isotopic profile x(z) is mathematically represented as a convolution of the profile with a Gaussian filter h(z),

$$\begin{array}{}\text{(2)}& h\left(z\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{\frac{-z^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}.\end{array}$$

Here, σ represents the diffusion length, in units equal to the depth z. The final measured record y(z) including some measurement noise term n(z) is given by,

$$\begin{array}{}\text{(3)}& y\left(z\right)=x\left(z\right)*h\left(z\right)+n\left(z\right),\end{array}$$

where * is the convolution operator. Due to the convolution theorem (Cochran et al., 1967), it is easier to work in the Fourier domain,

$$\begin{array}{}\text{(4)}& Y\left(f\right)=X\left(f\right)H\left(f\right)+N,\end{array}$$

where capitalised variables represent the corresponding Fourier transforms. Taking the absolute square of the Fourier transforms gives us the power spectra,

$$\begin{array}{}\text{(5)}& P\left(f\right)=P_{\mathrm{0}}\left(f\right)e^{(-\mathrm{2}\mathit{\pi }f\mathit{\sigma }{)}^{\mathrm{2}}}+P_{\mathrm{n}},\end{array}$$

where P(f) and P₀(f) are the power spectrum of y(z) and x(z) respectively, f represents frequency in the reciprocal units of z and P_n(f) is the power spectrum of the measurement noise. To recover the best estimate, $\widehat{x}\left(z\right)$, of the original isotopic signal x(z), we want to find some function g(z) that can reverse the effects of diffusion without inflating the measurement noise. After convolving our measured record y(z) with this function,

$$\begin{array}{}\text{(6)}& \widehat{x}\left(z\right)=g\left(z\right)*y\left(z\right),\end{array}$$

we want the mean square error between x(z) and our estimate $\widehat{x}\left(z\right)$ to be minimised. The function g(z) that provides such an estimate is the Wiener filter (Wiener, 1949), which has the Fourier transform,

$$\begin{array}{}\text{(7)}& G\left(f\right)={\displaystyle \frac{P_{\mathrm{0}}\left(f\right)H^{*}\left(f\right)}{P_{\mathrm{0}}\left(f\right)\left|H\right(f){|}^{\mathrm{2}}+P_{\mathrm{n}}(f)}}.\end{array}$$

Using the convolution theorem and the inverse Fourier transform ℱ⁻¹, $\widehat{x}\left(z\right)$ can be calculated from

$$\begin{array}{}\text{(8)}& \widehat{x}\left(z\right)={\mathcal{F}}^{-\mathrm{1}}\left(G\right(f\left)Y\right(f\left)\right).\end{array}$$

From Eq. (7), if P_n(f)=0 then the Wiener filter simplifies to $G\left(f\right)=\frac{\mathrm{1}}{H\left(f\right)}$, and completely reverses the smoothing effect, restoring the exact original record. In reality, the noise term reduces the amplification of the diffused frequencies, with a greater reduction the more the noise dominates the diffused signal. Wiener deconvolution balances the recovery of the desired signal and the suppression of over-amplifying undesired noise, such that the mean square error is minimised.

3.2 Highest resolvable frequency

To determine what effect measurement precision has on the recoverable resolution, we first need to define the highest resolvable frequency f_max. Intuitively, we can define a frequency as resolvable if its amplitude is above some fraction of its initial, undiffused amplitude. We call this the relative amplitude, A. For a diffused isotopic record, the relative amplitude A_diff is simply the diffusion transfer function,

$$\begin{array}{}\text{(9)}& A_{\mathrm{diff}}=H\left(f\right)=e^{(-\mathrm{2}\mathit{\pi }f\mathit{\sigma }{)}^{\mathrm{2}}}\end{array}$$

while for a deconvolved isotopic record, the relative amplitude A_decon will also depend on the Wiener filter (see Fig. 3),

$$\begin{array}{}\text{(10)}& A_{\mathrm{decon}}=H\left(f\right)G\left(f\right)={\displaystyle \frac{P_{\mathrm{0}}\left(f\right)\left|H\right(f){|}^{\mathrm{2}}}{P_{\mathrm{0}}\left(f\right)\left|H\right(f){|}^{\mathrm{2}}+P_{\mathrm{n}}(f)}}\end{array}$$

The corresponding highest resolvable frequencies (f_{max_diff}, f_{max_decon}) are thus the frequencies at which the relative amplitude drops to some minimum value (Fig. 3),

$$\begin{array}{}\text{(11)}& {\displaystyle }{A}_{\mathrm{diff}\mathrm{\_}\min }=H\left(f_{\max \mathrm{\_}\mathrm{diff}}\right)\text{(12)}& {\displaystyle }{A}_{\mathrm{decon}\mathrm{\_}\min }=H\left(f_{\max \mathrm{\_}\mathrm{decon}}\right)G\left(f_{\max \mathrm{\_}\mathrm{decon}}\right)\end{array}$$

Equation (11) can be rearranged for f_{max_diff},

$$\begin{array}{}\text{(13)}& f_{\max \mathrm{\_}\mathrm{diff}}={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }\mathit{\sigma }}}\sqrt{-{\displaystyle \frac{\ln A_{\mathrm{diff}\mathrm{\_}\min }}{\mathrm{2}}}}\end{array}$$

with dependency only on diffusion length and our chosen minimum relative amplitude.

For the deconvolved record, there are additional dependencies on the initial signal P₀(f) and the measurement noise P_n(f). In this study, we take P_n to be constant with no frequency dependence, a reasonable assumption for discrete measurements as the noise between measurements should be uncorrelated. For the initial signal, we look at two cases. Firstly, we consider firn records or ice cores retrieved from areas with high stratigraphic noise, for which we can also approximate P₀ as white-noise with no frequency dependence (Johnsen et al., 2000; Shaw et al., 2024). Given this assumption, we can rearrange Eq. (12) for the maximum resolvable frequency,

$$\begin{array}{}\text{(14)}& f_{\max \mathrm{\_}\mathrm{decon}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathit{\sigma }}}\sqrt{\ln {\displaystyle \frac{P_{\mathrm{0}}\left(\frac{\mathrm{1}}{A_{\mathrm{decon}\mathrm{\_}\min }}-\mathrm{1}\right)}{P_{\mathrm{n}}}}}\end{array}$$

Reducing the power spectrum of the noise by a factor η results in the same A_{decon_min} at some new frequency $f_{\max \mathrm{\_}\mathrm{decon}}^{\prime }$

$$\begin{array}{}\text{(15)}& f_{\max \mathrm{\_}\mathrm{decon}}^{\prime }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathit{\sigma }}}\sqrt{\ln {\displaystyle \frac{\mathit{\eta }{P}_{\mathrm{0}}\left(\frac{\mathrm{1}}{A_{\mathrm{decon}\mathrm{\_}\min }}-\mathrm{1}\right)}{P_{\mathrm{n}}}}}\end{array}$$

The relative gain in frequency, f_gain, from this noise reduction η is given by the ratio of these two frequencies,

$$\begin{array}{}\text{(16)}& f_{\mathrm{gain}}={\displaystyle \frac{f_{\max \mathrm{\_}\mathrm{decon}}^{\prime }}{f_{\max \mathrm{\_}\mathrm{decon}}}}=\sqrt{\mathrm{1}+{\displaystyle \frac{\ln \mathit{\eta }}{\ln \frac{P_{\mathrm{0}}\left(\frac{\mathrm{1}}{A_{\mathrm{decon}\mathrm{\_}\min }}-\mathrm{1}\right)}{P_{\mathrm{n}}}}}}\end{array}$$

For the long time scales of deep ice cores, the white-noise P₀(f) assumption is no longer applicable (Shaw et al., 2024). In this case, the climate signal will have some frequency dependence and thus rearranging for frequency is not analytically solvable. However, we can still rearrange Eq. (12) to give the measurement noise required for a specific frequency,

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& P_{\mathrm{n}}(f_{\max \mathrm{\_}\mathrm{decon}},A_{\mathrm{decon}\mathrm{\_}\min })\\ & =P_{\mathrm{0}}\left(f_{\max \mathrm{\_}\mathrm{decon}}\right)e^{-(\mathrm{2}\mathit{\pi }{f}_{\max \mathrm{\_}\mathrm{decon}}\mathit{\sigma }{)}^{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{1}}{A_{\mathrm{decon}\mathrm{\_}\min }}}-\mathrm{1}\right)\end{array}\end{array}$$

Notably, if we aim to restore 50 % of the amplitude of a given frequency, i.e. A_{decon_min}=0.5, the noise required would be

$$\begin{array}{}\text{(18)}& P_{\mathrm{n}}=P_{\mathrm{0}}\left(f_{\max \mathrm{\_}\mathrm{decon}}\right)e^{-(\mathrm{2}\mathit{\pi }{f}_{\max \mathrm{\_}\mathrm{decon}}\mathit{\sigma }{)}^{\mathrm{2}}}\end{array}$$

which is also the noise required for a SNR of 1 (first RHS term = second RHS term from Eq. 5). In the following, we take $A_{\mathrm{diff}\mathrm{\_}\min }=A_{\mathrm{decon}\mathrm{\_}\min }=\mathrm{0.25}$ as the lowest amplitude which we consider recovered. This means at our highest resolvable frequency there is a 75 % loss of the initial signal amplitude.

Figure 3 Concept of diffusion and Wiener deconvolution in the frequency domain, and the definition of f_max. (a) The power spectra of the isotopic signal before (grey) and after (green) diffusion, and after deconvolution (blue). The frequency at which the diffused signal intersects the measurement noise (i.e. SNR = 1) is shown by the green dot. (b) The Gaussian diffusion filter (orange) and Wiener filter (purple) in the Fourier domain. The highest resolvable frequency of the diffused record f_{max_diff} (dashed orange) and deconvolved record f_{max_decon} (dashed blue) occur where the corresponding relative amplitudes A_diff and A_decon intersect the A_min=0.25 line (dashed black).

4 Results

4.1 Generalised white-noise case

We first investigate how a reduction in measurement noise influences the maximum resolvable frequency for the case where P₀(f) is approximated as white-noise, and thus in this section will be notated as P₀. For this, we adopt a value of A_{min_decon}=0.25 and three different $\frac{P_{\mathrm{0}}}{P_{\mathrm{n}}}$ ratios (hereafter SNR₀) of 5, 10 and 20, and analyse the gain in resolvable frequency, f_gain, following Eq. (16) (Fig. 4). We see a non-linear relationship between f_gain and the measurement noise reduction, with a sharper gain for small reductions which diminishes the further the noise is decreased. Similarly, SNR₀ affects the gain, with noisier records benefiting the most. As a test of our analytical solution (Eq. 16), we also compute f_gain numerically and obtain results matching the analytical method. This validates our numerical approach and establishes confidence in using it for the more complex scenario where P₀ has some frequency dependence, and therefore no analytical solution for frequency is possible.

Figure 4 Relative gain in frequency f_gain as a function of measurement noise reduction factor ($\sqrt{\mathit{\eta }}$, thus measured in standard deviation), estimated for the white noise case assuming three different initial signal-to-noise ratios (SNR₀, i.e. the ratio between the initial signal P₀ and the measurement noise level P_n). The analytical results (dashed lines) match the corresponding numerical results (solid black lines), corroborating the numerical method.

4.2 Application to Beyond EPICA – Oldest Ice Core

Deeper ice cores spanning longer timescales store isotopic data that are more complex than white-noise. In these cases, too many parameters are involved to generalise the measurement noise and highest resolvable frequency relationship, so we instead find estimates through direct application to a specific ice core. Taking the example of the BE-OIC, we calculate absolute f_{max_diff} values throughout the record using our diffusion length profile, transformed into time units using the age-depth model, and A_{diff_min}=0.25 (dotted lines in Fig. 5a). Similarly, we also calculate f_{max_decon} values at the same ages for various measurement noise levels, giving a measurement noise and f_{max_decon} relationship predicted for different ages of ice in the BE-OIC (solid lines in Fig. 5a). For a closer inspection, we rescale the x axis by dividing all f_{max_decon} values by the f_{max_diff} at the corresponding ages. This presents the relative improvement offered by Wiener deconvolution for different measurement noise levels (Fig. 5b). We see the resolvable timescales for a standard 0.1 ‰ of δ¹⁸O measurement noise vary from 5 years (5 kyr old ice) to longer than 10 000 years (1.5 Myr old ice). Throughout the ice core, Wiener deconvolution increases f_max by a factor of 1.0 to 1.6. Further improvement is seen if the measurement noise is reduced to 0.01 ‰, with frequencies resolvable up to twice as high as the non-deconvolved record.

Figure 5 (a) Maximum permitted measurement noise to recover the highest resolvable frequency f_{max_decon} for selected ages of the BE-OIC after deconvolution, using A_min=0.25. Also shown are the highest resolvable frequencies with no deconvolution f_{max_diff} (vertical dotted lines). For improved visibility, the points at which the f_{max_diff} lines intersect the corresponding f_{max_decon} lines are marked by dots. The horizontal dashed lines represent a typical measurement noise of 0.1 ‰ (orange) and a 10 times reduction (blue). (b) Same as (a) but showing the frequency gain relative to the corresponding f_{max_diff} (no deconvolution).

4.3 Surrogate ice core record examples

While the highest resolvable frequency is a useful quantification parameter, visualising the direct improvement to a deconvolved depth profile can provide a more intuitive understanding of the significance. We generate stochastic depth series by first creating white-noise, scaling its Fourier transform to match the target power spectrum, and then applying the inverse Fourier transform to obtain the depth-domain signal. As a shallow ice case, we use the white-noise section of our P₀ model in Fig. 2 as our target power spectrum, i.e. P₀=ϵ from Eq. (1). We simulate this shallow ice depth profile at a 5 cm resolution, and smooth it with a diffusion length of σ=0.08 m, taken at an age of 1 kya from the modelled profile (Fig. B1), following Eq. (5). Duplicating the record, we add normally distributed random variations with a standard deviation of 0.1 ‰ and 0.01 ‰, representing measurement noise. We compare the original climate signal input to these two surrogate isotope records (Fig. 6a) and their Wiener deconvolutions (Fig. 6b). This process is repeated for a deep ice case, where the power spectrum is rescaled to the power-law regime of Fig. 2 (i.e. $P_{\mathrm{0}}\left(f\right)=\mathit{\alpha }{f}^{-\mathit{\beta }}$) and the diffusion length is σ=0.21 m, taken at an age of 1.2 Mya, (Fig. 7a, b). In both the shallow and deep ice cases, the surrogate records look similar despite differing noise levels. However, the improvement in the deconvolution is more apparent, with the less noisy record better representing the true, original signal. Some features occurring on timescales only recoverable with the lower measurement noise are restored much more accurately for the less noisy records (e.g. 56–57.5 m in Fig. 6b, 2489–2491 m in Fig. 7b). In these examples, the residuals between the original and deconvolved signals (Figs. 6c, 7c) show a mean reduction of 33 % and 49 % respectively for the less noisy record.

Figure 6 Example of the effect of the measurement noise level on the deconvolution in the depth domain for shallow ice. (a) The true, undiffused record (grey) and two diffused records (σ=0.08 m) with measurement noise standard deviations of 0.1 ‰ (orange) and 0.01 ‰ (blue). (b) Deconvolutions of the two diffused records in (a). (c) Residuals between the original and deconvolved records.

Figure 7 Example of the effect of the measurement noise level on the deconvolution in the depth domain for deep ice. (a) The true, undiffused record (grey) and two diffused records (σ=0.21 m) with measurement noise standard deviations of 0.1 ‰ (orange) and 0.01 ‰ (blue). (b) Deconvolutions of the two diffused records in (a). (c) Residuals between the original and deconvolved records.

5 Discussion

This study explores the role of measurement precision in restoring high frequency variability of diffused water isotope signals from ice cores. We achieve this by studying how the measurement noise level affects the maximum recoverable frequency after deconvolution. For this, we use the criterion where Wiener deconvolution is able to restore to at least 25 % of a given frequency's amplitude before diffusion. For records spanning shorter time periods, where the undiffused isotopic record can be approximated as white-noise, we derive a closed form expression for the dependency on the measurement noise. We find a non-linear relationship between measurement noise and the highest resolvable frequency f_{max_decon}. This relationship also depends on the ratio between the power spectra of the initial signal and the measurement noise. Over the longer timescales covered in deep ice, where the white noise assumption for the initial signal does not hold, we instead preform numerical calculations through direct application to a specific ice core. We use the BE-OIC as an example, estimating absolute frequency limits for different measurement noise levels. Additionally, we show the resolution improvement offered by Wiener deconvolution relative to the diffused, non-deconvolved case, finding similar gains from reducing measurement noise as the generalised white-noise case throughout the ice core.

5.1 Dependency on the relative amplitude

Our paper describes a method of quantifying the resolution limit f_max imposed by measurement noise. For this aim, a definition is needed for what amplitude loss is still acceptable to the researcher. Here, we choose a relative amplitude A_min=0.25, i.e. 25 % of the initial amplitude remaining, since we are interested in the restoration of the climate signal. One may require restoration closer to the initial amplitude, for which a higher minimum relative amplitude can be chosen, say A_min=0.9. Conversely, one could select a lower relative amplitude for a more optimistic f_max, although the extra frequencies will have minimal restoration and contain mostly measurement noise. For different relative amplitudes, the response of the corresponding f_max to changes in measurement noise varies (see Appendix C). Predictably, choosing A_min=0.05 slightly increases f_max values, but also reduces the relative frequency gain to approximately 1.4× for 0.01 ‰ noise, and even results in a smaller f_{max_decon} than f_{max_diff} for a noise 0.1 ‰ (Fig. C1). In contrast, setting A_min=0.9 demonstrates a much larger frequency gain after Wiener deconvolution, between 1.5× and 4× f_{max_diff} for 0.1 ‰ noise and from 4.5× to 6× for noise of 0.01 ‰ (Fig. C2). This shows that the benefit of reducing measurement noise is greater when striving for more accurate restorations of a signal's original amplitude.

5.2 Alternative deconvolution methods

Our results depend on the choice of deconvolution method, for which we have used the Wiener filter. One might ask if there are other deconvolution methods that can recover higher frequencies to a greater extent. Any linear response function will damp or amplify signal and noise equally, and thus not change the SNR. In principle, one could amplify higher frequencies than the Wiener filter allows, and thus get a deconvolved record with more high frequency variability that will be much noisier. Alternatively, a more conservative reconstruction can be made to keep the average SNR down, but in turn more high frequencies will be lost. The Wiener deconvolution method was selected as it is optimal in terms of minimising the root mean square error, and also due to its prominence in previous water isotope research (Johnsen, 1977; Johnsen et al., 2000; Gkinis et al., 2011) and other paleoclimatic proxy reconstructions (Liu et al., 2021).

5.3 Assumptions and limitations

Due to the complexity of the many parameters involved in this study, some assumptions were made for simplicity.

Throughout this study, measurement noise was approximated as white noise, which allowed for greater generalisation and simplified certain calculations. This assumption is appropriate for discretely sampled data such as the Dome C water isotope record, because consecutive measurements are statistically independent. However, when working with isotopic data from continuous flow analysis, additional complications from instrumental memory effects should be considered (Kahle et al., 2018).

Given the limited data on the physical properties of the borehole at the time of writing, the diffusion length model for the oldest ice core incorporates some parameter choices, which introduce uncertainty into the absolute values of f_max across depth and age. While these uncertainties affect the specific estimates of f_max, the relative gain from reducing measurement noise remains largely unaffected. Consequently, the impact on the core objective of this manuscript, exploring the relationship between frequency resolution and measurement noise, is small.

While our results are presented only for the oldest ice core case study, the method introduced in this paper can be applied to other specific ice cores to determine the level of measurement noise required to resolve the desired temporal resolution. In order to do so, some representation of P₀(f) must be available, such as in the oldest ice study where P₀(f) was estimated from the Dome C water isotope record. However, the convenience of an equivalent ice core from which P₀(f) can be estimated before measuring the new core will not always be available, especially for deep ice cores, and as such the P₀(f) estimation could be more challenging in future applications of this method.

5.4 Implications for the BE-OIC project

Our study shows that the relative frequency gain through Wiener deconvolution for a measurement noise of 0.1 ‰ is between 1.0 and 1.6, and increases to 1.6–2.0 for 0.01 ‰ of noise. The surrogate data examples illustrate how this difference offered by the noise reduction could be crucial for more accurate representation of deep ice records especially. If features of interest are expected to have high variability over timescales which benefit most from reduced measurement noise, then it is necessary to measure with enough precision to fully capture the important features. For example, one aim of the Beyond EPICA project is to recover glacial-interglacial patterns before the Mid-Pleistocene Transition, which have a periodicity of 41 kyr. Our estimate in Fig. 5 shows a measurement noise of 0.1 ‰ has an f_max of <0.1 kyr⁻¹ for 1.5 Myr old ice. To recover features like the asymmetry of glacial-interglacial cycles, or overshoots after deglaciation, information on frequencies >0.1 kyr⁻¹ could be necessary, and thus reducing the noise may prove very beneficial.

It is also worth comparing our maximum resolvable frequency of 0.1 kyr⁻¹ for the oldest ice with previous predictions. Fischer et al. (2013) considers ice with age densities >20 kyr m⁻¹ to no longer contain any measurable paleoclimate signal, justified as the maximum density after which greenhouse gases have potentially been homogenised by diffusion. Our limit defined and discussed in this paper concerns the diffusion and recovery of stable water isotopes, for which we predict a more optimistic threshold of up to 10 kyr cycles in the deepest ice. We therefore argue that while there are clear limitations on the gain achieved by reducing measurement noise, given the value and difficulty of retrieving such an old proxy record, acquiring high measurement precision is still of great importance to achieve the most paleoclimatic information possible from the oldest ice core, especially for the previously unextracted and precious >1 Myr old ice.

5.5 Future directions – How to decrease measurement noise

Reducing the δ¹⁸O measurement noise to the 0.01 ‰ suggested in this paper, or even further, is already achievable through multiple advances in the area. Firstly, the L2140-i Cavity Ring-Down Spectroscopy water isotope analyser from PICARRO now offers a δ¹⁸O measurement noise level of 0.025 ‰. Additionally, PICARRO claim using the high precision mode, with long integration time measurements of 54 min per sample, the measurement noise can be reduced to 0.01 ‰. Indeed, Steen-Larsen and Zannoni (2024) demonstrates even better measurement precision is possible, with Allan deviations as low as 0.005 ‰ for δ¹⁸O after measurement times on the order of 2–3 h per sample. While this is undoubtedly a long time to measure each sample, it may still be considered appropriate for the deepest ice in the BE-OIC. Since the highly sought-after and previously unrecovered time period of 800–1500 kya spans only the bottom ∼75 m (Chung et al., 2023), there are very few samples which would require such precision, and therefore the time commitment is less significant.

Another factor to consider is how the sampling resolution of discrete records affects the highest resolvable frequency. In strongly diffused deep ice, the additional high frequencies acquired through finer sampling are irrelevant for climate reconstructions as the measurement noise will dominate the signal. However, the data can be averaged in time to the desired frequency. If the measurement noise is uncorrelated, this reduces the measurement noise by a factor of $\sqrt{N}$ where N is the number of data points per bin. This in turn has an effect on the highest frequency for which the climate signal is recoverable. However, given the greatly increased time commitment necessary for cutting and measuring at a reduced discrete sampling size, an improvement of $\mathrm{1}/\sqrt{N}$ is not beneficial enough to advocate over alternative, more direct methods of reducing measurement noise. Alternatively, taking duplicate cores also allows any white measurement noise to be reduced by a factor of $\sqrt{N}$, where N is the number of cores. This in turn improves the SNR and enables higher frequencies to be recovered using Wiener deconvolution, but can be a difficult task for larger, deep ice core projects. Fundamentally, given the extra effort required, the ideal precision comes down to a trade-off between the improvement of the deconvolved record and the additional investment needed to obtain it.

6 Conclusions

In conclusion, we define the highest resolvable frequency f_max of water isotope ice core records, both before and after deconvolution, and analyse the limitation introduced by measurement noise for the deconvolved records. We find that enhancing measurement precision leads to a non-linear increase in f_max, with a similar gain across various isotopic signals expected in Antarctic ice cores. For records spanning up to centennial timescales, a general increase in f_max of approximately 1.5× can be achieved by decreasing the measurement noise by an order of magnitude. While difficult to generalise for ice cores covering millennial or longer timescales, due to more complex spectral properties, we apply our approach to the Beyond EPICA Oldest Ice Core (BE-OIC), whose deep ice is yet to be measured. Under our definitions, we see an f_max before diffusion for 1.5 Mya ice of approximately 0.05 kyr⁻¹, with measurement noise below 0.02 ‰ required to recover sub-10 000 year variability after Wiener deconvolution. Throughout the ice core, we find Wiener deconvolution with a standard δ¹⁸O measurement noise of 0.1 ‰ offers an f_max increase of about 1.4× the uncorrected record, which increases up to 2× for a measurement noise of 0.01 ‰. Ultimately, judging the value of improved measurement precision comes down to a case-by-case basis, weighing up the frequency gain against the extra effort required. For example, in the case of the BE-OIC, where layer thinning could result in orbital frequencies approaching f_max, it is arguably worth the extra precision, as the additional effort would only need applying to the few samples spanning >800 kya ice. With the advancement of measuring capabilities, the ice core community is already progressing towards improved precision, which will benefit the BE-OIC and future deep ice core projects.

Appendix A: Composite Dome C Spectra and Diffusion Lengths

The EDC δ¹⁸O data is shown in Fig. A1 and the corresponding power spectra in Fig. A2. The mean spectrum was estimated by interpolating each power spectrum to the highest frequency spacing (the green spectrum) and then averaging for each frequency.

The upper frequency cut-off for the fit was manually selected at 10 kyr⁻¹, as the firn diffusion is evident from the power spectrum. It can be justified mathematically by taking the diffusion length and layer thickness in the firn and showing the loss of amplitude at frequencies above 10 kyr⁻¹. At a depth of 30 m we can expect a diffusion length of σ≈0.08 m. For the Dome C site, the layer thickness at 30 m depth is 42 m kyr⁻¹. Thus the frequency in m⁻¹ corresponding to 10 kyr⁻¹ at this depth is

$$\begin{array}{}\text{(A1)}& f_{\mathrm{depth}}={\displaystyle \frac{\mathrm{10}{\mathrm{kyr}}^{-\mathrm{1}}}{\mathrm{42}{\mathrm{m}\mathrm{kyr}}^{-\mathrm{1}}}}=\mathrm{0.24}{\mathrm{m}}^{-\mathrm{1}}\end{array}$$

Figure A1 Water isotope data from the EDC core used to estimate the local climate spectrum. The data was binned at different timescales, from a depth of 0 m until the depth at which the sample spacing was larger than the given bin size. The vertical lines mark the deepest data point of the corresponding bin size.

Figure A2 Power spectra of the timeseries shown in Fig. A1. The mean of all the power spectra is shown in black.

Therefore, the remaining power of the 10 kyr⁻¹ frequency is equal to

$$\begin{array}{}\text{(A2)}& P_{\mathrm{remaining}}=e^{-(\mathrm{2}\mathit{\pi }\mathit{\sigma }{f}_{\mathrm{depth}}{)}^{\mathrm{2}}}=e^{-(\mathrm{2}\mathit{\pi }\cdot \mathrm{0.08}\cdot \mathrm{0.24}{)}^{\mathrm{2}}}=\mathrm{0.986}\end{array}$$

Or a power loss of 1.4 %. Any frequencies at a faster timescale than this will have lost even more power to firn diffusion, and are therefore not appropriate for our analogue of the isotopic record before diffusion, P₀(f).

Appendix B: Diffusion length calculation

Here we provide a first order calculation of the diffusion length for the BE-OIC. The latter is the result of the diffusion in the firn pore space and the solid diffusion within the ice grain boundaries (Gkinis et al., 2014; Holme et al., 2018). The first process ceases at the pore close-off while the latter persists down to the bottom of the core. The total signal as a function of depth z is described by (Gkinis et al., 2014)

$$\begin{array}{}\text{(B1)}& {\mathit{\sigma }}^{\mathrm{2}}\left(z\right)={\left[S\left(z\right){\mathit{\sigma }}_{\text{firn}}\left(z\right)\right]}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{ice}}^{\mathrm{2}}\left(z\right)\end{array}$$

where S(z) is the ice flow thinning.

B1 The firn diffusion calculation

We use a depth domain spanning from the surface (z=0) to a depth of 2579 m based on the results from Lilien et al. (2021). The latter refers to the depth at which the flowing ice slab meets the stagnant basal unit at the LDC site. For the calculation of the σ_firn quantity we use the analytical solutions described in (Gkinis et al., 2021b) using:

• Surface temperature T₀=213.15 K

• Accumulation $A=\mathrm{0.016}{\mathrm{ma}}^{-\mathrm{1}}$ ice equivalent

• Surface Density ${\mathit{\rho }}_{\mathrm{0}}=\mathrm{330}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$

• Close-off density ${\mathit{\rho }}_{\mathrm{co}}=\mathrm{804}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$

• Atmospheric pressure P=0.6 Atm

B2 The thinning function

For the calculation of the thinning function we use the Lliboutry model (Lliboutry, 1979) for the annual layer thickness described by:

$$\begin{array}{}\text{(B2)}& \mathit{\lambda }\left(z\right)=A\left[\mathrm{1}-{\displaystyle \frac{p+\mathrm{2}}{p+\mathrm{1}}}\left({\displaystyle \frac{z}{H}}\right)+{\displaystyle \frac{\mathrm{1}}{p+\mathrm{1}}}{\left({\displaystyle \frac{z}{H}}\right)}^{p+\mathrm{2}}\right]\end{array}$$

where the shape factor p=5.5 and the ice thickness H=2579 (Lilien et al., 2021). A steady state time scale is then calculated from Eq. (B2) as:

$$\begin{array}{}\text{(B3)}& \overline{t}=\underset{\mathrm{0}}{\overset{z}{\int }}{\left(\mathit{\lambda }\left(z^{\prime }\right)\right)}^{-\mathrm{1}}\mathrm{d}{z}^{\prime }\end{array}$$

We follow the approach shown in Parrenin et al. (2017), Lilien et al. (2021) introducing a change of variable in a pseudo-steady scheme as:

$$\begin{array}{}\text{(B4)}& \overline{t}=\underset{\mathrm{0}}{\overset{t}{\int }}{\displaystyle \frac{A_{\mathrm{edc}}\left(t^{\prime }\right)}{{\overline{A}}_{\mathrm{edc}}\left(t^{\prime }\right)}}\mathrm{d}{t}^{\prime },\end{array}$$

where A_edc is the accumulation at the Dome C site as a function of time (Bazin et al., 2013). We solve Eq. (B4) for t using the Brent algorithm (Brent, 1973) implemented in the Scipy toolbox (Jones et al., 2001) and obtain a timescale adjusted for a time dependent accumulation. For ages older than the bottom of EDC at 803 ka the normalized accumulation is assumed to be equal to 1.

B3 The ice diffusion calculation

For the calculation of ${\mathit{\sigma }}_{\text{ice}}^{\mathrm{2}}\left(z\right)$ we integrate the ice diffusivity in the time domain taking into account the effect of ice flow thinning (Gkinis et al., 2014; Holme et al., 2018):

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

We use the Ramsaeir parametrization for the diffusivity coefficient in solid ice (Ramseier, 1967):

$$\begin{array}{}\text{(B6)}& D_{\mathrm{ice}}=\mathrm{9.2}\times {\mathrm{10}}^{-\mathrm{4}}\cdot \exp \left(-{\displaystyle \frac{\mathrm{7186}}{T}}\right){\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}\end{array}$$

Figure B1 Modelled diffusion length profile for the BE-OIC (black). The diffusion lengths used in the simulated isotopic records in Sect. 4.3 were selected from the profile (white-noise case is red, power-law case is blue).

B4 Temperature profile calculation

We calculate the temperature profile by solving the parabolic differential equation accounting for advection using an inverted depth domain ie distance from top of the stagnant ice slab:

$$\begin{array}{}\text{(B7)}& {\displaystyle \frac{\partial T}{\partial t}}=D_{\mathrm{T}}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial h^{\mathrm{2}}}}+w{\displaystyle \frac{\partial T}{\partial h}}\end{array}$$

where D_T is the thermal diffusivity coefficient ($\mathrm{1.13}\times {\mathrm{10}}^{-\mathrm{6}}$) and the advection coefficient $w=-\mathit{\lambda }$. We use a automated finite difference scheme using the method of lines utilizing the Scikit-FDiff package (Cellier and Ruyer-Quil, 2019). The model run spans 10⁶ a with a time step of 100 a and a spatial resolution of 10 m.

We use a Dirichlet boundary condition for the surface with T₀=213.15 K and a Neumann boundary condition for the boundary between the stratified ice and the stagnant ice. The initial profile is linear with T₀ at the top and a bottom temperature at 272 K. Choosing the value of the Geothermal Heat Flux is a question that deserves a dedicated study. For this tentative calculation we have chosen $G=\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. This value lies in the low range of the spectrum indicated by Van Liefferinge et al. (2018), assuming though that we are omitting the stagnant ice slab with a thickness of ≈200 m. This value is a reasonable choice for a first order calculation of the temperature profile and the ice diffusion length thereafter. Further consideration of the temperature profile when the borehole is logged will obviously allow for a better estimation of the diffusion at the bottom of the BE-OIC.

Appendix C: Sensitivity to A_min

Figure C1 Same as Fig. 5 but with A_min=0.05.

Figure C2 Same as Fig. 5 but with A_min=0.9.