What climate signal is contained in decadal- to centennial-scale isotope variations from Antarctic ice cores?

Ice-core-based records of isotopic composition are a proxy for past temperatures and can thus provide information on polar climate variability over a large range of timescales. However, individual isotope records are affected by a multitude of processes that may mask the true temperature variability. The relative magnitude of climate and non-climate contributions is expected to vary as a function of timescale, and thus it is crucial to determine those temporal scales on which the actual signal dominates the noise. At present, there are no reliable estimates of this timescale dependence of the signal-to-noise ratio (SNR). Here, we present a simple method that applies spectral analyses to stable-isotope data from multiple cores to estimate the SNR, and the signal and noise variability, as a function of timescale. The method builds on separating the contributions from a common signal and from local variations and includes a correction for the effects of diffusion and time uncertainty. We apply our approach to firn-core arrays from Dronning Maud Land (DML) in East Antarctica and from the West Antarctic Ice Sheet (WAIS). For DML and decadal to multi-centennial timescales, we find an increase in the SNR by nearly 1 order of magnitude (∼0.2 at decadal and ∼1.0 at multi-centennial scales). The estimated spectrum of climate variability also shows increasing variability towards longer timescales, contrary to what is traditionally inferred from single records in this region. In contrast, the inferred variability spectrum for WAIS stays close to constant over decadal to centennial timescales, and the results even suggest a decrease in SNR over this range of timescales. We speculate that these differences between DML and WAIS are related to differences in the spatial and temporal scales of the isotope signal, highlighting the potentially more homogeneous atmospheric conditions on the Antarctic Plateau in contrast to the marine-influenced conditions on WAIS. In general, our approach provides a methodological basis for separating local proxy variability from coherent climate variations, which is applicable to a large set of palaeoclimate records.

1 Introduction

Ice cores represent key archives for studying polar climate variability on timescales beyond instrumental observations. The isotopic composition of water stored in the ice serves as a proxy for reconstructing past temperature variations (Dansgaard, 1964; Jouzel and Merlivat, 1984; Jouzel et al., 2003) over a wide range of timescales ranging from subannual to glacial–interglacial variations. This can provide valuable insights into the timescale dependence of climate variability, which is thought to be an inherent property of the climate system (Hasselmann, 1976; North et al., 2011; Lovejoy et al., 2013; Rypdal et al., 2015).

However, the interpretation of isotope records in terms of local atmospheric temperatures is complicated by a multitude of processes that distort the original relationship present in precipitation (Casado et al., 2018; Fujita and Abe, 2006; Sjolte et al., 2011; Steen-Larsen et al., 2011; Touzeau et al., 2016). To a first approximation, changes in isotopic composition are only recorded in the ice if there is snowfall, while the role of water vapour exchange processes in between precipitation events is still debated (Steen-Larsen et al., 2011; Stenni et al., 2016; Casado et al., 2018; Ritter et al., 2016; Münch et al., 2017a). Both the seasonality and interannual variability in the seasonality of precipitation events is known to strongly affect the isotopic composition of snow layers by introducing a bias to the mean values, e.g. the annual average (Laepple et al., 2011; Sime et al., 2009; Steig et al., 1994), and by adding variability to the signal (Laepple et al., 2018; Persson et al., 2011). Furthermore, uneven deposition of snow, in combination with steady and strong surface winds, lead to constant erosion, drift and vertical mixing of the surface snow, giving rise to strong spatial variability (Fisher et al., 1985; Münch et al., 2016, 2017a). Finally, once the snow is deposited, the diffusion of vapour through the firn column smoothes the isotope variations (Johnsen, 1977; Whillans and Grootes, 1985). This significantly reduces the high-frequency power spectral density of the variations (Johnsen et al., 2000; van der Wel et al., 2015) and substantially shapes the isotope variability (Laepple et al., 2018). Overall, these processes are to a first approximation not directly linked to variations in temperature and therefore add a significant amount of noise to the isotope records, especially in low-accumulation regions on the Antarctic Plateau. This has been demonstrated by the low representativity of individual ice-core measurements (Karlöf et al., 2006; Münch et al., 2016) and questions the usability of a direct interpretation of high-resolution isotope variations in terms of temperature variability (Laepple et al., 2018), particularly when the climate signal itself is only relatively weak. While deep ice cores show a consistent picture for the strong glacial–interglacial variations (Jouzel et al., 2007), it may be questionable whether the Holocene variability recorded in individual cores depicts the true temperature variability (Kobashi et al., 2011). Spatial or temporal averaging of isotope records thus provide important tools to reduce overall noise (Fisher et al., 1996; Münch et al., 2016; Stenni et al., 2017).

Previous studies provided first insights into the relationship between climate signal and noise for short spatial and temporal scales (Fisher et al., 1985; Münch et al., 2016). However, an extension to longer scales, which is important for the interpretation of Holocene temperature reconstructions, is still missing. Furthermore, in order to correct the isotope-inferred variability estimates for the noise contribution, it is necessary to know the variance of the noise across timescales, i.e. its spectral shape. Here, we present a simple spectral method to separate the local noise from spatially coherent signals using information from several isotope records, including a correction for diffusion and time uncertainty. We apply our model to two spatial arrays of firn cores: (1) from Dronning Maud Land in East Antarctica and (2) from the West Antarctic Ice Sheet. Our objective is to derive an improved estimate of the temperature variability in both regions and to learn about the timescale dependence of the signal-to-noise ratio in ice-core isotope data. For Dronning Maud Land, our results confirm the noise levels inferred in previous studies on short temporal scales (Münch et al., 2016, 2017a) and also suggest white noise on longer timescales, which results in an increase in the isotopic signal-to-noise ratio (SNR). Unexpectedly, the SNR inferred from the West Antarctic data is found to show the opposite timescale dependence. These results may point towards marked differences in the spatial and temporal scales of the isotope signals and reveal gaps in our current understanding of the underlying processes.

2 Data and methods

2.1 Isotope records from Dronning Maud Land and West Antarctica

We analyse published oxygen isotope records of annually dated firn cores from two contrasting Antarctic regions (Table 1): (1) Dronning Maud Land (DML) on the Antarctic Plateau and (2) the West Antarctic Ice Sheet (WAIS). While the DML core sites are located in a rather flat region relatively isolated from the coast (average elevation 2900 m), the WAIS core sites are lower in elevation (1600 m on average) and therefore potentially more influenced by marine conditions.

For DML, we use a total of 15 records, which were collected during the EPICA (European Project for Ice Coring in Antarctica) pre-site survey (Oerter et al., 2000) and published in Graf et al. (2002a). All records cover at least the last ∼ 200 years and form our first data set, DML1. Three of these records (B31–B33) span the last ∼ 1000 years and are therefore included in a second separate data set for this time span (DML2). Core B32 is in close proximity (∼1 km) to the deep EPICA DML (EDML) ice-core site at Kohnen Station (EPICA community members, 2006). The record chronologies were established from seasonal layer counting of chemical impurity records constrained by tie points from the dating of volcanic ash layers (Graf et al., 2002a). The resulting uncertainty of the chronologies has been reported to be ∼2 % of the time interval to the nearest tie point (Graf et al., 2002a), which translates to a maximum time uncertainty of ∼ 1.2 years for the short records and of ∼ 3.5 years for the long records. Since our method (Sect. 2.2) relies on all records having equal lengths, we restrict the time span for the DML1 data set to the period from 1801 to 1994 CE and to the period from 1000 to 1994 CE for the DML2 data set.

