Regional climate signal vs. local noise: a two-dimensional view of water isotopes in Antarctic firn at Kohnen Station, Dronning Maud Land

Abstract

In low-accumulation regions, the reliability of${\mathit{\delta }}^{18}\mathrm{O}$-derived temperature signals from ice cores within the Holocene is unclear, primarily due to the small climate changes relative to the intrinsic noise of the isotopic signal. In order to learn about the representativity of single ice cores and to optimise future ice-core-based climate reconstructions, we studied the stable-water isotope composition of firn at Kohnen Station, Dronning Maud Land, Antarctica. Analysing${\mathit{\delta }}^{18}\mathrm{O}$ in two 50 $\mathrm{m}$ long snow trenches allowed us to create an unprecedented, two-dimensional image characterising the isotopic variations from the centimetre to the 100-metre scale. Our results show seasonal layering of the isotopic composition but also high horizontal isotopic variability caused by local stratigraphic noise. Based on the horizontal and vertical structure of the isotopic variations, we derive a statistical noise model which successfully explains the trench data. The model further allows one to determine an upper bound for the reliability of climate reconstructions conducted in our study region at seasonal to annual resolution, depending on the number and the spacing of the cores taken.

Introduction

Ice cores obtained from continental ice sheets and glaciers are a key climate archive. They store information on past changes in temperature in the form of stable water isotopes , in greenhouse gas concentrations via trapped air and in many other parameters such as accumulation rates e.g. or aerosols e.g.. Analysis of the isotope ratios recorded in single deep ice cores provided milestones in the palaeo-climate research, including the investigation of glacial–interglacial climate changes and the existence of rapid climate variations within glacial periods .

In contrast to the coherent view established from polar ice cores for millennial and longer timescales, the reliability of single ice cores as archives of the Holocene climate evolution is less clear . The small amplitude of changes and the aim to reconstruct climate parameters at high temporal resolution poses a challenge to the interpretation of ice-core signals. This is especially true for low-accumulation sites, defined here for accumulation rates below $100{\text{mm w.e. yr}}^{-\mathrm{1}}$, which holds for large parts of the East Antarctic Plateau. There, the non-climate noise – that part of the isotopic record which cannot be interpreted in terms of temperature variations on regional or larger scales, hence including any meteorological, pre- and post-depositional effects that additionally influence the isotopic composition – may often be too high to accurately extract a climatic temperature signal . Despite the challenges, quantifying the Holocene polar climate variability is the key foundation to determine the range of possible future climate changes e.g.and references therein as well as to test the ability of climate models in simulating natural climate variability .

The quantitative estimation of climate variability from proxy data therefore requires an understanding of the non-climate noise in order to separate it from the climate signal e.g.. Several mechanisms influence the isotopic composition of snow prior to and after its deposition onto the ice sheet. On larger spatial scales, non-climate variability may be introduced by different moisture pathways and source regions e.g. as well as spatial and temporal precipitation intermittency . Irregular deposition caused by wind and surface roughness along with spatial redistribution and erosion of snow is a major contribution to non-climate variance on smaller spatial scales (“stratigraphic noise”, ). Wind scouring can additionally remove entire seasons from the isotopic record . Vapour exchange with the atmosphere by sublimation–condensation processes can influence the isotopic composition of the surface layers; diffusion of vapour into or out of the firn driven by forced ventilation may represent an additional component of post-depositional change. Finally, diffusion of water vapour through the porous firn smoothes isotopic variations from seasonal to inter-annual or longer timescales, depending on the accumulation rate .

In the last two decades, a number of studies analysed the representativity of single ice cores in recording sub-millennial climate changes. One well-studied region is the low-accumulation ($\sim 40\text{–}90{\text{mm w.e. yr}}^{-\mathrm{1}}$, ) Dronning Maud Land (DML) on the East Antarctic Plateau. Here, found low signal-to-noise variance ratios ($F$) in 200-$\mathrm{year}$-long firn-core records for oxygen isotopes ($F=0.14$) and accumulation rates ($F=0.04$), implying that the climate signal content of a single core is much smaller than the noise level ($14$ and $\mathrm{4}\mathit{\%}$, respectively). On a similar timescale, detected no relationship in electrical properties apart from volcanic imprints between firn cores. Similarly, high-resolution records of chemical trace species from three shallow ice cores showed a lack of inter-site correlation on decadal timescales. These results were supported by process studies comparing observed and simulated snow-pit isotope data . Whereas the model–data comparison was successful for coastal high-accumulation regions of DML ($400{\text{mm w.e. yr}}^{-\mathrm{1}}$), it largely failed on the dryer East Antarctic Plateau ($70{\text{mm w.e. yr}}^{-\mathrm{1}}$). This dependency between accumulation rate and signal-to-noise ratio was further demonstrated in studies across the Antarctic continent . From high-accumulation ($140\text{–}520{\text{mm w.e. yr}}^{-\mathrm{1}}$) Greenland ice cores, estimated signal-to-noise ratios clearly larger than 1.

Despite this large body of literature, quantitative information about the signal-to-noise ratios and the noise itself is mainly limited to correlation statistics of nearby cores. While a relatively good understanding of stratigraphic noise exists in Arctic records , this does not apply to low-accumulation regions of Antarctica where the accumulated snow is considerably reworked in and between storms .

Here, we provide a direct visualisation and analysis of the signal and noise in an East Antarctic low-accumulation region ($64{\text{mm w.e. yr}}^{-\mathrm{1}}$) by an extensive two-dimensional sampling of the firn column in two 50 $\mathrm{m}$ long snow trenches. Our approach, for the first time, offers a detailed quantitative analysis of the spatial structure of isotope variability on a centimetre to 100-metre scale. This is a first step towards a signal and noise model to enable quantitative reconstructions of the climate signal and their uncertainties from ice cores.

Data and methods

Near Kohnen Station in the interior of Dronning Maud Land, close to the EPICA deep ice core drilling site (EDML, $-75.0$${}^{\circ }$ S, $0.1$${}^{\circ }$ E; altitude $2892\mathrm{m}$ a.s.l.; mean annual temperature $-44.5{}^{\circ }\mathrm{C}$; mean annual accumulation rate $64{\text{mm w.e. yr}}^{-\mathrm{1}}$; ), two $1.2\mathrm{m}$ deep and approximately $45\mathrm{m}$ long trenches in the firn, named T1 and T2, were excavated during the austral-summer field season $2012/2013$ using a snow blower. Each trench was aligned perpendicularly to the local snow-dune direction. The horizontal distance between the starting points of T1 and T2 was $415\mathrm{m}$.

An absolute height reference was established using bamboo poles by adjusting their heights above ground with a spirit level. A control measurement with a laser level yielded in each snow trench a vertical accuracy better than $\mathrm{2}\mathrm{cm}$. No absolute height reference between the two trenches could be established, but, based on a stacked laser level measurement, the vertical difference between the trenches was estimated to be less than $20\mathrm{cm}$.

Both trenches were sampled for stable-water-isotope analysis with a vertical resolution of $\mathrm{3}\mathrm{cm}$. In T1, 38 profiles were taken at variable horizontal spacings between $0.1$ and $\sim \mathrm{\hspace{0.17em}2.5}\mathrm{m}$. In T2, due to time constraints during the field campaign, only four profiles at positions of 0.3, 10, 30 and 40 m from the trench starting point were realised. The sampling of each trench was completed within 24 h. All firn samples ($N=1507$) were stored in plastic bags and transported to Germany in frozen state. Stable isotope ratios were analysed using cavity ring-down spectrometers (L2120i and L2130i, Picarro Inc.) in the isotope laboratories of the Alfred Wegener Institute (AWI) in Potsdam and Bremerhaven. The isotope ratios are reported in the usual delta notation in per mil (‰) as $$\mathit{\delta }=\left({\displaystyle \frac{R_{\text{sample}}}{R_{\text{reference}}}},-,\mathrm{1}\right)\times {10}^{\mathrm{3}},$$ where $R_{\text{sample}}$ is the isotopic ratio of the sample (${}^{18}\mathrm{O}{/}^{16}\mathrm{O}$) and $R_{\text{reference}}$ that of a reference. The isotopic ratios are calibrated with a linear three-point regression analysis using in-house standards at the beginning of each measurement sequence, where each standard has been calibrated to the international V-SMOW/SLAP scale. Additionally, a linear drift-correction scheme and a memory-correction scheme (adapted from ) is applied, using three repeated measurements per sample. The analytical precision of the calibrated and corrected ${\mathit{\delta }}^{18}\mathrm{O}$ measurements is assessed by evaluating standards in the middle of each measurement sequence. This yields a mean combined measurement uncertainty of $0.09\text{‰}$ (RMSD). The ${\mathit{\delta }}^{18}\mathrm{O}$ trench data are archived under 10.1594/PANGAEA.861675 .

