Estimation and calibration of the water isotope differential diffusion length in ice core records

Abstract

Palaeoclimatic information can be retrieved from the diffusion of the stable water isotope signal during firnification of snow. The diffusion length, a measure for the amount of diffusion a layer has experienced, depends on the firn temperature and the accumulation rate. We show that the estimation of the diffusion length using power spectral densities (PSDs) of the record of a single isotope species can be biased by uncertainties in spectral properties of the isotope signal prior to diffusion. By using a second water isotope and calculating the difference in diffusion lengths between the two isotopes, this problem is circumvented. We study the PSD method applied to two isotopes in detail and additionally present a new forward diffusion method for retrieving the differential diffusion length based on the Pearson correlation between the two isotope signals. The two methods are discussed and extensively tested on synthetic data which are generated in a Monte Carlo manner. We show that calibration of the PSD method with this synthetic data is necessary to be able to objectively determine the differential diffusion length. The correlation-based method proves to be a good alternative for the PSD method as it yields precision equal to or somewhat higher than the PSD method. The use of synthetic data also allows us to estimate the accuracy and precision of the two methods and to choose the best sampling strategy to obtain past temperatures with the required precision. In addition to application to synthetic data the two methods are tested on stable-isotope records from the EPICA (European Project for Ice Coring in Antarctica) ice core drilled in Dronning Maud Land, Antarctica, showing that reliable firn temperatures can be reconstructed with a typical uncertainty of 1.5 and 2 ${}^{\circ }$C for the Holocene period and 2 and 2.5 ${}^{\circ }$C for the last glacial period for the correlation and PSD method, respectively.

Introduction

The stable water isotope ratio in precipitation is known to be related to the local atmospheric temperature e.g.. During evaporation of source water and condensation of water vapour during transport from the source region to the precipitation site, isotopic fractionation occurs. As a result the precipitation water is depleted in the heavy isotopes (${}^{18}\mathrm{O}$ and ${}^{\mathrm{2}}\mathrm{H}$) with respect to the source water. As this effect depends on the temperature difference between source region and precipitation location, the isotope concentration of precipitation water can be generally related to the temperature at which condensation occurred in the atmosphere. This relation has been heavily used in the past on ice cores in polar or alpine regions, providing valuable information on past climate. The interpretation of the isotopic signal in ice cores in terms of temperature is, however, complicated as the spatial gradient between temperature and isotope concentration that has been measured for Greenland and Antarctica is not necessarily applicable for temporal changes . For example, borehole thermometry for Greenland sites points to a much larger temperature difference between the Holocene and the last glacial period than predicted by the water isotope signal using the spatial gradient . Additional evidence for this larger change in temperature is found in the isotope ratios of nitrogen and argon in the air trapped in the ice , which provide independent temperature information during rapid climate variations but do not record long-term changes in local temperature. The difference between the spatial and temporal gradient may be caused by changes in seasonality of precipitation or by a changing source region or transport pathway .

Another method for extracting the temperature signal from the stable water isotope signal is based on the diffusion of the isotope signal in the firn stage. After snow falls onto the surface of an ice sheet, its isotope signal is subject to diffusional smoothing . This smoothing occurs mainly in the firn, where water molecules are in continuous exchange between the vapour and solid phase. showed that this diffusion process is related to the firn temperature and the local accumulation rate. It is therefore possible to reconstruct past climate by estimating the extent to which a layer has been subject to diffusion. This is quantified in terms of the squared diffusion length, which is the average squared displacement of a molecule due to diffusion. This quantity can be obtained from the power spectral density (PSD) of the isotope data by performing a linear fit to part of the spectrum . also suggested to use the difference in diffusion length between two water isotopes ($\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$) for palaeoclimate reconstruction. This difference in diffusion length can be obtained with the PSD method when the ratio of the PSDs of the two isotope signals is taken. The main advantage of this method is that factors other than diffusion that may influence the individual diffusion length are often common to both isotopes. These common factors cancel out to a large extent in the ratio of the PSDs as we will show in section .

One of the main drawbacks of the PSD technique is that it is dependent on parameter choices. In and the PSD is calculated using an autoregressive (AR) model. The resulting PSD spectra are sensitive to the order of autoregression that needs to be chosen in such a model. Furthermore, the slope found by linear regression depends on which part of the spectrum is included in the fit. As we will show, both parameter choices have a significant influence on the estimated diffusion length and can cause a deviation in the reconstructed temperature of several degrees Celsius. A calibration procedure is therefore essential for this method. In this paper we investigate the sensitivity of the determination of the differential diffusion length to these two parameters by applying the method to a large number of synthetic data sets generated for given climatic conditions. This allows us to objectively choose the best set of PSD parameters for a certain ice core section and thereby calibrate the method. Additionally, an estimate of the uncertainty associated with the determination of the diffusion length for the ice core section using this method is obtained.

To circumvent the issues related to the PSD method, we propose an alternative method for finding the differential diffusion length based on the correlation between the $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal. The correlation between the two isotopes of the precipitating snow is very high but decreases in the firn as diffusion progresses. This is due to the higher diffusion rate for $\mathit{\delta }{}^{18}\mathrm{O}$ compared to $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$. Compensating for this effect by diffusing the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal numerically forward in time leads to an optimum in correlation when the diffusion lengths of both isotopes are equal. This method is tested and compared to the PSD method using synthetic data sets. Finally, both methods are applied to an Antarctic ice core section for which climatic conditions are approximately known. To this end the differential diffusion length is translated to firn temperature using a density/water vapour diffusion model similar to that of . Here we also study the effect of impurity-dependent firnification on the derived diffusion temperatures.

Diffusion theory

In the firn stage the isotopic profile changes with time due to the diffusion process and by thinning of the layers due to densification and ice flow. The change in the isotopic concentration can be expressed as $$\mathrm{d}\mathit{\delta }=\frac{\partial \mathit{\delta }}{\partial t}\mathrm{d}t+\frac{\partial \mathit{\delta }}{\partial z}\mathrm{d}z,$$ where $z$ is the vertical coordinate in a coordinate system with its origin fixed to a horizontal layer as it moves downward into the firn and $t$ is time. Using Fick's second law, $$\frac{\partial \mathit{\delta }}{\partial t}={\mathrm{\Omega }}_{\text{f}}\frac{{\partial }^{\mathrm{2}}\mathit{\delta }}{\partial z^{\mathrm{2}}},$$ and the following expression for the vertical strain rate, $$\dot{{\mathit{\epsilon }}_z}=-\frac{\mathrm{1}}{z}\frac{\mathrm{d}z}{\mathrm{d}t},$$ the evolution of the isotope profile ${\mathit{\delta }}_{\mathrm{i}}$ is given by $$\frac{\partial {\mathit{\delta }}_{\mathrm{i}}}{\partial t}={\mathrm{\Omega }}_{\text{fi}}\frac{{\partial }^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{i}}}{\partial z^{\mathrm{2}}}-\dot{{\mathit{\epsilon }}_z}z\frac{\partial {\mathit{\delta }}_{\mathrm{i}}}{\partial z}$$ . Here, ${\mathrm{\Omega }}_{\text{fi}}$ is the firn diffusivity and $\dot{{\mathit{\epsilon }}_z}$ the vertical strain rate. The subscript i indicates the dependency on the isotopologue, where 2 is used for ${}^{\mathrm{2}}{\mathrm{H}}^{16}{\mathrm{O}}^{\mathrm{1}}\mathrm{H}$ and 18 for ${}^{\mathrm{1}}{\mathrm{H}}^{18}{\mathrm{O}}^{\mathrm{1}}\mathrm{H}$. The first term on the right-hand side describes the effect of firn diffusion, whereas the second term accounts for the thinning of the layers. Expressions for the diffusivity were developed by and . Adopting the expression from for the firn diffusivity ${\mathrm{\Omega }}_{\text{fi}}$, we can write $${\mathrm{\Omega }}_{\text{fi}}=\frac{m{p}_{\text{sat}}{\mathrm{\Omega }}_{\text{ai}}}{RT\mathit{\tau }{\mathit{\alpha }}_{\mathrm{i}}}\left(\frac{\mathrm{1}}{\mathit{\rho }},-,\frac{\mathrm{1}}{{\mathit{\rho }}_{\text{ice}}}\right),$$ where $m$ is the molar mass of water; $T$ the firn temperature; $p_{\text{sat}}$ the saturation vapour pressure over ice; $R$ the gas constant; and $\mathit{\rho }$ and ${\mathit{\rho }}_{\text{ice}}$ the density of firn and ice, respectively. Furthermore, the tortuosity $\mathit{\tau }$ accounts for the shapes of the interconnecting pores, and ${\mathrm{\Omega }}_{\text{ai}}$ and ${\mathit{\alpha }}_{\mathrm{i}}$ are the diffusivity of water vapour in air and the ice–vapour fractionation factor, respectively.