The WAIS data set selected for this study consists of five isotope records (Steig et al., 2013) collected during the US ITASE (International Trans-Antarctic Scientific Expedition) project (Mayewski et al., 2005), including the core WDC2005A from the WAIS Divide ice core (WDC; WAIS Divide Project Members, 2013), and cover the time interval from 1800 to 2000 CE. The oxygen isotope data of cores ITASE-2000-4 and ITASE-2000-5, previously published at subannual resolution (Steig et al., 2013), have been block-averaged in this study to obtain annual resolution data. The core selection was based on the constraint of covering a similarly large area and sufficiently long time period to allow a meaningful comparison with the DML1 records. The relative time uncertainty between the WAIS records, based on dating through annual layer counting of chemical trace species validated by identification of volcanic marker horizons, was reported to be no more than 1 year (Steig et al., 2005).

Graf et al. (2002b)Graf et al. (2002b)Steig (2013)

Table 1 Overview of the firn cores (sorted into three data sets) used in this study. Listed are the covered time spans of each core array (in yr CE), the number of records in each array (n), the region of origin (latitude–longitude range), the range of site elevations, local accumulation rates ($\dot{b}$) and 10 m firn temperatures (T_firn) reported in the original publications, and the range of intercore distances (d). The range of WAIS firn temperatures is based on ERA-Interim annual mean anomalies (Dee et al., 2011) with respect to the WDC2005A site.

^a Firn cores FB9804, FB9805, FB9807–FB9811, FB9813–FB9817, B31–B33. ^b Firn cores B31–B33. ^c Firn cores WDC2005A, ITASE-1999-1, ITASE-2000-1, ITASE-2000-4, ITASE-2000-5. ^d Oerter et al. (2000). ^e Kaspari et al. (2004). ^f WAIS Divide Project Members (2013). ^g Dee et al. (2011).

2.2 Model for the separation of signal and noise in the spectral domain

Our approach in general assumes that individual isotope records from a certain region contain two contributions: (i) a signal common to all cores from that region and (ii) independent noise components, which are, for example, related to spatial variability from redistribution of snow (stratigraphic noise) or to precipitation intermittency. By utilising several records we can disentangle both contributions and estimate the underlying common and noise signals. The approach is similar to the analysis of variance (Fisher et al., 1985) but uses the power spectra of the time series to derive timescale-dependent estimates.

More formally, given a core array of n isotope records, the power spectral density (PSD) of an individual record from site i, 𝒳_i(f), where f denotes frequency, is the sum ${\mathcal{X}}_i\left(f\right)=\mathcal{C}\left(f\right)+{\mathcal{N}}_i\left(f\right)$, where 𝒞(f) and 𝒩_i(f) are the original spectra of the common signal and the noise component prior to the smoothing by molecular diffusion of water vapour within the porous firn (Johnsen et al., 2000; van der Wel et al., 2015). To relate this with the actual measured signal, ${\widehat{\mathcal{X}}}_i\left(f\right)$, we additionally have to account for the measurement process, which adds additional noise to the diffused record. ${\widehat{\mathcal{X}}}_i\left(f\right)$ is thus given by

$$\begin{array}{}\text{(1)}& {\widehat{\mathcal{X}}}_i\left(f\right)={\mathcal{G}}_i\left(f\right)\left[\mathcal{C}\left(f\right)+{\mathcal{N}}_i\left(f\right)\right]+\mathrm{\Sigma }.\end{array}$$

Here, 𝒢_i(f) is a linear transfer function that acts as a low-pass filter and accounts for the diffusion process (Appendix A), and Σ is the measurement error. At annual resolution, Σ is typically at least 1 order of magnitude smaller than the stratigraphic noise level (Münch et al., 2016) and is therefore neglected in the following. We now calculate the mean spectrum, ℳ(f), of all individual spectra ${\widehat{\mathcal{X}}}_i\left(f\right)$. Assuming that the statistical properties of the individual noise terms are identical for all n records, we obtain

$$\begin{array}{}\text{(2)}& \mathcal{M}\left(f\right)={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\widehat{\mathcal{X}}}_i\left(f\right)=\stackrel{\mathrm{‾}}{\mathcal{G}}\left(f\right)\left[\mathcal{C}\left(f\right)+\mathcal{N}\left(f\right)\right],\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathcal{G}}\left(f\right)=\mathrm{1}/n{\sum }_{i=\mathrm{1}}^n{\mathcal{G}}_i\left(f\right)$ is the average diffusion transfer function. In contrast to forming the spectral mean, we can also average the n isotope records in the time domain and then calculate the PSD of this “stacked” record. Here, the noise component will be reduced by a factor of n compared to the mean spectrum in Eq. (2), since we assume independent noise between the sites. Additionally, we have to take into account the time uncertainty of the dated records. This does not affect the overall shape of broadband spectra (Rhines and Huybers, 2011) but diminishes the correlation between the records (Haam and Huybers, 2010) and thus their spectral coherence, which reduces the variance of the common signal in the stacked record. The PSD of the stacked record is thus

$$\begin{array}{}\text{(3)}& \mathcal{S}\left(f\right)\approx \stackrel{\mathrm{‾}}{\mathcal{G}}\left(f\right)\left[\mathrm{\Phi }\left(f\right)\mathcal{C}\left(f\right)+{\displaystyle \frac{\mathrm{1}}{n}}\mathcal{N}\left(f\right)\right],\end{array}$$

where Φ(f) is a linear transfer function accounting for the effect of time uncertainty. Applying the average diffusion transfer function, $\stackrel{\mathrm{‾}}{\mathcal{G}}\left(f\right)$, also to the spectrum of the stacked record is a valid approximation in the case of similar 𝒢_i(f), which we confirmed for our records (Appendix A). From the expressions for the mean spectrum (ℳ, Eq. 2) and the spectrum of the stacked record (𝒮, Eq. 3) we can now derive expressions for the spectra of the common signal 𝒞 and the noise 𝒩 (omitting the explicit frequency dependence in the notation here),

$$\begin{array}{}\text{(4a)}& {\displaystyle }\mathcal{C}& {\displaystyle }={\displaystyle \frac{n}{n-{\mathrm{\Phi }}^{-\mathrm{1}}}}{\mathrm{\Phi }}^{-\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathcal{G}}}^{-\mathrm{1}}[\mathcal{S}-n^{-\mathrm{1}}\mathcal{M}],\text{(4b)}& {\displaystyle }\mathcal{N}& {\displaystyle }={\displaystyle \frac{n}{n-{\mathrm{\Phi }}^{-\mathrm{1}}}}{\stackrel{\mathrm{‾}}{\mathcal{G}}}^{-\mathrm{1}}[\mathcal{M}-{\mathrm{\Phi }}^{-\mathrm{1}}\mathcal{S}].\end{array}$$

2.3 Transfer functions for vapour diffusion and time uncertainty