For the analysis of the measurements, we set up two coordinate systems for each trench (Fig. ). Surface coordinates refer to a local, curvilinear system with the horizontal axis along the surface height profile and the vertical axis denoting the firn depth below the local surface. Absolute coordinates adopt the mean surface height as a reference for a straight horizontal axis, completed by an absolute depth scale.

Coordinate systems used for the analysis of the trench isotope data: (1) a curvilinear coordinate system $(\mathit{\xi },\mathit{\zeta })$ (blue dashed lines, surface coordinates) with horizontal axis along the surface height profile and vertical axis denoting the depth below the local surface; (2) a Cartesian system $(x,z)$ (black lines, absolute coordinates) defined by the mean surface height.

[Figure omitted. See PDF]

Results

Trench isotope records

The firn samples obtained from trench T1 provide a two-dimensional image of the ${\mathit{\delta }}^{18}\mathrm{O}$ structure of the upper $\sim \mathrm{1}\mathrm{m}$ of firn on a horizontal scale of $\sim \mathrm{\hspace{0.17em}50}\mathrm{m}$ (Fig. a).

(a) The two-dimensional ${\mathit{\delta }}^{18}\mathrm{O}$ data set of trench T1 displayed on absolute coordinates. The solid black line shows the surface height profile, the long-dashed black line the mean surface height. Sampling positions are marked by the black dots above. White gaps indicate missing data. (b) The stratigraphy of trench T1 expressed as the seasonal layer profiles tracking the local ${\mathit{\delta }}^{18}\mathrm{O}$ extrema as explained in the text.

[Figure omitted. See PDF]

The surface height profile of the trench reflects the typical snow topography of the sampling region characterised by small-scale dunes with their main ridges elongated parallel to the mean wind direction . Trench T1 features one prominent dune located between 25 and 40 m, accompanied by a dune valley between 8 and 18 m, and some smaller-scale height variations. The peak-to-peak amplitude of the large dune undulation is $\sim 10\mathrm{cm}$; the entire height variations exhibit a standard deviation (SD) of $2.9\mathrm{cm}$.

Overall, the trench ${\mathit{\delta }}^{18}\mathrm{O}$ record shows a diverse picture. The delta values in T1 (Fig. a) span a range from $-54$ to $-34\text{‰}$ with a mean of $-44.4\text{‰}$ (SD $3.1\text{‰}$). A similar range of $-50$ to $-38\text{‰}$ is observed in T2 (Fig. ) with a mean of $-44.0\text{‰}$ (SD $2.7\text{‰}$). We can identify 8 to 10 alternating layers of enriched and depleted isotopic composition in the T1 record. The uppermost layer (first $\mathrm{6}\mathrm{cm}$ relative to the surface) essentially shows enriched (mean of $-42.7\text{‰}$) but also strongly variable ${\mathit{\delta }}^{18}\mathrm{O}$ values between $-54$ and $-34\text{‰}$ (SD $4.4\text{‰}$), thus already covering the range of the entire trench record. Stronger enrichment tends to be found in the valleys; however, the limited data do not allow one to conclude whether this is a general feature. In an absolute depth of $\mathrm{5}\text{–}20\mathrm{cm}$, a band of generally more depleted ${\mathit{\delta }}^{18}\mathrm{O}$ values is found exhibiting less horizontal variability compared to the first layer with a range of $-54$ to $-45\text{‰}$ (mean $-48.5\text{‰}$, SD $1.9\text{‰}$). The layering appears strongly perturbed in the depth of $\sim 60\text{–}100\mathrm{cm}$ for profile positions $<\mathrm{\hspace{0.17em}30}\mathrm{m}$. Here, a broad and diffuse region of rather constant ${\mathit{\delta }}^{18}\mathrm{O}$ values around $-40\text{‰}$ is observed, together with a prominent, $20\mathrm{cm}$ thick feature of high delta values between $18$ and $28\mathrm{m}$.

The four ${\mathit{\delta }}^{18}\mathrm{O}$ profiles obtained from trench T2 (Fig. ) show similar features as trench T1. We can identify roughly five cycles in each profile. However, the profiles diverge considerably at depths of $50\text{–}90\mathrm{cm}$, which coincides with the region of strong perturbations identified in T1.

To further analyse the isotopic layering, we determine the pronounced local maxima and minima of each T1 ${\mathit{\delta }}^{18}\mathrm{O}$ profile and visually assign summer and winter to the depths of these extrema. This results in consecutive horizontal curves tracing the vertical positions of seasonal extrema along the trench (seasonal layer profiles, Fig. b). Assuming that respective isotopic extrema occur at the same point in time (summer/winter), the seasonal layer profiles reflect the surface height profile for a given season. However, considering the highly variable isotopic composition observed at the current trench surface (Fig. a), this is a rough approximation and the seasonal layer profiles will likely overestimate the past surface height profiles. Nevertheless, the vertical undulations of the layer profiles show peak-to-peak amplitudes of $\mathrm{6}\text{–}24\mathrm{cm}$ (average SD $3.7\mathrm{cm}$), comparable to the present surface undulations, and the layers are vertically separated by approximately $20\mathrm{cm}$, in accord with the local mean annual accumulation rate of snow ($64{\text{mm w.e. yr}}^{-\mathrm{1}}$) and the mean firn density measured in trench T1 (${\mathit{\rho }}_{\text{firn}}=340{\text{kg m}}^{-\mathrm{3}}$). To study the similarity between the seasonal layer profiles and the present surface height profile, we calculate the standard deviation of their height differences (${\text{SD}}_{\text{surface}}$, hence, the SD of each layer profile on surface coordinates). This is compared to the standard deviation of the layer profiles on absolute coordinates (${\text{SD}}_{\text{horiz}}$). We find that the first layer profile closely follows the present surface (${\text{SD}}_{\text{horiz}}-{\text{SD}}_{\text{surface}}=1.8\mathrm{cm}$). For the second layer profile, the link with the surface is weaker (${\text{SD}}_{\text{horiz}}-{\text{SD}}_{\text{surface}}=1.5\mathrm{cm}$), and the layer profiles below $20\mathrm{cm}$ are on average horizontally aligned (${\text{SD}}_{\text{horiz}}-{\text{SD}}_{\text{surface}}=-0.6\mathrm{cm}$). This can be explained by an annual reorganisation of the stratigraphy so that aligning the isotopic variations on absolute coordinates is on average more appropriate than the alignment according to a specific surface height profile. The positive autocorrelation with a decorrelation length of $\sim \mathrm{\hspace{0.17em}6}\mathrm{cm}$ that is found from the vertical T1 ${\mathit{\delta }}^{18}\mathrm{O}$ variations after subtraction of the mean trench profile is consistent with this hypothesis.

Due to the on average horizontal stratigraphy of the isotopic composition in the larger part of the trench record all further plots and calculations will be evaluated on absolute coordinates.

The four ${\mathit{\delta }}^{18}\mathrm{O}$ profiles obtained from trench T2 displayed on absolute coordinates.

[Figure omitted. See PDF]

Histogram of all possible pairwise correlations ($N=152$) between single profiles of trench T1 and single profiles of trench T2. Displayed are the maximum correlations allowing vertical shifts of the T2 profiles of up to $\pm 12\mathrm{cm}$. Shown in red is the correlation between the mean profiles of T1 and T2 (Fig. ).

[Figure omitted. See PDF]

Single-profile representativity

