Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core

. In order to reconstruct the temperature of the North Greenland Ice Core Project (NGRIP) site, new measurements of δ 15 N have been performed covering the time period from the beginning of the Holocene to Dansgaard– Oeschger (DO) event 8. Together with previously measured and mostly published δ 15 N data, we present for the ﬁrst time a NGRIP temperature reconstruction for the whole last glacial period from 10 to 120 kyr b2k (thousand years before 2000 AD) including every DO event based on δ 15 N isotope measurements combined with a ﬁrn densiﬁcation and heat diffusion model. The detected temperature rises at the onset of DO events range from 5 ◦ C (DO 25) up to 16.5 ◦ C (DO 11) with an uncertainty of ± 3 ◦ C. To bring measured and modelled data into agreement, we had to reduce the accumulation rate given by the NGRIP ss09sea06bm timescale in some periods by 30 to 35 %, especially during the last glacial maximum. A comparison between reconstructed temperature and δ 18 O ice data conﬁrms that the isotopic composition of the stadial was strongly inﬂuenced by seasonality. We evidence an anticorrelation between the variations of the δ 18 O ice sensitivity to temperature (referred to as α ) and obliquity in agreement with a simple Rayleigh distillation model. Finally, we suggest that α might be inﬂuenced by the Northern Hemisphere ice sheet volume.


Introduction
Deep ice core drilling in Greenland revealed large millennial scale variability in the water isotopic composition δ 18 O ice during the last glacial, reflecting local temperature variations (Johnsen et al., 1972;Dansgaard et al., 1982), known today as Dansgaard-Oeschger (DO) events. The 25 DO events identified in the North Greenland Ice Core Project (NGRIP) ice core (NGRIP members, 2004) are characterised by a rapid temperature increase followed by a gradual cooling back to stadial conditions. These rapid temperature variations are at least of northern hemispheric extent (Voelker, 2002) and can be clearly traced in different ice core climate proxies, for example, δ 18 O ice (Dansgaard et al., 1993), dust content (Ruth et al., 2003) and other aerosol contents (Mayewski et al., 1997), greenhouse gas concentrations (Huber et al., 2006b;Schilt et al., 2010), as well as in other climate proxies such as sea sediments Deplazes et al., 2013) and speleothems (Fleitmann et al., 2009;Kanner et al., 2012;Wang et al., 2001).
It is widely assumed that DO events are linked to reorganisations and/or variations in the strength of the Atlantic Meridional Overturning Circulation (AMOC) (Broecker et al., 1990). In addition, a southward shift of the Intertropical Convergence Zone (ITCZ) during Greenland cold (stadials and Heinrich) events is suggested to explain the similar DO pattern seen in Greenland ice cores and in tropical records (Chiang and Bitz, 2005). While the actual trigger and the detailed mechanisms involved in DO events remain unclear,  Capron et al. (2012) recent work underlines the role of sea ice cover in the Arctic and especially in the Nordic Seas (Li et al., 2010;Dokken et al., 2013;Petersen et al., 2013) where its rapid retreat and regrowth would be able to alter Greenland's temperature on decadal or even shorter timescale. Most of the DO events have an Antarctic analogue called Antarctic Isotopic Maximum (AIM) (EPICA community members, 2006;Blunier and Brook, 2001;Capron et al., 2010a;Wolff et al., 2010). The slow and quite smooth AIM temperature increases precede the rapid Greenland warming by several hundreds to thousands of years (Capron et al., 2010b). These findings are in line with the concept of the "thermal bipolar seesaw" where the Southern Ocean is considered as a heat reservoir which delivers heat via the AMOC to the North Atlantic and the Northern Sea (Stocker and Johnsen, 2003).
A first qualitative Greenland temperature reconstruction for the last glacial has been obtained using δ 18 O ice and the present-day spatial relationship between temperature and water isotopes of 0.67 ‰ • C −1 (Johnsen et al., 1989), leading to DO temperature amplitudes of ∼ 7 • C (Johnsen et al., 1992). However, temperature reconstructions based on the deconvolution of the borehole temperature (Johnsen et al., 1995;Cuffey et al., 1995) made clear that the present-day spatial gradient cannot be applied to the past. An alternative quantitative method is provided by nitrogen isotopes of the entrapped air. This proxy is sensitive to local firn processes only: firn thickness, temperature and temperature gradient as well as air entrapping processes (Schwander, 1989;Severinghaus et al., 1998;Severinghaus and Battle, 2006;Huber et al., 2006a). The obtained temperature variations of 3 to 16 • C within decades Huber et al., 2006b;Landais et al., 2005;Lang et al., 1999;Severinghaus et al., 1998) are in general larger than the ones obtained by climate model simulations (up to 7.5 • C, e.g. review by Kageyama et al., 2013). Kageyama et al. (2013) underline the varying response of the oceanic circulation and the sea ice cover to the fresh water hosing between models used to simulate a stadial period, explaining the diversity of modelled temperature difference between stadials and interstadials. A complete Greenland temperature record covering the last glacial may help to investigate this model-data discrepancy. Moreover, while a clear relationship between the Greenland stadial durations and the AIM warming events has been established (e.g. EPICA community members, 2006), a possible relationship between Antarctic and subsequent Greenland warming amplitudes has not been tested for the full glacial period yet.
The aim of this paper is first to quantitatively reconstruct the NGRIP temperature over the full last glacial based on δ 15 N data. To do so, we present the first NGRIP δ 15 N data covering the period from the beginning of the Holocene back to DO 8, completing the already existing δ 15 N records (Table 1). The combined δ 15 N records together with firn modelling are used to reconstruct the past NGRIP temperature and accumulation rate from 10 to 120 kyr. Then, we investigate the sensitivity of Greenland water isotopes to the reconstructed temperature changes at millennial and orbital timescales and its possible relation to the ice sheet size.