Inspection of the terms in Eq. () reveals that firn diffusivity is essentially a function of firn temperature and firn density. The diffusivity in air, the fractionation factor and the saturation vapour pressure can all be parametrised as a function of firn temperature (see and for ${\mathrm{\Omega }}_{\text{ai}}$, and for ${\mathit{\alpha }}_{\mathrm{i}}$, and for $p_{\text{sat}}$). The dominating factor of these temperature-dependent parameters is the saturation vapour pressure. A higher firn temperature leads to a larger number of water molecules that are available for diffusive transport in the vapour phase and thus to a higher firn diffusivity. The dependency on the firn density is clear from the term in brackets in Eq. (), which equals the pore space volume per unit mass. Lower-density firn has a larger pore space and therefore a higher diffusion rate than high-density firn. Additionally, for the tortuosity we adopt the parameterisation used in based on data from , where this quantity is described as a function of firn density. Note that this parameterisation is a simplification, as tortuosity is a complex parameter which is not determined by the density of the firn alone. As a result, the main uncertainty in ${\mathrm{\Omega }}_{\text{fi}}$ is caused by the uncertainty in the tortuosity.

The diffused isotope signal at any time $t$ is related to the initial signal by a convolution with a Gaussian distribution: $$\mathit{\delta }(z,t)=\frac{\mathrm{1}}{\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\mathit{\delta }}_{\mathrm{0}}\left(z^{\prime }\right)\cdot \exp \left(-,\frac{(z^{\prime }-z{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathrm{d}{z}^{\prime },$$ where ${\mathit{\delta }}_{\mathrm{0}}$ is the initial isotope signal corrected for the total strain since deposition of the layer and $\mathit{\sigma }$ the SD of the Gaussian distribution. This parameter $\mathit{\sigma }$ is termed the diffusion length as its square is the average squared displacement of a molecule as a result of diffusion. For a constant diffusivity $\mathrm{\Omega }$ the expectation value of the squared vertical displacement of the molecules evolves with time $t$ as $${\mathit{\sigma }}^{\mathrm{2}}=⟨z^{\mathrm{2}}⟩=\mathrm{2}\mathrm{\Omega }t.$$ In firn the diffusivity is not constant in time, and densification of firn and deformation of ice leads to a thinning of the layers and thereby to a reduction of the diffusion length. Combining these processes, we can write the change in the squared diffusion length $d{\mathit{\sigma }}^{\mathrm{2}}$ in a time interval $\mathrm{d}t$ as $$\mathrm{d}{\mathit{\sigma }}^{\mathrm{2}}=\mathrm{2}{\mathrm{\Omega }}_{\text{fi}}\mathrm{d}t-\mathrm{2}{\dot{\mathit{\epsilon }}}_z\mathrm{d}t.$$ Thus the diffusion length at any depth in the ice is proportional to the integrated diffusivity since the time of precipitation corrected for the vertical strain. When using the Herron–Langway (HL) densification model and a prescribed accumulation rate, the differential equation can be solved by integrating over the time since deposition seefor a full derivation. Evaluating the resulting expression at the pore close-off depth leads to the diffusion length after firn diffusion is completed. This length is a function of firn temperature through the firn diffusivity. Firn diffusion stops at pore close-off after which the diffusion length will decrease due to the compression and deformation of the ice. Diffusion in the ice matrix will cause the length to increase, but this process is much slower than the diffusion in firn and generally becomes significant only for very old or relatively warm ice (e.g. near the bedrock).

This illustrates that together with an independently derived accumulation rate and an estimate of the thinning function a temperature history can be reconstructed by determining the firn diffusion length in the ice. This can in principle be done for either of the water isotopes, but it is better constrained when the combination of two isotopes is used, as will be shown in Sect. . For this reason we use the differential diffusion length $\mathrm{\Delta }\mathit{\sigma }$, which was introduced by and is defined by $$\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}={\mathit{\sigma }}_{18}^{\mathrm{2}}-{\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}.$$

Synthetic data

To be able to verify the accuracy and precision of the reconstructed values for the differential diffusion length, it is necessary to know the true values. As this is impossible for data from existing ice core records, we created synthetic data sets for this purpose. To give statistically robust results, these data sets were created in a Monte Carlo (MC) routine in which several parameters characterising the resulting data set were randomised as is outlined below. In total 4000 data sets were created for each MC run. Within an MC run the parameters that determine the diffusion length (firn temperature, mean accumulation rate and thinning of the ice) were the same for each data set, but parameters such as the amount of precipitation per event and the amplitude of the isotope signal were varied. In this paper six different runs will be discussed. The first MC run assumes climate conditions similar to present-day climate at the European Project for Ice Coring in Antarctica (EPICA) drill site in Dronning Maud Land (EDML) and no vertical strain (thinning $=$ 1). This was chosen as we will apply the two methods to an ice core section from this location. MC runs II and III assume the same climatic conditions but with thinning factors of 0.88 and 0.76, respectively, to investigate the effect of the compression of the layers. Runs IV, V and VI were created to simulate a colder climate at the EDML site, similar to what is estimated for the Last Glacial Maximum (LGM). We choose the annual accumulation for the latter runs to be half of the present-day value and the temperatures to range from 6 to 12 ${}^{\circ }$C lower than present-day temperature. A thinning of 0.6 is used for these runs as this is approximately the thinning for the LGM at EDML. An overview of the values of the parameters and the corresponding diffusion lengths for the six runs is given in Table .

The creation of a synthetic data set consists of several steps. First a precipitation record for both $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ is created, after which the diffusion and densification processes are simulated. Then thinning of the record due to ice flow is taken into account, after which the isotope records are synthetically sampled. Finally, analytical measurement noise is added to account for any measurement uncertainty.

Climatological and glaciological parameters used in the different Monte Carlo runs. The firn temperature, annual accumulation and thinning are input values, leading to the (differential) diffusion lengths given in the last three columns.

MC	$T$	Ann. acc.	thinning	${\mathit{\sigma }}_{18}^{\mathrm{2}}$	${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$	$\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$
run	(${}^{\circ }$C)	(mm ice)		(cm${}^{\mathrm{2}}$)	(cm${}^{\mathrm{2}}$)	(cm${}^{\mathrm{2}}$)
I	$-$44.6	69.8	1.00	49.3	40.8	8.55
II	$-$44.6	69.8	0.88	38.2	31.6	6.62
III	$-$44.6	69.8	0.76	28.5	23.6	4.94
IV	$-$50.6	34.9	0.60	17.4	14.2	3.25
V	$-$53.6	34.9	0.60	13.3	10.8	2.57
VI	$-$56.6	34.9	0.60	10.1	8.10	2.02

In creating the precipitation signals it is assumed that on average the isotope concentration of precipitation water varies sinusoidally in time (with a period of 1 $\mathrm{year}$). Within a given data set the number of precipitation events per year and the timing of these events within a year were varied randomly. Also the amount of precipitation per event was varied randomly, allowing for variations in the annual layer thicknesses, with the restriction that over the full record the average annual accumulation was fixed to the set value for the specific MC run. The average number of precipitation events and amount of precipitation per event was varied between different data sets within an MC run in a wide range around observed data from automatic weather stations (AWSs) at the EDML site . In addition to the seasonal component described by a sine curve, the synthetic $\mathit{\delta }{}^{18}\mathrm{O}$ signal also exhibits a random component, which was achieved by drawing numbers from a Gaussian distribution. The $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal is derived from $\mathit{\delta }{}^{18}\mathrm{O}$, by assuming an 8-times-larger amplitude and a combination of seasonal and random variations in deuterium excess (d-excess), in line with those observed at Neumayer station, coastal Dronning Maud Land . The amplitude of the seasonal variation of $\mathit{\delta }{}^{18}\mathrm{O}$ and the Gaussian distributions for the random components in the signals vary randomly from one data set to the other.

The effects of diffusion and densification on the isotope signals are calculated numerically in discrete time steps. A finite-difference technique is used to calculate the effect of diffusion for each time step, with the firn diffusivity set by Eq. (). To obtain a stable solution for the differential equation, we use an implicit method . The finite-difference technique used here is also used in the correlation method that will be discussed in Sect. . Our analysis has shown that a second-order approximation for the spatial derivative may not be sufficient to obtain reliable results with this method. Depending on the sampling resolution it may be necessary to use a higher-order method. For this reason we use an eighth-order approximation providing precise results. Initially, the density $\mathit{\rho }$ is set to a surface density of 360 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, which is similar to what is observed at the drill site. As a layer moves down into the firn, with a rate determined by the annual accumulation, the density and the layer thicknesses change according to the Herron–Langway model . This process continues until the pore close-off density (${\mathit{\rho }}_{\text{pc}}=804.3$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$; ) is reached. During each time step the diffusion length in ice equivalents is calculated from its value of the previous time step ($t_{m-\mathrm{1}}$) as $${\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{2}}\left(t_m\right)={\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{2}}\left(t_{m-\mathrm{1}}\right)+\mathrm{2}{\mathrm{\Omega }}_{\text{fi}}\mathrm{d}t{\left(\frac{\mathit{\rho }}{{\mathit{\rho }}_{\text{ice}}}\right)}^{\mathrm{2}},$$ where $\mathrm{d}t$ is the time step and ${\mathit{\rho }}_{\text{ice}}$ the density of ice (917 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$). Note that thinning due to densification is accounted for (by multiplying with $(\mathit{\rho }/{\mathit{\rho }}_i{)}^{\mathrm{2}}$) as the diffusion length is calculated in ice equivalents. Furthermore, we assume that thinning due to ice flow is negligible during firnification. Therefore, the second term on the right-hand side of Eq. (), which contains the vertical strain rate, is zero. For deep cores such as EDML the thinning due to ice flow is small in the firn stage. It can be included in the numerical calculation of the diffusion length, but this requires more computational effort. Here, we choose to correct for the thinning afterwards, leading to very similar results. The correction is applied by multiplying the diffusion lengths by the thinning factor.

Finally, from each data set a 20 $\mathrm{m}$ long section is extracted and synthetically sampled. This is done by taking the average of the layers within the sampling interval, weighted with the layer thickness. A random number drawn from a Gaussian distribution is then added to account for measurement uncertainty.

The $\mathit{\delta }{}^{18}\mathrm{O}$ record of one of the 4000 synthetic data sets from MC run I is given in Fig.  together with a 20 $\mathrm{m}$ section of the EDML ice core record. The two signals have very similar spectral properties as is evident from the power spectra given in Fig. b. The light grey lines in this figure refer to power spectra of 200 other randomly selected synthetic data sets out of the total of 4000 data sets from MC run I, which serve to illustrate the variability within the run. At low frequencies the range in the PSDs of the synthetic data sets is relatively large, reflecting the large range of amplitudes used in the creation of the data. The decrease in PSDs with frequency is similar for all data sets, which is a consequence of the fixed climatic conditions and thereby diffusion length within the run. Finally, at high frequencies the PSDs tend to stabilise at a constant value. This is the white-noise baseline caused by the measurement uncertainty, which is assumed to be 0.06 ‰, corresponding to a baseline value in the power spectrum of $1.8\times {10}^{-\mathrm{4}}{\mathrm{‰}}^{\mathrm{2}}\mathrm{m}$.

$\mathit{\delta }{}^{18}\mathrm{O}$ records of a synthetic data set (blue) and of the EDML ice core (red). For both records the sample size is 5 $\mathrm{cm}$. (b) shows the power spectral densities of each isotope record together with those of 200 other synthetic data sets of MC run I. The average PSD of the 200 runs is shown as the dark grey solid line. The dashed line is the average PSD of the initial signals (before diffusion) of these runs.

[Figure omitted. See PDF]

Reconstruction of the differential diffusion length

The extent to which a stable water isotope record has been subject to diffusion can be quantified with the diffusion length $\mathit{\sigma }$. Traditionally this quantity is obtained from the PSD of the isotope record . Determination of the differential diffusion length is done similarly using the ratio of the PSDs of $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ and $\mathit{\delta }{}^{18}\mathrm{O}$. One of the disadvantages of the PSD method is that the resulting values for the diffusion lengths are sensitive to the choice of parameters such as the order of autoregression used in a maximum-entropy calculation of the PSD and the cut-off frequency, where the PSD becomes dominated by the white measurement noise. For this reason we introduce a new method based on the Pearson correlation coefficient between the $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal. This method cannot reconstruct the diffusion length for the individual isotopes but is suitable for, and potentially more robust in, estimating differential diffusion lengths. A detailed description of each method is outlined in the following sections. Both methods are applied to a large number of synthetic data sets to assess their capacity in reconstructing the differential diffusion length accurately and precisely.

PSD method

By taking the square of the Fourier transform of Eq. (), the following expression for the PSD as a function of wave number $k$ is obtained: $$P\left(k\right)=P_{\mathrm{0}}\left(k\right)\cdot \exp (-k^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}),$$ where $P_{\mathrm{0}}\left(k\right)$ is the PSD of the initial signal corrected for compression and $k=\mathrm{2}\mathit{\pi }f$, with $f$ denoting frequency in ${\mathrm{m}}^{-\mathrm{1}}$. From this relation it is clear that diffusion acts as a Gaussian low-pass filter attenuating the high-frequency components more than the low-frequency components. In practice, the PSD of a measured ice core record will not become negligible at high frequencies but remains at relatively low values due to measurement uncertainty (see also Fig. ). For discrete samples measured in randomised order we can assume that the measurement error is uncorrelated, leading to a constant baseline in the frequency spectrum. This baseline needs to be determined either from the frequency spectrum of real ice core measurements or from the known precision of the analytical instrument. After this baseline is subtracted, the net PSD signal can be rewritten as $$\ln \left(P\right(k\left)\right)=\ln \left(P_{\mathrm{0}}\right(k\left)\right)-k^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}.$$ When the initial spectrum $P_{\mathrm{0}}\left(k\right)$ is white (i.e. $P_{\mathrm{0}}$ is independent of $k$), as is argued by , the above equation states that there is a linear relationship in the $\ln \left(P\right)$ vs. $k^{\mathrm{2}}$ space. If this is the case, the squared diffusion length can be obtained from a linear regression of the red part of the spectrum for each isotopic species. This method was recently applied to the North Greenland Ice Core Project (NorthGRIP) ice core in the study of .

To verify the assumption that the initial signal is independent of frequency, we have stacked precipitation data from the Global Network of Isotope in Precipitation (GNIP) database . The database contains the $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ isotope ratios and the amount of precipitation for many stations worldwide. For most stations the data are available in monthly resolution, but for a few stations sampling is performed at higher resolution. In Fig.  the PSD of stacked precipitation data is shown for Debrecen and Vernadsky. Debrecen was chosen as this is one of the few stations with a long record of higher than monthly resolution. Sampling was done on a weekly basis from January 2001 to November 2004 and on a daily basis from November 2004 to December 2009 . The station is located in eastern Hungary, where the annual accumulation is 59 $\mathrm{cm}$. The length of the stacked record used is 5.25 $\mathrm{m}$. Vernadsky (formerly named Faraday) was chosen as it is located on the Antarctic Peninsula and thereby one of the closest stations to Antarctic ice core drill sites. Samples collected at this location are monthly averages. The full Vernadsky record extends from 1964 to 2008, but as this record contains gaps, we use the longest continuous record of 8.7 $\mathrm{years}$ (from 1992 to 2000), for which the average annual accumulation was 48 $\mathrm{cm}$.

The power spectral densities of two stacked precipitation records retrieved from the GNIP database. In light colours the Debrecen data are given; dark colours refer to data from Vernadsky. The figure on the left shows the PSD as a function of frequency, whereas the figure on the right shows the natural logarithm of the PSD as a function of the wave number squared (bottom $x$ axis).

[Figure omitted. See PDF]

The PSD of the Debrecen record shows an annual peak at 1.7 $\mathrm{cycles}{\mathrm{m}}^{-\mathrm{1}}$ and a decreasing trend in the frequency range of 0–5 $\mathrm{cycles}{\mathrm{m}}^{-\mathrm{1}}$. For higher frequencies the PSD is relatively constant. If a linear regression was performed in the first half of the frequency spectrum of Fig. b, the resulting slope would significantly differ from zero. Also the Vernadsky record shows a trend in the PSD, and together with the large annual peak at 2 $\mathrm{cycles}{\mathrm{m}}^{-\mathrm{1}}$ a linear regression would result in a non-zero slope. This implies that the individual diffusion length derived for each isotopic species separately will be biased. The green lines in Fig. b show the ratio of the PSDs of $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ . For this we can write $$\ln \left(\frac{P_{\mathrm{2}}\left(k\right)}{P_{18}\left(k\right)}\right)=\ln \left(\frac{P_{\mathrm{0},\mathrm{\hspace{0.17em}2}}\left(k\right)}{P_{\mathrm{0},\mathrm{\hspace{0.17em}18}}\left(k\right)}\right)+k^{\mathrm{2}}\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}},$$ from which we can see that a linear regression would allow us to retrieve the differential diffusion length $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ if the first term on the right-hand side of this equation is independent of $k$. Since the initial signals $P_{\mathrm{0},\mathrm{\hspace{0.17em}2}}$ and $P_{\mathrm{0},\mathrm{\hspace{0.17em}18}}$ exhibit very similar patterns, their ratio is close to constant. This implies that the characteristics of the individual PSDs are the same for each of the two isotopes and due to the same effects during deposition but that they cancel out when taking the ratio of the PSDs. It is therefore not necessary to have any knowledge about the initial signal, in contrast to the single-isotope diffusion method. For this reason we prefer the differential diffusion method. In the remainder of this section we will therefore focus on reconstructing the differential diffusion length. When the methods are applied to real ice core data, we will include results from the single-isotope method for comparison.