For estimating the transfer functions for diffusion and time uncertainty, the inverses of which are used to correct the signal and noise spectra (Eq. 4), we use numerical simulations, since no closed form expressions are available. For this, we create surrogate records mimicking the individual core arrays and simulate the effects of diffusion and time uncertainty on the surrogate spectra of the stacked records. We model the firn diffusion as the convolution of the original time series with a Gaussian kernel (Johnsen et al., 2000). The width of the kernel is set by the diffusion length σ, which is site-specific and a function of depth. For modelling the time uncertainty we use the approach of Comboul et al. (2014). Model parameters are the rates of missing and double-counted annual layers, which we set to match the reported time uncertainties of the isotope records. Appendix B gives a detailed description of the simulations including the estimated transfer functions (Fig. B1). Because the diffusion correction (${\stackrel{\mathrm{‾}}{\mathcal{G}}}^{-\mathrm{1}}$ in Eq. 4) strongly amplifies both the raw spectra as well as their uncertainties on the fast frequencies, we restrict our analyses to the frequency region where the estimated transfer function, $\stackrel{\mathrm{‾}}{\mathcal{G}}$, assumes values ≥0.5, equivalent to a correction factor ≤2 in Eq. (4). This avoids large uncertainties and translates to a maximum frequency that is used for the analyses (hereafter the cut-off frequency) of 1∕5 yr⁻¹ for DML1, 1∕8.5 yr⁻¹ for DML2 and 1∕2.8 yr⁻¹ for WAIS (Fig. B1).

2.4 Timescale-dependent estimate of the signal-to-noise ratio

The frequency-dependent SNR, F(f), is defined as the ratio of the signal and noise spectra,

$$\begin{array}{}\text{(5)}& F\left(f\right)={\displaystyle \frac{\mathcal{C}\left(f\right)}{\mathcal{N}\left(f\right)}}.\end{array}$$

As we explicitly neglect measurement noise, this quantity is unaffected by the effect of diffusion or its correction (Eq. 4) but directly influenced by time uncertainty, which biases the uncorrected SNR towards zero. Typically, firn or ice-core isotope records are averaged onto a fixed temporal resolution Δt (the “averaging period”) during preprocessing. The signal-to-noise variance ratio after this averaging step is given by

$$\begin{array}{}\text{(6)}& \stackrel{\mathrm{‾}}{F}\left(f_{\mathrm{Nyq}}\right)={\displaystyle \frac{{\int }_{f_{\mathrm{0}}}^{f_{\mathrm{Nyq}}}\mathcal{C}\left(f\right)\mathrm{d}f}{{\int }_{f_{\mathrm{0}}}^{f_{\mathrm{Nyq}}}\mathcal{N}\left(f\right)\mathrm{d}f}},\end{array}$$

where f₀ is the lowest frequency of the spectral estimates and f_Nyq is the Nyquist frequency for the chosen averaging period, i.e. 1∕(2Δt). Since our records have different lengths, we choose a common minimum value for f₀ of 1∕100 yr⁻¹ for the integrations in Eq. (6). The value of $\stackrel{\mathrm{‾}}{F}$ can then be used to obtain the correlation between the time series of the common signal c and a stacked record $\stackrel{\mathrm{‾}}{x}$ built from the spatial average of n individual records (Fisher et al., 1985):

$$\begin{array}{}\text{(7)}& r(c,\stackrel{\mathrm{‾}}{x})\left(f_{\mathrm{Nyq}}\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+{\left(n\stackrel{\mathrm{‾}}{F}\left(f_{\mathrm{Nyq}}\right)\right)}^{-\mathrm{1}}}}}.\end{array}$$

2.5 Spectral analysis

Missing years in the published records (∼1.6 % of all data points) are linearly interpolated. Power spectra are then estimated using Thomson's multitaper method with three windows (Percival and Walden, 1993) with linear detrending before analysis. This approach is known to introduce a small bias at the lowest frequencies, and we omit the lowest frequency from all plots and calculations. In general, spectral smoothing is necessary to improve the quality of the estimates from short time series. Here, we use a Gaussian smoothing kernel with varying width that is proportional to the applied frequency and thus constant in logarithmic frequency space (Kirchner, 2005). The smoothing width in logarithmic units is chosen to be 0.1 for the WAIS data and 0.15 for the DML data. To avoid biased estimates at the low- and high-frequency boundaries, the kernel is truncated on both sides to maintain its symmetry. Logarithmic smoothing preserves power-law scaling of spectra and is thus useful for climate spectra.

3 Results

3.1 Illustrating the methodological approach

In order to illustrate our method (Sect. 2.2), we first use the DML1 data set to demonstrate the individual steps involved in the analysis.

Each power spectrum derived from an individual record of the DML1 firn-core array provides a timescale-dependent representation of the isotope variations in the study region (Fig. 1, thin grey lines). The differences between the individual spectra are not only due to differences in the isotope variability between the records, but also caused by the spectral uncertainty from the finite length of each record. The mean spectrum, ℳ, of all individual spectra reduces this spectral uncertainty and thus provides a more robust estimate of the single records' spectrum (Fig. 1, black line). It shows a two-part structure with a near-constant PSD above decadal timescales (i.e. is “white”) and a strong decrease in spectral power towards shorter timescales, which is expected from vapour diffusion in the firn (Eq. 2 and Appendix B).

Figure 1 Detailed results of estimated PSD for the DML1 data set. Thin grey lines show the individual power spectra for each record with the mean spectrum indicated by the thick black line. The dashed black line shows the null hypothesis according to which all isotope variations are noise; the brown line depicts the spectrum from averaging all records in the time domain (the stacked record). The extent of the blue (red) shading is proportional to the uncorrected PSD of the signal (noise) (Eq. 4). The vertical dashed line denotes the cut-off period for the diffusion correction (see Methods).

The mean spectrum divided by the number of records (ℳ∕n; here, n=15; Fig. 1, dashed line) is the expected spectrum if all isotope variations present in the firn-core records were independent noise. In comparison, averaging in the time domain across records that contain noise and additionally a common (i.e. spatially coherent) signal will result in a spectrum 𝒮, where the noise component is also reduced by 1∕n but with the common signal left unaltered (Eq. 3, neglecting time uncertainty), and which is thus located between the mean spectrum and the mean spectrum divided by n. This spectrum 𝒮 for the DML1 stack (Fig. 1, brown line) is for short timescales (periods from ∼ 3 to 5 years) nearly identical to the mean spectrum divided by n, which suggests that here the variability of the individual records is consistent with the null hypothesis of independent noise. The divergence from the white-noise level close to the 2-year Nyquist period is likely an artefact from the jitter of the annual averages as a result of the uncertainty in defining annual depth increments upon dating. In contrast to the short periods below 10 years, the individual records clearly contain a common isotope signal on longer timescales that “survives” the averaging process when building the stacked record and which increases in power with increasing timescale (Fig. 1, brown line). Hence, the differences between 𝒮 and ℳ∕n and between 𝒮 and ℳ inform us about the average signal and noise content of the individual isotope records, but they need to be corrected for the residual amount of noise contained in 𝒮 and for the loss of variance from diffusion and time uncertainty (Eq. 4).

Figure 2 Estimated signal (a, c) and noise (b, d) spectra of Antarctic isotope records. Results are shown based on the spectral correction model for the records from Dronning Maud Land (DML, a, b) and from the West Antarctic Ice Sheet (WAIS, c, d). The raw estimates (dotted lines) show the spectra prior to any correction, while corrected estimates differentiate between the correction for time uncertainty only (dashed lines) and the full correction including time uncertainty and diffusion (solid lines). For DML, results include both the short 200-year-long records (DML1; blue lines in a, red lines in b) and the longer records covering the last 1000 years (DML2; black lines).

3.2 Timescale dependence of DML and WAIS signal and noise variability