The isotope record of trench T1 (Fig. a) allows the quantification of the horizontal isotopic variability of the snow and firn column in our study region. We observe considerable horizontal variability with a mean variance of ${\mathit{\sigma }}_{\text{h,T1}}^{\mathrm{2}}\simeq 5.9(\text{‰}{)}^{\mathrm{2}}$, directly affecting the representativity of single trench profiles. To mimic the potential result obtained from correlating two snow pits taken at a distance of $500\mathrm{m}$, similarly done in many firn-core studies e.g., we calculate the pairwise Pearson correlation coefficient between single profiles of T1 and single profiles of T2. We account for potential surface undulations between the trenches by allowing bin-wise vertical shifts of $\pm$12 cm between the T1 and T2 profiles to maximise their correlation. The estimated correlations (Fig. ) are substantially scattered around a mean correlation of $\sim \mathrm{\hspace{0.17em}0.50}$ ($\text{SD}=0.13$). The relative majority ($\sim 43\mathrm{\%}$) of all possible profile pairs ($N=152$) shows a maximum correlation at a shift of $+\mathrm{3}\mathrm{cm}$, which is well below the estimated upper vertical height difference of the trenches.

Mean trench profiles

Spatial averaging is expected to improve the correlation between the trenches compared to the single profiles. We therefore correlate the mean trench profiles of T1 and T2, allowing again for bin-wise vertical shifts of the T2 profile to maximise the correlation.

The mean trench profiles (Fig. ) are highly correlated ($r_{\text{T1,T2}}=0.81$ for an optimal shift of $+\mathrm{3}\mathrm{cm}$; $p=0.01$, accounting for the full autocorrelation structure and allowing for vertical shifting), indicating a common isotopic signal reproducible over a spatial scale of at least $500\mathrm{m}$. It is interesting to note that this value is above most of the single inter-trench correlations (Fig. ).

In both mean profiles, we observe five seasonal cycles spanning a range of $\sim \mathrm{6}\text{–}\mathrm{7}\text{‰}$ at the surface, but being attenuated further down and exhibiting no clear sinusoidal shape in the depth range of $65\text{–}90\mathrm{cm}$. Interestingly, this obscured part without clearly depleted ${\mathit{\delta }}^{18}\mathrm{O}$ “winter” values is found in both trenches, indicating that this feature persists over at least $500\mathrm{m}$ and is thus likely of climatic origin, e.g. a winter with unusually low precipitation. Despite the statistically significant correlation, pronounced differences between the mean profiles are present, such as a significantly more depleted, and partially more enriched, isotopic composition of the T2 mean between 50 and 80 cm.

Comparison of the mean ${\mathit{\delta }}^{18}\mathrm{O}$ profiles (lines: seasonal, points: annual mean) from T1 (black) and T2 (red). To maximise the seasonal correlation ($r_{\text{T1,T2}}=0.81$), trench T2 was shifted by $+\mathrm{3}\mathrm{cm}$. For the first three depth bins, the number of existing observations varies on absolute coordinates between the trench profiles. To obtain non-biased seasonal mean profiles only the depth range covered by all profiles is used. Shading represents the range of the approximate annual-mean profiles due to different binning definitions. Note that their first and last value are biased since the trench data are incomplete here. Vertical dashed grey lines mark the six local maxima of the average of both seasonal mean profiles.

[Figure omitted. See PDF]

Observed and modelled inter-profile correlation as a function of profile spacing for T1. Observations for a given spacing are the mean across all possible profile pairs. Shading denotes the standard error of the mean assuming maximum degrees of freedom (DOF) of $N=12$ (estimated from the effective DOF of the horizontal trench data accounting for autocorrelation).

[Figure omitted. See PDF]

To analyse annual-mean ${\mathit{\delta }}^{18}\mathrm{O}$ time series we use different binning methods to average the seasonal trench data with bins defined by (1) the six local maxima determined from the average of the two mean trench profiles, (2) the five local minima, (3) the midpoints of the ascending slopes flanking the maxima and (4) the midpoints of the descending slopes. To display the data on an absolute time axis we assign the year 2012 to the first annual bin. The annual-mean time series derived from the four possible binning sets are averaged to obtain a single time series for each trench (Fig. ). The correlation of the average annual-mean ${\mathit{\delta }}^{18}\mathrm{O}$ time series of $r_{\stackrel{\mathrm{‾}}{\text{T1}},\stackrel{\mathrm{‾}}{\text{T2}}}={0.87}_{-0.20}^{+0.07}$ (range represents the four binning methods) is comparable to that of the mean seasonal profiles (0.81). However, five observations of annual means are too short to reliably estimate the correlation and its significance.

Spatial correlation structure

We have shown that spatial averaging significantly increases the correlation between the trenches. To learn more about the spatial correlation structure of the trench isotope record, we investigate (1) the inter-profile correlation as a function of profile spacing for T1 and (2) the inter-trench correlation between different sets of mean profiles from T1 and the mean T2 profile.

The inter-profile correlation is estimated as the mean of the correlations obtained from all possible T1 profile pairs separated by a given spacing, allowing a tolerance in the horizontal position of $\mathrm{5}\mathrm{\%}$. For the inter-trench correlation, we define a T1 profile stack as the spatial average across a certain number of T1 profiles separated by a given distance, and determine all possible equivalent stacks. The inter-trench correlation with the mean T2 profile is then recorded as the mean across the correlations between the mean T2 profile and all possible T1 stacks.

The inter-profile correlation approaches 1 for nearest neighbours and rapidly drops with increasing inter-profile spacing before it stabilises around a value of $\sim \mathrm{\hspace{0.17em}0.5}$ for spacings $\mathit{\gtrsim }10\mathrm{m}$ (Fig. ). For the inter-trench correlation, we find a steady increase in the correlation with the T2 reference with increasing number of profiles used in the T1 stacks (Fig. ). Additionally, the correlation increases with a wider spacing between the individual profiles of a stack.

The observed decrease of the inter-profile correlation with distance suggests a horizontal autocorrelation of the isotopic composition. Indeed, a positive autocorrelation of the horizontal ${\mathit{\delta }}^{18}\mathrm{O}$ variations of T1 with a decorrelation length of $\mathit{\lambda }\simeq 1.5\mathrm{m}$ is found by applying a Gaussian kernel correlation which accounts for the irregular horizontal sampling. As we do not expect any climate-related part of the isotopic record to vary on such small spatial scales, we attribute the observed autocorrelation to noise features.

Observed and modelled correlations between T1 profile stacks and the mean of all T2 profiles depending on the number of profiles in the T1 stack for three selected inter-profile spacings. Observed results for given spacing and number of profiles are the mean across the correlations obtained for all possible unique stacks and only calculated when at least $15$ stacks are available.

[Figure omitted. See PDF]

Statistical noise model

The inter-profile correlation $r_{XY}$ provides an estimate of the signal-to-noise variance ratio $F$ of single profiles , $$F=r_{XY}/(\mathrm{1}-r_{XY}).$$ Neglecting the small-scale correlation, we estimate $F$ from the data using the mean inter-profile correlation for the profile spacings between 10 and 35 m and find $F=0.9\pm 0.1$.

Based on our findings, we develop a simple statistical model: we assume that each trench profile consists of the sum of a common climate signal $S$ and a noise component $w$ independent of the signal. The noise component is modelled as a first-order autoregressive process (AR(1)) in the horizontal direction. Then, the inter-profile correlation coefficient between profiles $X$ and $Y$ becomes a function of their spacing $d$ (see Appendix ), $$r_{XY}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+F^{-\mathrm{1}}}}\left\{\mathrm{1},+,F^{-\mathrm{1}},\exp ,\left(-,{\displaystyle \frac{d}{\mathit{\lambda }}}\right)\right\}.$$ Here, $F^{-\mathrm{1}}=\text{var}\left(w\right)/\text{var}\left(S\right)$ is the inverse of the signal-to-noise variance ratio. Using our estimate for $F$ and the value for $\mathit{\lambda }$ obtained in the previous section, the model reproduces the observed inter-profile correlations (Fig. ). Applying the same parameter values, the theoretical inter-trench correlation (Eq. ) is also in good agreement with the empirical results (Fig. ). This validates the model and the parameter values ($F,\mathit{\lambda }$) from the intra-trench ($\sim 10\mathrm{m}$) to the inter-trench spatial scale ($\sim 500\mathrm{m}$).

Discussion

Our trench data confirm earlier results that individual ${\mathit{\delta }}^{18}\mathrm{O}$ firn-core records from low-accumulation regions are strongly influenced by local noise . Going beyond this finding, our two-dimensional ${\mathit{\delta }}^{18}\mathrm{O}$ data set also allows one to determine the spatial structure and to learn about the causes of the noise. In this section, we discuss our findings in the context of the possible noise sources and derive implications for seasonal to inter-annual climate reconstructions based on firn cores.

Local stratigraphic noise and regional climate signal