The PSDs of the isotope data can be calculated using the maximum-entropy method (MEM; ), a parametric technique proposed by . The underlying model is an AR model . The order $M$ of this AR process can have a large influence on the resulting frequency spectrum. In Fig. a two spectra calculated from the same synthetic data set but with different orders $M$ are depicted. A higher order leads to a higher frequency resolution, and possibly also to the appearance of spurious peaks. On the other hand, a lower order may lead to too smooth a spectrum. The choice of order $M$ in the calculation of the PSD has an influence on the differential diffusion length obtained by a linear regression. This is clearly visible in Fig. c, which shows the dependance of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ on the order $M$ for the data of Fig. a.

Figure a shows the power spectral densities of the $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signals of a synthetic data set. The dark colours in (a), (b) and (d) refer to PSDs calculated with MEM technique with order $M=20$. The light colours in these figures correspond to PSDs calculated with order 50. In (b) the natural logarithm of the ratio of the PSDs is given as a function of squared wave number. (c) and (d) show the obtained values for the differential diffusion length as a function of the order of AR and cut-off frequency, respectively. The red dashed lines in these figures indicate the true value of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ for this data set.

[Figure omitted. See PDF]

A second parameter that needs to be chosen in the PSD method is the cut-off frequency ($f_{\text{co}}$). When performing the fit of the PSD, only the red part of the frequency spectrum can be considered. For higher frequencies the PSD is dominated by the measurement noise of the isotope signal.