Paleothermometery method based on air δ 15 N measurements
To reconstruct the surface temperature evolution at the NGRIP site, we combine air δ 15 N measurements with simulations performed with a firn densification and heat diffusion model (Severinghaus et al., 1998;Leuenberger et al., 1999;Lang et al., 1999;Huber et al., 2006b;Schwander et al., 1997). The upper most 50 to 100 m of an ice sheet, where the air in the still open pores can exchange with the atmosphere, are called firn. In this layer the snow is gradually transformed to ice and at its bottom, the air is enclosed into bubbles. The processes relevant for the results presented here occur in the firn. With respect to the air molecule mobility, one can divide the firn into three parts: convective, diffusive and non-diffusive zones. In the generally small upper convective zone the composition of the air remains nearly unchanged compared to ambient air, whereas a gradual enrichment occurs in heavy molecules within the diffusive zone due to gravitation (i.e. 15 N/ 14 N ratio is increased). The enrichment stops at the lower end of the diffusive zone, at the lock in depth (LID). With the barometric formula (Craig et al., 1988;Schwander, 1989) one finds (in delta-notation): where z is the depth of the diffusive zone, g the acceleration constant, m the mass difference between the isotopes, R the ideal gas constant and T the mean firn temperature. Present-day δ 15 N measurements in firns show that no further gravitational enrichment occurs below the LID (Landais et al., 2006;Buizert et al., 2012). A temperature gradient in the firn causes a second important process called thermal diffusion, which leads, for example during a surface warming, to an additional enrichment of the 15 N/ 14 N ratio towards the bottom of the firn. Because gas diffuses about ten times faster than heat through the firn (Paterson, 1994), the thermal diffusion signal of a rapid temperature increase at the beginning of DO events can fully develop before the warmer temperatures propagate through the firn column and therefore slowly diminish the thermal signal. Thermal diffusion can be characterised as (Severinghaus et al., 1998) where T t and T b stand for the temperature at the top and the bottom of the firn, α T for the thermal diffusion constant (α T = 4.61198 × 10 −3 · ln(T /113.65 K), T stands for the average firn temperature ), for the thermal diffusion sensitivity and T for the temperature difference (T t − T b ). After several hundreds of years, when the temperature distribution in the firn is uniform again (neglecting geothermal influences), the thermal diffusion effect disappears but the signal of rapid warming is preserved by the 15 N/ 14 N ratio in air enclosed in bubbles.
As the air is able to diffuse down into the firn, a climatic signal will be registered in the gas phase deeper down in the ice core than in the ice phase. The resulting difference in depth is called depth and the age difference of the ice and the gas at a certain depth called age.

New and published δ 15 N measurements
The δ 15 N data set used for this work is a composite of measurements of two laboratories, Laboratoire des Sciences du Climat et de l'Environnement (LSCE), Gif-sur-Yvette and the Climate and Environmental Physics Division (KUP) of the Physics Institute of the University of Bern (Table 1). The data from KUP (Holocene to DO 17) has been measured with an online set-up with an uncertainty of ±0.02 ‰ in δ 15 N Huber and Leuenberger, 2004). For the new data set, 384 ice samples (up to 45 cm long) have been measured covering the period from the beginning of the Holocene to DO 8. Eleven samples have been rejected due to a failure of the system. The remaining 373 measurements have been divided into 600 data points, each representing 10 to 25 cm of ice depending on the original sample length which is influenced by the ice availability and ice quality. The ice samples investigated at LSCE (DO 18 to 25) have been measured with a melt-refreeze technique with a precision of ±0.006 ‰ (Landais et al., 2003(Landais et al., , 2004c. Additional LSCE δ 15 N data between DO 21 and DO 23 is published here for the first time (Table 1).