After this detailed description for DML1, we now turn to the results of the estimated signal and noise spectra for all three data sets. In general, the shape of the signal spectra is, as a result of the corrections, clearly distinct from the mean spectrum of the individual isotope variations. This is seen, for example, in the corrected DML1 signal spectrum which indicates a much more steady increase in PSD from short to long timescales (Fig. 2a, solid blue line) compared to the mean spectrum (Fig. 1, black line). This is confirmed by the three 1000-year-long DML2 records, which have a common signal that exhibits a roughly similar power spectrum in the range of timescales that overlap (∼ 10- to 100-year period). We partly expect this, since the longer cores are also part of the DML1 data set, but the increase in signal power in this frequency band seems to be a general feature over the last 1000 years. In addition, the DML2 signal spectrum shows a further and similar increase for timescales beyond the 100-year period. This change in spectral shape from the mean to the signal spectrum results from two contributions in the correction process: firstly, the average isotope variability is corrected for the noise contribution (“uncorrected signal”, dotted lines in Fig. 2a). This correction is mostly apparent on the longer timescales, leading to a higher slope in PSD of the signal spectrum compared to the mean spectrum (for DML1, the signal increases from the 10- to 100-year period by a factor of ∼2.6, the mean spectrum by a factor of ∼1.4). Secondly, the corrections for the loss in spectral power by the effects of diffusion and time uncertainty (here important for time periods from ∼ 20- to 5-years) lead to a smaller increase in PSD of the signal spectra with increasing timescale compared to the raw estimates (compare dotted, dashed and solid lines in Fig. 2a). In sharp contrast to DML, the corrections for the WAIS records yield a signal spectrum (Fig. 2c, solid line) that exhibits a rather flat appearance throughout and indicates decreasing PSD towards centennial timescales. Much of the flat character is caused by the diffusion and time uncertainty corrections, which strongly amplify the raw signal power on the subdecadal timescales; the decrease in PSD on the long timescales is a result of the correction for the noise contribution. We note, however, that the spectral uncertainty at the lowest frequency is high, since fewer data points contribute to the estimated PSD for lower frequencies when using logarithmic smoothing.

Figure 3 Estimated timescale dependence of ice-core isotope signal-to-noise ratios. Results are shown for the DML (black) and WAIS (blue) isotope records. The results for DML are based on combining the spectra from DML1 and DML2 (see text).

A difference between the regions can also be seen in the noise spectra (Fig. 2b and d). Prior to any correction, the DML1 and DML2 noise spectra are different on shorter timescales (≲ 20-year period), but become consistent with each other after the correction, showing essentially white PSD (Fig. 2b). In comparison, the corrected WAIS noise spectrum shows an increase in spatially incoherent isotope variations towards longer timescales (Fig. 2d), which correlates with the decrease in signal power on centennial timescales (Fig. 2c).

The corrected signal and noise spectra directly provide an estimate of the timescale dependence of the SNR (Eq. 5). To derive a single estimate for DML, we combine the DML2 spectra for timescales above the decadal period with the respective DML1 spectra from the subdecadal frequency band. Again, the results are very different between DML and WAIS (Fig. 3). For DML, the SNR shows a continuous increase with timescale, with values ranging from ∼ 0.1 at the 5-year period to ≲ 0.2 at the decadal timescale and increasing further to ∼1 at multi-centennial timescales. The estimate for WAIS exhibits the opposite timescale dependence, with SNR values of ∼0.5–0.7 for interannual timescales that continuously decrease to ∼0.1 at centennial timescales.

Figure 4 Estimated correlation of stacked isotope records from (a) DML and (b) WAIS with their common signal as a function of the averaging period and the number of firn cores included in the average. The correlations are based on Eq. (7) with the same signal and noise spectra that were used for Fig. 3.

A complementary picture is obtained from the expected correlation between the time series of a stacked record and the underlying common signal as a function of both the number of averaged records and the temporal averaging period (Eq. 7, Fig. 4). For DML, the correlation for an individual record at interannual resolution is rather low, with roughly 0.4–0.5 (Fig. 4a); for WAIS, this correlation is slightly higher (∼0.5–0.6, Fig. 4b). For longer averaging periods, the correlation for DML shows a steady increase with the averaging period, in line with the increase in SNR, reaching values around 0.6 for single records and multi-decadal averaging periods (Fig. 4a). Naturally, the correlation further increases with the number of averaged cores. For WAIS, the correlation with the underlying signal decreases with the averaging period and only increases when more records are averaged (Fig. 4b).

4 Discussion

4.1 Physical interpretation of the analysis

We presented a method and the results of separating the variability recorded by Antarctic isotope records into two contributions: local variations (noise, Eq. 4b) and spatially coherent variations (signal, Eq. 4a). We now assess the physical meaning of these terms by discussing the relevant spatial scales of the major processes that influence the isotopic composition of the records: (i) temperature variations, (ii) atmospheric circulation, (iii) precipitation intermittency and (iv) stratigraphic noise.

Classically, the isotopic composition of firn and ice cores is interpreted as being related to variations in local air temperature (Dansgaard, 1964; Jouzel and Merlivat, 1984; Jouzel et al., 2003). The extent to which spatially distributed isotope records are influenced by a common temperature signal then depends on the decorrelation scale of the temperature anomalies, which is, for annual mean temperatures, typically of the order of hundreds to thousands of kilometres and increases with timescale (Jones et al., 1997). For our study regions, annual mean temperature variations from the ERA-Interim reanalysis (Dee et al., 2011) exhibit decorrelation scales of ∼1200 km (Appendix C), which are much larger than the individual intercore distances (<400 km, Table 1). This implies that, if the temperature variations were fully recorded by the array of firn-core records, they would to a large extent contribute to the estimated signal spectrum. In fact, for the temperature reanalysis fields in the study regions, the average of all correlations between sites separated by less than the maximum intercore distances (Fig. C1) suggests that the estimated signal spectra for DML (WAIS) would contain 87 % (94 %) of the total temperature variability, while the remaining fraction would be interpreted as noise.

However, changes in large-scale atmospheric circulation can lead to variations in the source and the pathways of the moisture, which can affect the isotopic composition of the precipitation that is formed, independently of local temperature changes (Schlosser et al., 2004). In general, the spatial scales of such variations are unclear. For the studied DML core array, which is located on the rather flat and remote Antarctic Plateau, one could expect that the effect is small and possibly spatially coherent, thus contributing to the estimated signal term. For WAIS, isotope data have been interpreted to also reflect changes in atmospheric circulation and sea-ice cover of the adjacent oceans (Küttel et al., 2012; Steig et al., 2013). These processes might exhibit, through the ice-sheet topography, a regional expression at the individual firn-core sites, which are generally lower in elevation and located near the ice divide, as suggested by the significant spatial differences in observed intercore correlations (Küttel et al., 2012). Such effects would affect the estimated noise component.

Major additional contributions to the overall variability of isotope data arise from precipitation intermittency, i.e. interannual variations in the seasonality of precipitation events (Laepple et al., 2018; Persson et al., 2011), and from stratigraphic noise (Fisher et al., 1985; Karlöf et al., 2006; Münch et al., 2016) created during deposition of the snow. Both processes exhibit significantly different decorrelation scales. The decorrelation scale of precipitation intermittency is related to the decorrelation scale of the precipitation itself, which is expected to be coherent on a local scale (i.e. ∼100 m) and to decorrelate on scales generally smaller than the temperature decorrelation scales. In both regions, the similar analyses of ERA-Interim data as for the temperature variations suggest decorrelation scales of the precipitation of 300–500 km (Appendix C). This is supported by an analysis of measurements from autonomic weather stations (Reijmer and van den Broeke, 2003) that indicate that accumulation in DML arises from many small (1–2 cm) and a few larger events (≳5 cm) of snowfall, where the major events occur only a few times per year without any clear seasonality but often over large areas. In contrast, stratigraphic noise is a short-scale phenomenon. Its generation is related to the uneven deposition (Fisher et al., 1985) and the constant erosion and redeposition of the surface snow by wind (Münch et al., 2016). Both processes are connected to the local surface undulations, as these directly interact with the snow deposition but also strongly shape the near-surface wind field. This is suggested by the observed decorrelation scale of the stratigraphic noise at EDML (<5 m, Münch et al., 2016), which is similar to the typical scale of the surface undulations.