The choice of the cut-off frequency also has a significant influence on the resulting value for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$. This is illustrated in Fig. d. For this data set (from MC run I) a firn temperature of $-$44.6 ${}^{\circ }\mathrm{C}$ and accumulation rate of 69.8 $\mathrm{mm}\mathrm{ice}$ per year was used. This leads to a theoretical squared differential diffusion length of 8.55 ${\mathrm{cm}}^{\mathrm{2}}$, indicated with the red dashed line in Fig. c and d. For these conditions a deviation in $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ of 1 ${\mathrm{cm}}^{\mathrm{2}}$ is equivalent to a deviation in temperature of 1.6 ${}^{\circ }\mathrm{C}$. This illustrates both the importance of using the right parameter values for a certain data set and the large uncertainty associated with this method. In the following we will describe a method to chose the parameter values objectively and thereby reduce the uncertainty of this method.

In order to find the optimum values for these two parameters, the PSD method was applied to the 4000 synthetic data sets within each Monte Carlo run. As the synthetic data sets within an MC run all refer to the same climatological and glaciological conditions, this approach will also give insight into how large the stochastic variability is in the method. In Fig.  the deviation of the mean of these synthetic data sets from the theoretical value (${\mathrm{\Delta }}_{\text{th}}$) and the SD ($\mathit{\varsigma }$) are depicted as a function of AR order $M$ and cut-off frequency. An upper limit for the cut-off frequency for a certain data set is found at the frequency at which one of the PSDs reaches a negative value after subtracting the baseline. The grey shading indicates the percentage of data sets for which the differential diffusion length can be calculated, i.e. where the cut-off frequency was lower than this upper limit.

Contour plots showing the average deviation from the theoretical value (left) and the SD of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ for all ensemble members of Monte Carlo run I as a function of the order of autoregression and the cut-off frequency. The grey shading indicates the percentage of data sets for which the differential diffusion length could be evaluated. The position in the two-parameter space with the lowest total error (as defined in the main text) is indicated by the red star.

[Figure omitted. See PDF]

From these figures it becomes clear that for most solutions the accuracy (the deviation from the theoretical value) is much higher than the precision (the SD). The large SD shows how much the obtained value for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ can vary for different realisations (different ice core records) under the same climatological conditions. In practice, when applying the method to a real ice core record, we only have one realisation for which the method will then generate a value for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ within the range determined by the SD. From Fig.  we also note that the highest precision and accuracy are not necessarily at the same location in the two-parameter space. To find the optimum choice of both parameters, we therefore define a total error (TE) as a quadratic sum of the two uncertainties (TE $=$ ${\left({\mathrm{\Delta }}_{\text{th}}^{\mathrm{2}},+,{\mathit{\varsigma }}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}$). The optimum values for the different Monte Carlo runs are listed in Table . For all runs a low AR order is preferred. This is most likely due to less scattering in the PSD, which leads to a higher precision than those PSDs created with a higher AR order. For runs with a low $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ the diffusion has a weaker effect on the high-frequency part of the spectrum than those with high $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$. As a consequence a larger part of the spectrum can be included in the fit, which explains the larger cut-off frequencies for these runs. Using the optimum values, the offset to the theoretical value is less than 0.3 ${}^{\circ }\mathrm{C}$ for all MC runs, which can be neglected compared to the total error, which is typically 2–2.5 ${}^{\circ }\mathrm{C}$.

The optimum parameters and the resulting mean values, deviations from theoretical values and SDs for the differential diffusion lengths of the Monte Carlo runs obtained with the PSD method for 5 $\mathrm{cm}$ samples. TE is the total error as described in the main text. The total error and the offsets from the real values are also converted to the corresponding errors in firn temperature for each run.

MC run	$\mathrm{\Delta }{\mathit{\sigma }}_{\text{th}}^{\mathrm{2}}$	$f_{\text{co}}$	AR	$\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$	Offset	SD	TE	offset	TE
	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{m}}^{-\mathrm{1}}$)	order	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${}^{\circ }\mathrm{C}$)	(${}^{\circ }\mathrm{C}$)
I	8.55	4.7	20	8.39	$-$0.16	1.14	1.16	$-$0.3	1.9
II	6.62	5.4	20	6.52	$-$0.10	0.87	0.87	$-$0.2	1.9
III	4.94	5.9	20	4.86	$-$0.08	0.67	0.67	$-$0.2	1.9
IV	3.25	6.5	20	3.24	$-$0.01	0.65	0.65	0.0	2.6
V	2.57	7.1	20	2.57	0.00	0.51	0.51	0.0	2.6
VI	2.02	8.5	20	2.08	0.05	0.37	0.37	0.3	2.2

Correlation method