Strategy for surface temperature reconstruction
The temperature reconstruction was done with the firn densification and heat diffusion model from Schwander et al. (1997) by using the NGRIP ss09sea06bm age scale (in years before 2000 AD or yr b2k) (NGRIP members, 2004;Johnsen et al., 2001;Andersen et al., 2006). The ss09sea06bm timescale is the most appropriate since it is the only age scale with accumulation rate reconstructions over the entire studied time period. Transferring our final model input to the GICC05 age scale leads only to minor changes within the uncertainties in the model output compared to the used ss09sea06bm age scale, so we are confident that our results are valid for both age scales.
The input parameters of the model are age, accumulation rate and temperature. The temperature is calculated according to which expresses the qualitative similarity between the temperature and δ 18 O ice . 35.1 ‰ and 241.6 K stand for NGRIP Holocene values (NGRIP members, 2004), α fit and β fit denote scaling parameters to take into account the changing sensitivity between temperature and δ 18 O ice as outlined in detail below. The used accumulation and δ 18 O ice data from the NGRIP ss09sea06bm timescale has been splined with a cut off period (COP) of 200 years in order to match the observed variability in δ 15 N with the model. The basic idea is to vary the temperature (α fit and β fit ) and the accumulation rate in such a way that the model is able to reproduce the measured δ 15 N data. Note that the use of Eq.
(3) to find this best fit is not compelling but it is technically meaningful because it expresses a first order relation between temperature and isotopes. The fitted values of α fit and β fit are not unique and α fit is only in rough agreement with the α sensitivity discussed later. Three steps are followed to infer the NGRIP surface temperature: (i) rough adjustment of the temperatureδ 18 O ice sensitivity (α fit ) to approximate the shape of the modelled δ 15 N, (ii) refinement of the adjustment by varying the accumulation rate and a temperature shift (β fit ) for best fitting 890 P. Kindler et al.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core the age and (iii) final manual tuning where needed. The effect of each step on the modelled δ 15 N and age can be seen for the time period DO 4 to 7 in Fig. 1.
In the first step, 19 different α fit scenarios have been calculated by varying the α fit value in 0.02 steps from 0.24 to 0.60 ‰ • C −1 . As in Huber et al. (2006b) the initial accumulation rate was reduced by 20 % and a constant temperature shift of β fit = +4 K was introduced to match approximately the measured age. Without these corrections the modelled age is underestimated by up to 400 years in some parts, especially in the time period from 60 to 10 kyr (Fig. 7). In order to find the best α fit in a 2 kyr time window, the sum of the squared differences between the δ 15 N model output of a scenario and the spline through the measured data was calculated and the scenario with the smallest sum was determined. Thereafter the time window was shifted continuously by 250 years and the same procedure was applied until all data points were covered. At the end of step one, we obtain a best α fit value every 250 years, these values are splined for further calculations with a COP of 2 kyr. With these values, which vary roughly between 0.28 and 0.42 ‰ • C −1 , the model is able to coarsely reproduce the measured δ 15 N data ( Fig. 1, blue lines).
To further refine the adjustment in a second step, the accumulation rate is varied from 60 to 100 % in 5 % steps and the temperature offset β fit from +0 to +8 K in 1 K steps while retaining the α fit values from step one. This yields 81 different combinations and we proceed in the same way as in step one to find the best accumulation rate and temperature shift to match the measured δ 15 N data. The time dependent temperature shift β fit varies within 2 K around +4 K. With these three tuned parameters (α fit , β fit and accumulation) the model is able to reproduce the measured δ 15 N data generally well, both the amplitudes as well as the timing of the DO events. The age output of the model and the calculated depth are also in agreement with the measured data ( Fig. 1, green lines).
The parts of the reconstruction which do not yet match are now tuned manually (using parameters β fit and accumulation) on short periods (centuries) in a third step so that the model reproduces the measured data within their uncertainty ( Fig. 1, red lines). Manual adjustment was necessary especially in the Bølling Allerød-Younger Dryas transition and in the periods of DO 16 and 17, DO 19 and 20 and DO 23 to 25. The tuning was done by a linear adjustment of the temperature evolution between two selected time points in a given time period (e.g. a DO event). In addition, the accumulation rate was enhanced or lowered slightly in parts where the timing was not yet satisfactory. After step three, modelled δ 15 N as well as age and depth are in agreement with the measured data.
The 2σ uncertainty of the temperature amplitudes of DO events amounts to ±3 • C (Huber et al., 2006b). Deviations (e.g. amplitude and timing) of our temperature reconstruction compared to other studies using the δ 15 N method can be explained by the fact that other investigators have used different models (e.g. Goujon model, Goujon et al., 2003) and data-adjustment methods or by using a different approach (e.g. the combined use of δ 15 N and δ 40 Ar (Severinghaus et al., 1998;Landais et al., 2004b;Kobashi et al., 2011)).

