Past surface temperatures at the NorthGRIP drill site from the difference in firn diffusion of water isotopes

A new ice core paleothermometer is introduced based on the temperature dependent diffusion of the stable water isotopes in the firn. A new parameter called differential diffusion length is defined as the difference between the diffusion length of the two stable water isotopologues 2H1H16O andH2O. A model treatment of the diffusion process of the firn and the ice is presented along with a method of retrieving the diffusion signal from the ice core record of water isotopes using spectral methods. The model shows how the diffusion process is highly dependent on the inter-annual variations in the surface temperatures. It results in a diffusion length longer than if the firn was isothermal. The longer diffusion length can be explained by the strong nonlinearly behaviour of the saturation pressure over ice in the range of the surface temperature fluctuations. The method has been tested on δ18O and δD measurements, spanning the transition from the last glacial to the holocene, from the NorthGRIP ice core. The surface temperature reconstruction based on the differential diffusion resembles other temperature reconstructions for the NorthGRIP ice core. However, the Allerød warming is seen to be significantly warmer than observed in other ice core based temperature reconstructions. The mechanisms behind this behaviour are not fully understood. The method shows the need of an expansion of high resolution stable water isotope datasets from ice cores. However, the new ice core paleothermometer presented here will give valuable insight into past climate, through the physical process of isotope diffusion in the firn column of ice sheets. Correspondence to: S. B. Simonsen (sbs@nbi.ku.dk)