A horizontally stratified trench without horizontal isotopic variations would yield perfectly correlated single profiles. Opposed to that, our records (Table ) show a significant variability in horizontal direction with mean variances $\left({\mathit{\sigma }}_{\text{h,T1}}^{\mathrm{2}},\simeq ,5.9,,(,\text{‰},{)}^{\mathrm{2}}\right.$, $\left..{\mathit{\sigma }}_{\text{h,T2}}^{\mathrm{2}},\simeq ,5.3,,(,\text{‰},{)}^{\mathrm{2}}\right)$ that are smaller but of the same order of magnitude as the mean down-core variances $\left({\mathit{\sigma }}_{\text{v,T1}}^{\mathrm{2}},\simeq ,9.5,,(,\text{‰},{)}^{\mathrm{2}}\right.$, $\left..{\mathit{\sigma }}_{\text{v,T2}}^{\mathrm{2}},\simeq ,7.3,,(,\text{‰},{)}^{\mathrm{2}}\right)$. In consequence, coherent isotopic features between single profiles separated by the trench distance are only found by chance (Fig. : the median correlation is $0.49$, only for two pairs ($\sim 1.3\mathrm{\%}$) the correlation is $>\mathrm{\hspace{0.17em}0.8}$). Thus, single firn profiles from our study region are no representative recorders of climatic isotope signals on the vertical scales analysed here.

On the horizontal scale of the trenches ($\sim 10\text{–}500\mathrm{m}$), we expect that stratigraphic noise dominates the isotopic variations . The observed length scale of the horizontal decorrelation of the noise ($\mathit{\lambda }\sim 1.5\mathrm{m}$) is similar in magnitude as that on which the local small-scale surface height variations occur, indicating that stratigraphic noise is in fact the prominent noise component in our data.

Despite the low single-profile representativity, the trench record contains a climate signal becoming apparent through the inter-profile correlation of $\sim \mathrm{\hspace{0.17em}0.5}$ remaining on scales on which the stratigraphic noise is decorrelated ($\mathit{\gtrsim }\mathrm{\hspace{0.17em}10}\mathrm{m}$). It appears to be regionally ($\mathit{\lesssim }\mathrm{\hspace{0.17em}1}\text{km}$) coherent as suggested firstly by the comparable values of the inter-profile correlation for spacings $\mathit{\gtrsim }\mathrm{\hspace{0.17em}10}\mathrm{m}$ and the mean correlation between single T1–T2 records (Fig. ), and secondly by the common seasonal signal observed in the mean trench profiles (Fig. ).

Noise is always reduced by averaging profiles; here, the autocorrelation causes nearby profiles to share more common noise variance than profiles at a larger spacing. Therefore, albeit the same number of profiles is averaged, stacks using a larger profile spacing will exhibit less common noise variance and hence a larger proportion of the underlying signal (Fig. ). Our results show a minimum profile spacing of $\sim 10\mathrm{m}$ to be optimal.

Representativity of a ${\mathit{\delta }}^{18}\mathrm{O}$ firn profile stack expressed as the correlation with a hypothetical climate signal depending on the number of profiles averaged and their inter-profile spacing. For annual resolution, the two limiting cases discussed in the text are displayed (a best-case scenario, b worst-case scenario), each for $\mathrm{2}\mathrm{m}$ (black) as well as $10\mathrm{m}$ (blue) inter-profile spacing. As a reference, in each case the seasonal representativity is shown in red for $10\mathrm{m}$ inter-profile spacing.

[Figure omitted. See PDF]

Representativity of isotope signals on seasonal to inter-annual timescales

For quantitative climate reconstructions from proxy data, a robust estimate of the climate signal is necessary. Based on our statistical noise model, we can estimate the isotopic climate signal content of a profile stack for our study region depending on the number of averaged profiles and their spacing.

To this end, we define the climate representativity of a trench profile stack as the correlation between the stack and a common climate signal (Eq. ). This signal is identified with the coherent isotope signal of the trench records. A physical interpretation of the climate representativity is then the upper bound of the correlation with a local temperature record, for example from a weather station. However, bearing in mind other influences such as meteorology (variable storm tracks, changing moisture source regions, precipitation-weighting), the true correlation will be lower. In the limit of independent noise our definition of climate representativity is equivalent to the expression derived by .

In general, climate signals are timescale-dependent. For example, the seasonal amplitude of the isotopic signal is much larger than any variations between the years. On the other hand, one expects larger changes of the climate signal on longer timescales, such as glacial–interglacial cycles. Moreover, not only the climate signal but also the noise can be a function of the timescale. One extreme example for this is the non-climate oscillations of the isotopic composition on up to centennial timescales which have been indicated by snow-pit studies around Vostok station and linked to the movement of accumulation waves on various scales . Since the climate representativity (Eq. ) depends on the ratio $F$ of signal and noise variance, it is in consequence also a function of the timescale.

Here, we assess the climate representativity of firn isotope profiles from our study region for two specific timescales: (1) the original resolution of the trench data and (2) an annual resolution based on binning the trench data.

Analysing the original data, which is dominated by variations on seasonal timescales, the climate representativity can be readily calculated with the model parameters obtained in Sect. . For the analysis at annual resolution, estimates of both annual signal and noise variance are necessary to assess the variance ratio $F$. However, the shortness of our trench data only allows heuristic estimates (see Appendix  for details). Specifically, for the annual noise variance we discuss two limiting cases: for case I we assume that the vertical noise is white (best-case scenario), for case II that the vertical noise shows complete inter-dependence on the sub-annual timescale (worst case). The inverse of the annual signal-to-noise variance ratio, $F_{\text{annual}}^{-\mathrm{1}}=\text{var}(w{)}_{\text{annual}}/\text{var}(S{)}_{\text{annual}}$, used in the model is then $\sim \mathrm{\hspace{0.17em}1.2}$ for case I and $\sim \mathrm{\hspace{0.17em}8.7}$ for case II. A summary of the noise levels is given in Table .

Variance levels of the two trenches: the horizontal variance is the mean variance of all depth layers on absolute coordinates; the down-core variance is the mean vertical variance of all respective trench profiles. The seasonal as well as the annual variance levels denote the variances of the respective mean seasonal and annual ${\mathit{\delta }}^{18}\mathrm{O}$ time series of the two trenches (Fig. ). All numbers are in units of (‰)${}^{\mathrm{2}}$.

Trench	Horizontal ${\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{2}}$	Down-core ${\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}$	Seasonal ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{v}}^{\mathrm{2}}$	Annual ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{a}}^{\mathrm{2}}$
T1	5.9	9.5	5.1	1.15
T2	5.3	7.3	3.3	0.21

For single profiles, the climate representativity estimated at seasonal resolution is $0.69$ (Fig. ). At annual resolution, single profiles show a representativity of $0.67$ in the best-case scenario (Fig. a) and of $0.32$ in the worst-case scenario (Fig. b).

Noise variance and standard deviation (SD) of the trench data together with the ratio of measurement uncertainty (${\mathit{\sigma }}_{\text{CRDS}}=0.09\text{‰}$) and respective noise SD, given for different resolutions and for the two limiting cases of the annual noise variance. The decadal noise level estimates are calculated from the annual noise variances accounting for full forward diffusion.

Resolution	Variance in (‰)${}^{\mathrm{2}}$	SD in ‰	${\mathit{\sigma }}_{\text{CRDS}}/\text{SD}$
seasonal	5.9	2.43	4 $\mathrm{\%}$
annual: case I	0.84	0.92	10 $\mathrm{\%}$
annual: case II	5.9	2.43	4 $\mathrm{\%}$
decadal: case I	0.08	0.28	32 $\mathrm{\%}$
decadal: case II	0.56	0.75	12 $\mathrm{\%}$

Similar to the correlation between the trenches (Fig. ), the representativity increases with the number of profiles averaged with a stronger increase for larger inter-profile spacings. However, spacings above $10\mathrm{m}$ do not yield a further increase as the stratigraphic noise is largely decorrelated. To obtain a climate representativity of $0.8$ at annual resolution with profiles separated by $10\mathrm{m}$, a minimum of 3–16 cores are needed (from best to worst case). Demanding a representativity of $0.9$, the number of cores required increases to 6–37.

Probability of detecting a linear temperature trend of $0.5{}^{\circ }\mathrm{C}$ (50 yr)${}^{-\mathrm{1}}$ ($p=0.05$) (solid lines) and of determining the strength of the trend with an accuracy of $25\mathrm{\%}$ (dashed lines) as a function of the number of firn cores averaged and for the two scenarios of the annual noise variance discussed in the text (black: best case, blue: worst case).