Damping of the δ 15 N signal in the firn
High-resolution (20 mm) NGRIP δ 18 O ice data has shown that stadial-interstadial transitions can occur within a few years, e.g. the transition at DO 8 occurred within 21 years (Thomas et al., 2009;Steffensen et al., 2008). We cannot observe such a rapid increase at the onset of DO 8 in the NGRIP δ 15 N record (1.1 m data resolution) nor in the NEEM record (0.55 m data resolution, Fig. 2a). Also the duration of the temperature increase at DO 8 (roughly 200 years) of our reconstructed temperature is in contrast to the duration suggested by the high-resolution δ 18 O ice data. In addition to the lower δ 15 N resolution, two reasons could explain the discrepancy of the longer duration observed in the measured δ 15 N data as well as in the reconstructed temperature: (i) during extremely rapid events (decades), the δ 15 N (especially the δ 15 N therm ) formed at the LID is damped by gas diffusion in the firn and mostly during the bubble closeoff process that occurs over a depth interval of more than 10 m (Spahni et al., 2003), resulting in a smoother δ 15 N signal in the ice. The attenuation effect on the gas signal during the enclosure process at the bottom of the firn is included neither in the Schwander nor in the Goujon model (Schwander et al., 1997;Goujon et al., 2003). To our knowledge, only Grachev and Severinghaus (2005) have studied such an effect on δ 15 N variations and thus on surface temperature reconstruction. (ii) As it is explained in Sect. 2.3, the δ 18 O ice data, which is used in the model input for a first temperature estimate, has a resolution of 55 cm (corresponding to 10 to 100 years) and is splined with a COP of 200 years in order to reduce the signal variability in the firn densification and heat diffusion model to a level, which is in agreement with the measured δ 15 N profile.
To investigate possibility (i), we model a highly simplified scenario to mimic DO 8, where the temperature is increased by +10 • C (from −46 to −36 • C) with a concomitant accumulation increase from 0.06 to 0.10 m ice equivalent (m i.e.) within 20 years, according to the timing proposed by Thomas et al. (2009) for DO 8. Hereby we first calculate with the Schwander model the corresponding δ 15 N signal (without the proposed damping, dark blue dotted line in Fig. 2a) which is then used as an input for the firn model of Spahni et al. (2003). The Spahni model calculates the smoothing of an initial surface signal at the LID with respect to gas diffusion and the additional gradual bubble enclosure. The dashed (effect of diffusion) and the solid line (diffusion and gradual gas enclosure) scenarios in Fig. 2a are calculated by folding the original δ 15 N signal (dotted line) with the characteristic age distribution of the initial cold state. This approach is physically not fully correct. On the one hand, the used δ 15 N signal input for the Spahni model is not a surface signal but is generated in the firn itself. On the other hand, using a constant age distribution is not correct either: after a certain time period the surface warming penetrates to depths where the bubbles are enclosed (Fig. 2b, dashed line), altering the snow structure and modifying the age distribution. To allow for a small temperature rise at the LID, we use an age distribution corresponding to −45 • C, slightly higher than the initial temperature of −46 • C.
In Fig. 2a, the results of our simplified model calculations are compared to the measured δ 15 N values of DO 8 from NGRIP (green diamonds) and NEEM (violet diamonds), where a temperature increase of 8.8 ± 1.2 • C has been estimated . As the shape of the measured data and the modelled dark blue solid line agrees well, our simplified calculations indeed suggest for the considered scenario (+10 • C in 20 years) a significant damping of the signal amplitude of roughly −30 % due to gradual bubble enclosure compared to diffusion only. We refrain from calculating the damping of longer temperature increases because the LID-temperature would be affected significantly by the surface temperature increase, complicating the estimation of a unique and appropriate age distribution. In general, the damping decreases with increasing duration of a temperature rise as well as with a smaller temperature amplitude and vice versa (Spahni et al., 2003).
Following possibility (ii), we spline the initial temperature scenario (+10 • C in 20 years) with 200 years and use this data to calculate a δ 15 N signal with the Schwander model (Fig. 2a, grey line). Figure 2a shows that the effect of this initial smoothing of the input data is of similar magnitude regarding the signal amplitude as the effect arising from the signal damping calculated by the Spahni model.
While the duration of temperature increases may be overestimated in our reconstruction, we are confident that our calculated temperature amplitudes at DO events using www.clim-past.net/10/887/2014/ 892 P. Kindler et al.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core approach (ii) should be valid within the given uncertainty. Remaining local disagreements between modelled and measured data can be cleared by manual tuning as described in Sect. 2.3. These arguments also provide confidence in other studies which used δ 18 O ice data smoothed up to 70 years or more to estimate the corresponding site temperature (Capron et al., 2010a;Guillevic et al., 2013;Huber et al., 2006b;Landais et al., 2005).
For future work it is advisable to investigate the δ 15 N damping in the firn in more detail with transient models and a highly resolved measurement of an exemplary DO event (e.g. DO 1,5,8 or 19). Additionally, it would make sense to implement a gradual bubble enclosure process in firn densification and heat diffusion models which are currently used to reconstruct the temperature evolution of a specific site. Also the calculations for δ 15 N therm should be further refined in these models. Currently, an instantaneous signal during a surface warming is assumed at the bottom of the firn and the temperature gradient, which is important for the determination of the δ 15 N therm , depends only on two temperature points (top and bottom). As these two assumptions do not take into account the time lag of the signal due to diffusion nor the actual temperature profile in the firn, the model may slightly overestimate the δ 15 N therm at the beginning of a temperature increase.

Results and discussion
The results are divided into 3 sections. First we will discuss the new temperature reconstruction (Sect. 3.1), followed by a comparison between the NGRIP δ 18 O ice -temperature relationship on millennial, glacial and orbital timescales (Sect. 3.2). Finally, we will discuss the inferred accumulation rate (Sect. 3.3).

