Constraints on post-depositional isotope modifications in East Antarctic firn from analysing temporal changes of isotope profiles

Abstract

The isotopic composition of water in ice sheets is extensively used to infer past climate changes. In low-accumulation regions their interpretation is, however, challenged by poorly constrained effects that may influence the initial isotope signal during and after deposition of the snow. This is reflected in snow-pit isotope data from Kohnen Station, Antarctica, which exhibit a seasonal cycle but also strong interannual variations that contradict local temperature observations. These inconsistencies persist even after averaging many profiles and are thus not explained by local stratigraphic noise. Previous studies have suggested that post-depositional processes may significantly influence the isotopic composition of East Antarctic firn. Here, we investigate the importance of post-depositional processes within the open-porous firn ($\mathit{\gtrsim }$ 10 cm depth) at Kohnen Station by separating spatial from temporal variability. To this end, we analyse 22 isotope profiles obtained from two snow trenches and examine the temporal isotope modifications by comparing the new data with published trench data extracted 2 years earlier. The initial isotope profiles undergo changes over time due to downward advection, firn diffusion and densification in magnitudes consistent with independent estimates. Beyond that, we find further modifications of the original isotope record to be unlikely or small in magnitude ($\ll$ 1 ‰ RMSD). These results show that the discrepancy between local temperatures and isotopes most likely originates from spatially coherent processes prior to or during deposition, such as precipitation intermittency or systematic isotope modifications acting on drifting or loose surface snow.

Introduction

The isotopic composition of water measured in firn and ice cores is an important climate proxy. The abundance ratios of the stable water isotopologues in falling snow are shaped by different fractionation processes in between the moisture source and the precipitation site, including evaporation , air-mass advection and Rayleigh distillation and snow formation . Hence, isotope ratios can be linked to the climatic conditions at the local or moisture source site. For instance, physical modelling of the large-scale hydrological cycle and the fractionation processes has validated the link between the isotopic composition of precipitation and local temperature and references therein previously inferred for polar ice sheets, where observational evidence has suggested a robust relationship at large spatial scales (i.e. continental) between the isotopic composition of snow and annual-mean temperature at the sampling sites . Isotope data archived in polar ice cores have therefore become an invaluable means to infer past site temperature variations e.g. or changes in the moisture sources e.g. and show, at least qualitatively, a globally consistent picture of glacial–interglacial to millennial-scale climate changes . However, it is questioned whether the assumption holds in general that pre-depositional fractionation processes alone are the main influence on the isotopic composition of firn and ice, while seemingly fulfilled for large spatial and temporal scales. Particularly in low-accumulation areas for which the snow surface is exposed to the atmosphere for a substantial time, a variety of processes are thought to considerably modify the original atmospheric isotope signal during or after deposition of the snow, thus from seasonal to interannual timescales e.g..

For the Antarctic Plateau, modifications of the original isotope signal that is imprinted in precipitation are generally expected. In buried snow and firn, the isotopic composition is affected by diffusion of interstitial water vapour and by densification ; however, these processes do not lead to any net change in the isotopic composition. In contrast, the seasonal intermittency of precipitation and accumulation can bias the original signal, induce variability or lead to a combination of both . In combination with the low accumulation rates on the Antarctic Plateau, precipitation intermittency also increases the time the surface is exposed to the atmosphere . These conditions might favour fractionation, diffusive and advective processes that can considerably alter the snow's original isotopic composition, acting either post condensation (on falling or drifting snow) or post deposition on snow at the surface or within the open-porous firn column which is no longer subject to erosion but still in contact with the atmosphere. For instance, the exchange of water vapour between the first metre of firn and the overlying atmosphere through diffusion and wind-driven ventilation can introduce vapour with a different isotopic signature to the firn and significantly change the isotopic composition. Isotopic exchanges between the top layer of snow and the lower atmosphere have been observed on daily scales at the NEEM site in Greenland and on diurnal scales at Kohnen Station in East Antarctica . Isotopic fractionation associated with sublimation, condensation and recrystallisation processes within the near-surface firn might change the initial isotope signal, as indicated by observations and lab experiments . Since these post-depositional processes also depend, besides temperature, on other climatic variables such as wind speed and relative humidity, any seasonal or interannual variations in these variables would induce additional variability in the isotope record. However, for East Antarctica, a quantitative assessment of the individual processes based on firn-core data is still outstanding, and their importance for shaping the isotope signal in the near-surface firn remains poorly constrained.

An additional, important source of variability in low-accumulation firn-core records is the spatial variability from stratigraphic noise caused by uneven deposition and the constant wind-driven erosion, redistribution and vertical mixing of the snow surface. A previous study from Kohnen Station in Dronning Maud Land, East Antarctica, has shown that the spatial variability can be overcome by averaging across a suitable number of single profiles extracted from snow trenches . This yielded a spatially representative isotope signal on a horizontal scale of approximately 500 m. However, contrasting the isotope data with instrumental observations from a nearby automatic weather station AWS, suggests that this regional signal does not necessarily represent a regional temperature signal (Fig. ). Whereas the isotope record shows strong year-to-year variability, the observed temperature variations are characterised by a regular seasonal cycle and small interannual changes. This discrepancy stresses the importance of contributions other than regional temperature alone to the formation of the isotope signal, such as precipitation intermittency and changes during or after deposition. Since quantitative knowledge on seasonal and interannual variations of accumulation amounts is still sparse on the Antarctic Plateau , in this study we investigate whether post-depositional isotope modifications in the open-porous firn contribute to the observed discrepancy between the isotope data and local temperatures at Kohnen Station.

One way to address the question of post-depositional modification is to compare two firn-core isotope profiles obtained at different times and to measure the nature in which the first profile has been modified. However, due to stratigraphic noise, the comparison of two single records sampled at different times will always confound temporal isotope changes and spatial variability. Therefore, in this study we present and use data from a new extensive snow trench campaign yielding 22 profiles of isotopic composition from two trenches and compare these with the data of the previous trench campaign conducted 2 years earlier. By generating representative records from the spatial averaging of single profiles, together with the theoretical understanding of stratigraphic noise, our study allows us for the first time to quantitatively follow the isotopic changes over a time span of 2 years. We use independent knowledge on firn diffusion and densification to subtract these effects from the observed temporal modifications. Therefore, beyond simply stating the problem of local isotope–temperature discrepancy, we go further and can quantitatively estimate the influence of post-depositional change for our study site. This is an important step towards better constraining the isotope signal formation in East Antarctic firn.