Introduction
The temperature dependency of stable water isotopes in precipitation has been observed since the 1950s (Dansgaard, 1953(Dansgaard, , 1964;;Dansgaard et al., 1969).This dependency has been used to reconstruct past climate from ice cores (e.g.Johnsen et al., 1989).However, the isotope-temperature relationship does not consistently hold through time.The biases using stable water isotopes is the dependency on the source area, air-parcel trajectory, isotopic ratio of oceanic water, seasonality and cloud-surface temperature (Jouzel et al., 1997).Therefore, other methods have been used to reconstruct the past temperature record from ice cores.Over long time scales, borehole thermometry has been used to constrain the stable water isotope temperature proxy (Johnsen et al., 1995) giving an estimate for the past surface temperature at the precipitation site.During rapid climate changes, the fractionation of gas isotopes 15 N 14 N and 40 Ar 36 Ar (Severinghaus et al., 1998;Jouzel, 1999;Severinghaus and Brook, 1999), can be used as a paleothermometer.
A new method of estimating paleotemperatures from ice core records is presented in this paper to shed light on the surface temperature history of the ice sheet at the drill sites, especially at times of rapid climate changes.We will test the method implied by Johnsen et al. (2000Johnsen et al. ( , 2006) ) with available data from the NorthGRIP ice core.Thereby establishing a new ice core paleothermometer based on the difference in the temperature dependent diffusion of two water isotopes δ 18 O and δD.Here, the δ-value is defined as where R is the ratio of the species 18 O 16 O or D H and R ocean is the ratio of oceanic water.The present oceanic ratio is Published by Copernicus Publications on behalf of the European Geosciences Union.1328 S. B. Simonsen et al.: Past surface temperatures at the NorthGRIP drill site defined by "Standard Mean Ocean Water" (SMOW), which was defined by Craig (1961) with respect to the National Bureau of Standards reference sample 1.The amplitude of the original isotopic signal is smoothed by diffusion in the firn column 1 , mainly through the vapour phase.This diffusion in the vapour phase is highly temperature dependent and is relatively well understood (Rempel and Wettlaufer, 2003).Therefore, the diffusion signal can be recovered from the ice to give information on the temperature changes affecting the measured isotopic signal in the ice cores.

Modelling isotopic diffusion
Two models have been setup to describe the diffusion processes of stable water isotopes in dry firn; Whillans and Grootes (1985) and Johnsen (1977).Both models give similar results for the diffusion of the water isotopes (Johnsen et al., 2000).Following the approach from Johnsen (1977) and Johnsen et al. (2000), the origin of the coordinate system's z-axis follows the snow layer as it moves down through the firn and ice matrix.Accounting for the vertical strain rate and assuming a uniform strain rate εz the diffusion equation becomes where δ = δ(z ,t) is the smoothed and compressed isotope profile at time, z = zexp t 0 ε(t )t is the vertical compression of the original profile and D i is the diffusivity of the isotopic species i.In order to solve the diffusion equation, Johnsen (1977) defined the diffusion length σ given by 1 2 The solution to the diffusion equation is then derived by the convolution ( * ) of δ(z,0) with a Gaussian having the stan- (Johnsen, 1977;Johnsen et al., 2000).In the firn column, the strain rate is given by where u is the horizontal velocity and ρ is modelled by the Herron-Langway densification model (Herron and Langway, 1980).Having established the frame of the diffusion process, the diffusion in both firn and ice will now be discussed separately.
1 Firn can be defined as snow which has survived one melt season without being transformed into ice (Paterson, 2002).Here we define the firn column as the transformation stage of snow into ice, where an interconnected network of pores makes the diffusion of stable water isotopes possible throughout the vapour phase.

Diffusion in firn
In the firn the diffusion process is caused by an exchange of water species in the vapour phase, which smoothes the seasonal amplitude of the isotopic ratio of δ 18 O and δD.The diffusion is restricted by the shape of the interconnected pores in the firn given by the tortuosity factor τ , resulting in an effective diffusivity D eff = D a i τ of stable water isotopes inside the pore space (Johnsen et al., 2000).The diffusivity of air D a i has been measured in the lab (Merlivat, 1978;Cappa et al., 2003;Luz et al., 2009) and we use the value found by Merlivat (1978) in the following.Based on Schwander et al. (1988) and Jean-Baptiste et al. (1998), Johnsen et al. (2000) showed that the tortuosity can be assumed to be a function of the density ρ where b τ = 1.30 and is the effective pore close off density.Assuming the firn grains are well mixed and in isotopic equilibrium with the vapour, the diffusivity of firn can be expressed using the effective diffusivity and the firn density where m is the molar weight of water, p is the saturation vapour pressure over ice, R is the gas constant, T K is the temperature (K) and α i is the fractionation factor for isotopologue i.The saturation vapour pressure has been found empirically to be (Johnsen et al., 2000).To model the diffusion process in the firn column, the temperature at which the diffusion takes place is an important parameter.The firn temperature is modelled using the general transfer equation (Paterson, 2002, pp.224), driven by the surface temperature.The surface temperature is modelled by a superposition of two cosines (Benson, 1962, p. 44), giving the observed narrow summer maxima and a broader winter minima.We parameterize the surface temperature fluctuations as where t is the time and T mean is the mean annual surface temperature at the site.A = 16.5, B = 3.0 and b = 0.5236 are constants controlling the temperature amplitude and the narrowness of the summer temperature peak.The parametrization in Eq. ( 9) has been validated by comparison to automated weather station data from the GISP2 site during the period from June 1989 to November 1996.Based on Eq. ( 9) the monthly temperatures in the firn have been calculated and the temperature history of a layer moving down through the firn can be derived, seen as the dotted line in Fig. 1a. Figure 1b shows how the calculated diffusion length of δ 18 O becomes longer when deriving the diffusion length with a polythermal firn instead of only using the mean annual temperature.The main effect of the polythermal firn is in the top meters, as seen in Fig. 1c.This effect is due to the exponential nature of the saturation vapour pressure over ice, giving a stronger signal from the summer warming.The modelled diffusion lengths for δ 18 O and δD in the firn show the value of incorporating changes in the surface temperature.By allowing the surface temperature to vary according to Eq. ( 9), we can obtain the most accurate picture of the diffusion in the firn by numerical modelling.

Diffusion in ice
Below the effective pore close-off the diffusion of the isotopic ratio is expected to be single crystal self-diffusion with a diffusivity of ice D ice = 3.96 • 10 4 exp 7273 T m 2 yr (Ramseier, 1967), which is significantly less than the firn diffusivity for stable water isotopes.Assuming the diffusion length at the pore close off σ τ , Eq. ( 3) gives the self diffusion length of the ice σ ice where the index sec refers to the section of ice core where the samples are taken.It has been speculated that other processes are involved in the diffusion of the isotopic signal other than the single crystal self-diffusion, as the report of excess diffusion in the GRIP ice core by Johnsen et al. (2000).However, similar excess diffusion has not been found in other ice cores in Greenland.Hence, the diffusion process of ice is assumed here to only be the single crystal self-diffusion.The segregation of firn and ice diffusion results in additive properties of the diffusion length.When the ice diffusion is equivalent to the strain corrected firn diffusion, the ice has lost the information of what happened in the firn with respect to the temperature dependent diffusion.Therefore, this can be seen as the lower boundary of the paleotemperature reconstruction based on diffusion studies.Based on the modelled diffusion length in Fig. 3, this boundary is found in the NorthGRIP ice core at about Greenland interstadial 8 (GI-8, at depth 2060 m and dated around 38 kyr before year 2000 (b2k)).

Differential diffusion length
The differential diffusion length ( σ ) is defined from the diffusion lengths of the δ 18 O and δD Using the equations derived above (Eqs.6-10), the differential diffusion length at the base of the firn (at pore close-off) can be modelled as a function of the accumulation and the surface temperature, see Fig. 2.This model will be used in the following to estimate the temperature from the differential diffusion length and the accumulation record.
In order to estimate the past temperature from Fig. 2, the measured diffusion length of an ice section σ sec has to be corrected for the thinning occurring from pore close-off until the final depth.Assuming only self-diffusion below the pore close-off the differential diffusion length correction simplifies to Here σ 2 τ is the differential diffusion length at the effective pore close-off.
Based on the isotope borehole temperature reconstruction for NorthGRIP (Johnsen et al., 1995, personal communication) seen in the lower panel of Fig. 3, the diffusion length for NorthGRIP ice core has been modelled without the ice diffusion in Fig. 3.The figure also shows the diffusion length of δ 17 O, which is also being measured in ice cores (Landais et al., 2008).Using δ 17 O, a second differential diffusion length of stable water isotopes can be defined by This differential diffusion of δ 17 O and deuterium is predicted to have a stronger signal in the isotopic record and, thereby, be a better candidate for differential diffusion studies.However, 17 O data have not yet become available for studies of this differential diffusion.Hence, in the following, the differential diffusion will refer to Eq. ( 11).

Spectral methods and estimations
Having modelled the diffusion process of stable water isotopes, the diffusion signal can be retrieved from the measured profile using spectral methods.Using the convolution theorem, Eq. ( 4) can be transformed to the frequency domain to give where is the Fourier transform and the k is the wave number k = 2π ω.Expressed as the power spectral density (PSD) Eq. ( 14) gives where P is the PSD of the isotopic signal measured in the ice core and P 0 is the power density spectrum of the original isotopic signal as precipitated when the snow was deposited.
The diffusion length of an isotopic data point can be determined from Eq. ( 15).The differential diffusion length PSD is then defined as Snow precipitated on the surface of the ice sheet is assumed to have white-noise characteristics with a peak corresponding to the annual climate cycle, later the wind drift shifts the spectra to blue and finally diffusion shifts the spectra to a red-noise spectra (Fisher et al., 1985).Equation ( 15) shows this red-noise property of the PSD of isotopic data caused by the diffusion.According to Eq. ( 15), the diffusion length can be estimated from the PSD of the data by a linear fit to the red part of the spectra.However, the precision of the fit is highly dependant on the method of PSD estimation.The Maximum Entropy Method (MEM) (Andersen, 1974) has proven to be a reliable method used for PSD estimation of short isotopic datasets (Simonsen, 2008), down to about 200 sample points, compared to other methods, such as FFT or the Wiener-Khinthine's theorem.When interested in the diffusion estimates, great spectral resolution is not needed, but the order of the MEM has to be chosen to minimize the error in the linear fit of the red part of the spectra.The autoregressive nature of the MEM leads to correlated points in the PSD.The uncertainty of the fit can then be estimated from the one standard deviation confidence interval and the degrees of freedom N DOF = N poles − 2, where N poles is the number of poles in the fit and 2 accounts for the two estimated parameters in the fit.
In addition to the shape of the spectra described by Eq. ( 15), a white-noise tail is often observed in the PSD of the stable water isotope record.Assuming uniform white noise, the white-noise tail can be used to estimate the noise levels of the measurements and may be subtracted from the estimated PSD.
The discrete sampled isotopic record, with a sampling size of , gives rise to an aliasing effect in the spectra.The aliasing effect can be estimated by evaluating the Fourier transform of the sampling theorem (Press et al., 2007) and Eq. ( 15) at the nyquist frequency f nq = 1 2 , σ 2 aliasing is the aliasing effect and should be subtracted from the estimated diffusion length for δ 18 O and δD.As long as both δ 18 O and δD have been sampled with the same resolution, the aliasing effect vanishes for the differential diffusion length.

Differential diffusion as surface temperature estimator at the last glacial transition
As shown above, the differential diffusion is strongly temperature dependent and can be retrieved from the isotopic record using spectral methods.From the definition of the differential diffusion, the single crystal self-diffusion of the ice matrix can be neglected.The signal retrieved from the ice cores can be directly linked to the surface temperature at the period of the precipitation through the understanding of the diffusion process in the firn column.The modelled differential diffusion in Fig. 2 is the key to the interpretation of the retrieved differential diffusion length.By correcting the measured differential diffusion length for strain effects and using the accumulation record, the past surface temperatures can be derived.Ice core data have to be corrected for the general flow of the ice in the proximity of the drill site.In the last part of this study, we test the performance of the differential diffusion method with measurements from the NorthGRIP ice core (NGRIP Members, 2004).Recently a section of the NorthGRIP ice core has been measured for both δ 18 O and δD with a resolution of 5 cm.These measurements provide the ideal dataset for testing the spectral method of reconstructing past climate from the differential diffusion.
A number of ice flow models have been developed for the NorthGRIP site, among them are Grinsted and Dahl-Jensen (2002) (In the following, this model will be referred to as the GDJ-model) and Johnsen et al. (2001).Such flow models can provide both the time scale and strain history in relation to the depth of the ice core.For the NorthGRIP site, annual layer counting has provided the most accurate time scale, the GICC05 Greenland Ice Core Chronology 2005 (GICC05) (Andersen et al., 2006;Svensson et al., 2006;Svensson et al., 2008).The GICC05 dates the NorthGRIP ice core throughout the last 60 kyr.In the following, a version of ss09sea (Johnsen et al., 2001) corrected to fit GICC05 will provide the strain history for the NorthGRIP ice core.Applying the strain history from ss09sea to the observed layer thickness in GICC05 provides us with the annual accumulation record needed for estimating past temperatures from the diffusion record of the NorthGRIP isotope series.
Since the available dataset covers a time period known for its abrupt climate changes, it would not be adequate to only derive one PSD from the full dataset.The MEM estimation of the PSD is known to be capable of handling short datasets.Therefore, a running mean estimation of the differential diffusion length is proposed, with a window size of N, for each of the windows the differential diffusion is found along with the uncertainty.The differential diffusion of the window is converted into a surface temperature using the mean accumulation from the GICC05.To ensure reproducibility and to suppress outliers the temperature reconstruction is done both for the running mean, running forward and backwards in time, and the mean of the two is calculated.
Besides the uncertainty in the linear fit, two other sources or errors in the reconstruction are known.First, the annual layer thickness is well constrained in the GICC05 time scale with an error at the onset of the GI-1 of ±93 yr (Rasmussen et al., 2006) and at the onset of GI-3 ±416 yr (Andersen et al., 2006;Svensson et al., 2008, Table 1).The error of the annual layer thickness is much smaller than the reported cumulative counting error of the GICC05 and may be neglected in the error estimate.Second, the error from the strain rate may invoke greater uncertainties in the reconstructed temperatures.The model structure makes it hard to assess the errors induced by the uncertainty in accumulation and strain rate.Therefore, in the following temperature reconstruction of the transition in the NorthGRIP ice core, the only error presented is the error estimate of the linear fit to the power spectra.
To evaluate the sensitivity of the temperature reconstruction based on the differential diffusion of the stable water isotopes, we have to look into the nature of the modelled diffusion length at the pore close-off, as seen in Fig. 2. The contour plot clearly shows how an uncertainty in the accumulation/stain rate is more pronounced in a cold climate conditions than warm climate conditions giving a constant accumulation for the two periods.
The temperature reconstruction of the Bølling/Allerødwarming (GI-1, see Björck et al., 1998) and the Younger Dryas (YD) are shown in Fig. 4. The temperature reconstruction has not been corrected for the white measuring noise, since the limited sample size combined with the depth of the ice core samples makes it impossible to pinpoint the whitenoise tail.The typical noise level for δ 18 O measurements are 0.07 0 / 00 and 0.5 0 / 00 for δD, the measuring noise would give rise to an error in the estimated differential diffusion length of 0.18 0 / 00 using typical values of the parameters in the differential diffusion estimation at the time of the reconstruction.The error due to measuring uncertainties is well within the error of the fitted differential diffusion length and may be neglected in the temperature reconstruction as long as the white-noise tail is not part of the fitted area.
To assess the bias of assumptions made to produce the reconstructed surface temperature shown in Fig. 4, the reconstruction has also been done using different approaches for modelling the dependency of the differential diffusions on temperature and accumulation.Despite not all being fully realistic, three approaches are also shown in Fig. 4. The three approaches are; (1) assuming flow and accumulation history derived by the GDJ-model, (2) another parameterization of the Herron-Langway densification model based on an isothermal version of the Nabarro-Herring type creep parameterization by Arthern et al. (2010), and (3) assuming the pore close at 830 kg per m 3 .These three approaches cannot be seen as a real confidence interval for the reconstructed surface temperature by differential diffusion, but show the effect of the model assumptions.
Finally, the temperature profile (in black) in Fig. 4 is believed to be the surface temperature during the period as implied by the differential diffusion of the two stable water isotopes.

Discussion
The borehole isotope temperature reconstruction and the differential diffusion temperature are in general agreement of the surface temperature in the cold Younger Dryas and the warm Bølling (GI-1e).However, the reconstructed Allerød warming does not match the borehole isotope reconstruction, especially at the time around 13.5 kyr b2k.Here the surface temperatures are just as warm as in the Bølling.This warming in the Allerød cannot be seen directly in the δ 18 O record nor in the excess.The only indication in favour of the estimated warming is the lowering in the excess just at the 13.5 kyr b2k.Another test to verify this warming seen by the differential diffusion is to look at the accumulation record itself, since a warming normally results in an increased layer thickness.However, when looking at the GICC05 (Svensson et al., 2006), no increases in the accumulation are recorded for the period.The lack of an increase in the accumulation of NorthGRIP record contradicts the differential diffusion temperature reconstruction.The only temperature record reassembling the behaviour seen in the differential diffusion surface temperature of the Allerød period is the sea surface temperatures (SST) found in marine sediment cores from the Iberian Margin at the coast of Portugal (Bard et al., 2000) and the Cariaco Basin on the coast of Venezuela (Lea et al., 2003) for the same time period.The reason for the Allerød to show up as a warm anomaly in the differential diffusion surface temperature is not understood.However, the dataset of high resolution measurements of δ 18 O and δD have to be expanded to see if the Allerød warming is just an anomaly in the temperature reconstruction by differential diffusion or if such deviations from the isotope temperatures are common in the record when using differential diffusion to reconstruct surface temperatures.
Three major shifts in the deuterium excess can be recognized in the record shown in Fig. 4, are also clearly seen as shifts in differential diffusion temperature reconstruction.The deuterium excess is believed to hold information about the precipitation source area (Johnsen et al., 1989;Masson-Delmotte et al., 2005), and the shifts have been defined as the transition from one climate period to another (Steffensen et al., 2008).The major shifts in the deuterium excess are seen in the low frequency part of the PSD and this frequency is excluded from the temperature reconstruction.Therefore, it is reinsuring that the temperature reconstruction sees these major shifts in deuterium excess as climate transitions.
Our method for reconstructing past surface temperatures from the difference in the diffusion for stable water isotopes is highly dependent on accurate knowledge of past accumulation and strain histories.Therefore, the method is applicable in locations, such as Greenland, where interannual changes in impurity content makes it possible to date the ice core with great precision from annual layer counting.In locations with lower accumulation, the success of the method becomes reliable on the accuracy of the model strain and accumulation history.However, at low accumulation sites, interesting climate transitions may be located at shallower depths and may not be as vulnerable to flow features as the deep transition in the NorthGRIP ice core.

Conclusions
Based on the definition of the differential diffusion and the modelling of the diffusion processes in the firn column it is shown how the differential diffusion length can be theoretically recovered down to about GI-8 in the NorthGRIP ice core.The methods have been applied to a dataset of highresolution measurements of the two stable water isotopes δ 18 O and δD from the glacial transition.This results in a new surface temperature reconstruction of the last glacial transition, which mostly agrees with previous temperature reconstructions made for this site.However, the Allerød warming is seen to be much warmer in the differential diffusion reconstruction than otherwise seen in the ice cores and it more resembles the SST history of the marine sediment core from the Cariaco basin.The mechanisms behind this behaviour are not yet known.
There is a need for expanding the archive of high resolution measurement of stable water isotopes in ice cores for improving the use of the differential diffusion paleothermometry.Recently, online water isotope measurements of ice cores with the use of melter systems and infrared spectroscopy have been developed (Gkinis et al., 2011).These methods have the potential to yield measurements of very high resolution that can be beneficial for the type of temperature reconstructions we present here.Therefore, the method presented here will give valuable insight into past climate, based on the physical process of diffusion and directly coupled to the surface temperature at the time of firn densification.The method presented here is only applicable for reconstructing the temperature signal before the ice self-diffusion becomes a dominant factor and cannot be corrected for.

Fig. 1 .
Fig. 1.Temperature profiles and diffusion lengths for the NorthGRIP firn column modelled under present day conditions.(a) The temperature field to a given month of the year.(b) The model diffusion length of δ 18 O both under constant temperature conditions and using the temperature variations seen in (a) from the surface to the pore close off.(c) The diffusion length of the top 4 m for selected months.

Fig. 2 .
Fig. 2. Modelled differential diffusion length at the pore close-off depending on accumulation and surface temperature.

Fig. 3 .
Fig. 3. Modelled diffusion lengths along the NorthGRIP ice core, based on NorthGRIP isotope borehole thermometry (shown in the lower panel).Also shown for reference is the δ 18 O record for NorthGRIP.As seen in the top panel, the modelled differential diffusion lengths are almost not seen in the interstadials in the Last Glacial.For the modelling of the diffusion length of δ 17 O the fractionation constants reported by Barkan and Luz (2007) are used.

Fig. 4 .
Fig. 4. Top panel: (Black curve) The temperature reconstruction based on the differential diffusion with a window of N = 200, which is moved N = 11.(Red shading) The confidence interval of the fitted PSD.(Green curve) Surface temperature reconstruction using flow parameters by DGJ.(Magenta curve) Surface temperature reconstruction using parameterization of the Herron-Langway firn densification model by Arthern et al. (2010).(Cyan curve) Surface temperature reconstruction using a close-off density of 830 kg per m 3 .For comparison the Isotope borehole thermometry are also show in blue.Middle panel: the deuterium excess record for the period.Major climatic shifts are marked with vertical green lines both on the top and middle panel.Lower panel: the δ 18 O record of the NorthGRIP ice core on the GICC05 timescale, for the last 60 kyr.All times are on the GICC05 time scale.