[Figure omitted. See PDF]

The modelled single-profile climate representativity at annual resolution appears consistent with previous findings from Dronning Maud Land. estimated a low signal-to-noise variance ratio of $F=0.14$ obtained from the cross-correlations of $16$ annually resolved ${\mathit{\delta }}^{18}\mathrm{O}$ records from an area of $500\text{km}\times \mathrm{\hspace{0.17em}200}\text{km}$. Due to the large inter-core spacings, the stratigraphic noise in the records is decorrelated and the variance ratio $F$ can be translated into a single-profile representativity of $r_{SX}=\mathrm{1}/\sqrt{\mathrm{1}+F^{-\mathrm{1}}}\simeq \mathrm{\hspace{0.17em}0.35}$, consistent with our findings for the worst-case scenario. However, the records analysed in are also subject to dating uncertainties, additional variability caused by spatially varying precipitation-weighting and other effects. Therefore, the similar representativities are not necessarily caused by the high stratigraphic noise level assumed in the worst-case scenario. In addition, our trench data indicate vertical autocorrelation of the noise (Fig. b and Sect. ). Thus, the true climate representativity for our study region will likely be in between our limiting estimates.

Stratigraphic noise affects not only isotopes but also other parameters measured in ice cores, such as aerosol-derived chemical constituents. investigated the seasonal to inter-annual representativity of ion records from five Greenland firn cores taken at varying distances from 7 to 10 m in the vicinity of the NEEM drilling site. Using the definition of representativity based on , they found inter-annual representativities of $\sim 0.55\text{–}0.95$, depending on the number of averaged cores and the ion species considered. These numbers are slightly higher than our best-case-scenario results for ${\mathit{\delta }}^{18}\mathrm{O}$, which is expected since the accumulation rate at the NEEM site is about 3 times higher than at Kohnen Station .

Our estimates for the climate representativity of firn cores hold as long as the signal-to-noise variance ratio $F$ does not change. Variance-affecting processes such as diffusion and densification have equal influence on signal and noise and thus do not alter the ratio $F$. On the other hand, only one component might change over time; for example, the noise variance might vary due to changing environmental conditions, or the variability of the climate could have been different in the past for certain time periods. Nevertheless, given the stability of the Holocene climate, we do not expect first-order changes of the signal and noise properties over time. However, we do expect a timescale dependency of the climate signal with more variance associated with longer timescales e.g.. The signal-to-noise variance ratio and the climate representativity of firn cores will improve considerably on these scales.

Implications

Our noise level and implied climate representativity estimates underline the challenge of firn-core-based climate reconstructions at seasonal to annual resolution in low-accumulation regions. For our study site, we now discuss implications of our noise model concerning (1) the required measurement precision of water isotopes in the case of classical isotope thermometry, (2) the potential noise fraction in isotope signals of the EDML ice core and (3) the detectability of a temperature trend.

Our estimates of the stratigraphic noise level are based on the upper 1 m of firn. Due to the shortness of the data our results are limited by insufficient knowledge of the vertical noise covariance structure for inter-annual and longer timescales for which we now assume white-noise behaviour. The noise of isotopic data obtained from deeper parts of the firn column is affected by diffusion and densification. Densification is only of importance when studying the isotopic time series in the depth domain since in that case constantly sampled data will contain noise levels on varying timescales. We estimate the effect of diffusion and find that at decadal resolution below the firn–ice transition the noise level at Kohnen Station is only $\mathrm{5}\mathit{\%}$ smaller as compared to the undiffused case (Appendix  and Table ).

The noise of an isotopic signal includes the stratigraphic noise as well as noise caused by the measurement process. Since the stratigraphic noise is a function of the number of analysed cores, and measurement precision is often related to measurement time, obtaining the best signal is a trade-off between measurement precision and the amount of analysed samples.

At seasonal as well as annual resolution, the measurement uncertainty of the trench data of ${\mathit{\sigma }}_{\text{CRDS}}=0.09\text{‰}$ is much lower ($\sim \mathrm{4}\text{–}10\mathrm{\%}$) than the standard deviation of the stratigraphic noise (Table ). This ratio is independent of the temporal resolution if a lower resolution is obtained by averaging annually resolved data as both contributions decrease by the same amount in the averaging process, assuming independence between the samples. In such a case, priority should be given to measuring and averaging across multiple cores in order to reduce the (stratigraphic) noise levels instead of performing high-precision measurements on single cores. As an example, with the cavity ring-down spectrometers used for this work faster measurements are possible by reducing the number of repeated measurements per sample and applying a memory correction . We explicitly note that this possibility is limited to classical single-isotope (${\mathit{\delta }}^{18}\mathrm{O}$) reconstructions as it can affect the data usability for diffusion- or deuterium-excess-based inferences.

If a lower temporal resolution is obtained by a coarser sampling of the cores, the measurement error to stratigraphic noise ratio will depend on the analysed resolution (Table ). For a resolution corresponding to 10 years, our measurement uncertainty might amount to up to $32\mathrm{\%}$ of the stratigraphic noise level, accounting for full diffusion. The noise level of single cores would become comparable to the measurement uncertainty for averages over $\sim \mathrm{\hspace{0.17em}104}$ or $\sim \mathrm{\hspace{0.17em}735}$ $\mathrm{years}$ (best- or worst-case scenario of annual noise level).

The deep EPICA DML ice core obtained in the vicinity of Kohnen Station reflects the climate evolution in Antarctica over the last $\mathrm{150\hspace{0.17em}000}$ $\mathrm{years}$ . studied a core section covering the last $6000$ $\mathrm{years}$ at decadal resolution. We find a ${\mathit{\delta }}^{18}\mathrm{O}$ variance for this section of $\sim \mathrm{\hspace{0.17em}0.57}(\text{‰}{)}^{\mathrm{2}}$. Using our diffusion-corrected stratigraphic noise variance estimates would imply that $\sim \mathrm{\hspace{0.17em}15}\text{–}100\mathrm{\%}$ (from best to worst case) of the observed decadal variance in the core might be noise (Table ), masking the underlying climate variability. We note that this is only a rough estimate as the shortness of the trench data does not allow one to fully assess the decadal noise covariance. In any case, averaging across multiple cores seems necessary in low-accumulation regions to reconstruct the climate variability of the last millennium. Alternatively, if only the magnitude of variability is of interest, the proxy variability has to be corrected for the noise contribution e.g..

As a final example of applying our noise model, we test the influence of stratigraphic noise on the detectability of a linear trend at Kohnen Station. This is motivated by the finding of that in the last 50 $\mathrm{years}$ the surface temperature over East Antarctica has warmed by about half a degree. While both the climate signal as well as the relationship between local temperature and isotopic signal are complex, we assess the detectability with a simplified model experiment. For this, we assume the climate signal to be a purely linear trend ($0.5{}^{\circ }\mathrm{C}$ (50 yr)${}^{-\mathrm{1}}$) and a linear isotope-to-temperature relationship ($\mathrm{1}{\text{‰ K}}^{-\mathrm{1}}$), further influenced only by post-depositional noise. In a Monte Carlo approach repeated ${10}^{\mathrm{5}}$ times, we create stacks from $50\mathrm{yr}$ long ${\mathit{\delta }}^{18}\mathrm{O}$ profiles with post-depositional noise variances based on our two limiting cases (Table ), accounting for an average effect of diffusion on the annual noise level over the 50 $\mathrm{years}$ (Appendix ) and assuming independent noise between the profiles (inter-profile spacings $\mathit{\gtrsim }\mathrm{\hspace{0.17em}10}\mathrm{m}$), and vary the number of averaged profiles. A trend in the stacked profile is successfully detected for an estimated trend that is significantly larger than zero ($p=0.05$); the estimated slope is defined to be correct when it lies in a range of $25\mathrm{\%}$ around the true slope. The probability of trend detection/slope determination is then the ratio of successful reconstructions to total number of realisations.