Comparison of oxygen isotope data and 2 m air temperature at Kohnen Station, Antarctica. (a) Mean ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ profiles of trenches T13–1 (black) and T13–2 (red) (modified after ) on original 3 cm (lines) as well as annual resolution (points with uncertainty shading from shifting the range of the annual bins). (b) 2 m air temperature (black lines: monthly means; blue points: annual means) recorded by the automatic weather station AWS9 located at Kohnen Station $<\mathrm{1}$ km from the trench excavation sites (Fig. ). Note the different timescales (a is based on counting and binning the extrema of the isotope data, b shows true calendar dates). The mean of the 2 m air temperature shown here lies about 3.5 ${}^{\circ }$C above the published local 10 m firn temperature (Table ).

[Figure omitted. See PDF]

Data and methods

Sampling and measurements

A pair of firn trenches, each with a horizontal length of 50 m and a depth of 3.4 m, was excavated using a snow blower in the austral summer field season 2014/2015 near Kohnen Station , the location of the EPICA Dronning Maud Land deep ice core drilling site (Fig.  and Table ). This campaign extends the published oxygen isotope data set obtained from two shallower ($\sim \mathrm{1}$ m) trenches in 2012/2013 . From the new trenches, we present the top 1.75 m of the data which are expected to cover the period imprinted in the trenches of the first campaign. To avoid contamination, the new trench positions were shifted relative to the previous ones by 160 and 300 m, respectively, and are separated by 550 m (Fig. ). In the remaining part of the manuscript, “T13” will refer to the pair of previous trenches from 2012/2013, “T15” to the pair of new trenches from 2014/2015.

Information on the EPICA Dronning Maud Land (EDML) drilling site at Kohnen Station, Antarctica. Listed are approximate position (latitude, longitude), elevation, 10 m firn temperature ${\stackrel{\mathrm{‾}}{T}}_{\text{firn}}$, mean annual accumulation rate of snow $\stackrel{\mathrm{‾}}{b}$ and mean daily wind speed ${\stackrel{\mathrm{‾}}{v}}_{\text{wind}}$ ($\pm \mathrm{1}$ standard deviation).

Drilling site	Latitude	Longitude	Elevation	${\stackrel{\mathrm{‾}}{T}}_{\text{firn}}$	$\stackrel{\mathrm{‾}}{b}$	${\stackrel{\mathrm{‾}}{v}}_{\text{wind}}$
	${}^{\circ }$ N	${}^{\circ }$ E	m a.s.l.	${}^{\circ }$C	mm w.e. yr${}^{-\mathrm{1}}$	m s${}^{-\mathrm{1}}$
EDML	$-$75.0${}^{\text{a}}$	${\mathrm{0.1}}^{\text{a}}$	${\mathrm{2892}}^{\text{a}}$	$-$44.5${}^{\text{a}}$	${\mathrm{64}}^{\text{a}}$/${\mathrm{82.5}}^{\text{b}}$	$\mathrm{4.4}\pm {\mathrm{2.3}}^{\text{c}}$

${}^{\text{a}}$ . ${}^{\text{b}}$ Mean of snow stake measurements 2013–2015. ${}^{\text{c}}$ AWS9 data 1998–2013 .

Map of the study area at Kohnen Station. Snow trenches are shown as black and red lines, firn-core sites as green-filled circles. The drilling site of the EPICA Dronning Maud Land (EDML) ice core is shown as a black star, the position of the automatic weather station (AWS9) as a blue-filled diamond. The main wind direction 57${}^{\circ }$ from geographic north, is indicated with a black arrow. The trenches were aligned perpendicularly to the local snow-dune direction.

[Figure omitted. See PDF]