In this section a different method for estimating the differential diffusion length from the stable water isotope records is presented. The main goal is to come up with a method that is not sensitive to parameter choices and therefore gives a more robust estimate of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$. The technique is based on the Pearson correlation coefficient $R$ between the $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal. In precipitation records this correlation coefficient is very close to 1, as the processes that are responsible for the isotope signals are the same for both isotopes. For example, the correlation coefficients for the precipitation data of Debrecen and Vernadsky used in Fig.  are 0.95 and 0.96, respectively. After deposition the different diffusion rates of the two isotopic molecules lead to a decorrelation in the firn with time. We will show that this effect can be used to derive the differential diffusion length. To illustrate this, we first consider an ideal synthetic data set in which the correlation of the signals before diffusion is perfect ($R=\mathrm{1}$). This is accomplished by deriving a  $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ record from the synthetic $\mathit{\delta }{}^{18}\mathrm{O}$ signal using the meteoric water line : $$\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}=\mathrm{8}\cdot \mathit{\delta }{}^{18}\mathrm{O}+10\mathrm{‰}.$$ The d-excess signal is then by definition constant at 10 ‰. These isotope signals are then diffused numerically, simulating firn diffusion at conditions similar to those currently at EDML (as in MC run I). Due to the isotope-dependent diffusivity the resulting $\mathit{\delta }{}^{18}\mathrm{O}$ profile is smoothed more than the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ profile. This is also reflected in the diffusion lengths at pore close-off depth of 7.0 $\mathrm{cm}$ ice equivalent (i.eq.) for $\mathit{\delta }{}^{18}\mathrm{O}$ and 6.4 $\mathrm{cm}$ i.eq. for $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ . This difference in diffusion rate also leads to a d-excess signal that is no longer constant with depth, as can be seen in Fig. a, and to a lowering of the correlation coefficient to a value of 0.9947. Note that the d-excess signal in Fig.  is now entirely created by the diffusion process and is not an atmospheric signal. This also implies that seasonal d-excess variations in real firn and ice cores are at least partly created by diffusion . This explains the phenomenon observed in ice cores that annual layer counting in $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ is limited at low-accumulation sites due to the diffusive loss of signal, while (diffusion created) seasonal variation in d-excess can still be counted (see e.g. ). If the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal is artificially diffused further relative to $\mathit{\delta }{}^{18}\mathrm{O}$ using the numerical diffusion method (see above), the correlation between the $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal increases again until the diffusion lengths of the two signals are equal: $${\mathit{\sigma }}_{18}^{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{2},\text{add}}^{\mathrm{2}},$$ where ${\mathit{\sigma }}_{\mathrm{2},\text{add}}^{\mathrm{2}}$ is the additional numerical diffusion length. When the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal is diffused beyond this point, the correlation drops, as is evident in Fig. b. Thus a maximum in correlation coefficient occurs when the total diffusion length (firn diffusion and numerical diffusion) of $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ equals the firn diffusion length of $\mathit{\delta }{}^{18}\mathrm{O}$. The additional diffusion necessary to reach the maximum correlation is equivalent to the differential diffusion value that is sought. That is, $${\mathit{\sigma }}_{\mathrm{2},\text{add}}^{\mathrm{2}}=\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}.$$ Note that at the point of maximum correlation the d-excess signal is also constant again in the idealised synthetic data set (see dashed line in Fig. a).

Illustration of the correlation method for ideal conditions in which the correlation coefficient of the isotope data before diffusion is 1. The $\mathit{\delta }{}^{18}\mathrm{O}$, $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ and d-excess signals after firn diffusion has taken place are depicted by the solid lines in the left figure. The dashed lines show the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ and d-excess profiles after the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ profile is diffused further until the total $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ diffusion length equals that of $\mathit{\delta }{}^{18}\mathrm{O}$. In the figure on the right the correlation coefficient between $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ is given as a function of the additional diffusion length as the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ profile is diffused further.

[Figure omitted. See PDF]

In reality the initial d-excess signal will not be constant, and therefore the initial correlation will be below 1. However, for more realistic initial signals a decrease in correlation as a result of firn diffusion is also to be expected. To study the correlation method for more realistic conditions, it was applied to the synthetic data sets, which have varying d-excess with both a seasonal component and a random component. For all these data sets the maximum in correlation coefficient is found as a function of the additional diffusion length in the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal. This leads to an approximate Gaussian distribution of obtained diffusion lengths with an average of 8.4 ${\mathrm{cm}}^{\mathrm{2}}$ and SD of 0.7 ${\mathrm{cm}}^{\mathrm{2}}$. From this, it is clear that the introduction of varying d-excess in the initial signal causes a spread in the resulting values of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$, but the mean value is still close to the theoretical value (8.55 ${\mathrm{cm}}^{\mathrm{2}}$).

For the data sets discussed above the measurement uncertainty has not been included in the synthetic data yet. The effect of adding measurement noise is shown in Fig. a and b, where the mean and SD of the distributions for a certain run are plotted as a function of the measurement uncertainty. As the uncertainty in the measurement increases, the distribution of the obtained diffusion lengths of each data set becomes wider, indicating that the precision decreases. For $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ this is accompanied by a deviation towards higher values. This occurs only for $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ because only the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ record is numerically diffused forward in time to maximise the correlation. Accordingly, the method loses accuracy as the uncertainty in the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ measurement increases.

The mean and SD of the differential diffusion length of the 4000 data sets of Monte Carlo run I as a function of the uncertainty in the isotopic measurement are depicted in the top figures. The uncertainty is given as the SD of a Gaussian distribution centred around zero from which random numbers are drawn. The figure on the bottom shows the absolute measurement deviation (as described in the main text) for two different sampling resolutions. The data points given here are obtained from all six Monte Carlo runs and for several assigned values for the uncertainty in the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ measurement. The solid lines represent linear regressions to the data points which can be used to correct $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ for the deviation that arises as a result of measurement uncertainty.

[Figure omitted. See PDF]