Drilling a single core, the probability of detecting the trend or reconstructing its slope is around $25\mathrm{\%}$ in the best-case and below $10\mathrm{\%}$ in the worst-case scenario (Fig. ). To reliably ($>80\mathrm{\%}$ of the cases) detect the warming over the East Antarctic Plateau, our results suggest that averaging across at least $\sim \mathrm{5}\text{–}35$ firn cores taken at spacings of $10\mathrm{m}$ (Fig. ) is needed, depending on the scenario for the annual noise variance. Inferring the right slope would need 3 times that number of cores. We note that more realistic assumptions about the isotopic signal (natural climate and atmospheric variability, varying isotope–temperature relationship, etc.) further complicate the trend detectability.

Conclusions

We presented extensive oxygen stable water isotope data derived from two snow trenches excavated at Kohnen Station in the interior of Dronning Maud Land, Antarctica. The two-dimensional approach allowed a thorough investigation of the representativity of single firn-core isotope profiles, as well as of the spatial structure of the signal and noise over spatial scales of up to $500\mathrm{m}$ and a time span of approximately 5 years.

The trench data confirm previous studies that single low-accumulation ($\le 100{\text{mm w.e. yr}}^{-\mathrm{1}}$) isotope profiles only show a weak coherent signal at least on sub-decadal timescales. We also demonstrate that the spatial average of a sufficient number of profiles provides representative isotopic signals, consistent with our finding that the local noise has a small horizontal decorrelation length ($\sim 1.5\mathrm{m}$). This also suggests stratigraphic noise to be the major contribution to the horizontal isotopic variability. A statistical noise model based on a first-order autoregressive process successfully explains the observed covariance structure and allows one to reproduce the correlation statistics between the trenches.

Based on these results we infer appropriate sampling strategies. At our low-accumulation ($64{\text{mm w.e. yr}}^{-\mathrm{1}}$) site an optimal spacing of about $10\mathrm{m}$ is necessary for a sufficient decorrelation of the stratigraphic noise. We estimate the climate representativity of isotope profiles depending on the number of averaged firn cores and the inter-core spacing. Our estimates show that at seasonal resolution five cores at the optimal spacing are necessary to obtain representative ($r>0.9$) isotope signals; at annual resolution up to $\sim \mathrm{8}$ times as many cores are needed. As climate variations are typically stronger on longer timescales than analysed here, the climate representativity of firn- and ice-core reconstructions for slower climate changes will likely be higher.

We present two examples of how the stratigraphic noise might hamper the quantitative interpretation of isotope in terms of climate variations at our study site. Our data suggest that at least $15\mathit{\%}$ of the decadal variations seen in the EPICA DML ice core over the last 6000 $\mathrm{years}$ might be post-depositional noise, but the climate signal might also be masked by a much higher decadal noise level. A simplified model experiment shows that the faithful reconstruction of the recent positive temperature trend observed over the East Antarctic Plateau likely requires averaging across at least 5–35 firn cores. For single-proxy (${\mathit{\delta }}^{18}\mathrm{O}$) reconstructions this task could be rendered more easily by the fact that the annual noise level is substantially larger than typical measurement uncertainties. Thus, monitoring the measurement error depending on sample throughput could allow fast measurements for the benefit of analysing many cores. Alternatively, using indirect methods based on diffusion or gas isotope ratios might circumvent the problem of stratigraphic noise.

Since the stratigraphic noise is related to irregular re-deposition and erosion of snow and the formation of surface dunes, it primarily depends on the local accumulation rate, besides further factors such as wind strength, temperature, seasonal timing of the precipitation and snow properties. Therefore, we expect that our representativity results improve (worsen) for regions with higher (lower) accumulation rates. In effect, our results are likely applicable for large parts of the East Antarctic Plateau, but similar studies in West Antarctica and Greenland – regions with considerably higher accumulation rates – are needed. In addition, studies with deeper trenches that cover a longer time period, complemented by spectral analyses of nearby firn cores, are necessary to enhance our knowledge of the vertical noise covariance structure. This is crucial to determine the climate representativity on longer timescales. Deeper trenches would also allow one to link our representativity results to actual correlations with temperature time series derived from weather stations. The latter is part of ongoing work at Kohnen Station.

Data availability

The trench oxygen isotope data presented in this study are archived at the PANGAEA database (http://www.pangaea.de) under 10.1594/PANGAEA.861675. PANGAEA is hosted by the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI), Bremerhaven, and the Center for Marine Environmental Sciences (MARUM), Bremen, Germany.

Derivation of noise model

Definitions

We consider isotope profiles $X_i\left(z\right)$ at equidistant spacings $\mathrm{\Delta }\mathrm{\ell }$ where $z$ is depth on absolute coordinates and $i$ refers to the profile's horizontal position along a snow trench, ${\mathrm{\ell }}_i={\mathrm{\ell }}_{\mathrm{0}}+i\times \mathrm{\Delta }\mathrm{\ell }$, with some arbitrary starting position ${\mathrm{\ell }}_{\mathrm{0}}$ (Fig. ). This and all subsequent nomenclature is summarised in Table .

We assume each $X_i\left(z\right)$ as a sum of a common signal $S\left(z\right)$ and a noise term $w_i\left(z\right)$ independent of $S$, $$X_i\left(z\right)=S\left(z\right)+w_i\left(z\right).$$

The noise $w_i\left(z\right)$ is modelled as an AR(1) process in the horizontal direction, $$w_i\left(z\right)=a{w}_{i-\mathrm{1}}\left(z\right)+{\mathit{\epsilon }}_i\left(z\right),$$ where $a$ is the autocorrelation parameter with $\mathrm{0}\le a\le \mathrm{1}$ and ${\mathit{\epsilon }}_i\left(z\right)$ are independent random normal variables (white noise). We assume the same variance $\text{var}\left(w\right)$ of the noise in both the horizontal and the vertical direction.

The mean of a set of $N$ trench isotope profiles ${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}$ (profile stack) is defined by the indices $\mathit{\left\{}i\mathit{\right\}}=\mathit{\{}{i}_{\mathrm{1}},i_{\mathrm{1}}+i_{\mathrm{2}},i_{\mathrm{1}}+i_{\mathrm{2}}+i_{\mathrm{3}},\mathrm{\dots },i_{\mathrm{1}}+i_{\mathrm{2}}+\mathrm{\cdots }+i_N\mathit{\}}$. This nomenclature of incremental steps simplifies the expressions obtained later. ${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\left(z\right)$ is given by the signal $S\left(z\right)$ and the mean of the noise terms, $${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\left(z\right)=S\left(z\right)+{\displaystyle \frac{\mathrm{1}}{N}}\left(w_{i_{\mathrm{1}}},+,w_{i_{\mathrm{1}}+i_{\mathrm{2}}},+,\mathrm{\cdots },+,w_{i_{\mathrm{1}}+i_{\mathrm{2}}+\mathrm{\cdots }+i_N}\right)\left(z\right).$$

The Pearson correlation of two single profiles $X_i$ and $X_{i+j}$ is $$\begin{array}{cc}\text{corr}\left(X_i,,,X_{i+j}\right)& ={\displaystyle \frac{\text{cov}\left(X_i,,,X_{i+j}\right)}{\sqrt{\text{var}\left(X_i\right)\text{var}\left(X_{i+j}\right)}}}\\ & & ={\displaystyle \frac{\text{var}\left(S\right)+\text{cov}\left(w_i,,,w_{i+j}\right)}{\text{var}\left(S\right)+\text{var}\left(w\right)}},\end{array}$$ using the independence of signal and noise and the stationarity of $w$.

Summary of the nomenclature used for the statistical noise model.

Symbol	Description
$z$	absolute depth below mean snow height
$X_i$	trench isotope profile at position ${\mathrm{\ell }}_i$
$\mathrm{\Delta }\mathrm{\ell }$	spacing of adjacent profiles
$S$	climate signal contained in $X_i$
$w_i$	noise contained in $X_i$
${\mathit{\epsilon }}_i$	white-noise component of $w_i$
${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}$	profile stack
$a$	autocorrelation parameter; $a=\exp \left(-,\mathrm{\Delta },\mathrm{\ell },/,\mathit{\lambda }\right)$
$\mathit{\lambda }$	horizontal noise decorrelation length
$d$	inter-profile distance
$N$	number of profiles
${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}$	relative effective noise variance of stack ${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}$
$F$	signal-to-noise variance ratio

Sketch of a snow trench used for the derivation of the statistical noise model. Vertical isotope profiles $X_i$ are spaced at constant intervals of $\mathrm{\Delta }\mathrm{\ell }$ at locations $l_i=l_{\mathrm{0}}+i\mathrm{\Delta }\mathrm{\ell }$. The horizontal distance $d$ between two profiles $X_{i_{\mathrm{1}}}$ and $X_{i_{\mathrm{1}}+i_{\mathrm{2}}}$ is defined by the incremental index $i_{\mathrm{2}}$, $d=i_{\mathrm{2}}\mathrm{\Delta }\mathrm{\ell }$.

[Figure omitted. See PDF]

(a) Relative effective noise variance ${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}$ of a profile stack ${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}$ as a function of the number of profiles averaged for a profile spacing of $\mathrm{\Delta }\mathrm{\ell }=\mathrm{1}\mathrm{m}$ and for different values of the autocorrelation parameter $a$. The limiting case of white noise ($a=\mathrm{0}$) is indicated by a dashed line. (b) ${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}$ as a function of the autocorrelation parameter $a$ for different numbers of averaged profiles and profile spacings.

[Figure omitted. See PDF]

The correlation of a profile stack ${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}$ and the signal is given by $$\begin{array}{cc}\text{corr}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},,,S\right)& ={\displaystyle \frac{\text{cov}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},,,S\right)}{\sqrt{\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)\text{var}\left(S\right)}}}\\ & & ={\displaystyle \frac{\text{var}\left(S\right)}{\sqrt{\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)\text{var}\left(S\right)}}}.\end{array}$$