Fieldwork for the new T15 trench campaign was conducted as follows: horizontal profiles of the surface height variations were obtained along each trench using a levelling instrument. The uncertainty of these profiles is estimated from the reading accuracy of the levelling rod of 0.5 cm. The windward walls of the trenches were cleaned after excavation by slicing off a thin firn layer. Firn profiles were then sampled directly off the wall with a vertical resolution of 3 cm and a horizontal spacing of 5 m, yielding 11 profiles in each trench. The vertical resolution is small enough to evaluate the seasonal cycle of the isotope data of $\sim \mathrm{20}$ cm ; the interprofile distance of 5 m corresponds approximately to three times the decorrelation length of the stratigraphic noise observed in the T13 record . At both trenches, excavation and subsequent sampling of the profiles was conducted in two consecutive stages (2 times $\sim \mathrm{1}$ m depth); each stage was completed within 24 h. All firn samples ($N=\mathrm{1214}$) were stored in plastic bags, tightly packed, transported to Germany in frozen state and analysed for oxygen (${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$) and hydrogen (${}^{\mathrm{2}}\mathrm{H}{/}^{\mathrm{1}}\mathrm{H}$) isotope ratios at the isotope laboratory of the Alfred Wegener Institute (AWI) in Potsdam, using a cavity ring-down spectrometer (L2130i, Picarro Inc.). The results are reported in the usual delta notation (oxygen isotopes: ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$, hydrogen isotopes: $\mathit{\delta }\mathrm{D}$) in per mil (‰) relative to the international VSMOW/SLAP

VSMOW (Vienna Standard Mean Ocean Water) and SLAP (Standard Light Antarctic Precipitation) are reference waters with VSMOW being defined as 0 ‰ both for oxygen and hydrogen isotopes, and, by convention, SLAP being fixed to $-$55.5 and $-$428 ‰ for ${\mathit{\delta }}^{\mathrm{18}}$O and $\mathit{\delta }$D, respectively.

scale. Calibration and correction of the raw measurements was performed as described in . The mean combined measurement uncertainty is 0.08 ‰ for ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (root mean square deviation, RMSD) and 0.8 ‰ for $\mathit{\delta }\mathrm{D}$, assessed by evaluating a standard not used in the calibration and correction procedure. All data are archived under 10.1594/PANGAEA.876639 .

Trench depth scale

Following , we record and display the trench isotope data with respect to an absolute height reference given by the respective maximum of the surface height profile of each trench. Note that the surface layer of the trench records is incomplete on this depth scale (up to $\sim$ 10 cm for T13 and $\sim \mathrm{18}$ cm for T15) due to the surface undulations. Averaging of trench profiles is performed relative to the absolute height reference. Therefore, the number of data points contributing to a mean profile is lower and varies in the surface layer. This part is marked by dashed lines for all mean profiles and is excluded from all quantitative calculations. Our conclusions are therefore limited to firn depths below $\sim \mathrm{10}$ cm but are, however, as will be shown, not essentially changed when including the surface layer.

Spatial variability of average trench profiles

Spatial variability arising from stratigraphic noise is a major contribution to the overall variability of individual trench isotope profiles . Its magnitude $\mathit{\omega }$ can be estimated from the horizontal variability of the trench isotope record. Averaging across individual trench profiles reduces the total noise variability. Specifically, stratigraphic noise can be modelled by a first-order autoregressive process with a horizontal decorrelation length for the study region of $\mathit{\lambda }\simeq \mathrm{1.5}$ m . Then, the residual noise variability of a mean profile built by averaging across $N$ individual records is $${\mathit{\epsilon }}_{\text{res}}={\displaystyle \frac{\mathit{\omega }}{N^{\mathrm{2}}}}(N+f(N,d,\mathit{\lambda }\left)\right)\equiv {\displaystyle \frac{\mathit{\omega }}{N_{\text{eff}}}},$$ where $f(N,d,\mathit{\lambda })$ is a function of $N$, $\mathit{\lambda }$ and of the interprofile distances $d$. Equation () can equivalently be expressed through the effective number of records, $N_{\text{eff}}$. For independent noise (zero autocorrelation, $\mathit{\lambda }\to \mathrm{0}$), $f(N,d,\mathit{\lambda }\to \mathrm{0})\to \mathrm{0}$ and thus $N_{\text{eff}}\to N$.

Sampling and statistical properties of the trench ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ records from the field seasons 2012/2013 T13, and 2014/2015 (T15, this study). Listed are number and distance of sampled profiles, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ values and variance, correlation of mean trench profiles and estimated signal-to-noise variance ratios (SNR) after . Correlations are maximised through allowing relative vertical shifts (optimal shift given in brackets). 67 % confidence intervals (CI) for the variance estimates account for autocorrelation of the data. Average signal-to-noise ratios are given with an uncertainty of 1 standard error (SE).

Trench record	T13–1	T13–2	T15–1	T15–2
Number of profiles:	38	4	11	11
Profile distances (m):	$\sim$ 0.1–2.5	10, 20	5	5
${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (‰):
range: min/max	$-$54/$-$34	$-$50/$-$38	$-$56/$-$32	$-$55/$-$33
mean (SD)	$-$44.4 (3.1)	$-$44.0 (2.7)	$-$44.7 (3.8)	$-$44.5 (3.8)
${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ variance ((‰)${}^{\mathrm{2}}$):
mean horizontal (67 % CI)	5.9 (5.2–7.0)	5.3 (4.2–7.0)	7.0 (6.1–8.3)	6.6 (5.7–7.7)
mean vertical (67 % CI)	9.5 (8.3–11.1)	7.3 (5.9–9.6)	13.8 (12.0–16.3)	14.2 (12.3–16.8)
Mean profile correlation (optimal shift)	0.81 ($+\mathrm{3}$ cm)		0.91 ($-$0.5 cm)
SNR ($\pm$1 SE)	$\mathrm{0.9}\pm \mathrm{0.4}$	$\mathrm{0.5}\pm \mathrm{0.5}$	$\mathrm{1.0}\pm \mathrm{0.3}$	$\mathrm{1.5}\pm \mathrm{0.5}$

Quantification of downward advection, firn densification and firn diffusion

We expect that within 2 years the original T13 isotope profiles have been compressed through densification of the firn, advected downwards due to new snow fall and affected by firn diffusion. To quantify these effects, certain site-specific parameters have to be known.

Firn densities are provided independently of the trench data by high-resolution X-ray Computer Tomography data of the firn cores B41 and B42 (core distance $\sim \mathrm{10}$ m, ) drilled in vicinity to the trenches ($\sim \mathrm{1}$ km, Fig. ). The average firn density in the first metre is $\sim \mathrm{330}$ kg m${}^{-\mathrm{3}}$. The densification rate relative to the surface is $\sim \mathrm{2}$ % m${}^{-\mathrm{1}}$ when regressing density against depth over the first 2 m, $\sim \mathrm{7}$ % m${}^{-\mathrm{1}}$ when regressing over the first 5 m.

The local annual accumulation rate of snow was 28.8 $\pm$ 0.4 cm ($\pm$1 standard error) in the year 2013 and 20.8 $\pm$ 0.3 cm in 2014, which was estimated from an array of snow stake measurements conducted near the trench excavation sites. In general, the recent local accumulation rate strongly depends on the measurement site, with values ranging from 20–30 cm of snow per year which is up to 50 % larger than the published longtime mean (Table ).

In case of isothermal firn, diffusion of interstitial water vapour leads to overall smoothing of an isotope profile which can be described as the convolution with a Gaussian kernel . The amount of smoothing (the width of the Gaussian convolution kernel) is controlled by the diffusion length $\mathit{\sigma }$ which increases monotonically in the upper firn layer . We model $\mathit{\sigma }$ according to with diffusivity after . Firn density is the main driver for the depth dependency of the diffusion length. For the calculations we smooth the stacked B41/B42 density data by fitting a quadratic polynomial in the square root of the depth. For the concept of differential diffusion, we consider a firn layer which is located at the average depth $z_{\mathrm{1}}$ and has thickness $\mathrm{\Delta }z$ over which the increase in diffusion length ($\mathrm{\Delta }\mathit{\sigma }$) is small compared to the layer thickness, $\mathrm{\Delta }\mathit{\sigma }/\mathrm{\Delta }z\ll \mathrm{1}$. Now the firn layer is advected downwards to the depth $z_{\mathrm{2}}$. The total amount of diffusion that acted since the layer has been at the surface is the convolution of the layer's initial isotope profile at the surface, ${\mathit{\delta }}_{\mathrm{0}}$, with a diffusion length $\mathit{\sigma }\left(z_{\mathrm{2}}\right)$. Equivalently, since the Gaussian convolution is a linear operation, we can express this as the diffusion of ${\mathit{\delta }}_{\mathrm{0}}$ with $\mathit{\sigma }\left(z_{\mathrm{1}}\right)$, followed by diffusion of the resulting profile with the differential diffusion length:

The new T15 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ data set. Displayed are the isotope records of trench T15–1 (a) and trench T15–2 (b) as two-dimensional colour images, and the mean profiles from averaging across the individual profiles of each trench (c), displayed for the optimal vertical shift of the T15–2 mean profile (see text). The trench surface height profiles are given by solid black lines, the near-surface part of each mean profile is marked by dashed lines since the trench data are incomplete there (see Sect. ). The vertical scale in (a, b) is strongly exaggerated.

[Figure omitted. See PDF]

$$\stackrel{\mathrm{̃}}{\mathit{\sigma }}=\sqrt{{\mathit{\sigma }}^{\mathrm{2}}\left(z_{\mathrm{2}}\right)-{\mathit{\sigma }}^{\mathrm{2}}\left(z_{\mathrm{1}}\right)}.$$ For the T13 isotope profiles, we account for an approximate average effect of differential diffusion over 2 years by considering the average diffusion lengths calculated over the depth of the T13 profiles before and after downward advection, neglecting the seasonal variations in firn temperature.

Statistical tests

We use the Kolmogorov–Smirnov (KS) test to assess whether distributions of differences between mean trench profiles vary. Autocorrelation of the data is accounted for with a modified version of the standard test adopting effective degrees of freedom of $n(\mathrm{1}-a)$ . Here, $n$ is the total number of data points for each profile and $a$ the estimated autocorrelation parameter at lag 1. The KS test compares the empirical cumulative distribution functions of the data and is thus sensitive to differences in both mean and variance.

Results

New T15 isotope data and qualitative comparison with T13

The two new T15 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ trench records measured in 2015/2016 (Fig. a, b) are qualitatively consistent with the T13 data measured 2 years earlier. The isotopic variability within the first metres of firn is characterised by roughly horizontal, alternating layers of enriched and depleted isotopic composition (Fig. a, b) which are separated on average by the annual layer thickness of snow (20–30 cm) and thus likely indicative of the climatic seasonal cycle. In addition, stratigraphic noise leads to significant horizontal variability, visible through discontinuous and inhomogeneous layering as well as patchy features, for example at the surface of trench T15–2 (Fig. b).

We find similar statistical properties for the data of each trench campaign (Table ). The higher variances in the vertical direction of the T15 records are partly expected for autocorrelated data in combination with a larger record length, in addition to the contribution by the strongly enriched layer around 170 cm depth.

Averaging across all individual profiles of each T15 trench reduces the noise level and yields mean profiles that are highly correlated (correlation $r=\mathrm{0.91}$, RMSD $\sim \mathrm{1.2}$ ‰, Fig. c) and thus spatially representative. We maximised this match by allowing vertical shifts of the T15–2 mean profile. Using linearly interpolated data on a resolution of 0.5 cm, we find an optimal shift of $-$0.5 cm. We note the exceptionally high delta values at the top of the T5–2 mean profile which stem from a prominent dune at the trench surface (Fig. ). However, on the absolute depth scale this top part has no overlap with the T15–1 mean profile and therefore does not contribute to the total T15 mean profile discussed below. Despite their representativity, the T15–1 and T15–2 mean profiles show strong year-to-year variability confirming the discrepancy to local temperature previously found for T13 (Fig. ). This also becomes apparent through the increase in average T15 summer maxima (Fig. c) which is statistically significant ($p<\mathrm{0.01}$) but not captured by the evolution of local summer temperatures (Fig. ).

Our first findings show that at our study site both the nature of the regional isotope signal and the stratigraphic noise are comparable between the two trench campaigns. In the following sections we quantitatively assess to what extent the original T13 signal can be recovered with the T15 trenches obtained 2 years later. For this task, we use a single data set for T13 and T15 from averaging across each pair of mean profiles (Fig. ), accounting for the optimal vertical shifts that maximise each intertrench correlation (Table ).

Expected isotope profile changes between 2013 and 2015

We analyse to what extent the T13 record can be recovered from the new T15 data, and which changes have modified the original record. Within the 2 years, we expect that the T13 isotope profiles are advected downwards, compressed by densification and smoothed by firn diffusion. Testing for additional isotope modifications hence requires first estimating the magnitudes of those expected processes. We do this in two ways: firstly, we use data that is independent of the trench records. Secondly, to check consistency with the first estimate, we determine the optimal parameter set that minimises the difference between the T13 and T15 mean profiles.

The mean oxygen isotope profiles of the T15 (this study) and T13 () trenches on their original depth scale. The incomplete surface layer of the trenches is marked by dashed lines.

[Figure omitted. See PDF]

Using the available independent snow stake and density data, we obtain the following estimates: the annual accumulation rates suggest a downward advection of the T13 profiles after 2 years of $\sim \mathrm{50}$ cm. Further, we expect additional diffusional smoothing of the T13 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ profiles according to a differential diffusion length (Eq. ) of $\stackrel{\mathrm{̃}}{\mathit{\sigma }}\sim \mathrm{1.9}$ cm. The estimated densification rate at the study site of $\sim \mathrm{2}$–7 % m${}^{-\mathrm{1}}$ implies a compression of the T13 profiles after 2 years of approximately 1–4 cm.

For the second estimate, we vary the three parameters (downward advection $\mathrm{\Delta }$, differential diffusion length $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, compression $\mathit{\gamma }$) in order to minimise the root mean square deviation between the T15 and T13 mean profiles. To avoid an influence on our results, we choose the range of tested parameter values independently of the trench data: for the downward advection, we apply vertical shifts between 40 and 60 cm, comprising the snow-stake-based range of the recent annual accumulation rates. We vary the differential diffusion length from $\mathrm{0}$ to 8 cm, which is equivalent to additional diffusional smoothing of the original T13 mean profile from zero to the maximum possible amount at the firn–ice transition. Finally, compression is applied for values between $\mathrm{0}$ and 10 cm (equivalent to $\mathrm{0}$ to $\sim \mathrm{5}$ times the observed average densification rate). We obtain the best agreement ($\text{RMSD}=\mathrm{0.92}$ ‰, Fig. ; $r=\mathrm{0.93}$) between the T15 and the modified T13 mean profile ($=$ T13${}^{\ast }$) for the optimal parameters ${\mathrm{\Delta }}_{\text{opt}}=\mathrm{50.5}$ cm, ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\text{opt}}=\mathrm{2.3}$ cm and ${\mathit{\gamma }}_{\text{opt}}=\mathrm{3.5}$ cm (Fig. ). These trench-based parameter estimates agree reasonably well with the independent estimates from above, showing that the trench data are compatible with our assumptions and parameterisations for downward advection, densification and diffusion. Indeed, using the independent parameter estimates (${\mathrm{\Delta }}_{\text{ind}}=\mathrm{50}$ cm, ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\text{ind}}=\mathrm{1.9}$ cm, ${\mathit{\gamma }}_{\text{ind}}=\mathrm{2.2}$ cm from mean over estimated densification rate) to modify the original T13 mean profile ($=$ T13${}^{\ast \ast }$) results in a deviation from T15 ($\text{RMSD}=\mathrm{0.94}$ ‰, $r=\mathrm{0.93}$) that is only slightly higher compared to T13${}^{\ast }$.

Effect of downward advection, firn diffusion and linear compression due to densification on the misfit (root mean square deviation, RMSD) between the T15 and the modified T13 mean profile. We record the RMSD for each point in the three-dimensional parameter space of downward advection, compression and diffusion. For each diffusion–compression pair, the figure shows the local minimum in RMSD (contour lines) from varying across the range of advection values, hence the RMSD for the optimal downward advection (colour scale). The global minimum in RMSD is marked with a black dot. Varying the downward advection has in fact the largest influence on the RMSD.

[Figure omitted. See PDF]

We note that the largest portion of optimising the fit between T15 and T13${}^{\ast }$ is accounted for by the downward advection. This is obvious from only shifting the T13 mean profile vertically to find the maximum correlation with T15, without accounting for diffusion and densification. We find an optimal shift of 48.5 cm ($r=\mathrm{0.88}$) with a minimum misfit of $\text{RMSD}=\mathrm{1.07}$ ‰ (Fig. ). Thus, the gain in RMSD is only small when adding diffusion and densification according to T13${}^{\ast }$ (black dot in Fig. ) or T13${}^{\ast \ast }$, but still appears significant given the above-found consistency in the magnitude of the trench-based and independent estimates. This is further supported by the fact that no second minimum in RMSD exists outside the region bounded by the contour line of only downward advection ($\text{RMSD}=\mathrm{1.07}$ ‰, Fig. ) where the magnitudes of diffusion and densification are unrealistically high.

Expected changes of the T13 and comparison to the T15 mean profile. (a) Upper panel: the original T13 mean profile (blue) and its modification by diffusion (black: 2-year diffusion with differential diffusion length $\stackrel{\mathrm{̃}}{\mathit{\sigma }}=\mathrm{2.3}$ cm) as well as densification (red: linear compression of $\mathit{\gamma }=\mathrm{3.5}$ cm). Lower panel: the original T13 mean profile (blue) compared to the joint effect of 2-year diffusion and linear compression (green, T13${}^{\ast }$). (b) The T15 mean profile (black) in comparison to the T13 mean profile after modifying the latter according to (1) the optimal parameters for downward advection, incremental diffusion and linear densification (green, T13${}^{\ast }$) and to (2) the corresponding parameters obtained independently from the trench records (orange, T13${}^{\ast \ast }$). Additionally, the difference between T15 and T13${}^{\ast \ast }$ is shown (violet lines, axis to the right). For comparison, the grey dotted line marks the difference between T15 and T13 only shifted optimally ($\mathrm{\Delta }=\mathrm{48.5}$ cm). Vertical dashed lines indicate the isotopic summer maxima which are not in phase with the difference curve.

[Figure omitted. See PDF]

The visual agreement of the trench mean profiles after modifying T13 according to downward advection, diffusion and densification is remarkable regarding cyclicity and, to a lesser extent, the amplitude of the isotope variations (Fig. b). However, deviations occur throughout most of the record's overlap (Fig. b) and are even amplified where the amplitude of the T13 profile prior to diffusion was smaller than for T15 (depths of $\sim \mathrm{70}$ cm and around $\sim \mathrm{125}$–140 cm). Here, locally additional diffusion does not lead to an improved match, although overall it reduces the mismatch between the profiles (Fig. ). In general, the profile deviations are relatively large compared to the influence of firn diffusion and densification on the original T13 profile (Fig. a), which calls for studying further processes in order to explain them.

Do the remaining differences represent temporal or spatial variability?

We have shown that downward advection, firn diffusion and densification contribute to the temporal modification of the original T13 profiles as expected from independent data and theoretical considerations. Taking these processes into account leads to a good match of the trench mean profiles (Fig. b). Most of the match is achieved by accounting for the downward advection; adding the effects of diffusion and densification yields a slightly greater improvement (gain in RMSD of $\sim \mathrm{0.15}$ ‰). However, deviations between the profiles on the order of $\sim \mathrm{1}$ ‰ RMSD still remain. These can have two causes: firstly, additional temporal changes driven by unaccounted post-depositional processes such as firn ventilation or sublimation; secondly, remaining spatial variability, since we average a large but finite number of records which do not originate from the exact same position. We can thus deduce the importance of additional post-depositional change for our study site if we quantify the contribution of spatial variability. In the following, this is done in two ways: (1) by using the statistical model for stratigraphic noise and (2) by analysing the distributions of the profile differences.

Variability of the trench data sets. The histograms depict (a) the distribution of the spatial differences between the two mean profiles of the T13 (T13–1 vs. T13–2, blue) and the T15 trenches (T15–1 vs. T15–2, black), and (b) the combined distribution from (a) (grey) compared to the distribution of the temporal differences between the T15 and the T13${}^{\ast \ast }$ mean profiles (red). All distributions' mean values are not significantly different from zero (all $p\ge \mathrm{0.4}$, accounting for autocorrelation).

[Figure omitted. See PDF]

According to the statistical noise model, the effective number of profiles that contribute to the T13 and T15 mean profiles (Eq. ) is $N_{\text{eff}}=\mathrm{13}$ for T13 and $N_{\text{eff}}=\mathrm{20}$ for T15. The residual noise of the mean profiles arising from spatial variability is thus the noise level before averaging ($\mathit{\omega }\sim \mathrm{5}$–7 (67 % CI: 4–8) (‰)${}^{\mathrm{2}}$, Table ) divided by $\mathrm{13}$ and $\mathrm{20}$, respectively. We assume that the residual noise terms are independent of each other. Therefore, the uncertainty of the difference between the T13 and T15 mean profiles due to stratigraphic noise is the sum of each residual spatial variability or $\sim \mathrm{0.6}$–0.9 (0.5–1.0) (‰)${}^{\mathrm{2}}$. For comparison, the square of the RMSD between the T13${}^{\ast }$ (T13${}^{\ast \ast }$) and T15 mean profile (the temporal variability) is 0.85 (0.88) (‰)${}^{\mathrm{2}}$. This agreement between RMSD and estimated residual spatial variability indicates that the remaining profile differences between the modified T13 mean profile and T15 (Fig. b) are likely consistent with stratigraphic noise. We note, however, that the squared RMSD lies towards the upper end of the estimated range of residual stratigraphic noise. This also applies to the RMSD between the T15–1 and T15–2 mean profiles (square of RMSD of 1.44 (‰)${}^{\mathrm{2}}$ vs. uncertainty from residual stratigraphic noise of $\sim \mathrm{1.0}$–1.4 (0.8–1.6) (‰)${}^{\mathrm{2}}$). This could indicate that part of the spatial variability on the scale of the intertrench distances ($\sim \mathrm{500}$ m) is not explained by our stratigraphic noise model.

We therefore make a formal statistical test comparing spatial and temporal variability, which accounts for the full extent of spatial uncertainty and makes no assumption about the covariance of the noise, by analysing the deviations between the mean trench profiles. We find that the spatial differences between the mean profiles of each trench campaign (T13–1 vs. T13–2 and T15–1 vs. T15–2, Fig. a) are statistically indistinguishable ($p>\mathrm{0.5}$ from modified KS test, combining all possible directions of calculating the differences), meaning that the statistical distributions of the intertrench-scale spatial variability are not significantly different between the years 2013 and 2015. This suggests that we can combine these spatial differences as a joint measure of spatial variability and compare them to the temporal differences between the T15 and the modified T13 mean profile (Fig. b). By applying the modified KS test once more, again the null hypothesis that both differences follow the same distribution cannot be rejected (all $p>\mathrm{0.5}$ for using T13${}^{\ast \ast }$ to avoid overfitting). Thus, the temporal differences between T13${}^{\ast \ast }$ and T15 likely just arise from the fact that the trenches have different locations and can be therefore explained by spatial variability alone.

In summary, both methods show no evidence for any temporal changes of the trench record over the course of 2 years apart from downward advection accompanied by firn diffusion and densification. The remaining deviations that are observed between the mean profiles of the two trench campaigns can be entirely explained by residual spatial variability arising from stratigraphic noise.

Discussion

We presented and analysed a new extensive data set of 22 oxygen isotope profiles obtained at Kohnen Station from two 50 m long and $\sim \mathrm{180}$ cm deep snow trenches. The new trench campaign was designed such that it allows for a direct comparison with a trench data set obtained from the same site 2 years earlier in order to test for post-depositional effects. In the following, we first discuss our results concerning the expected processes that have influenced the trench isotope profiles over the observed time period, then we discuss our findings regarding the possibility of additional post-depositional changes.

Densification, diffusion and stratigraphic noise

We found a strong resemblance between the mean oxygen isotope profiles from the trench field campaigns of 2013 and 2015 (Fig. b), achieved mostly by accounting for the downward advection and further improved by adding the effects of water vapour diffusion within the firn and firn densification that occurred over the 2 years (Fig. ). The estimated magnitudes of these processes obtained from matching both records are consistent with independent estimates from snow stakes, diffusion theory and independent density profiles.

The estimated small compression of the T13 profiles is reasonable given the low densification rate observed in the top metres of nearby firn cores. However, our assumption of a linear profile compression with depth is certainly a rough approximation given the actually observed seasonal firn density variation , which might indicate a stronger density change with depth of summer compared to winter layers. However, in general the seasonality of densification in Antarctic firn is largely unclear and references therein.

Our data-based estimate of differential firn diffusion agrees with theoretical expectations and in total leads to a further reduction in RMSD between the T13 and T15 mean profiles compared to the case of downward advection and densification alone (Fig. ). More specifically, the diffusion correction improves the match of the trench mean profiles in the medium depth range but also results in higher deviations of the profile minima at the top and bottom parts of the overlap (Fig. ), where the amplitude of the T13 profile had been smaller than for T15 already prior to diffusion. Part of this mismatch might be reduced by accounting for the seasonally varying firn temperature resulting in stronger (weaker) attenuation of summer (winter) layers caused by the seasonal difference in diffusion length which is largest close to the surface . In general, firn diffusion is still an active area of research , and progress in this field could conceivably result in an improved understanding of our data.

Stratigraphic noise is a major contribution to the overall variability of isotope profiles . Our large trench data set allows a significant reduction of the noise level by averaging across the single profiles. This is done in two steps: first we average across the local (intra-trench) scale, then we average the resulting mean profiles to account for potential uncertainties on the 500 m (intertrench) scale. Furthermore, we can estimate the remaining uncertainty of the trench mean profiles based on our theoretical understanding of stratigraphic noise. As a result, we found that the difference of the T13 and T15 mean profiles still exhibits an uncertainty of $\sim \mathrm{0.77}$–0.95 ‰ (SD). Thus, the trench data allow us to detect any additional post-depositional changes of the T13 profiles that exceed a detection limit of $\sim \mathrm{1}$ ‰ RMSD. Obviously, a lower detection limit would be beneficial but is in practice constrained by the amount of field work, given the high local stratigraphic noise level as observed from the mean horizontal isotope variability (Table ).

Additional post-depositional modifications

Based on the above results we have shown that the remaining differences between the 2013 and 2015 data sets are, after accounting for downward advection, firn diffusion and densification, likely consistent with spatial variability from stratigraphic noise. In other words, we conclude that at our study site the impact of any additional post-depositional changes of the isotopic composition of the firn below $\sim \mathrm{10}$ cm is on average below the residual stratigraphic noise level, thus $\ll \mathrm{1}$ ‰ RMSD. We limited our conclusion to this depth range due to the applied absolute depth scale resulting in a lower and varying number of available data points in the surface layer. However, looking at this part of the modified T13 mean profile (dashed lines of T13${}^{\ast }$ or T13${}^{\ast \ast }$, Fig. b) also does not show any solid evidence of additional post-depositional changes.

Our conclusion is also supported by comparing the nature of the differences between the mean profiles (Fig. b) with the expected effect of post-depositional modification processes. Studied processes all point to isotopic enrichment, such as sublimation (; ) and wind-driven firn ventilation . Specifically, the latter modelling study showed that firn ventilation can result in isotopic annual-mean enrichment from the strong enrichment of isotopic winter layers, compensating an observed slight depletion of summer layers. For South Pole conditions (annual-mean temperature $-$50 ${}^{\circ }$C, accumulation rate 84 mm w.e. yr${}^{-\mathrm{1}}$, mean surface wind speed 5 m s${}^{-\mathrm{1}}$), the effect amounts to $\sim \mathrm{3}$ ‰ for firn ventilation until the layers are advected below the influence of the atmosphere, thus after several years . The environmental conditions at the South Pole are comparable to Kohnen Station (Table ), suggesting a similar influence of ventilation on the isotopic composition of the firn. The higher temperatures at Kohnen Station would even imply a slightly stronger enrichment . However, if we analyse the difference curve of the T15 and T13${}^{\ast \ast }$ mean profiles (Fig. b) we find no evidence for firn ventilation. Comparing the direct seasonal counterparts, the first winter layer, which was closest to the surface at the time of excavation of T13 and thus presumably being under strongest influence of the atmosphere, is more depleted in isotopic composition in T15 than in T13${}^{\ast \ast }$, in contrast with the expectation from firn ventilation. Moreover, despite the fact that the first three summer layers exhibit more depleted values, which would be in line with ventilation, the remaining summer layers do not confirm this finding, and none of the average annual differences show enrichment: for the first annual cycle, T15 exhibits an average difference from T13${}^{\ast \ast }$ of $-$1.6 ‰ ($-$1.3 ‰ including the surface region), for the other annual cycles the averages are $-$0.4, $\pm$0 and $-$0.1 ‰. Also in general, the difference curve (Fig. b) does not show any clear seasonal timing which might be expected for a systematic post-depositional modification. Instead, minimum and maximum differences appear rather randomly across the seasons. In addition, the global average difference of about $-$0.45 ‰ is not significantly different from zero ($p=\mathrm{0.4}$, accounting for autocorrelation). We nevertheless note that the RMSD of the first overlapping annual cycle is above our stated detection limit for post-depositional change. However, this limit applies to the average over the record's entire overlap and does not account for the possibility of autocorrelated differences. Finally, we note the seeming increase with depth of the annual-mean differences towards more positive values (Fig. b), which is also indicated by the slight skewness of the corresponding histogram (Fig. b). However, the trend is not strongly significant ($p=\mathrm{0.12}$, accounting for autocorrelation), and the KS test of the distribution of the differences showed that mean and variability of the residual temporal differences are likely explained by the spatial distribution alone. In addition, we obtain similar results (not shown) when we apply our analysis to the trench d-excess ($d:=\mathit{\delta }\mathrm{D}-\mathrm{8}\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$) data, a second-order parameter that is potentially more sensitive to post-depositional fractionation processes . The spatial and residual temporal differences between the corresponding d-excess mean profiles follow the same distribution ($p>\mathrm{0.5}$), and the histogram of the temporal differences is even more symmetric than for ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$.

In summary, all evidence suggests that post-depositional modifications from firn ventilation or sublimation are unlikely to contribute to the deviations between the T15 and the modified T13 mean profiles and that the shape of the difference curve only arises from the statistical nature of stratigraphic noise, smoothed by diffusion. We nevertheless note the possibility that additional post-depositional changes are still present at Kohnen Station but not detectable in our analysis. Wind-driven firn ventilation might exist with an effect that is, however, much weaker than expected and thus masked by the stratigraphic noise level. One possible explanation for the discrepancy between the firn ventilation model results and our data could be that the model misrepresents the isotopic signature of the surface vapour advected into the firn. Another possibility is weaker firn temperature gradients at Kohnen Station compared to the South Pole, preventing significant vapour deposition. Assessing these possibilities in detail is, however, beyond the scope of our study. Interestingly, if the seeming trend in difference values was significant, it would suggest an oriented post-deposition process that is yet unknown. In any case, the stronger profile differences for the first overlapping annual cycle might indicate modification processes that are constrained to the very surface layer. In addition, the RMSD between T15 and T13${}^{\ast }$ can be further minimised if one allows shifts in the mean value of the T13 profile (new minimum RMSD of $-$0.82 ‰ for a shift in mean of $-$0.4 ‰) which is an interesting observation, yet without any obvious explanation. However, based on the presented evidence, these possibilities are speculative and further field studies are needed to test them.

Our study underlined the pronounced discrepancy at Kohnen Station between interannual variations of isotope ratios in the firn and local temperatures and showed that this feature is not only spatially (over distances of $\sim \mathrm{500}$ m) but also temporally representative over a period of 2 years. Furthermore, given the sum of our above findings, it is unlikely that post-depositional modifications of the isotopic composition of the open-porous firn (below depths of $\sim \mathrm{10}$ cm and probably also not in shallower depths) are the cause of the observed discrepancy. Since a strong relationship between isotopes in precipitation samples and local temperature has been observed at different sites of the Antarctic Plateau , this cause must instead be sought in processes working directly at or above the firn surface. At least two explanations for this seem possible. (1) Seasonal variation and intermittency of precipitation cause the discrepancy between isotope and local temperature data . At Kohnen Station, a large part of the annual accumulation is assumed to occur in winter since little or no precipitation is observed in the summer field seasons. However, the exact seasonal and interannual variation of accumulation is still unclear due to the lack of sufficiently precise, year-round observations . The available surface height changes derived from sonic altimeters of automatic weather stations are difficult to separate into events of drifting snow and true snowfall . (2) Isotope modification occurring directly at or above the surface is the key driver for shaping the interannual isotope variations. Such processes might be acting on falling, loose or drifting snow, or on the top layer (first few centimetres) of deposited snow . The fact that our trench records are reproducible on spatial scales of at least 500 m implies that the atmospheric parameters and conditions controlling potential processes would also need to be spatially coherent.

Conclusions

Many studies, including our present one, show that interannual isotope records from the dry Antarctic Plateau are inconsistent with local temperature variations. However, beyond simply stating the problem, we take two steps further: (1) we use the average over $\mathrm{2}\times \mathrm{11}$ isotope profiles to obtain a spatially representative record. (2) We designed our study such that it allows for testing for post-depositional effects over a time span of 2 years.

Our results provide important constraints on the formation of the stable water isotope signal and its propagation with depth in East Antarctic firn: the trench records show a pure downward advection of the isotope signal within the open-porous firn ($\mathit{\gtrsim }$ 10 cm depth), further influenced only by firn diffusion and densification, with no evidence for substantial additional post-depositional modification. Hence, once the signal is archived at this stage, we do not expect any significant change in the mean values deeper down, reinforcing the credibility of palaeoclimate studies using ice core isotope data. However, from our analysis we can constrain post-depositional changes only down to the level of stratigraphic noise. Therefore, qualitatively, firn ventilation and sublimation might still be present with effects that are, however, very small or constrained merely to the surface layer for which the lower number of data points in our study prevents quantitative analyses. These constraints lead us to conclude that the observed discrepant isotope–temperature relationship on the interannual timescale must be caused either by processes prior to or during deposition.

To improve our understanding of the interannual isotope signal, we suggest a mixture of field and modelling efforts. Year-long isotope studies (e.g. in seasonal intervals) with a focus on the near surface would help to constrain isotope modifications at the interface of surface snow and atmosphere. Further, the role of precipitation and accumulation intermittency has to be clarified, e.g. through measuring wet-deposition tracers and improved accumulation measurements. These studies should optimally be accompanied by monitoring and modelling the atmospheric water vapour isotopic composition as well as modelling the potential exchange and fractionation processes between the loose or deposited snow at the surface and the overlying atmosphere.

Our results again underline the role of stratigraphic noise for the total variability of isotope records. Spatial averaging is thus essential for improving the signal-to-noise ratio and thereby separating spatial from temporal variability. Alternatively, single records can only be compared faithfully for temporal changes when their spatial separation is well below the spatial decorrelation length of the stratigraphic noise, which minimises the amount of spatial variability between the records. The effects of potential isotope modifications depend substantially on the time the surface is exposed to the atmosphere, thus on accumulation rate and seasonal timing of precipitation. Comparable recovering efforts at other ice-coring sites are hence highly needed. Our data indicate that present models might overestimate the expected influence of wind-driven firn ventilation; however, regions with higher wind speeds and lower accumulation rates might still be susceptible towards post-depositional changes within the open-porous firn. A deeper understanding of the isotope signal formation in Antarctic firn is, beyond holding intrinsic interest, essential for deciphering the temperature signal archived in ice core records and thus crucial for their palaeoclimatic interpretation.

The trench stable water isotopologue data presented in this study are archived at the PANGAEA database (https://www.pangaea.de) under 10.1594/PANGAEA.876639 (Münch et al., 2017). 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.

Software to reproduce all plots and the main analyses presented in this paper is available as an R package under https://bitbucket.org/ecus/trenchr.

TM, TL, SK and JF designed the trench campaign, SK and TM led the field work. TM and TL designed the analysis. TM performed the research and wrote the manuscript. HM supervised the isotope measurements. All authors contributed significantly to the discussion of the results and the refinement of the manuscript.

The authors declare that they have no conflict of interest.

Acknowledgements

We would like to thank all scientists, technicians and the logistic support who worked at Kohnen Station in the 2014/2015 austral summer, especially Holger Schubert and Tobias Binder, for assistance in creating the trench data set. Thanks are also due to Christoph Manthey of the isotope laboratory in Potsdam for assistance in measuring the isotope data. We are grateful to Gerit Birnbaum and Carleen Tijm-Reijmer for providing the AWS9 data. This work significantly benefited from discussions with Mathieu Casado and from proofreading and many valuable comments by Andrew Dolman. All plots and numerical calculations were carried out using the software R: A Language and Environment for Statistical Computing. We acknowledge the support of this work by the Initiative and Networking Fund of the Helmholtz Association Grant VG-NH900. We thank Joël Savarino for the kind handling of the manuscript and two anonymous reviewers for their constructive comments that helped to improve a first version.The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association. Edited by: Joel Savarino Reviewed by: two anonymous referees