To study this effect in detail, we arbitrarily selected 20 data sets from each Monte Carlo run. For these data sets the differential diffusion length is first found for the records without measurement uncertainty, indicated by $(\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}{)}_{\mathrm{0}}$ in the following. Next, measurement uncertainty in $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ is added to this data set by adding a number drawn from a Gaussian distribution centred around zero to each isotopic value. The correlation method is then applied again, leading to a new value for the differential diffusion length. This is repeated 100 times for each data set and each uncertainty to ensure statistically robust results, leading to a mean differential diffusion length including measurement noise indicated by $(\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}{)}_m$. The results of this analysis are given in Fig. c. The $y$ axis of this figure gives the absolute measurement deviation (AMD), which we define as $$\begin{array}{cc}& \text{AMD}=\text{sgn}\left((,\mathrm{\Delta },{\mathit{\sigma }}^{\mathrm{2}},{)}_m,-,(,\mathrm{\Delta },{\mathit{\sigma }}^{\mathrm{2}},{)}_{\mathrm{0}}\right)\\ & & \cdot {\left(\left|(,\mathrm{\Delta },{\mathit{\sigma }}^{\mathrm{2}},{)}_m,-,(,\mathrm{\Delta },{\mathit{\sigma }}^{\mathrm{2}},{)}_{\mathrm{0}}\right|\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$ The quantity on the $x$ scale can be regarded as a noise-to-signal ratio of the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ signal. The uncertainty in the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ measurement $s_{\mathrm{2}}$, as given by the SD of the Gaussian distribution from which errors are drawn, is divided by the average change in $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ between neighbouring points $d_{\mathrm{2}}$: $$d_{\mathrm{2}}=\frac{\mathrm{1}}{N-\mathrm{1}}\sum _{n=\mathrm{2}}^N\left|\mathit{\delta },{{}^{\mathrm{2}}\mathrm{H}}_n,-,\mathit{\delta },{{}^{\mathrm{2}}\mathrm{H}}_{n-\mathrm{1}}\right|,$$ with $N$ being the number of data points in the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ record. Using these scales a linear relationship appears for each sampling resolution. This shows that the AMD is a function of the sampling resolution and the measurement uncertainty. A linear regression through all data points gives a relation that can be used to correct the offset caused by the measurement uncertainty. Note that this correction function is independent of the MC run. In Fig. b the open symbols represent the mean and SD of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ after this correction. Applying the method to each of the Monte Carlo runs leads to the results given in Table . As was done in the PSD method tests, realistic values for measurement uncertainty of 0.06 and 0.40 ‰ are used for $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ , respectively. Compared to the PSD method, the accuracy of the correlation method after correction is lower, with a systematic offset of up to 1 ${}^{\circ }\mathrm{C}$ for LGM conditions. The total error, however, is still dominated by the SD and is slightly smaller than that of the PSD method. When the correlation method is used, the firn temperature can be reconstructed with a total uncertainty of 1.5–2.0 ${}^{\circ }\mathrm{C}$.

The mean values and deviation from theory and SD for the differential diffusion lengths of the Monte Carlo runs obtained with the correlation method for 5 $\mathrm{cm}$ samples. Measurement uncertainties of 0.06 for $\mathit{\delta }{}^{18}\mathrm{O}$ and 0.40 for $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ are used. Column 6 and 8 give the total error in the squared differential diffusion length and in the resulting firn temperature, respectively.

MC run	$\mathrm{\Delta }{\mathit{\sigma }}_{\text{th}}^{\mathrm{2}}$	$\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$	Offset	SD	TE	Offset	TE
	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${}^{\circ }$C)	(${}^{\circ }$C)
I	8.55	8.59	0.04	0.90	0.90	0.1	1.5
II	6.62	6.68	0.05	0.68	0.68	0.1	1.4
III	4.94	5.01	0.07	0.52	0.52	0.2	1.5
IV	3.25	3.40	0.15	0.46	0.49	0.6	1.9
V	2.57	2.72	0.14	0.38	0.41	0.7	1.9
VI	2.02	2.19	0.16	0.31	0.35	1.0	2.0

Comparison of the methods

In the following section we will compare the performance of the two methods discussed above in terms of reliably reconstructing the differential diffusion length. To assess the accuracy of each method, we compare the mean value for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ of the 4000 data sets for each run with the theoretical value set by the climatological and glaciological parameters for the run. The precision is quantified by the SD of the 4000 values for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$. Values for accuracy and precision for all the MC runs are given in Tables  and  as absolute errors. In the discussion below, however, we will use relative numbers (expressed as percentage of the theoretical value) rather than absolute ones. This is also useful in interpreting the deviations and uncertainties in terms of the firn temperature that can be derived from the differential diffusion length. The relation between temperature and $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ is non-linear, but for these data sets a deviation of 10 % of the true value in $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ corresponds to an offset of roughly 1.3 ${}^{\circ }\mathrm{C}$ in the reconstructed temperature.

In general the total error in each of these MC runs is dominated by the SD as this value is significantly larger than the deviation from the theoretical value. However, we will discuss both accuracy and precision separately because the importance of each quantity might vary with application. For example, in quantifying the temperature difference between the Holocene and LGM the accuracy becomes more important. In contrast, in situations where temperature changes are small (for example within the Holocene or around the LGM), the precision is the crucial factor. As will be shown below, the precision can also be influenced by the sampling strategy or replicate analysis, whereas the accuracy is not affected by this.

In terms of accuracy, the PSD method performs slightly better than the correlation method with a difference between the mean and theoretical value of less than 2.5 % for all runs, i.e. a systematic error in the reconstructed temperature of less than 0.3 ${}^{\circ }\mathrm{C}$. Such accuracy is also obtained with the correlation method for the Monte Carlo runs I, II and III, but for the other runs higher deviations up to 7 % (equivalent to up to 1.0 ${}^{\circ }\mathrm{C}$) are found. On the other hand, the precision for the PSD method is slightly lower than for the correlation method. The one SD using the PSD method for Monte Carlo runs I, II and III is approximately 13 % of the mean value, whereas it is up to 20 % for Monte Carlo runs IV, V and VI (corresponding to 1.9 and 2.6 ${}^{\circ }\mathrm{C}$, respectively). The correlation method results in better precision: an uncertainty of approximately 10 % of the theoretical value for runs I, II and III and about 15 % for runs IV, V and VI (corresponding to 1.5 and 1.9 ${}^{\circ }\mathrm{C}$, respectively).

The numbers given in the discussion above are based on analysis of 20 $\mathrm{m}$ long ice core sections sampled at a resolution of 5 $\mathrm{cm}$. Table  shows how the precision changes for MC runs I and V when section length and sampling resolution are varied. For both methods we see an increase in uncertainty of about 50 % when section length is halved (keeping the same resolution). Sampling the ice core at higher resolution results in an increase in precision. The uncertainty of the PSD method is reduced by 30 to 43 % when sampling is performed at 2.5 $\mathrm{cm}$ instead of 5 $\mathrm{cm}$ resolution. Reducing the sample size to 1 $\mathrm{cm}$ instead of 5 $\mathrm{cm}$ leads to a decrease in uncertainty of 50 to 65 % for these MC runs. Although it is not as large as for the PSD method, also for the correlation method a significantly higher precision is obtained with higher sampling resolution. Sampling at 2.5 $\mathrm{cm}$ leads to a reduction in uncertainty of 15 to 43 % relative to 5 $\mathrm{cm}$ sampling. For 1 $\mathrm{cm}$ sampling this reduction amounts to 23 to 52 %. The precision at which the firn temperature can be reconstructed with diffusion thermometry thus strongly depends on the sampling strategy. Analysis of synthetic data can provide essential information in assessing the optimal sampling strategy for a certain ice core section.

Sensitivity of the two methods to sample size $\mathrm{\Delta }z$ and section length $L$ given for Holocene (MC run I) and glacial (MC run V) conditions at the EDML drill site location.

			PSD method				correlation method
	$\mathrm{\Delta }z$	$L$	$\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$	Offset	SD	TE	$\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$	Offset	SD	TE
	(cm)	(m)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)	(${\mathrm{cm}}^{\mathrm{2}}$)
MC run I	1.0	10	8.36	$-$0.20	0.90	0.92	8.44	$-$0.12	1.03	1.04
	1.0	20	8.40	$-$0.15	0.61	0.63	8.43	$-$0.13	0.71	0.72
	1.0	40	8.39	$-$0.17	0.43	0.46	8.40	$-$0.15	0.50	0.53
	2.5	10	8.36	$-$0.19	1.19	1.21	8.38	$-$0.18	1.09	1.10
	2.5	20	8.39	$-$0.17	0.80	0.82	8.37	$-$0.19	0.76	0.78
	2.5	40	8.39	$-$0.16	0.56	0.59	8.35	$-$0.20	0.54	0.58
	5.0	10	8.41	$-$0.15	1.76	1.77	8.63	0.07	1.34	1.34
	5.0	20	8.39	$-$0.16	1.14	1.16	8.59	0.04	0.90	0.90
	5.0	40	8.40	$-$0.16	0.77	0.79	8.56	0.01	0.63	0.63
MC run V	1.0	10	2.55	$-$0.02	0.27	0.28	2.62	0.05	0.27	0.27
	1.0	20	2.55	$-$0.02	0.19	0.19	2.61	0.04	0.18	0.19
	1.0	40	2.56	$-$0.01	0.14	0.14	2.62	0.05	0.13	0.14
	2.5	10	2.56	$-$0.02	0.44	0.44	2.67	0.09	0.31	0.33
	2.5	20	2.56	$-$0.02	0.29	0.29	2.65	0.08	0.22	0.23
	2.5	40	2.56	$-$0.01	0.20	0.20	2.66	0.08	0.16	0.18
	5.0	10	2.62	0.04	0.74	0.74	2.73	0.16	0.55	0.57
	5.0	20	2.57	0.00	0.51	0.51	2.72	0.14	0.38	0.41
	5.0	40	2.58	0.01	0.33	0.33	2.71	0.14	0.26	0.30

From the analysis of synthetic data it follows that the influence of measurement uncertainty on the precision of the two methods is small compared to the influence of section length and sampling resolution. For the PSD method, measurement uncertainty causes a baseline in the frequency spectra, which needs to be subtracted before the ratio of the two spectra is calculated. But as the influence of this subtraction is very small in the low frequency part of the spectrum, it has negligible influence as long as the cut-off frequency is chosen correctly. In the correlation method the measurement uncertainty initially has a significant influence on the resulting differential diffusion length as is evident from Fig. . However, after applying the signal-to-noise correction discussed in the previous section, the increase in uncertainty in $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ with increasing isotopic measurement uncertainty is very small. It is thus clear that for both methods knowledge of the measurement uncertainty is required. In general, this knowledge is available from routine measurements of standard waters with the same analytical instruments. Additionally, if the sampling resolution is high enough, the measurement uncertainty can also be estimated from the baseline in the PSD spectrum.

Conversion of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ to firn temperature

To convert the obtained differential diffusion length to firn temperature, we use a model similar to those described in and . Input parameters for this model include the accumulation rate and altitude of the precipitation site as well as the thinning of the ice since deposition. Reconstructions for accumulation rate and total thinning are available from the Antarctic ice core chronology AICC2012, . The altitude of the precipitation site is important as it is used to estimate the air pressure, which influences the diffusivity in the pore space. Estimates of the altitude of the site at the time of precipitation are obtained from ice flow models reconstructing the past ice flow. For the EDML core the 3-D ice flow model of is used. For each of these parameters we calculate the 20 $\mathrm{m}$ running means, which serve as input to the model. The modelled diffusion lengths are calculated for a range of mean surface temperatures with a resolution of 0.1 ${}^{\circ }\mathrm{C}$. The temperature corresponding to each section is then found by linearly interpolating the modelled values to the value obtained from the isotope data.