NGRIP surface temperature reconstruction
The reconstructed NGRIP temperature evolution from 10 to 120 kyr b2k, together with the used δ 18 O ice data and the measured and modelled δ 15 N data, is shown in Fig. 3. The temperature evolution for the period Younger Dryas-Holocene (DO 0) to DO 8 is presented here the first time whereas the temperature evolution for DO 9 to DO 25 is a reanalysis of existing data.
To define the temperature amplitude of a DO event, we specify the onset and the end of the event based on the following criteria. The DO event onset corresponds to the difference quotient exceeding 0.25 • C/50 years (0.05 • C decade −1 ), which is equivalent to about one-tenth of the increase rates for DO 9 to 17 found by Huber et al. (2006b). The end of the event was determined likewise (i.e. when the temperature increase rate became smaller than 0.05 • C decade −1 ). The 50 years time interval was chosen because it is long enough to overcome small-scale temperature variations. When one would use a longer time interval of 100 years or longer, the durations of the temperature increases tend to be prolonged in an unrealistic way and the corresponding temperature amplitudes are in general slightly reduced.
Generally, the DO events feature temperature amplitudes between 9 and 13 • C which are independent of the background of the climate state (different marine isotope stages (MIS)). The smallest temperature amplitude of 5 • C is obtained for DO 25. This is 2 • more than the result from Capron et al. (2012) but within their uncertainty of ±2.5 • C. On the other hand the largest amplitude of our temperature reconstruction is DO 11 with 16.5 • C, where Huber et al. (2006b) found a slightly smaller rise of 15 • C but still within the given uncertainty.
DO events 2, 7 (Fig. 3c), 11 and 18 feature a temperature increase in two steps which is different from the three DO types found by Capron et al. (2010a). The duration of the small drop between the two increases in temperature generally lasts for less than 100 years with amplitudes of less than 1 • C. As both sections of such a split temperature increase exhibit about the same temperature increase rate they were treated as one increase. Also DO 9 features basically a stop (Fig. 3d). As the temperature plateau in this reconstruction is of longer duration (about 300 years) than the following temperature increase (140 years), the temperature plateau was not taken into account which leads to T = 6.5 • C. This is lower in amplitude than the findings from Huber et al. (2006b) ( T = 9 • C), which would be better in line when one includes both temperature rises resulting in T = 8 • C. DO 5 (Fig. 3b) may also be considered to have a small temperature reversal at the beginning. When the 2 • C temperature drop just before the onset of the DO event is included, we obtain T = 14.5 • C. As the second (main) temperature increase starts more gradually, which is not the case for the DO events which feature a two-step temperature increase mentioned above, we exclude the small temperature drop before DO 5 from our temperature increase calculations.
Another type of temperature increase can be observed for DO 8, 12, 22 and maybe 17 where a slight long-term warming occurs before the start of the rapid temperature increase. They are roughly 3.5 • C in 950 years (DO 8), 3 • C in 1300 years (DO 12), 1 • C in 3500 years (DO 17) and 1 • C in 2800 years (DO 22). The long-term warming before DO 8 is comparable to the findings of Huber et al. (2006b) whereas they find a more pronounced long-term warming of about 3 to 4 • C for DO 17. Looking at MIS 3, each stadial where a Heinrich (H) event took place (Heinrich stadial) exhibits a long-term warming (Fig. 3a, orange shaded stadials). The yellow shaded stadials in Fig. 3a, situated around the last glacial maximum, experience H events but do not manifest any obvious long-term warming before a DO event. The warming at the NGRIP site during H stadial 4 (GS-9) could be of larger spatial extent as Guillevic et al. (2013) found also a slight long-term warming at the NEEM site.  Based on pollen records, Genty et al. (2010) and Sánchez Goñi et al. (2000) demonstrated that Europe featured a colder climate during periods of H events. In contrast, our NGRIP temperature reconstruction does not show particularly cold temperatures during these stadials.
The finding of the linear relationship between the duration of the Greenland stadial and the corresponding AIM temperature amplitude is explained by the bipolar seesaw mechanism (EPICA community members, 2006;Capron et al., 2010a;Vallelonga et al., 2012;Stocker and Johnsen, 2003). We see no correlation between the NGRIP DO temperature amplitudes and the corresponding Greenland stadial durations (Fig. 4a) nor the AIM temperature amplitudes (Fig. 4b). Therefore, the different temperature amplitudes at the onset of DO events may not be primarily influenced by the seesaw mechanism nor by different AMOC strengths but rather by regional amplification mechanisms such as sea ice changes (Gildor and Tziperman, 2003;Li et al., 2010;Petersen et al., 2013;Dokken et al., 2013). This is in line with a comparison of model simulations of different complexities (Kageyama et al., 2013), where the largest Arctic warming is produced by the AMOC restart associated with the largest sea ice retreat.