These differences in decorrelation scales also become apparent when analysing the relative contributions of both processes to the total isotope variability. At EDML, stratigraphic noise provides ∼50 % of the total variance on the seasonal timescale, as suggested from the average correlation of individual shallow isotope profiles separated above (>5 m) the decorrelation scale of the stratigraphic noise but below 1 km (Münch et al., 2016, 2017a). A much higher noise level of at least 80 % of the total variance (i.e. higher by a factor of ≥1.6) has been independently inferred from comparing the observed seasonal isotope variability to the expectation from a profile that is a mixture of a deterministic seasonal cycle and diffused noise (Laepple et al., 2018). The apparent mismatch with the noise level from stratigraphic noise can be reconciled by asserting that a significant part of the common isotope signal at EDML on the local scale is coherent noise from precipitation intermittency, leading to an underestimation of the total noise level when analysing only the interprofile correlations. On the larger spatial scales of the firn-core arrays (array scale) analysed here, the effect of precipitation intermittency should then appear, at least partly, as a noise term, since the spatial precipitation patterns are expected to decorrelate on these scales. Therefore, in general, the variability contribution by stratigraphic noise will fully appear in the noise spectra, while precipitation intermittency is expected to partly contribute to the noise spectra but also partly appear in the signal spectra.

In summary, given the large decorrelation scales of atmospheric temperature variations and the generally smaller scales of the other terms, one could interpret the estimated isotope signal spectra to a first approximation as temperature signals. However, this clearly will be an upper bound of the true temperature signal, since other processes can also lead to spatially coherent isotope signals. Furthermore, we have neglected the transfer function from isotopic ratios to temperatures, and other less constrained effects that affect the isotope signal from the atmosphere to the snow (Casado et al., 2018).

4.2 Interpretation of the estimated signal and noise spectra

The raw noise spectra, i.e. prior to correction, derived from the two DML data sets, exhibit a clear imprint from the diffusional smoothing in the firn. This is suggested by their common decrease in PSD towards shorter time periods (Fig. 2b), since the effect of diffusion is stronger at higher frequencies. It is corroborated by comparing the loss in PSD between the two data sets, which for DML2 is stronger towards the high-frequency end and also extends further towards longer time periods. This is due to the stronger diffusional smoothing in the older sections of the cores that are only contained in the DML2 records, since the diffusion process had more time to act there since the deposition of the snow. The applied correction reconciles both noise estimates, as it takes these differences into account, and reveals a close-to-constant noise level (white noise) across the range from subdecadal to multi-centennial timescales.

This near constancy of the noise level suggests that the noise-creating processes are independent of the timescale. This seems plausible for stratigraphic noise as it is also indicated by the observed small memory in the interannual variations of the surface topography at EDML (Laepple et al., 2016; Münch et al., 2016). Furthermore, we would not expect strong changes in surface wind regimes or depositional characteristics for the rather stable climatic conditions over the studied time periods. To test whether the effect of precipitation intermittency additionally contributes to the noise spectrum, as suggested from the decorrelation scales of the ERA-Interim precipitation fields, we apply our spectral analysis to the available isotope profiles on the local scale (∼100 m) at EDML (Münch et al., 2017a, b) and compare the estimated noise spectra between the local scale and the array scale (Fig. 5). Although the spectral estimates do not overlap, the results indicate an offset between the noise levels at both spatial scales. For the 5-year period, the PSD of the noise on the array scale is about 1.7 times higher than at the local scale. Some increase in the noise level is expected from the small contribution of uncorrelated temperature variability, but the ∼1.7-fold increase suggests that most of the additional noise arises from spatially incoherent precipitation intermittency.

Figure 5 Comparison of DML noise spectra as a function of intersite distance. The red curve (array scale, ∼100 km) shows the section for periods ≤ 50 years of the composite noise spectrum from the DML1 and DML2 firn-core records (Fig. 2b and text). The blue curve (local scale, ∼100 m) depicts the noise spectrum estimated from 22 shallow (∼3.5 m depth) trench profiles from EDML (Münch et al., 2017a, b). The depth scale of the trench data has been converted into time units using a constant accumulation rate of 25 cm of snow per year, neglecting the small compression by densification (Münch et al., 2017a). Blue shading denotes the spectral range from an assumed 20 % uncertainty of the accumulation rate. Note that the trench noise spectrum is not corrected for the diffusional smoothing, which has a negligible effect on the trench records for the 5-year period and thus does not affect the comparison of the noise levels.

This supports the interpretation of the estimated DML signal spectra (Fig. 2a) as temperature variability, assuming that the influence of atmospheric circulation changes is small. Fitting a power law of the form f^−β, where f denotes frequency, yields a slope of the DML2 signal spectrum of roughly β∼0.6 for variations above the 20-year period. This is higher than the slopes of the mean spectra of the single records (β∼0.2 for DML2, β∼0 for DML1), which underlines again that the high noise level in individual records strongly masks the true spectral shape of the signal. Such a power-law increase in temperature variability is not unexpected since it was also observed in spectra from marine-proxy-inferred sea surface temperature variations (Laepple and Huybers, 2014) and from other proxy and instrumental temperature records (Huybers and Curry, 2006; Pelletier, 1998).

The resulting SNR (Fig. 3) of the DML data illustrates on which timescales the signal dominates the isotope variability recorded in individual records, and how this translates into the representativity of isotope variations when averaging records in space and time (Fig. 4). We can test the validity of our spectral approach by comparing the estimated SNR to published estimates. Graf et al. (2002a) analysed the same data set as the present study and found an SNR at annual resolution of $\stackrel{\mathrm{‾}}{F}=\mathrm{0.14}$ based on the average intercore correlations (r=0.35). Similar values have been obtained using a larger set of cores from the same region that cover a shorter time period (Altnau et al., 2015). From integrating the estimated signal and noise spectra up to the minimum averaging period constrained by the diffusion correction (∼ 2.5 years), we find an SNR of $\stackrel{\mathrm{‾}}{F}\sim \mathrm{0.2}$ (correlation r∼0.4, Fig. 4a), which is, as expected, slightly higher than the published estimates, since we corrected for the effect of time uncertainty, as well as due to the effect of the slightly different underlying averaging periods. Our results further explain the strong agreement between individual ice cores on glacial–interglacial timescales. Averaging to multi-decadal resolution results in a correlation between a single core and the signal of around 0.6 (Fig. 4a), which rises to ∼0.7 for the available centennial averaging periods. Moreover, the steady increase in SNR on the analysed timescales (Fig. 3) might indicate a further increase in SNR towards even longer timescales. However, clearly our results also underline again that for assessing Holocene temperature variability on timescales shorter than multi-decadal, the averaging of records is essential to increasing their representativity (Fig. 4a).