showed that using seasonal variations in temperature that are present in the firn to a depth of approximately 15 $\mathrm{m}$ leads to a small deviation compared to simply taking the annual mean surface temperature. For this reason the temperature variation is included in the model. The amplitude of the assumed sinusoidal temperature signal in the model is 13.5 ${}^{\circ }\mathrm{C}$, in line with 2 $\mathrm{m}$ temperature measurements from AWSs at EDML . Below the surface, temperature of the firn is calculated using heat conduction p. 206.

A very sensitive parameter in the calculation of the diffusion length is the density of the firn as this determines the volume of the pore space through which most of the diffusion takes place. In the calculation of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ at pore close-off the HL model is used to describe density as a function of depth, with temperature and accumulation rate as fixed parameters. Recently have shown that densification rate may also be influenced by the concentration of impurities in the firn. A larger concentration of impurities would lead to a lowering of the activation energy, thereby increasing the densification rate. As the densification process is faster in that case, the time for diffusion in the firn stage is shorter, resulting in shorter diffusion lengths. Especially when studying glacial time periods for which the impurity concentration can be 1 to 2 orders of magnitude larger than present day this would have a large effect on the modelled diffusion lengths. However, also for the data we present here, which only span part of the Holocene, this effect would not be negligible. Although this impurity effect is not seen in some other ice cores (see e.g. ), and therefore still under debate, we will include it in our diffusion model to investigate its possible effect on the reconstructed temperature.

Figure  shows the density–depth profile for the B32 core drilled at the same site as the deep EDML ice core. The density is measured using three different techniques: (1) volumetric method where density is calculated from measurement of volume and mass of core pieces, (2) gamma ray attenuation and (3) radioscopic imaging. The density profiles obtained with the volumetric and gamma ray methods are almost identical, whereas the radioscopic method results in density values that are on average approximately 10 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ lower. In all cases, however, the measured densities in the depth range 30–90 $\mathrm{m}$ are systematically lower than predicted by the classical HL model without impurity effect. To incorporate the impurity effect in the HL model, use the calcium (${\mathrm{Ca}}^{\mathrm{2}+}$) concentration as a parameter to reflect the total impurity concentration. With this parameter they then modify the activation energy. used the high-resolution radioscopic imaging data together with a high-resolution ${\mathrm{Ca}}^{\mathrm{2}+}$ record to determine the best parameters in their calcium-adjusted HL model. We will use the same parameters, but with the average ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration instead of the high-resolution record as in general a full seasonally resolved record may not be available. The resulting profile is shown in Fig. a (dashed red line) together with the classical HL density profile (dashed blue line). The new density profile is lower than the classical HL profile, which indicates that the impurity concentration for this location is lower than the mean impurity effect that is already implicitly included in the classical HL model. It is clear that the modified HL model agrees much better with the data than the classical HL model.

Density–depth profiles for the B32 ice core drilled close to the main EDML ice core (a). Density was measured volumetrically, using gamma ray attenuation and with radioscopic imaging. The dashed lines show the density profiles calculated with the classical Herron–Langway model and with the Herron–Langway model adjusted using the average ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration of 1.7 $\mathrm{ppb}$ . (c) shows the differential diffusion lengths calculated using the classical model and the ${\mathrm{Ca}}^{\mathrm{2}+}$-adjusted model for two different temperatures. The visible decreasing trend is mostly due to the progressive thinning of the ice. The ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration used in the calculations is shown in (b).

[Figure omitted. See PDF]

A 20 $\mathrm{m}$ running mean of the ${\mathrm{Ca}}^{\mathrm{2}+}$ record is calculated (Fig. b) and used as input for the isotope diffusion model. Figure c shows the model calculations for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ for two different temperatures. Blue lines refer to the model calculations with the classical HL density profile, and red lines refer to the model that includes the impurity effect. It is clear that higher values of the differential diffusion length are obtained with the ${\mathrm{Ca}}^{\mathrm{2}+}$-adjusted HL model, which agrees with the observation of a lower-density profile. For the period in the EDML section analysed here, whether or not the impurity effect is included in the model leads to a difference of more than $\mathrm{1}{}^{\circ }\mathrm{C}$ in the temperature reconstruction. For glacial periods this would be several degrees Celsius as the ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration is at least 1 order of magnitude higher. This illustrates the importance of a good understanding of the densification process for past-temperature reconstructions using isotope diffusion thermometry.

Figure  shows the main result of the model calculations. Contour lines for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ are depicted as a function of the average annual temperature and accumulation rate for both models. This figure can be used as a calibration field for the Holocene part of the EDML ice core record. With the differential diffusion length determined by either of the two methods described in this paper and an independent estimate for the accumulation rate, the temperature of the firn can be determined from this plot. Note that this calibration field is calculated using fixed values for the altitude of the site and ${\mathrm{Ca}}^{\mathrm{2}+}$ concentrations and is therefore not valid for other time periods or other ice core records.

Contour lines of the squared differential diffusion length ($\mathrm{in}{\mathrm{cm}}^{\mathrm{2}}$) as a function of firn temperature and accumulation rate for the EDML ice core site. Values for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ at pore close-off are calculated with the classical HL model (blue lines) and with the modified HL model. For the latter a ${\mathrm{Ca}}^{\mathrm{2}+}$ concentration of 2.2 $\mathrm{ppb}$ is used, which corresponds to the average concentration in the 123–289 $\mathrm{m}$ depth interval of the EDML ice core. $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ is given in ice equivalents assuming a thinning of 1. The star indicates the present-day temperature and accumulation rate at the drill site. Note that the accumulation rate is plotted on a log scale.

[Figure omitted. See PDF]

Application to real ice core records

In Sect.  the accuracy and precision of both methods for finding the differential diffusion lengths were assessed for synthetic data with spectral properties very close to real ice core conditions at EDML. In the following, we will apply the methods to part of the EDML ice core record to investigate how the methods perform when applied to real ice core data. The data, in part published previously in , spans a section of 166 $\mathrm{m}$ from the Holocene period (depth: 123–289 $\mathrm{m}$). This section is selected as this period is believed to represent a relatively stable climatological period with a temperature and accumulation rate similar to present-day conditions. The methods are applied using a 20 $\mathrm{m}$ running window. For this part of the ice core record this window length corresponds to about 300 $\mathrm{years}$.

A first requirement for both methods to be applied is that both isotope records ($\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$) are continuous over the whole section that is analysed. In practice, however, sample loss can occur, which leads to gaps in one or both of the records. The most straightforward approach to fill these data gaps is to interpolate between the neighbouring points in the individual isotope record. However, this can lead to unrealistic values for the deuterium-excess signal. From our experience with the EDML record, the best method for gap filling is a spline interpolation of the d-excess signal and calculating the missing isotope value from the interpolated d-excess and the other isotope species. Using this method prevents unrealistic outliers in both the d-excess and the primary isotopic record. If a data gap occurs in both isotopic records, gap filling can only be done using the primary records. In these cases we have applied interpolations using the spline method. The EDML record we present here consists of 3320 data points, with 14 data gaps in the $\mathit{\delta }{}^{18}\mathrm{O}$ record. The $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ record has six data gaps, all coinciding with a data gap in $\mathit{\delta }{}^{18}\mathrm{O}$.

After applying this gap-filling procedure, the two methods for determining the differential diffusion length can be applied. One of the first steps in the PSD method is to determine the noise level that needs to be subtracted from the PSD. The determination of the base line level from the PSD spectrum is, however, very sensitive to outliers in the isotope record. In Fig.  a section of the EDML $\mathit{\delta }{}^{18}\mathrm{O}$ record is depicted. The PSD spectra are calculated for the isotope data in a running window of 20 $\mathrm{m}$ length. The baseline level is calculated for each of the PSDs by calculating the average of the PSD for frequencies above 8 $\mathrm{cycles}{\mathrm{m}}^{-\mathrm{1}}$. When this level is plotted as a function of depth (the blue line in the bottom panel of Fig. ), two sudden step changes in the baseline level appear. These step changes are in the opposite direction and exactly 20 $\mathrm{m}$ apart, suggesting that a single value in the $\mathit{\delta }{}^{18}\mathrm{O}$ record is causing them. Deleting this particular data point and replacing it with a value obtained with a spline interpolation between the neighbouring data points causes the step changes to disappear, as indicated by the red curve in the same figure. This illustrates how sensitive this method of estimating the measurement uncertainty is to outliers. For both methods to determine the differential diffusion length the measurement uncertainty is an important parameter. A single outlier in either or both of the isotope records can, thus, have a significant effect on the resulting value for $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$. Careful sample handling and accurate and precise isotope measurements are therefore essential for a reliable reconstruction of past temperatures using diffusion thermometry.