Similarly, the correlation of two profile stacks with indices $\mathit{\left\{}i\mathit{\right\}}$ and $\mathit{\left\{}j\mathit{\right\}}$, assuming independent noise between the sets, is obtained from $$\begin{array}{cc}\text{corr}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},,,{\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}j\mathit{\right\}}}\right)& ={\displaystyle \frac{\text{cov}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},,,{\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}j\mathit{\right\}}}\right)}{\sqrt{\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}j\mathit{\right\}}}\right)}}}\\ & & ={\displaystyle \frac{\text{var}\left(S\right)}{\sqrt{\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}j\mathit{\right\}}}\right)}}}.\end{array}$$

Derivation of model correlations

To derive the explicit correlations ()–() for the AR(1) noise model, we need expressions for the noise variance, $\text{var}\left(w\right)$, the noise covariance in horizontal direction, $\text{cov}(w_i,w_{i+j})$, and the variance of a profile stack, $\text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)$.

The former two are given by $$\begin{array}{}& & \text{var}\left(w\right)={\displaystyle \frac{\text{var}\left(\mathit{\epsilon }\right)}{\mathrm{1}-a^{\mathrm{2}}}},& & \text{cov}(w_i,w_{i+j})={\displaystyle \frac{\text{var}\left(\mathit{\epsilon }\right)}{\mathrm{1}-a^{\mathrm{2}}}}{a}^j=\text{var}\left(w\right)a^j.\end{array}$$ The index $j$ can be expressed here by the distance $d={\mathrm{\ell }}_{i+j}-{\mathrm{\ell }}_i$ between the profiles $X_i$ and $X_{i+j}$ and the spacing of adjacent profiles $\mathrm{\Delta }\mathrm{\ell }$ as $j=d/\mathrm{\Delta }\mathrm{\ell }$. Further, for an AR(1) process the lag one autocorrelation is given by $a=\exp \left(-,\mathrm{\Delta },\mathrm{\ell },/,\mathit{\lambda }\right)$ with the decorrelation scale $\mathit{\lambda }$. It follows from Eq. () that the horizontal noise covariance decreases exponentially with distance $d$ as $$\text{cov}(w_i,w_{i+j})=\text{var}\left(w\right)\exp \left(-,{\displaystyle \frac{d}{\mathit{\lambda }}}\right).$$

The variance of a profile stack ${\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}$ is calculated according to $$\begin{array}{cc}& \text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)=⟨{\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}^{\mathrm{2}},\left(,z,\right)⟩-{⟨{\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},\left(,z,\right)⟩}^{\mathrm{2}}=\text{var}\left(S\right)\\ & & +{\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}⟨{\left(w_{i_{\mathrm{1}}},+,w_{i_{\mathrm{1}}+i_{\mathrm{2}}},+,\mathrm{\cdots },+,w_{i_{\mathrm{1}}+i_{\mathrm{2}}+\mathrm{\cdots }+i_N}\right)}^{\mathrm{2}},\left(,z,\right)⟩,\end{array}$$ where $⟨\cdot ⟩$ denotes the expected value. Using the multinomial identity $({\mathit{\xi }}_{\mathrm{1}}+{\mathit{\xi }}_{\mathrm{2}}+\mathrm{\cdots }+{\mathit{\xi }}_N{)}^{\mathrm{2}}={\sum }_{i=\mathrm{1}}^N{\mathit{\xi }}_i^{\mathrm{2}}+\mathrm{2}{\sum }_{i=\mathrm{1}}^{N-\mathrm{1}}{\sum }_{j>i}{\mathit{\xi }}_i{\mathit{\xi }}_j$ yields $$\begin{array}{cc}& \text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)=\text{var}\left(S\right)+{\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}\mathit{\left\{}N\text{var}\right(w)+\\ & \mathrm{2}(⟨w_{i_{\mathrm{1}}},w_{i_{\mathrm{1}}+i_{\mathrm{2}}}⟩+⟨w_{i_{\mathrm{1}}},w_{i_{\mathrm{1}}+i_{\mathrm{2}}+i_{\mathrm{3}}}⟩+\mathrm{\cdots }+.\\ & .⟨w_{i_{\mathrm{1}}},w_{i_{\mathrm{1}}+i_{\mathrm{2}}+i_{\mathrm{3}}+\mathrm{\cdots }+i_N}⟩+⟨w_{i_{\mathrm{2}}},w_{i_{\mathrm{2}}+i_{\mathrm{3}}}⟩+\mathrm{\cdots }+.\\ & & .⟨w_{i_{\mathrm{2}}},w_{i_{\mathrm{2}}+i_{\mathrm{3}}+\mathrm{\cdots }+i_N}⟩+\mathrm{\cdots }+⟨w_{i_{N-\mathrm{1}}},w_{i_{N-\mathrm{1}}+i_N}⟩\left)\mathit{\right\}}.\end{array}$$ By applying Eq. () for the horizontal covariance of the noise we obtain $$\begin{array}{cc}& \text{var}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}}\right)=\text{var}\left(S\right)+\text{var}\left(w\right)\times \\ & & {\underbrace{{\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}}}\mathit{\{}N+\mathrm{2}\left(a^{i_{\mathrm{2}}},+,a^{i_{\mathrm{2}}+i_{\mathrm{3}}},+,\mathrm{\cdots },+,a^{i_{\mathrm{2}}+\mathrm{\cdots }+i_N},+,a^{i_{\mathrm{3}}},+,\mathrm{\cdots },+,a^{i_{\mathrm{3}}+\mathrm{\cdots }+i_N},+,\mathrm{\cdots },+,a^{i_N}\right)\mathit{\}}}}_{{\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}},\end{array}$$ where we define ${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}$ as the relative effective noise variance of the profile stack. In the limiting case of $a=\mathrm{0}$ (zero autocorrelation) ${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}=\mathrm{1}/N$, in the limit of $a=\mathrm{1}$ (perfect autocorrelation) ${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}=\mathrm{1}$. In general, ${\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}$ is a function of both $N$ and the spacing of the profiles averaged (Fig. ).