δ 18 O ice -temperature relationship
In this section we will investigate the relationship between δ 18 O ice and temperature, referred to as α, on different timescales. In general δ 18 O ice follows temperature (Fig. 3a) but the sensitivity is changing with time (Landais et al., 2004a). While a rather similar evolution of variations between temperature and δ 18 O ice can be observed during DO 20 to 14, a significant deviation occurs, for example, during DO 8. The interstadials of DO 21 (80 to 85 kyr) and 23 (98 to 105 kyr) feature even a slight decoupling of temperature and δ 18 O ice ; the temperature seems to be constant while the δ 18 O ice has a decreasing trend. The modelling of the temperature in these periods was more delicate due to rare age match points (Fig. 7). If the reconstructed temperature evolution for DO 21 and 23 is correct, the observed decoupling between temperature and δ 18 O ice may indicate a changing atmospheric circulation during longer interstadials. On the other hand, we observe also a divergence during the stadial preceding DO 24: temperature has a decreasing trend, δ 18 O ice an increasing trend. As the measured age could not be matched well by the model, caution should be taken in interpreting this period.
In the following two sections, the δ 18 O ice -temperature relationship will be discussed on glacial and orbital timescales. Note, that α is defined here as the actual slope between the reconstructed temperature and δ 18 O ice within a certain time period, as outlined below, and not in relation to the Holocene condition as expressed in Eq. (3). To calculate α, we use the 200 years splined δ 18 O ice data which served as an initial temperature estimation (Sect. 2.3).