For WAIS, the higher SNR at interannual timescales compared to DML (average SNR between 5 and 10-year periods of 0.43 for WAIS compared to 0.13 for DML; Fig. 3) is consistent with the general notion that higher accumulation rates (on average, ∼3 times as high as at DML; Table 1) result in higher SNR (Fisher et al., 1985; Steen-Larsen et al., 2011). However, one would expect, to a first approximation, a similar increase in SNR with timescale in both regions. In contrast to this expectation, the WAIS results show strongly different timescale dependencies with the tendency of a decrease in signal power and an increase in noise level on longer timescales (Fig. 2), resulting in an overall decrease in SNR (Fig. 3). This explains the only slightly higher correlations at interannual timescales of a stacked isotope record with the common signal at WAIS compared to DML (Fig. 4), since the correlation is based on the integrated value of the SNR (Eqs. 6 and 7). If they are correct, these findings suggest that there is, in general, no simple scaling of the SNR with accumulation rate but that additional effects need to be involved.

The shape of the signal and noise spectra at subdecadal timescales is sensitive to the diffusion and time uncertainty corrections (Fig. 2), which could thus contribute to the discrepant SNR estimates. While the diffusion correction for DML led to a consistent white-noise level of both noise estimates (Fig. 2b) at these timescales, lending support to our approach, the WAIS noise spectrum keeps decreasing towards higher frequencies, even after applying the diffusion correction (Fig. 2d). This result could be caused either by a too weak diffusion correction or by an overestimation of the time uncertainty leading to an excessive reduction in high-frequency noise levels that cannot be compensated by the low diffusion correction. To test both hypotheses, we varied the strengths of the diffusion and time uncertainty corrections. This has different impacts on the estimated spectra: while the diffusion correction equally applies to both raw signal and noise spectra, the time uncertainty correction has a proportional influence on the signal spectrum but only an indirect influence on the noise spectrum (Eq. 4). Indeed, halving the time uncertainty would increase the noise spectrum close to the cut-off frequency only by a maximum factor of ∼1.4, and this therefore cannot reconcile the remaining decrease in the corrected noise spectrum. By contrast, an overall doubling of the diffusion length for all WAIS records would amplify the noise spectrum much more strongly with a maximum factor of ∼4 at the cut-off frequency, leading to interannual noise levels similar to those observed for centennial periods. However, a much stronger diffusional smoothing at WAIS compared to the expectation seems implausible given that the same physical mechanisms of the diffusion process should be valid for East Antarctica as well as West Antarctica, and no anomalously high diffusional smoothing has been observed for the upper part of the WAIS Divide ice core (Jones et al., 2017). Additionally, a much stronger diffusion correction would also imply a much steeper decrease in the signal spectrum towards longer timescales (Fig. 2c), contrary to the expectation. We conclude that there is no obvious reason for the applied correction approach to strongly over- or underestimate the WAIS signal and noise spectra on the interannual timescales. Moreover, by also taking into account the shape of the spectra for the signal and noise on the longer timescales (periods from ∼30 to 100 years), our results might therefore indicate a true increase in local variability at the WAIS sites across the timescales studied and a close-to-constant, or even decreasing, coherent signal variability.

These results suggest firstly that an additional noise process, apart from stratigraphic noise and precipitation intermittency, must contribute to the noise spectrum towards longer timescales. Secondly, the shape of the signal spectrum either implies a nearly white-noise temperature signal or some process that destroys the coherence of the large-scale temperature field on longer timescales upon recording by the firn-core isotope records. Since there is no obvious reason for a fundamentally different Holocene climate variability in West compared to East Antarctica, the second possibility seems more likely.

WAIS isotope data have been reported to covary with local temperatures, but also with the large-scale atmospheric circulation and the sea-ice cover of the adjacent oceans (Küttel et al., 2012; Noone and Simmonds, 2004; Steig et al., 2013). A pronounced regional imprint on the recorded isotope variability is also suggested from the spatial differences in intercore correlations for an extended set of the available WAIS firn cores (Küttel et al., 2012). In particular, cores east of the WAIS Divide display strong differences in the recorded isotope variability despite their rather small spatial separation (Gregory and Noone, 2008), which is suggested to be related to slow variations in the dominant circulation patterns. While such differences have not explicitly been reported for the core sites west of the divide studied here, we hypothesise that signatures of sea-ice variations or meridional inflow could affect the isotopic composition at the individual firn-core sites to different extents. This is motivated by model-based results (Noone and Simmonds, 2004) linking an increase in sea-ice cover to an earlier ascent of long-range transported air masses to the continent, reducing the influence of local air from the surface, while less sea ice inhibits this early ascent and allows more mixing of surface air. Variations in sea ice could thus influence the isotopic composition of air masses by controlling the influence of local moisture, which might only be recorded by certain WAIS cores depending on the elevation of the air mass transport and characteristics of the core positions, such as their surrounding topography, elevation, or distance to the coast. Decadal trends or slower changes in these variations could then destroy the recording of the large-scale coherent temperature field by the firn cores and thus cause a loss in spectral signal power in the isotope record towards longer timescales and an increase in the noise level (Fig. 2c, d). In contrast, the East Antarctic Plateau including the DML firn-core sites is higher in elevation and might be more shielded from marine influences by the steep topographic slopes leading to a more coherent signal. This might also explain, besides differences in core quality, the rather low agreement among deep West Antarctic cores on millennial timescales compared to East Antarctic cores (WAIS Divide Project Members, 2013). However, since our inferences are speculative, further studies are needed to help disentangle the role of variations in West Antarctic temperature, atmospheric circulation and sea ice on the recorded isotope variability.

5 Conclusions

We presented a simple spectral method to separate the variations recorded by isotope records into a local (noise) and a common (i.e. spatially coherent) signal component. We applied this method to firn cores from the East Antarctic Dronning Maud Land (DML) and the West Antarctic Ice Sheet (WAIS) to estimate, for the first time, the isotopic signal-to-noise ratio (SNR) as a function of timescale. This is of fundamental interest for interpreting isotope records obtained from individual ice cores, since it provides an upper limit on the SNR of the temperature signal recorded by the cores. For DML, the SNR at the interannual timescale is very low, but it steadily increases with timescale reaching values above 0.5 for multi-decadal and slower variations. Therefore, only on these timescales should isotope changes from individual cores be interpreted in terms of regional climate variations. For WAIS, the results are counterintuitive. On interannual to decadal timescales, the estimated SNR is higher than 0.5, which would support the regional climate interpretation of the isotope records. For longer timescales, however, the estimated SNR decreases, suggesting that local variations start to dominate the recorded isotope variability.

Our method further allows the estimation of the power spectra of the coherent isotope signal. For DML, the spectra of single cores largely resemble white noise. In contrast, the derived signal spectrum shows increasing variability towards longer timescales. Such an increase is also observed in instrumental temperature records and other climate proxies. The marked difference between the raw interpretation of single cores – as it is usually done – and the signal spectra derived from the core array demonstrates the relevance of the signal and noise separation. The interpretation of the WAIS isotope signal is more challenging, since the signal shows a close-to-constant spectral power even after applying our method. We speculate that this might be due to atmospheric circulation variations that create a local imprint at the different firn-core sites. This might prevent a coherent recording of the large-scale atmospheric temperature field. To test this hypothesis, we suggest to analyse firn-core arrays as a function of the average separation distance between the individual core sites within each array. This could allow us to investigate whether the stable-isotope data record a stronger coherent signal on a regional scale (e.g. ∼1–10 km) than on the larger scales analysed here. A similar approach in DML would also help to better constrain the spatial correlation scales of the precipitation intermittency.