Part of the $\mathit{\delta }{}^{18}\mathrm{O}$ record (top figure) and the corresponding noise level (bottom figure) that is determined as the average PSD level for frequencies higher than 8 $\mathrm{cycles}{\mathrm{m}}^{-\mathrm{1}}$. The noise level is calculated with a window of 20 $\mathrm{m}$ length. The sudden steps in the blue profile are caused by a single data point. Replacing this data point (see inset) by a value obtained by interpolating the neighbouring points leads to a much lower noise level (red curve).

[Figure omitted. See PDF]

In our EDML record consisting of 3320 samples, we identified 11 samples as outliers in the $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$ record based on the noise level in the PSD spectra of 20 $\mathrm{m}$ sections. For 8 of these 11 samples the $\mathit{\delta }{}^{18}\mathrm{O}$ value was also considered an outlier. The fact that most outliers were found in both isotopes suggests that it is not the measurement itself that led to these outliers, but rather errors in handling of samples (contamination, water vapour exchange or misnumbering) prior to or during the measurement.

Instead of using the noise levels obtained from the PSD, it is also possible to use a fixed value for the measurement uncertainty based on repeated measurements of standards. This may lead to more robust results and is necessary in case the signal at higher frequencies is of the same or higher magnitude as the noise level. This will be the case for short diffusion lengths or relatively low resolution measurements. In Fig.  we show the differential diffusion lengths for the EDML data using both varying noise levels and a fixed noise level. The chosen fixed noise levels correspond to an uncertainty of 0.04 and 0.45 (1$\mathit{\sigma }$ SD) for $\mathit{\delta }{}^{18}\mathrm{O}$ and $\mathit{\delta }{}^{\mathrm{2}}\mathrm{H}$, respectively, which agrees well with what we expected from repeated standard measurements. For the PSD method a cut-off frequency of 5.0 $\mathrm{cycles}{\mathrm{m}}^{-\mathrm{1}}$ and AR order of 20 was chosen in line with the optimal parameters found by synthetic data sets (MC runs I and II). Here we apply the same cut-off frequency for the calculation with fixed and variable noise level to enable comparison. However, a higher noise level would naturally lead to a lower cut-off frequency as the signal-to-noise ratio is lower in such a case. Comparing the results with variable and constant noise level, we can conclude that for most of the record the effect of a changing noise level is negligible. Only in the depth section 250–270 $\mathrm{m}$, where the uncertainty in the measurement as determined by the PSD is close to double the assumed uncertainty in the fixed noise level scenario, do the two curves deviate significantly from each other.

The squared differential diffusion lengths obtained with the PSD and correlation method (a). For both methods a calculation is done with a fixed noise level and with a noise level determined by the baseline in the PSD spectra (b and c).

[Figure omitted. See PDF]

Finally, local surface temperatures for the EDML record are reconstructed by combining the values of $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ obtained from the isotope records with those obtained from the densification/diffusion model described in Sect. . Figure  shows this temperature record, with 1$\mathit{\sigma }$ and 2$\mathit{\sigma }$ uncertainty bands for both methods based on the uncertainty in $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ in the MC runs (see Tables  and ). For most of the record the two methods agree within their 1$\mathit{\sigma }$ uncertainty interval, but for the first 60 $\mathrm{m}$ (130–190 $\mathrm{m}$ depth) the difference between them is larger. Here, the temperature reconstructed with the correlation method is on average 3 ${}^{\circ }\mathrm{C}$ higher than with the PSD method. Note that the recent mean annual temperature measured at the EDML site of $-$44.6 ${}^{\circ }\mathrm{C}$ is always within 2$\mathit{\sigma }$ error of the diffusion temperature reconstruction for both methods. We attribute the wiggles in the reconstructed temperature over the entire 166 $\mathrm{m}$ ice core section and the differences between the results of the PSD and the correlation method to the sensitivity of the methods to imperfections in the data quality. Taking the average reconstructed temperature over the entire 166 $\mathrm{m}$ ice core section leads to a mean $\pm$ SD of the reconstructed temperature of $-$45.3 $\pm$ 1.3 ${}^{\circ }\mathrm{C}$ for the PSD and $-$43.5 $\pm$ 1.6 ${}^{\circ }\mathrm{C}$ for the correlation method, which agrees well with the recent mean annual temperature at EDML. Note that the SD over the entire record also agrees favourably with the uncertainty determined by our calibration with synthetic data for each 20 $\mathrm{m}$ section, providing an additional consistency check of our approach.

Local surface temperature reconstructions based on the differential diffusion lengths given in Fig.  and the model described in the main text. The shaded bands give the 1$\mathit{\sigma }$ and 2$\mathit{\sigma }$ uncertainties and are based on the uncertainty in the differential diffusion lengths. Also given are temperature reconstructions obtained with the PSD method applied to single isotopes (green lines). The recent measured annual temperature at EDML is $-$44.6 ${}^{\circ }\mathrm{C}$. The age scale on the top $x$ axis is the AICC2012 age scale.

[Figure omitted. See PDF]

For comparison we also show the temperature records that would be obtained if the PSD method were applied to single-isotope records. Although we obtain on average higher temperatures with both ${\mathit{\sigma }}_{18}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}$ compared to $\mathrm{\Delta }{\mathit{\sigma }}^{\mathrm{2}}$ derived with the PSD method, within the uncertainty the methods agree for this data set. This suggests that the PSD spectrum of the initial record (before diffusion took place) did not have a significant trend. However, this may not be the case in general as deposition regimes vary in time and space.

Conclusions

In this paper two methods for retrieving the differential diffusion length of the stable water isotope signal in ice cores were investigated. We have shown that temperature reconstruction using the PSD method on a single isotope can lead to significant biases because the PSD of the initial signal in precipitation may influence the obtained diffusion length. Therefore, unless whiteness of the initial spectrum can be proven, the uncertainty in the reconstructed temperature is larger for the single-isotope method than for the double-isotope method. Furthermore, before the PSD method is applied to ice core data, it needs to be calibrated by synthetic data. This will allow for an objective determination of the order of autoregression that is used to compute the PSD and the cut-off frequency. The quality of such a calibration depends on how well the synthetic data approximates the real data and thus on the deposition regime and the climatic conditions assumed. For the correlation method a calibration is not required as it is not dependent on the choice of parameters, making this a more robust method. However, also for this method the use of synthetic data is required to be able to constrain the uncertainty with which the differential diffusion length can be obtained. Doing such an analysis before ice core samples are cut is also useful in determining the best sampling strategy (length of the section and sampling resolution) to obtain the required accuracy and precision. Application to part of the EDML record from the Holocene period has shown that high data quality is required for diffusion thermometry. Identification of outliers in the isotope records is possible from the baseline of the PSD spectra. When the existing data for the EDML ice core (20 $\mathrm{m}$ long sections at 5 $\mathrm{cm}$ resolution) are used, the uncertainty (1$\mathit{\sigma }$ SD) in the reconstructed firn temperature for the Holocene period is about 2 ${}^{\circ }$C with the PSD method and 1.6 ${}^{\circ }$C with the correlation method. For samples from the glacial period this uncertainty increases by up to 0.6 ${}^{\circ }$C. As such, in addition to borehole thermometry and noble-gas measurements, diffusion thermometry is a useful independent tool for inferring differences in temperature between glacial and interglacial periods and thus for estimating the temporal gradient between water isotope concentration and local temperature.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Advanced Grant Agreement no. 226172 (MATRICs) awarded to H. Fischer.

This work is also a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Science Foundation–European Commission (EC) scientific programme, funded by the EC and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the UK. The main logistic support was provided by IPEV and PNRA (at Dome C) and AWI (at Dronning Maud Land). This is EPICA publication no. 300. Edited by: M. Schneebeli