δ 18 O ice -temperature relationship of the last glacial
Before comparing the δ 18 O ice data with our reconstructed temperature, the δ 18 O ice data has been corrected for the influence of the increased isotopic composition of the ocean during the glacial stage after Jouzel et al. (2003): δ 18 O sea denotes the isotopic difference of the sea water which is due to the increased ice masses at the poles compared to present conditions (Bintanja et al., 2005;Lisiecki and Raymo, 2005). Figure 5 shows that the δ 18 O icecorr -temperature relationship during the last glacial period (15 to 119 kyr) has a mean value of 0.52 ‰ • C −1 (black line). As there are uncertainties in both variables, the temperature and the δ 18 O icecorr , we used the geometric mean regression to determine α. The grey line indicates the modern Greenland spatial dependency of α = 0.80 ‰ • C −1 , the green and the red diamond represent glacial and modern mean values, respectively.
How can the observed data be explained? As described in Huber et al. (2006b) and as we suggest for interstadials, a general precipitation source temperature cooling during the glacial would shift the grey line into the grey dashed line (shown by the blue arrow), which means that less depleted precipitation is expected. This can be understood by means of a simple Rayleigh model: as the source temperature gets colder, the temperature gradient between source and site becomes smaller and therefore the isotopic composition is less depleted. Huber et al. (2006b) explained the smaller α (during MIS 3) compared to the present-day value by changes in annual distribution of precipitation: less winter precipitation in the stadial, more winter precipitation in an interstadial, respectively (Fig. 5, black solid arrows). The fact that the δ 18 O ice is not only dependent on the mean annual temperature but also amongst others on the seasonality of the precipitation has been discussed in other publications (Steig et al., 1994; Fawcett et al., 1997;Krinner et al., 1997;Masson-Delmotte et al., 2005;Jouzel et al., 2007). For the NGRIP site in particular, at present-day large spatial gradients of precipitation seasonality are observed (Steen-Larsen et al., 2011). Small changes in circulation patterns could therefore be responsible for significant seasonality changes in the vicinity of the NGRIP site.
Compared to Huber et al. (2006b), we suggest an additional explanation for the values of the isotopic precipitation during stadial conditions. It has been proposed (Gildor and Tziperman, 2003;Kaspi et al., 2004) and shown that sea ice changes around Greenland can have a significant influence on Greenland climate (Li et al., 2005(Li et al., , 2010. As the sea ice extent is enlarged during stadial conditions (Broecker, 2000), the precipitation source location for Greenland precipitation is pushed southwards to warmer ocean conditions which can be traced in the d-excess . In addition to the expected low stadial δ 18 O ice values (grey dashed line) the isotopy should therefore be further depleted because of the increased temperature gradient between precipitation source and the NGRIP site (Rayleigh). This suggested behaviour is shown in Fig. 5 by the qualitative red line. The discrepancy between observed and expected δ 18 O ice values is even enlarged compared to the interpretation from Huber et al. (2006b), suggesting that the seasonality effect (black dashed arrow) which is used to explain the mismatch between measured data and the enhanced isotopic sensitivity (red line) is even more pronounced during stadials.

δ 18 O ice -temperature relationship on orbital time scales
Here we investigate the δ 18 O ice -temperature relationship on orbital timescales. For this, we calculated α on a 10 kyr time window, every 2 kyr (Fig. 6a), together with the corresponding correlation coefficient (R 2 ). Blue dots correspond to time windows with a robust correlation (0.8 < R 2 ) and black dots to those with a weaker correlation (0.6 < R 2 < 0.8), occurring generally in time periods with rare and small abrupt climatic changes. The calculated α variations seem to follow obliquity (Fig. 6a, green line). When looking at single DO events (1 kyr time window ending at the DO peak, small grey diamonds in Fig. 6a), the variations are more scattered and the α-obliquity relationship is almost impossible to detect. The imprint of obliquity in the source-site temperature gradient has been previously evidenced in the d-excess measurements performed in both Antarctica and Greenland ice cores (Vimeux et al., 1999;Masson-Delmotte et al., 2005). Indeed, a low obliquity implies an increase of the source temperature and a cooling of the northern latitudes. The NGRIP δ 18 O ice record and our reconstructed temperature profile (10 kyr spline, Fig. 6a (NGRIP members, 2004;Bintanja et al., 2005) splined with a COP of 200 years (upper black line); 10 kyr temperature spline (red line); obliquity (green line, Berger and Loutre, 1991); δ 18 O sea as a proxy of ice sheet volume (lower black line, Bintanja et al., 2005); α values calculated for single DO events in a 1 kyr window (grey diamonds); α calculated over a time period of 10 kyr with a correlation of 0.6 < R 2 < 0.8 and 0.8 < R 2 , respectively (black and blue dots). α was calculated using the temperature record as shown in Fig. 3  probably further reduced by positive feedbacks such as extending sea ice and albedo effects resulting in a further increase of the source-site temperature gradient. Our α reconstruction covering the last glacial period shows that also the sensitivity of NGRIP water isotopes is affected by obliquity (Fig. 6a); α minima of around 0.3 ‰ • C −1 correspond to obliquity maxima while α values reach up to 0.6 ‰ • C −1 during obliquity minima. This is in qualitative agreement with a schematic Rayleigh distillation model (Fig. 6b), where the δ 18 O of precipitation is mainly explained by the source-site temperature gradient. In the case of an obliquity decrease, we move on the Rayleigh precipitation curve from the right black dot to the left one, which exhibits a steeper gradient (Fig. 6b, blue line). This steeper gradient would produce an enhanced isotopic effect for the same temperature variation in Greenland, therefore increasing α which is in agreement with our observations. While α and obliquity are roughly anticorrelated from 112 to 70 kyr, the α maximum at 20 kyr appears to be concomitant with the maximum of the ice sheet volume (Fig. 6a, black line, Bintanja et al., 2005). Note that the time uncertainty related to ice sheet size data is 4 kyr (Lisiecki and Raymo, 2005). Therefore it seems that both obliquity and Northern Hemisphere ice sheet volume have an impact on the sensitivity of Greenland water isotopes. Masson-Delmotte et al. (2005) and Landais et al. (2004a) have already proposed that larger ice sheets (favoured by smaller obliquity) would strongly reduce the arrival of precipitation in winter and hence increase the seasonality effect, thus decreasing α. While this suggestion remains valid, our data suggest a counter-acting mechanism; in the period 70 to 115 kyr a concomitant increase of α and ice sheet size can be observed. Bromwich et al. (2004) suggested that the Laurentide ice sheet modified the air mass trajectories, possibly influencing Greenland by very depleted precipitation coming from Pacific sources. Again considering a simple Rayleigh model, these precipitations would (i) lower the δ 18 O level and (ii) be more likely to be sensitive to temperature fluctuations in Greenland. As a consequence, α would increase with the Northern Hemisphere ice sheet volume, which is what we observe here. Our first continuous record of α over a complete glacial-interglacial cycle therefore supports that both obliquity and probably Northern Hemisphere ice sheet volume influences the hydrological cycle and therefore the δ 18 O of precipitation.  Fig. 7. (a) Accumulation rate from the ss09sea06bm timescale (blue line; NGRIP members, 2004;Johnsen et al., 2001), reduced accumulation rate used for the temperature reconstruction (red line). (b) Modelled age with original accumulation rate (blue) and reduced accumulation rate (red), depth measurements (green dots).

Accumulation rate
As mentioned in Sect. 2.3 and illustrated in Fig. 7, we reduced the accumulation rate in some parts substantially to adjust the modelled δ 15 N as well as depth and age to match the measured data. From 12 to 64 kyr the accumulation was lowered by 20 to 30 % with a mean value of 24 % and from 64 to 120 kyr by 0 to 20 % with a mean value of 4 %. During some few periods the accumulation rate was even reduced by 30 to 35 %.
A comparison of our temperature reconstruction with reduced accumulation to a scenario with unchanged accumulation is shown in Fig. 7. The red lines show the used reduced accumulation rate and the corresponding modelled age values. The blue lines are obtained with an unchanged accumulation rate as described above in Sect. 2.3 but without a final manual adjustment which would lead only to minor changes regarding age values. One can clearly see that in the period from 10 to 60 kyr the 100 %-accumulation scenario underestimates the age up to 400 years. Most significant deviations occur during the last glacial maximum. Not shown is the modelled δ 15 N data for the 100 %-accumulation scenario, for which a substantial disagreement in the timing between measured and modelled data occurs in this period (modelled δ 15 N values are shifted to younger ages). According to our data adjustment, the accumulation reduction seems to be smaller in the older half of the record, nevertheless also here our model suggests in some periods (around 99 and 110 kyr) a reduction of 10 to 20 %.
A general test to verify the Schwander model (Schwander et al., 1997) regarding accumulation rates and corresponding LID (which influence the age and depth) can be performed with modern data from different ice core sites. With the present-day NGRIP accumulation rate of 0.19 m i.e. yr −1 and a mean surface temperature of −31.5 • C (NGRIP members, 2004), the Schwander model (Schwander et al., 1997) calculates a LID of 66.8 m, which is in excellent agreement to measurements of 66 to 68 m (Huber et al., 2006a).
Consequently, in absence of any warming event, the model is also able to calculate the δ 15 N grav . According to Figs. 3 and 7, the coldest temperatures of the reconstruction are around −55 • C and the lowest accumulation rates are in the order of 4 cm i.e. yr −1 . These values are similar for Dome C present-day surface climatic conditions: accumulation 2.7 cm i.e. yr −1 and a temperature of −54.5 • C (Landais et al., 2006). So we can consider the Dome C characteristics as a test for the Schwander model for NGRIP lower boundary stadial conditions. When one runs the Schwander model with Dome C present-day surface conditions it is able to reproduce the LID (94.6 m) and the δ 15 N grav (0.51 ‰) in close agreement with observations of LID = 98 m and δ 15 N grav = 0.52 ‰ (Landais et al., 2006). Therefore we are confident that the Schwander model is capable of modelling NGRIP glacial climatic conditions. The observed necessity for a partial reduction in the accumulation rate in order to reconstruct a δ 15 N-based temperature with firn densification and heat diffusion models was recognised also in previous studies and therefore supports our findings. To reconstruct the temperature evolution from DO 9 to 17 (38 to 64 kyr), Huber et al. (2006b) reduced the accumulation rate from the ss09sea06bm timescale constantly by 20 %. The average of our time dependent reduced accumulation rate over the same period is similar (22.5 %). To test the validity of their model, Huber et al. (2006b) also modelled δ 15 N values of DO 18 to 20 (64 to 79 kyr) measured by Landais et al. (2004a) and found that in this period one can use mostly the unchanged accumulation rate, which is also in line with our findings. Also Guillevic et al. (2013) who reconstructed the NGRIP temperature over DO 8,9 and 10 (36.5 to 43 kyr on the GICC05 timescale), using the same δ 15 N data but the Goujon firnification model (Goujon et al., 2003), had to reduce the accumulation rate by 26 % in this section. Again, when averaged we find for the same period a similar accumulation reduction of 29 %. Interestingly, Landais et al. (2004aLandais et al. ( , 2005 and Capron et al. (2010aCapron et al. ( , 2012 worked with the original accumulation rate to reconstruct the temperature evolution over DO 18 to 25. These findings are mostly confirmed by our work where we use accumulation rates of approximately 100 % around the DO events mentioned above. An exception is DO 24 for which we had to lower the accumulation rate by 10 to 20 %. As the ss09sea06bm timescale has been validated by the GICC05 annual layer counting timescale  back to 64 kyr, the accumulation reduction in this period cannot be due to an inaccuracy in the timescale. We therefore confirm that the Dansgaard-Johnsen ice flow model (Dansgaard and Johnsen, 1969;NGRIP members, 2004) which was used to establish the NGRIP accumulation rate probably provides generally too high glacial accumulation rates, especially during the coldest periods. This tendency is likely due to a too simplistic ice flow model .

Conclusions
We present for the first time a continuous temperature reconstruction for the whole glacial period (10 to 120 kyr) based on new and published δ 15 N measurements performed on the trapped air of the NGRIP ice core. In line with previous studies, we find surface temperature rises from +5 to +16.5 • C (±3 • C, 2σ ) at the onset of abrupt events (Huber et al., 2006b;Landais et al., 2005;Capron et al., 2010aCapron et al., , 2012. Most of the DO events show a one-step temperature increase. In contrast, DO 2, 7, 11 and 18 feature a two-step temperature rise with a small interruption of generally less than 100 years between stadial and interstadial conditions. No particular cold temperatures characterise stadials associated with H events and a long-term warming of about one to three degrees is observed during the Heinrich stadials 4, 5 and 6 of MIS 3.
With a highly simplified model calculation, we tried to assess the δ 15 N smoothing in the firn due to the gradual bubble enclosure and found a damping of roughly −30 % for a +10 • C temperature increase within 20 years. Our smoothing of the model input data has an effect which is in the same range. To be able to better quantify the damping in the firn in the future, it would be advantageous to implement this effect in the used firn densification and heat diffusion models.
When comparing the δ 18 O icecorr data with the corresponding temperatures during the last glacial and taking into account a warmer source temperature for the precipitation during stadials, we propose that the seasonality effect of precipitation during stadials is even more pronounced than previously assumed. Our first reconstruction of NGRIP α during the full glacial period shows (i) an anticorrelation between obliquity and a long-term (10 kyr) α, supported qualitatively by a simple Rayleigh distillation model, and (ii) a possible influence of the Northern Hemisphere ice sheet size on α which may be explained by enhanced Pacific vapour sources at the NGRIP site when the Laurentide ice sheet volume increases.
Associated with the NGRIP δ 18 O ice profile, our reconstructed temperature provide useful constraints for future investigations on the parameter α based on water isotope modelling aiming at better quantifying in particular, the respective influences of obliquity on the source-site temperature gradient, and the ice volume.
In order to match measured and modelled δ 15 N data as well as age, we had to reduce the accumulation rate given by the ss09sea06bm age scale significantly, particularly in the second half of the last glacial from 12 to 64 kyr where a mean reduction of 24 % was applied. This further emphasises the fact that the Dansgaard-Johnsen ice flow model partly overestimates the NGRIP accumulation rate (Dansgaard and Johnsen, 1969;NGRIP members, 2004;Huber et al., 2006b;Guillevic et al., 2013) Supplementary material related to this article is available online at http://www.clim-past.net/10/887/2014/ cp-10-887-2014-supplement.zip.