For final expressions of the correlation functions ()–(), we define the signal-to-noise variance ratio $F:=\frac{\text{var}\left(S\right)}{\text{var}\left(w\right)}$ and use Eqs. () and () to obtain $$\begin{array}{cc}& \text{inter-profile corr.:}\\ & & \text{corr}\left(X_i,,,X_{i+j}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+F^{-\mathrm{1}}}}\left\{\mathrm{1},+,F^{-\mathrm{1}},\exp ,\left(-,{\displaystyle \frac{d}{\mathit{\lambda }}}\right)\right\},& \text{stack-signal corr.:}\\ & & \text{corr}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},,,S\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+F^{-\mathrm{1}}{\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}}}},& \text{stack-stack corr.:}\\ & & \mathrm{\text{corr}}\left({\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}i\mathit{\right\}}},,,{\stackrel{\mathrm{‾}}{X}}_{\mathit{\left\{}j\mathit{\right\}}}\right)={\displaystyle \frac{\mathrm{1}}{{\left\{\left(\mathrm{1},+,F^{-\mathrm{1}},{\mathit{\sigma }}_{\mathit{\left\{}i\mathit{\right\}}}^{*\mathrm{2}}\right),\left(\mathrm{1},+,F^{-\mathrm{1}},{\mathit{\sigma }}_{\mathit{\left\{}j\mathit{\right\}}}^{*\mathrm{2}}\right)\right\}}^{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

Estimation of parameters

To evaluate the correlation functions ()–() we need estimates of the decorrelation length $\mathit{\lambda }$ and of the timescale-dependent variance ratio $F^{-\mathrm{1}}$.

For the trench data at seasonal resolution, we obtain a variance ratio of $F^{-\mathrm{1}}\simeq 1.1\pm 0.1$ from the observed inter-profile correlations of T1 (Fig. ) for profile spacings $\ge 10\mathrm{m}$, and an estimate of the decorrelation length of $\mathit{\lambda }\simeq 1.5\mathrm{m}$ from the horizontal autocorrelation of the T1 ${\mathit{\delta }}^{18}\mathrm{O}$ data. We validate the parameters by comparing the predicted (Eq. ) and observed correlations between profile stacks derived from T1 and T2 (Fig. ). This assumes independent noise between T1 and T2, a valid approximation given that the trench distance ($\sim 500\mathrm{m}$) is much larger than $\mathit{\lambda }$. Relying on the assumption of equal noise variance in the horizontal and vertical direction, a second estimate of $F^{-\mathrm{1}}\sim 1.6$ can be obtained from the observed mean T1 down-core variance (identified with signal and noise) subtracted by the observed mean T1 horizontal variance (i.e. noise, see Table ).

Going from the original seasonal resolution of the trench data to an explicit annual resolution, the short data sets only allow limited estimations. We thus make use of the following simple heuristic arguments. The annual signal variance is estimated from the mean annual ${\mathit{\delta }}^{18}\mathrm{O}$ time series of each trench neglecting the residual noise contributions and averaging both variance estimates to obtain $\text{var}(S{)}_{\text{annual}}\simeq 0.68(\text{‰}{)}^{\mathrm{2}}$ (Table ). The annual noise variance, $\text{var}(w{)}_{\text{annual}}$, is calculated from the seasonal noise variance estimated by the mean horizontal T1 variance of $\text{var}\left(w\right)\simeq 5.9(\text{‰}{)}^{\mathrm{2}}$. Physically, we expect a vertical autocorrelation of the noise due to the underlying processes (stratigraphic noise, ; diffusion, ), which is also indicated by our data (Fig. 1b). However, due to the limited vertical trench data, the vertical noise autocorrelation cannot be reliably estimated and we discuss two limiting cases: in case I the vertical noise is independent (white noise) and the seasonal noise variance therefore reduced by the number of samples included in the annual average ($N\approx \mathrm{7}$); in case II the vertical noise shows complete inter-dependence on the sub-annual timescale and its variance is not reduced by taking annual means. The resulting annual variance ratios of noise over signal are $$\begin{array}{cc}{F}_{\text{annual}}^{-\mathrm{1}}& ={\displaystyle \frac{\text{var}(w{)}_{\text{annual}}}{\text{var}(S{)}_{\text{annual}}}}\\ & & \simeq {\displaystyle \frac{\mathrm{1}}{0.68}}\times \left\{\begin{array}{l}0.84,\\ 5.9\end{array}\right.=\left\{\begin{array}{ll}1.2,& \text{for case I},\\ 8.7,& \text{for case II}.\end{array}\right.\end{array}$$ For all longer timescales, we generally assume white-noise behaviour for the noise covariance.

Reduction of noise level by diffusion

The integral over the power spectrum $P\left(f\right)$ of a time series $X\left(t\right)$, where $f$ denotes frequency and $t$ time, is equal to the total variance of $X$ , $$\text{var}\left(X\right)=\mathrm{2}\underset{\mathrm{0}}{\overset{f_N}{\int }}P\left(f\right)\text{d}f.$$

Here, $f_N$ is the Nyquist frequency according to the sample resolution of $X$.

For white noise, the power spectrum is a constant, $P_{\mathrm{0}}\left(f\right)=P_{\mathrm{0}}=\text{const.}$ Evaluation of Eq. () then gives $\mathrm{2}{f}_N{P}_{\mathrm{0}}$.

In case of diffusion, the initial power spectrum $P_{\mathrm{0}}\left(f\right)$ is changed for a given diffusion length $\mathit{\sigma }$ and local annual layer thickness $\dot{b}$ according to $$P\left(f\right)=P_{\mathrm{0}}\left(f\right)\exp \left\{-,{\left(\mathrm{2},\mathit{\pi },\mathit{\sigma },{\dot{b}}^{-\mathrm{1}},,f\right)}^{\mathrm{2}}\right\},$$ where $\dot{b}$ is introduced in order to work in the temporal domain, thus with frequency measured in yr${}^{-\mathrm{1}}$. For white noise, the integral in Eq. () is straightforward to solve, yielding the noise variance at a given resolution $f_N$ accounting for diffusion: $$\begin{array}{cc}& \mathrm{2}{P}_{\mathrm{0}}\underset{\mathrm{0}}{\overset{f_N}{\int }}\exp \left\{-,{\left(\mathrm{2},\mathit{\pi },\mathit{\sigma },{\dot{b}}^{-\mathrm{1}},,f\right)}^{\mathrm{2}}\right\}\text{d}f=\\ & & P_{\mathrm{0}}\sqrt{\mathit{\pi }}/\left(\mathrm{2}\mathit{\pi }\mathit{\sigma }{\dot{b}}^{-\mathrm{1}}\right)\text{erf}\left(\mathrm{2}\mathit{\pi }\mathit{\sigma }{\dot{b}}^{-\mathrm{1}}{f}_N\right).\end{array}$$

To estimate how much diffusion has reduced the annual trench noise level at the firn–ice transition at decadal resolution ($f_N=0.05{\text{yr}}^{-\mathrm{1}}$), we assume $\dot{b}=\mathrm{7}{\text{cm yr}}^{-\mathrm{1}}$ and a constant diffusion length of $\mathit{\sigma }=\mathrm{8}\mathrm{cm}$ . Evaluation of Eq. () yields a reduction of the annual noise power of $\sim 0.095\left[{\text{yr}}^{-\mathrm{1}}\right]P_{\mathrm{0}}$. Comparing the diffused with the undiffused case shows that at decadal resolution, diffusion only has a relative effect of $\sim \mathrm{5}\mathit{\%}$ on the reduction of the annual noise level.

Thus, diffusion only has a minor influence on decadal and longer timescales at our study site. However, on shorter timescales it has to be accounted for. The annual noise levels given in Table  are therefore only valid for the uppermost part of the firn column where diffusion is almost negligible. In the deeper parts of the firn, they have been affected by diffusion. In our simplified trend detection experiment, we assume over the first $\mathrm{9}\mathrm{m}$ of firn (approximately the last 50 years) an average annual layer thickness of $15\mathrm{cm}$ and a mean diffusion length of $\mathit{\sigma }\sim \mathrm{5}\mathrm{cm}$ . Evaluation of Eq. () then gives an average reduction of the annual noise level to $\sim 73\mathit{\%}$ of the original value.

Thomas Laepple, Sepp Kipfstuhl and Johannes Freitag designed and led the field work. Thomas Münch and Thomas Laepple designed the analysis. Thomas Münch led the isotope measurements, performed the research, developed the noise model and wrote the manuscript. Hanno Meyer supervised the isotope measurements in Potsdam. All authors contributed significantly to the discussion of the results and the revision of the manuscript.

Acknowledgements

We thank all the scientists, technicians and the logistic support who worked at Kohnen Station in the 2012/2013 austral summer, especially Melanie Behrens, Tobias Binder, Andreas Frenzel, Katja Instenberg, Katharina Klein, Martin Schneebeli, Holger Schubert, Jan Tell and Stefanie Weissbach, for assistance in creating the trench data set. We further thank the technicians of the isotope laboratories in Bremerhaven and Potsdam, especially York Schlomann and Christoph Manthey. All plots and numerical calculations were carried out using the software R: A Language and Environment for Statistical Computing. This work was supported by the Initiative and Networking Fund of the Helmholtz Association Grant VG-NH900. We thank the editor and two anonymous reviewers for their constructive comments that helped to improve the manuscript. Edited by: E. Wolff Reviewed by: two anonymous referees