We conclude that the pronounced differences seen between East and West Antarctica could thus be related to the differences in the topographic settings and the different marine influences (Noone and Simmonds, 2004; WAIS Divide Project Members, 2013). Attempts to reconstruct the temperature variability in these regions based on firn and ice cores should thus not only aim to average as many records as possible, but also consider the spatial coherence of the circulation and precipitation patterns in order to establish an optimal strategy for averaging or obtaining new firn-core isotope records. Additionally, extending our analyses to data derived from non-isotope-based temperature proxies could give further insights into the true spectral shape of temperature variability in Antarctica.

Finally, our approach of separating signal and noise from a set of spatially distributed records is also applicable beyond Antarctic ice cores. The challenge of low and timescale-dependent SNR is common to many high-resolution climate archives, and the number of nearby core sites is continuously increasing. Therefore, we envision our approach to better constrain the reliability of proxy records as a function of timescale in general and to allow a significant improvement of our knowledge of past climate variability.

Appendix A: Theoretical spectra for site-specific diffusion

We derive the effect of site-specific diffusion on the estimates of the mean spectrum and the spectrum of the stack given a core array consisting of n isotope records.

The spectrum of a time series x(t) is given by the absolute square of its Fourier transformation, $\mathcal{X}=\mathcal{F}\left(x\right){\mathcal{F}}^{*}\left(x\right)$, where ℱ^*(x) denotes the complex conjugate of ℱ(x). Given an isotope record at site i, $x_i\left(t\right)=c\left(t\right)+{\mathit{\epsilon }}_i\left(t\right)$, where c(t) denotes the common signal, and ε_i(t) is noise, the time series after diffusion is $x_i^{\prime }\left(t\right)=g_i\star \left(c\left(t\right)+{\mathit{\epsilon }}_i\left(t\right)\right)$. Here, ⋆ denotes convolution and g_i is a Gaussian convolution kernel, the standard deviation of which is the diffusion length σ_i that is a function of depth (or, equivalently, time since deposition) and depends on site i due to its dependence on local temperature, atmospheric pressure and accumulation rate (Appendix B). The spectrum of x^′ is then

$$\begin{array}{ll}{\displaystyle }{\mathcal{X}}_i& {\displaystyle }=\mathcal{F}\left(g_i\right)\mathcal{F}(c+{\mathit{\epsilon }}_i)\cdot {\mathcal{F}}^{*}\left(g_i\right){\mathcal{F}}^{*}(c+{\mathit{\epsilon }}_i)\\ \text{(A1)}& {\displaystyle }& {\displaystyle }={\mathcal{G}}_i\left(\mathcal{C}+{\mathcal{N}}_i\right),\end{array}$$

with $\mathcal{C}=\mathcal{F}\left(c\right){\mathcal{F}}^{*}\left(c\right)$, $\mathcal{N}=\mathcal{F}\left({\mathit{\epsilon }}_i\right){\mathcal{F}}^{*}\left({\mathit{\epsilon }}_i\right)$ and where the linear site-specific diffusion transfer function was defined as ${\mathcal{G}}_i:=\mathcal{F}\left(g_i\right){\mathcal{F}}^{*}\left(g_i\right)$¹. Assuming that the statistical properties of the individual noise terms are the same for all sites, 𝒩_i≡𝒩, the mean across n such spectra is

$$\begin{array}{}\text{(A2)}& \mathcal{M}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathcal{X}}_i=\stackrel{\mathrm{‾}}{\mathcal{G}}\left(\mathcal{C}+\mathcal{N}\right)\end{array}$$

with the average diffusion transfer function $\stackrel{\mathrm{‾}}{\mathcal{G}}:=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n{\mathcal{G}}_i$.

The average in the time domain of n independent noise signals subject to site-specific diffusion, $y_{\mathrm{noise}}\left(t\right)=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n{g}_i\star {\mathit{\epsilon }}_i\left(t\right)$, has a spectrum of

$$\begin{array}{ll}{\displaystyle }{\mathcal{S}}_{\mathrm{noise}}& {\displaystyle }=\mathcal{F}\left(y_{\mathrm{noise}}\right){\mathcal{F}}^{*}\left(y_{\mathrm{noise}}\right)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n^{\mathrm{2}}}}\sum _i\sum _j\mathcal{F}\left(g_i\right){\mathcal{F}}^{*}\left(g_j\right)\mathcal{F}\left({\mathit{\epsilon }}_i\right){\mathcal{F}}^{*}\left({\mathit{\epsilon }}_j\right).\end{array}$$

We have $\mathcal{F}\left({\mathit{\epsilon }}_i\right){\mathcal{F}}^{*}\left({\mathit{\epsilon }}_j\right)\ne \mathrm{0}$ only for i=j, and hence

$$\begin{array}{}\text{(A3)}& {\mathcal{S}}_{\mathrm{noise}}={\displaystyle \frac{\mathrm{1}}{n}}\stackrel{\mathrm{‾}}{\mathcal{G}}\mathcal{N},\end{array}$$

which shows that averaging n independent noise records reduces the spectral power by 1∕n, as expected.

By contrast, the spectrum of the average of n coherent signals subject to site-specific diffusion, $y_{\mathrm{coh}}\left(t\right)=\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n{g}_i\star c\left(t\right)$, is

$$\begin{array}{ll}{\displaystyle }{\mathcal{S}}_{\mathrm{coh}}& {\displaystyle }=\mathcal{F}\left(y_{\mathrm{coh}}\right){\mathcal{F}}^{*}\left(y_{\mathrm{coh}}\right)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n^{\mathrm{2}}}}\mathcal{F}\left(c\right){\mathcal{F}}^{*}\left(c\right)\sum _i\sum _j\mathcal{F}\left(g_i\right){\mathcal{F}}^{*}\left(g_j\right)\\ \text{(A4)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n}}\mathcal{C}\left(\stackrel{\mathrm{‾}}{\mathcal{G}}+{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i\ne j}\mathcal{F}\left(g_i\right){\mathcal{F}}^{*}\left(g_j\right)\right).\end{array}$$

For small differences in the diffusion lengths between the sites, the transfer functions are approximately the same and we can simplify the second term in brackets in Eq. (A4),

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\mathrm{1}}{n}}\sum _{i\ne j}\mathcal{F}\left(g_i\right){\mathcal{F}}^{*}\left(g_j\right)\approx {\displaystyle \frac{\mathrm{1}}{n}}\sum _{i\ne j}\mathcal{F}\left(g_i\right){\mathcal{F}}^{*}\left(g_i\right).\end{array}$$

From rearranging the summation terms we find that this is identical to $(n-\mathrm{1})\stackrel{\mathrm{‾}}{\mathcal{G}}$, and Eq. (A4) becomes

$$\begin{array}{}\text{(A6)}& {\mathcal{S}}_{\mathrm{coh}}\approx \stackrel{\mathrm{‾}}{\mathcal{G}}\mathcal{C}.\end{array}$$

In this approximation, the average over diffused coherent signals does not reduce the spectral power, as one would expect. We tested the approximation for our core arrays by comparing simulation results with surrogate data between the cases of independent and coherent noise. We find that Eq. (A6) is a reasonable approximation for the full expression Eq. (A4). However, we note that Eq. (A4) has slightly less power on the high-frequency end compared to Eq. (A6), since the site-specific diffusion destroys part of the coherence between the sites.

Appendix B: Estimates of the transfer functions for diffusion and time uncertainty

To estimate the transfer functions for diffusion ($\stackrel{\mathrm{‾}}{\mathcal{G}}$) and time uncertainty (Φ), we simulate for each core array n surrogate records of the same time duration as the original isotope records. For each record, we simulate the effects of diffusion and time uncertainty and then calculate in each case and for each core array the spectrum of the average record (𝒮_diffusion or 𝒮_time). At the same time, the spectrum of the average record of the initial surrogate data (𝒮₀) is calculated without applying diffusion or time uncertainty. The entire process is repeated k times and the resulting spectra are averaged to reduce the spectral uncertainty. We report $\stackrel{\mathrm{‾}}{\mathcal{G}}=\langle {\mathcal{S}}_{\mathrm{diffusion}}\rangle /\langle {\mathcal{S}}_{\mathrm{0}}\rangle$ and $\mathrm{\Phi }=\langle {\mathcal{S}}_{\mathrm{time}}\rangle /\langle {\mathcal{S}}_{\mathrm{0}}\rangle$ as the transfer functions, where 〈⋅〉 denotes the average over the k simulations where we use k=10⁵.

For the diffusion simulations, each surrogate record is smoothed by convolution with a Gaussian kernel of width σ, where σ is the diffusion length, measured in years, sensitive to ambient temperature, pressure and the density of the firn (Whillans and Grootes, 1985) and therefore site-specific and a function of depth. We treat the dependency on density according to Gkinis et al. (2014) with diffusivity after Johnsen et al. (2000); firn density is modelled according to the Herron–Langway model (Herron and Langway, 1980). Surface air pressure is calculated from the barometric height formula using the local elevation and firn temperature. For the surface firn density we assume a constant value of 340 kg m⁻³ for all sites. The range of input parameters and their references are given in Table 1. Missing temperature information for the cores FB9807, FB9811 and FB9815 are filled with the temperatures from the nearby cores B32 (1 km), FB9812 (19 km) and FB9805 (28 km) respectively. We simulate diffusion lengths measured in centimetres for core lengths of 500 m at a resolution of 1 cm, convert the values into years and interpolate them onto a regular time axis at a resolution of 0.5 years. This higher temporal resolution compared to the annual target resolution of the isotope records is necessary for numerical stability of the diffusion results at the high-frequency end. The diffusion transfer functions are then interpolated in frequency space onto the frequency axis of the isotope records. Simulated diffusion lengths at a depth corresponding to 200 years after deposition of the snow range, across the firn-core sites, from 0.6 to 1.3 years (8.3–11.4 cm) for DML and from 0.3 to 0.6 years (8.9–11.4 cm) for WAIS; for DML and 1000 years after deposition of the snow, diffusion lengths range from 1.1 to 1.6 years (8.2–9.1 cm).

To simulate the effect of time uncertainty, we apply the modelling approach by Comboul et al. (2014). This model represents an unconstrained process, meaning that the time uncertainty monotonically increases with modelled age, thus deeper in the core. To account for volcanic tie points, where the chronologies of the isotope records are fixed (neglecting time uncertainties of the volcanic chronologies themselves), we modify the Comboul approach by forcing the model back to zero time uncertainty at the reported volcanic tie points (constrained process) in the form of a Brownian Bridge process. In contrast to diffusional smoothing, it is not a priori clear whether time uncertainty acts as a linear transfer function on the spectrum of the average record. We tested this by using different spectral characteristics of the surrogate data, adopting the cases of white noise as well as power laws with spectral slopes of β=0.5, 1, 1.5 and 2, and found that the modelled effect of time uncertainty is insensitive to the spectral shape of the input data and indeed acts as a linear transfer function. The Comboul model includes as parameters the rates of missing and double-counted annual layers. We assume equal rates for these processes, which we set so that the yielded maximum time uncertainties (maximum standard deviation over simulated chronologies) match the reported time uncertainties of the isotope records. The reported value of the DML records (∼2 % of the time interval to the nearest tie point, Graf et al., 2002a) erroneously implies that the uncertainty increases linearly with time (Comboul et al., 2014). Here, as a best guess, we choose a process rate that yields a maximum standard deviation of 2 % for a constrained process of duration equal to the mean interval length between the tie points. In summary, we use process rates of γ=0.013 for the DML and of γ=0.027 for the WAIS cores.

The resulting transfer functions used for the main results are shown in Fig. B1. Diffusion shows a stronger smoothing for the longer records of DML2 than for DML1 as it had more time to act. WAIS shows less smoothing compared to DML1 due to the higher accumulation rates. The influence of diffusion is negligible for frequencies below the decadal period for all core arrays. The effect of time uncertainty is negligible for frequencies below the decadal period for DML1 and WAIS, but only for frequencies below the 20-year period for DML2 due to the longer distances between the volcanic tie points.

Figure B1 Estimates of the spectral transfer functions for the effects of site-specific diffusion (a) and time uncertainty (b) for the three studied core arrays, DML1 (black), DML2 (red) and WAIS (blue). The estimates are based on simulations with surrogate data as explained in the text. Both effects are largely negligible beyond the decadal period. The horizontal dashed line in (a) at 0.5 denotes the transfer function value, which was chosen in order to constrain the diffusion correction (see main text). It corresponds to minimum time periods for the analysis of the respective spectral data that are indicated by the vertical dashed lines.

Appendix C: ERA-Interim temperature and precipitation field decorrelation scales

To estimate the present-day spatial decorrelation scales of the atmospheric temperature and the precipitation fields in our study regions, we resort to ERA-Interim reanalysis data (Dee et al., 2011; European Centre for Medium-Range Weather Forecasts, 2018), since direct observations are sparse. ERA-Interim is in general regarded to be the most reliable reanalysis product, both for temperature and precipitation in Antarctica; the reanalysis temperatures are, despite significant biases in mean values, well correlated with station observations at monthly to interannual timescales (Bracegirdle and Marshall, 2012; Jones and Lister, 2015), and the reanalysis produces the closest match to Antarctic snowfall, regarding both mean values and seasonal-to-interannual variability (Bromwich et al., 2011; Palerme et al., 2017).

To obtain decorrelation scales, we calculate the correlations between reanalysis time series from the grid box closest to the location of the EDML site (75^∘ S, 0^∘ W) with the respective time series of all other grid boxes of the data set south of 60^∘ S. Using an exponential fit to the correlations of the form $\exp (-d/\mathrm{\Delta })$, where d is the distance between EDML and a particular grid box, yields the decorrelation scale through the fit parameter Δ. The same approach is applied to the WDC site (79.5^∘ S, 112^∘ W).

Figure C1 Present-day temperature decorrelation across DML and WAIS. Shown are the correlations (grey dots) of ERA-Interim annual-mean temperatures (Dee et al., 2011) at (a) EDML and (b) WDC with all other grid boxes below 60^∘ S. Black lines show an exponential fit to the data; dashed lines indicate the point at which the correlations in the model have dropped to 1∕e – at 1202 (1232) km for EDML (WDC). Coloured shading indicates the range of site distances of the studied firn-core arrays (Table 1); the average correlation across each range is marked by the filled diamond – average correlation of 0.87 (0.94) for EDML (WDC).

For annual mean temperatures, the analysis yields decorrelation scales of around 1200 km for both regions (Fig. C1). The decorrelation scales of precipitation intermittency are expected to depend on the decorrelation scales of the seasonal precipitation distribution rather than on the total annual precipitation amount (Persson et al., 2011). However, our analysis suggests that the precipitation decorrelation scales are actually more or less insensitive to the chosen variable, with estimates between roughly 300 and 500 km for both regions for a variety of chosen fields (total precipitation amount, the difference between summer (months DJF) and winter (months JJA) precipitation amounts, the fraction of the summer amount over the winter (winter over summer) precipitation amount and the fraction of the summer amount (winter) over the total precipitation amount).