Spatial gradients of temperature, accumulation and δ 18 O-ice in Greenland over a series of Dansgaard–Oeschger events

. Air and water stable isotope measurements from four Greenland deep ice cores (GRIP, GISP2, NGRIP and NEEM) are investigated over a series of Dansgaard– Oeschger events (DO 8, 9 and 10), which are representative of glacial millennial scale variability. Combined with ﬁrn modeling, air isotope data allow us to quantify abrupt temperature increases for each drill site (1 σ = 0.6 ◦ C for NEEM, GRIP and GISP2, 1.5 ◦ C for NGRIP). Our data show that the magnitude of stadial–interstadial temperature increase is up to 2 ◦ C larger in central and North Greenland than in northwest Greenland: i.e., for DO 8, a magnitude of +8.8 ◦ C is inferred, which is signiﬁcantly the inferred at GISP2. spatial for accumulation is with


Introduction
The last glacial period is characterized by rapid climatic instabilities at the millennial timescale occurring in the Northern Hemisphere and recorded both in marine and terrestrial archives (Voelker, 2002;Bond et al., 1993).The NGRIP (North Greenland Ice Core Project) ice core, northern Greenland, offers a high resolution water isotopes record where 25 rapid events were identified and described with a precise timing (NGRIP members, 2004).These events consist of a cold phase or stadial (Greenland stadial, GS), followed by a sharp temperature increase of 9 to 16 • C at the NGRIP site as constrained by gas isotope measurements (Landais et al., 2004a(Landais et al., , 2005;;Huber et al., 2006).The warm phase is referred as Greenland interstadial (GI).Temperature then gradually cools down, sometimes with a small but abrupt cooling in the end, to the next stadial state.These temperature variations are associated with significant changes in accumulation rate, with annual layer thicknesses varying by a factor of two between GS and GI at NGRIP (Andersen et al., 2006;Svensson et al., 2008).
The identification of ice rafted debris horizons during GS in North Atlantic sediments (Heinrich, 1988;Bond et al., 1993;Elliot et al., 2001), together with proxy records pointing to changes in salinity (Elliot et al., 2001(Elliot et al., , 2002)),  (Rasmussen and Thomsen, 2004;Kissel et al., 2008) and Atlantic Meridional Overturning Circulation (AMOC) (McManus et al., 1994), had led to the theory that DO (Dansgaard-Oeschger) events are associated with large scale reorganizations in AMOC and interhemispheric heat transport (Blunier and Brook, 2001).The identification of a systematic Antarctic counterpart to each Greenland DO event (EPICA community members, 2006;Capron et al., 2010) supports this theory.This observation can be reproduced with a conceptual see-saw model using the Antarctic ocean as a heat reservoir and the AMOC as the way to exchange heat between Antarctica and Greenland (Stocker and Johnsen, 2003), as well as with climate models (e.g., Roche et al., 2010).Atmospheric teleconnections can also be at play (Chiang et al., 2008) and influence the timing of Antarctic warming with respect to Greenland cooling (Buiron et al., 2012).
Coupled atmosphere-ocean climate models are now able to reproduce the temperature pattern of DO events in Greenland in response to AMOC changes induced by freshwater forcing in the high latitudes of the Atlantic Ocean (Kageyama et al., 2010).However, modeled amplitudes of temperature changes are typically between 5 and 7 • C (Ganopolski and Rahmstorf, 2001;Li et al., 2005Li et al., , 2010;;Otto-Bliesner and Brady, 2010), significantly smaller than the temperature increase of 8-16 • C reconstructed based on ice core data (Landais et al., 2004a;Huber et al., 2006).The correct amplitude of temperature change over the Bølling-Allerød is only reproduced in a fully coupled and high resolution atmosphere-ocean global circulation model (Liu et al., 2009).However, a large part of the simulated warming is due to the simultaneous large changes in insolation (not at play for most DO events) and atmospheric CO 2 concentration.During DO events, limited CO 2 changes (20 ppm) are depicted by existing data (Indermühle et al., 2000;Bereiter et al., 2012;Ahn et al., 2012).This model-data mismatch motivates us to strengthen the description of the magnitude and spatial patterns of DO temperature changes, using different ice core sites.An improved regional description of past changes in Greenland climate will also be needed for the comparison with regional climate models recently equipped with water stable isotopes (Sjolte et al., 2011).So far, no systematic comparison of the signature of DO events in water stable isotopes and temperature (from gas isotopes) has been conducted over an array of drilling sites.This is the main target of this study.
In 2010, bedrock was reached at the North Eemian Ice Core Drilling (NEEM) site, northwest (NW) Greenland.A new deep ice core, 2.5 km long, is now available (Dahl-Jensen et al., 2013).In this paper, we present new data from the NEEM ice core, together with existing and new measurements conducted on the GISP2 (Greenland Ice Sheet Project 2) and NGRIP ice cores on DO events 8 to 10.The location of these drilling sites is depicted in Fig. 1 and their present-day characteristics are summarized in Table 1 (see also Johnsen et al., 2001).At present, the main source of NEEM precipitation is located in the North Atlantic, between 30 • N and 50 • N (Steen-Larsen et al., 2011).The recent interannual variability of water stable isotopes (δ 18 O, δD) shows similarities with the variability of the Baffin Bay sea-ice extent.Unlike central Greenland where snow falls year round, NW Greenland precipitation occur predominantly in summer according to snow pit studies (Shuman et al., 1995(Shuman et al., , 2001) ) and model simulations (Steen-Larsen et al., 2011;Sjolte et al., 2011;Persson et al., 2011).This specificity of the precipitation seasonality explains the particularly weak fingerprint of the North Atlantic Oscillation in NEEM shallow ice cores (Steen-Larsen et al., 2011) compared to GISP2 (Barlow et al., 1993).These regional peculiarities are of particular interest because past changes in precipitation seasonality are likely to affect water stable isotopes values.
Water isotopes are to a first order markers of local condensation temperature changes at the precipitation site (Dansgaard, 1964).However, they are also affected by evaporation conditions (temperature, relative humidity and wind regime, e.g., Merlivat and Jouzel, 1979;Johnsen et al., 1989), atmospheric transport and distillation, condensation conditions, as well as seasonality of precipitation (Werner et al., 2000;Werner et al., 2001;Masson-Delmotte et al., 2005): they are integrated tracers of the hydrological cycle and quantitative indicators of past site-temperature change, albeit with a timevarying relationship with local surface temperature.This temporal variability of the isotope-temperature relationship has been verified in Greenland thanks to independent constraints on past temperatures, either based on the inversion of borehole temperature data or derived from gas isotopes (e.g., Cuffey and Clow, 1997;Dahl-Jensen et al., 1998;Severinghaus and Brook, 1999;Lang et al., 1999;Johnsen et al., 2001;Landais et al., 2004a;Huber et al., 2006;Vinther et al., 2009).
Using the isotopic composition of nitrogen (δ 15 N) from N 2 trapped in the ice bubbles allows us to quantify the amplitude of past rapid-temperature changes (e.g., Severinghaus and Brook, 1999;Landais et al., 2004a;Huber et al., 2006;Grachev and Severinghaus, 2005;Kobashi et al., 2011).At the onset of a DO event, the firn surface warms rapidly but its base remains cold because of the slow diffusion of heat in snow and ice.The resulting temperature gradient in the firn leads to thermal fractionation of gases: the heavy nitrogen isotopes migrate towards the cold bottom of the firn, where air is progressively trapped into air bubbles.As a result, a sharp peak in δ 15 N is seen in the gas phase as a counterpart to the rapid increase in water stable isotopes in the ice phase.Using δ 15 N data and firn modeling, past surface temperature variations can be reconstructed (Schwander et al., 1997;Goujon et al., 2003).This method has already been applied to specific DO events on the NGRIP, GRIP and GISP2 ice cores (Lang et al., 1999;Huber et al., 2006;Landais et al., 2004aLandais et al., , 2005;;Goujon et al., 2003;Severinghaus and Brook, 1999;Capron et al., 2010) and will be applied here for the first time to the NEEM ice core.
For this first study of regional variability of temperature changes over DO events, we focus on the series of DO events 8, 9 and 10 during Marine Isotope Stage 3 (MIS3, 28-60 ka b2k, thousand years before 2000 AD).This period is indeed the most widely documented for the millennial scale variability in a variety of natural archives.It is characterized by a large terrestrial ice volume (Bintanja et al., 2005), low atmospheric greenhouse gas concentration (Schilt et al., 2010), decreasing obliquity, low eccentricity and therefore small fluctuations in Northern Hemisphere summer insolation (Laskar et al., 2004).During MIS3, iconic DO events are particularly frequent with short lived interstadials (Capron et al., 2010) and constitute a clear target for model-data comparisons.
In this study, we first present a new δ 15 N profile covering DO events 8 to 10 on the NEEM ice core.We then produce past temperature and accumulation reconstructions for NEEM and compare them with scenarios obtained with the same method for NGRIP, GISP2 and GRIP, investigating the water isotope-temperature relationships for these four different locations.We finally discuss the implications of our results in terms of regional climate variations., Grootes and Stuiver (1997).f Gundestrup et al. (1994).g Johnsen et al. (1992), 20-220 a b2k average of the GRIP ice core.

Nitrogen isotope data
The isotopic composition of nitrogen (δ 15 N) was measured on the NEEM core from 1746.8 to 1811.6 m depth at Laboratoire des Sciences du Climat et de l'Environnement (LSCE), France.We have a total of 84 data points with replicates, with an average depth resolution of 77 cm corresponding to an average temporal resolution of ∼ 70 a.
On the GRIP ice core, 18 data points covering DO 9 have been measured with an average resolution of 61 cm (equivalent to about 43 a).For these two data sets, we have used a melt-refreeze technique to extract the air from the ice (Sowers et al., 1989;Landais et al., 2004b).The collected air is then measured by dual inlet mass spectrometry (Delta V plus, Thermo Scientific).Data are corrected for mass interferences occurring in the mass spectrometer (Sowers et al., 1989;Bender et al., 1994b).Dry atmospheric air is used as a standard to express the results.The final pooled standard deviation over all duplicate samples is 0.007 ‰.
For the NGRIP core, δ 15 N was measured at the University of Bern (73 data points from Huber et al. (2006) and 36 new data points on DO 8).A continuous flow method was used for air extraction and mass spectrometry measurement (Huber and Leuenberger, 2004).The associated uncertainty is 0.02 ‰.
For the GISP2 ice core, nitrogen isotopes were measured at Scripps Institution of Oceanography, University of California, using the melt-refreeze technique from Sowers et al. (1989) with a pooled standard deviation of 0.0065 ‰ for the 70 δ 15 N data points (Orsi et al., 2013).In addition to these data, argon isotopes were also measured using the method from Severinghaus et al. (2003) (46 samples, pooled standard deviation of 0.013 ‰).We use the δ 18 O bag data (one data point corresponds to an average over 55 cm) from the NEEM ice core measured at the Centre for Ice and Climate (CIC), University of Copenhagen, with an analytical accuracy of 0.07 ‰.For the NGRIP and GRIP ice cores, we use the bag data previously measured at CIC, with the same precision (NGRIP members, 2004;Johnsen et al., 1992).The GISP2 δ 18 O data (Grootes and Stuiver, 1997, 20 cm resolution) are associated with a precision of 0.05 to 0.1 ‰.

Timescale
NEEM, NGRIP, GISP2 and GRIP ice cores are all dated according to the Greenland Ice Core Chronology 2005, GICC05 (Vinther et al., 2006;Rasmussen et al., 2006;Andersen et al., 2006;Svensson et al., 2008).This timescale has been produced based on annual layer counting of several parameters measured continuously on the NGRIP, GRIP and Dye3 ice cores and featuring a clear annual cycle, back to 60 ka b2k.The Maximum Counting Error (MCE), which can be regarded as a 2σ error estimate (Rasmussen et al., 2006) is 1439 a at 38 ka b2k.To transfer this timescale to the NEEM ice core, match points between peaks of electrical conductivity measurements and dielectrical properties, measured continuously on the ice cores, have been used.The obtained timescale for NEEM is called GICC05-NEEM-1 (S.O. Rasmussen, personal communication, 2010).The GISP2 core is matched to NGRIP using the same method as Rasmussen et al. (2008), with match points from I. Seierstad (personal communication, 2012).Using these match points, we scale the Meese et al. (1997) GISP2 timescale, based on annual layer counting, to the GICC05 timescale, in order to keep the information from the annual layer count while producing an age scale consistent with GICC05.The GICC05 age scale gives the age of the ice at each depth, and thus the annual layer thickness at each depth, but not the accumulation rate.This age scale is independent from estimation of thinning and past accumulation rate.

Accumulation rate
Accumulation rate histories for NEEM, NGRIP and GRIP are obtained using a Dansgaard-Johnsen (DJ) ice flow model (Dansgaard and Johnsen, 1969).The ice flow parameters of the model are tuned to obtain the best match between modeled and observed depth-age horizons in the ice cores.The thinning function calculated from the DJ model is then used to correct the observed annual layer thicknesses in the core for the effect of ice flow induced thinning, thereby producing an accumulation rate history.The NEEM version of the DJ model (Buchardt, 2009) is tuned in order to match the GICC05 timescale.For NGRIP, the accumulation rate was first calculated using the ss09sea06bm age scale (Johnsen et al., 2001;Grinsted and Dahl-Jensen, 2002;NGRIP members, 2004).We recalculated that accumulation according to the more accurate GICC05 timescale.We did the same for the GRIP ss09sea-accumulation rate from Johnsen et al. (2001).Note that the ss09sea, ss09sea06bm and the GICC05 timescales agree within the GICC05 uncertainty between 28 and 60 ka b2k.For GISP2, the accumulation rate was first estimated with a 1 m resolution based on the coupled heat and ice flow model from Cuffey and Clow (1997), with the layer counted timescale from Alley et al. (1993), Meese et al. (1994) and Bender et al. (1994b).This timescale has known issues in the vicinity of DO 8 (Orsi et al., 2013;Svensson et al., 2006), which causes the accumulation history derived from it to be also wrong.Orsi et al. (2013) used the layer thickness from the GICC05 timescale to recalculate the accumulation history.Cuffey and Clow (1997) suggested 3 accumulation scenarios and Orsi et al. (2013) use the "200 km margin retreat" scenario adapted to the GICC05 timescale, compatible with the firn thickness and age derived from δ 15 N data.This accumulation scenario has also been proved to best reproduce ice sheet thickness variations (Vinther et al., 2009).

Ice-gas depth data
Figure 2 presents the NEEM δ 15 N profile over the sequence DO 8-10.The peaks of δ 15 N at 1769.4,1787.5, and 1801.0 m are the result of the maximum temperature gradient in the firn corresponding to the abrupt temperature increases of DO 8, 9 and 10.We assume that δ 15 N peaks and δ 18 Oice peaks are synchronous (see Sect.A2) and thus relate the maximum firn temperature gradient to the peaks in δ 18 O-ice at 1758.1, 1776.8, and 1790.5 m.The depth differences between the temperature increases recorded in the gas and ice phases, named depth, can thus directly be inferred as 11.3,10.7,and 10.5 m over DO 8,9,and 10,respectively (Fig. 2,Table 2,points 1,5,and 7,respectively).We propose another match point between weaker peaks of δ 15 N and δ 18 O, see match point 3 in Table 2 and Fig. A3.
δ 15 N also increases with accumulation increase, which deepens the firn (see Sect. 2.2), and we believe that this effect explains the beginning of δ 15 N increase at the onset of each DO event.Several abrupt transitions (Bølling-Allerød and DO 8) have been investigated at high resolution (Steffensen et al., 2008;Thomas et al., 2008), also showing that the accumulation increases before the δ 18 O shifts with a time lead up to decades and ends after the completion of the δ 18 O increase.We observe the same feature for DO 8, 9, and 10 on the NEEM core.We thus match the onset of the δ 15 N increase at the beginning of DO events to the onset of accumulation increase, which occurs before the δ 18 O increase (Table 2 and Fig. A3, match points 2, 6, 8).Finally, match point 4 is a step in accumulation that we relate to the same step seen in δ 15 N variations.At NGRIP, DO 8, 9 and 10 are seen at 2086.6, 2116.2, and 2139.4 m in the gas phase and at 2068.0, 2099.1, and 2123.3 m in the δ 18 O from the ice phase (Fig. 2, Table 2 and Fig. D1, match points 2, 3, 4, respectively).For DO 8, δ 18 O shows a double peak and we use the middle depth for this match point.We propose another match point at the end of DO 8 between weaker peaks of δ 15 N and δ 18 O (match point 1).All these depth match points will be used in Sect.3.1, combined with firn modeling, to reconstruct past surface temperature and accumulation.

Model description
To reconstruct a surface temperature scenario from the δ 15 N profiles, we use a classical approach consisting of fitting the output of a firnification and heat diffusion model with the δ 15 N records (Schwander et al., 1997;Lang et al., 1999;Huber et al., 2006;Goujon et al., 2003;Landais et al., 2004a;Kobashi et al., 2011;Orsi et al., 2013).We use here the semi-empirical firnification model with heat diffusion by Goujon et al. (2003).Adapted to each ice core (see method in Appendix A), this model calculates for each ice age and hence for each corresponding depth level the initial firn depth (defined here as the depth where diffusion of gases stops i.e., lock-in-depth, LID), the age difference between ice and gas at the LID ( age), and the temperature gradient between the bottom and the top of the firn.It is then possible to calculate the δ 15 N as the sum of two effects: -gravitational effect (Craig et al., 1988;Schwander, 1989): the heavy isotopes preferentially migrate towards the bottom of the firn according to the barometric equation: with m being the mass difference between the light and heavy isotope, g the acceleration constant, z the firn depth, R the ideal gas constant, and T mean the mean firn temperature.An increase in accumulation rate increases the firn column depth and therefore increases δ 15 N grav ; on the other hand, a high temperature accelerates the densification processes and shallows the LID.
-thermal effect (Severinghaus et al., 1998): the cold part of the firn is enriched in heavy isotopes according to with T t and T b being the temperatures of the top and bottom parcel, respectively, α T the thermal diffusion constant, the thermal diffusion sensitivity (Grachev and Severinghaus, 2003), and T the temperature difference between top and bottom of the firn.A transient temperature increase after a stable cold period will create a transient peak in δ 15 N therm .
The model needs input temperature, accumulation and dating scenarios with a depth-age correspondence.In the standard version of the Goujon model, the temperature scenario is based on a tuned variable relationship between water isotopes and surface firn temperature, with where α and β can be variable over time.The reconstructed temperature has thus the shape of the water isotope profile but the temperature change amplitudes are constrained by tuning α and β in order for the modeled δ 15 N to match the measured δ 15 N. Several earlier studies have shown that the temporal values of α are lower than the present-day spatial slope for Greenland of 0.80 ‰ • C −1 (Sjolte et al., 2011;Masson-Delmotte et al., 2011), which can be used as a maximum value.

Temperature and accumulation reconstruction
To reconstruct continuous temperature and accumulation scenarios for DO 8 to 10, we run the firnification model from 60 to 30 ka b2k with a time step of one year and try to reproduce the δ 15 N data as well as the depth match points.Figure 3 shows the comparison between the measured and modeled (scenarios d1 to d3) δ 15 N over DO 8-10 at NEEM.First we try to reproduce the δ 15 N data by varying the temperature alone: the measured δ 15 N amplitudes of DO 8, 9, and 10 can be reproduced with temperature increases at the GS-GI transitions of 8.8, 6.0, and 7.7 • C, respectively (Fig. 3, reconstruction d1).This scenario nicely reproduces both the mean δ 15 N level and the amplitude of the δ 15 N peaks.However, the modeled δ 15 N peaks are systematically at a too shallow depth.To model a larger depth, we systematically lower the temperature scenario used in reconstruction d1 (Fig. 3) by 3.5 • C.This systematically deepens the LID, increasing both depth and δ 15 N (Fig. 3, reconstruction d2).The modeled depth is therefore closer to the measured one and the amplitude of the δ 15 N peaks is still correct but the mean δ 15 N level is systematically too high.From this experiment, we conclude that it is not possible to match both δ 15 N data and depth by tuning only the temperature scenario.
Several explanations can be proposed to explain the underestimation of the depth by the model: -the tuning of the Goujon model (LID density, vertical velocity field) is not appropriate for the NEEM site and predicts a too shallow LID.However, we show in Appendix A3 that different tuning strategies have no impact on the modeled LID; -the Goujon model is not appropriate for the NEEM site.However, this model is valid for present-day at NEEM (see Appendix A1) and has also been validated for a large range of temperature and accumulation rates covering the expected glacial climatic conditions at NEEM (Arnaud et al., 2000;Goujon et al., 2003;Landais et al., 2006).Moreover, using other firnification models (the Schwander model on NGRIP Huber et al. (2006) and a Herron Langway model on NEEM, see Appendix C) with similar forcings in temperature and accumulation rate does not reproduce the measured depth either; -fundamental parameters are missing in the description of current firnification models.A recent study has shown that the firn density profile could be strongly influenced by impurities, the density increasing with calcium and dust content in the ice (Hörhold et al., 2012).Calcium and dust in Greenland ice cores are both originating from low-latitude Asian deserts and their content is influenced by source strength and transport conditions (Svensson et al., 2000;Ruth et al., 2007).They covary in Greenland ice cores from seasonal to millennial timescales (Hörhold et al., 2012;Thomas et al., 2008;Steffensen et al., 2008).During cold periods (glacials, stadials), calcium and dust content in Greenland ice cores are strongly enhanced compared to warm periods (interglacials, interstadials) (Mayewski et al., 1997;Ruth et al., 2007;Wolff et al., 2009).Taking this effect into account, the modeled LID during the glacial period should be shallower than the one calculated with the current version of the firnification model calibrated on present-day observations.This would further en-hance the disagreement between modeled and observed depth; -the glacial firn layer at NEEM could be subject to variations in the extent of the convective zone due to katabatic winds and/or low accumulation rate.Considering present-day Antarctic sites as an analogue for the past NEEM firn, a convective zone of 0-3 m (like Dome C, Landais et al., 2005) to 12 m (like Vostok, Bender et al., 1994a) can be considered (Sect.A4 and Fig. A2), event though we suggest that the range 0-3 m is most likely (see discussion in Sect.A4).A convective zone has no direct impact on the depth but it lowers the δ 15 N level.Accounting for such convective zone in the firn model requires lower temperatures, which increases the depth.
-the forcing in accumulation of the firnification model is not correct.To match the observed depth with a correctly modeled δ 15 N, we need to significantly decrease the accumulation rate compared to the original DJ estimation.
With a 2 m convective zone, by adjusting changes in accumulation rate and the δ 18 O-temperature relationship (Fig. 3c, b), we manage to reproduce the δ 15 N profile as presented in Fig. 3, scenario d3.This best δ 15 N fit corresponds to a mean accumulation reduction of 34 % (30 to 40 %, depending on the DO event).Because the depth-age correspondence is imposed by the layer counting, this accumulation rate reduction by 34 % directly implies the same 34 % decrease in the ice thinning.If we use this accumulation scenario as input for the DJ model, with keeping the original DJ accumulation scenario in the remaining ice core sections, the output timescale is just at the limit of the age uncertainty estimated by annual layer counting.With a 12 m convective zone (Fig. 3, d4), the δ 15 N profile can be reproduced using the temperature scenario d3 systematically lowered by 2 • C and the DJ accumulation rate reduced by 28 %.
For NGRIP, the Goujon model can reproduce the measured δ 15 N profile with the correct depth when using a convective zone of 2 m and the DJ accumulation rate reduced by 26 % over the whole section (Fig. D1).Alternatively, we can use a 12 m convective zone with a 19 % reduction in accumulation; this impacts the mean temperature level, which has to be lowered by 2 • C (not shown).For GRIP, using a 2 m convective zone, we have to decrease the DJ accumulation rate by 40 %.Note that this reduced GRIP DJ accumulation rate is then very close to the GISP2 accumulation rate for the same period (Fig. 4).We further discuss past changes in accumulation rate in Sect.3.4.
Based on these calculations, we conclude that reducing the DJ accumulation scenario is necessary to match both δ 15 N data and depth with a firnification model over the sequence of DO 8-10, even when accounting for uncertainties linked with the presence of a convective zone.This reduction has no impact on the reconstructed rapid temperature variations but requires a lower mean temperature level (Fig. 3, d3, d4).Our 19 to 26 % accumulation reduction for NGRIP supports the findings by Huber et al. (2006) where the original accumulation scenario was reduced by 20 %, without convective zone.

Uncertainties quantification
Following the same method for the four cores, we estimate the uncertainty (1σ ) associated with the temperature increases T at the onset of the DO events to be ∼ 0.6 • C for NEEM, GRIP, and GISP2, and ∼ 1.5 • C for NGRIP.For the δ 18 O increases, δ 18 O, the uncertainty is estimated to be ∼ 0.05 ‰ for NEEM, 0.04 ‰ for NGRIP, 0.06 ‰ for GRIP and 0.02 ‰ for GISP2.The thermal sensitivity of δ 18 O, defined as α = δ 18 O/ T, is associated with an uncertainty of 0.05, 0.08, 0.04 and 0.02 ‰ • C −1 for NEEM, NGRIP, GRIP and GISP2, respectively.The detailed calculations are given in Appendix B. Note that for DO 8, 9 and 10, the GRIP and GISP2 ice cores depict δ 18 O increases significantly different from each other even though they are geographically very close to each other.Grootes et al. (1993) calculated an 89 % common variance between these two cores for the interval 9-104 ka b2k and suggested local variability to explain the remaining differences.

Regional δ 18 O and temperature patterns
Our best guess temperature and accumulation reconstructions for the four Greenland sites are displayed in Fig. 4 as a function of the GICC05 timescale.Our temperature re-construction for NGRIP is in good agreement with the one from Huber et al. (2006) where a different firnification model was used (see Fig. D1 in Appendix D).For the GISP2 core, the temperature reconstruction for DO 8 follows the same approach (Orsi et al., 2013): temperature and accumulation scenarios are used as inputs to the Goujon firnification model and constrained using δ 15 N and δ 40 Ar measurements.Four different accumulation scenarios were used, with a GS to GI increase of 2, 2.5, 3, and 3.5 times.The difference in temperature increase between these four scenarios is very small (1σ = 0.07 • C).We report in Table 3 the mean temperature increase for these four scenarios.
For a systematic comparison between the different ice core records, we have used a ramp-fitting approach (Mudelsee, 2000) to quantify the start, end and amplitude of DO increases in δ 18 O, temperature and accumulation: each parameter is assumed to change linearly between GS and GI states.The magnitude of DO increases are then estimated as the difference between the mean GS and GI values (Table 3).The time periods used on each DO event for this statistical analysis are shown in Figs.A3 and D1.

Temperature sensitivity of δ 18 O for present-day and glacial climate
For all four sites, the temporal sensitivity of water isotopes to temperature varies from 0.34 to 0.63 ‰ • C −1 , being therefore systematically smaller than the presentday spatial gradient of 0.80 ‰ • C −1 (Table 3 and Sjolte et al., 2011).This reduction can be explained by precipitation intermittency/seasonality effects (Steig et al., 1994;Jouzel et al., 1997): under glacial boundary conditions, atmospheric models depict a shift of Greenland precipitation towards summer; this has been linked to a southward shift of the winter storm tracks due to the position of the Laurentide ice sheet (Werner et al., 2000(Werner et al., , 2001;;Krinner et al., 1997;Fawcett et al., 1997;Kageyama and Valdes, 2000).During cold periods, summer snow may represent most of the annual accumulation, inducing a bias of the isotopic thermometer towards summer temperature and lowering α compared to the spatial gradient (associated with a classical Rayleigh distillation).So far, seasonality changes have not been systematically investigated in climate model simulations aiming to represent DO events such as driven by freshwater hosing.In reduced sea-ice experiments by Li et al. (2005) using an atmospheric general circulation model, a 7 • C temperature increase and a doubling of the accumulation rate are simulated in GI compared to GS, accompanied with a relatively higher winter snow contribution that could partly explain the low α that we observe here.
Another argument in favor of such seasonality change comes from observations in the NGRIP ice core: records of different ion and dust show synchronous annual peaks during stadials for these species, whereas peaks occur at different periods of the year during interstadials, as for present-day ( Andersen et al., 2006).A first explanation is that during GS the accumulation rate is so low that ion/dust income are dominated by dry deposition all around the year, producing an ion/dust rich layer at the snow surface.This layer is then separated from the one from the following year by the incoming summer snow, resulting in apparent annual synchronous peaks of all species.This supports the hypothesis of a dramatic decrease of winter precipitation during GS at NGRIP.Andersen et al. (2006) also suggested that changes in transport paths may account for the observed pattern.Indeed, the presence of the Laurentide ice sheet (LIS) has been suggested to allow a split jet stream (Andersen et al., 2006).A shift of the path from south to north of the LIS during GI-GS may explain their data.The GS-GI impurities patterns are therefore again in favor of different atmospheric circulation patterns between GS and GI.No such high resolution measurements are yet available for GRIP, GISP2 and NEEM.In the mean time, studies of the second order parameter deuterium excess suggest that the main source of water vapor is shifted southwards during GS (Johnsen et al., 1989;Masson-Delmotte et al., 2005;Jouzel et al., 2007;Ruth et al., 2003).The enhancement of the source-site temperature gradient enhances isotopic distillation and produces precipitation with lower δ 18 O levels during cold periods, increasing α.Contradicting earlier assumptions (Boyle, 1997), conceptual distillation models constrained by GRIP deuterium excess data suggest that this effect is most probably secondary (Masson-Delmotte et al., 2005).

Regional differences between the ice cores sites
The magnitude of GS-GI temperature rise is significantly increasing from NW Greenland to Summit for DO 8 and 10: +8.8 ±0.6 • C at NEEM, +10.4 ±1.5 • C at NGRIP, and +11.1 ±0.6 • C at GISP2 for DO 8; for DO 10, T is largest at NGRIP and smallest at NEEM.In the mean time, the amplitude of δ 18 O is decreasing from NW to central Greenland: for DO 8, δ 18 O is 5.6 ‰ for NEEM, 4.7 ‰ for NGRIP, 4.2 ‰ for GISP2, and 4.6 ‰ for GRIP (Table 3 and Appendix B for uncertainties estimation).As a result, the α coefficient decreases from NEEM to GISP2 for DO 8 and 10.For DO 9, no significant regional difference can be detected in reconstructed temperatures nor in δ 18 O.For this small event, the signal is small compared to the data uncertainties, therefore any spatial gradient becomes difficult to verify.For DO 8 and 10, the larger temporal values of α encountered at NEEM are probably explained by smaller precipitation seasonality effects for this site, which is already biased towards summer at present-day by a factor of 2 to 3.5 (Steen-Larsen et al., 2011;Sjolte et al., 2011;Persson et al., 2011).In other words, because warm periods already undersample the winter snow at NEEM, a winter snow reduction during cold periods at NEEM cannot have an effect as strong as for the NGRIP and the GISP2/GRIP sites, where precipitation are indeed distributed year-round for present-day.We note that α decreases with site elevation (Tables 1 and 3).Interestingly, for DO 8 and 10 the spatial pattern of DO α distribution appears consistent with the spatial patterns of present-day interannual slopes (for summer or winter months), which are also higher in the NW sector (Sjolte et al., 2011).We also note that our spatial patterns of temperature and accumulation increases for DO 8 and 10 are consistent with the pattern obtained by Li et al. (2010).In this study, a 5 • C temperature warming is simulated at Summit with an amplitude decreasing from south to north in response to sea-ice reduction in the Nordic seas.
In addition to differences in seasonality/precipitation intermittency, differences in moisture transportation paths may also modulate the spatial gradients of α over DO events.Although this effect cannot fully explain the fact that the slopes are systematically lower than at present, it could contribute to the difference between the slopes at NEEM, NGRIP and GISP2 (see Sect. 3.3.1).Several studies conducted with isotopic atmospheric general circulation models equipped with water tagging have indeed revealed different isotopic depletions related to the fraction of moisture transported from nearby or more distant moisture sources under glacial conditions (Werner et al., 2001;Charles et al., 2001).In particular, changes in storm tracks were simulated in response to the topographic effect of the Laurentide ice sheet, resulting in the advection of very depleted Pacific moisture towards North Greenland.Indeed, systematic offsets between water stable isotope records of GRIP and NGRIP have been documented during the last glacial period (NGRIP members, 2004).So far, we cannot rule out that changes in moisture origin may cause differences in δ 18 O variations between NEEM, NGRIP, and GISP2.Assessing the importance of source effects will require the combination of deuterium excess and 17 O excess data with regional isotopic modeling and remains beyond the scope of this study.

Past surface accumulation rate reconstruction and glaciological implications
For NEEM, NGRIP and GRIP (reduced DJ-accumulation) as well as GISP2 (original accumulation from Cuffey and Clow, 1997), accumulation variations follow annual layer thickness variations: for each ice core, the smallest accumulation increase is seen for DO 9 and the largest one where the temperature increase is largest (DO 8 for NEEM, DO 8 and 10 for NGRIP).Accumulation shifts therefore scale with temperature variations (Table 3).This is in agreement with the thermodynamic approximation considering the atmospheric vapor content, and thus the amount of precipitation, as an exponential function of the atmospheric temperature.Com-paring the four sites, NEEM and NGRIP show similar accumulation rates, whereas the accumulation is clearly higher at GISP2 and GRIP over the whole time period.
One important finding of our study is the requirement for a lower accumulation rate at NEEM, NGRIP and GRIP over DO 8-10, compared to the initial accumulation rate given by the DJ ice flow model.Taking into account the presence of a possible convective zone at NEEM, our reconstruction based on firn modeling needs to reduce the original DJ accumulation rate by 28 % (12 m convective zone) to 35 % (2 m convective zone).
Several lines of evidence point to an overestimation of the glacial accumulation rates given by the DJ model.First, in their temperature reconstruction for DO 9 to 17, Huber et al. (2006), using the firnification model of Schwander et al. (1997), also had to decrease the accumulation rate calculated by the DJ model by 20 % everywhere to fit the observed depth.We found similar results applying the Goujon model to the NGRIP ice core over DO 8-10.For DO 9 on GRIP, we need to reduce the DJ-accumulation rate by 40 %; the resulting accumulation rate is then very similar to that at GISP2 (Fig. 4).For GISP2 on DO 8, Orsi et al. (2013) used the Goujon model and an accumulation rate of 0.059 m.i.e.a −1 for the GS preceding DO 8, as calculated by the ice flow model from Cuffey and Clow (1997) adapted to the GICC05 timescale.
It is very unlikely that GISP2 and GRIP have significantly different accumulation histories.Indeed, at present GRIP is located 28 km east of GISP2 and the ice divide is not between them, therefore no foehn effect is expected (Buchardt et al., 2012).A small (8 %) accumulation difference is reported for the last 200 a (Meese et al., 1994;Johnsen et al., 1992).During the glacial period, the expansion of the ice sheet margins is expected to produce a flatter topography in central Greenland, further reducing a possible foehn effect.Therefore, large differences in past accumulation rates between GISP2 and GRIP are not climatically plausible; the observed discrepancy must be an artifact of the different methodologies deployed to estimate accumulation rates.Note that during the glacial inception, Landais et al. (2004aLandais et al. ( , 2005) ) were able to reproduce the measured δ 15 N at NGRIP with the original timescale (ss09sea06bm, NGRIP members, 2004) and accumulation values from the DJ model.In the climatic context of the glacial inception, marked by higher temperatures compared to DO 8-10, firnification model and DJ ice flow models seem to agree.
Altogether, these results suggest that the DJ model is consistent with firn constraints during interglacials and inceptions, but a mismatch is obvious during DO 8-10, likely representative of glacial conditions.We now summarize three potential causes that could produce an overestimation of glacial accumulation in the DJ model (for a detailed presentation of this model we refer to Dansgaard and Johnsen, 1969).a. Wrong age scale produced by the DJ model: the DJ model could underestimate the duration between two given depths in the ice core and thus overestimate the accumulation rate.However, for the NEEM ice core from present until 60 ka b2k, the DJ model is tuned in order to produce an age scale in agreement with the GICC05 timescale (Buchardt, 2009).For the NGRIP core, the ss09sea06bm timescale produced by the DJ model together with the accumulation rate has been validated back to 60 ka b2k by comparison to the GICC05 time sale (Svensson et al., 2008).The cumulated uncertainty associated with the GICC05 timescale at 60 ka b2k is 2600 a (Svensson et al., 2008), which could explain 5 % maximum of accumulation reduction, assuming a systematic undercounting of the annual layers.For the glacial inception at NGRIP, Svensson et al. (2011) have counted annual layers on particular sections during DO 25 and the glacial inception and confirm the durations proposed by the ss09sea06bm timescale.We therefore rule out a possible wrong timescale as the main cause for the disagreement on these particular periods.
b.The changing shape of the Greenland ice sheet over time (thickness and margin location) may affect the reconstructed accumulation rate, as suggested by model studies from Cutler et al. (1995) and Cuffey and Clow (1997).The DJ model assumes a constant ice sheet thickness over time for NGRIP (Grinsted and Dahl-Jensen, 2002) and a variable one for NEEM (Buchardt, 2009;Vinther et al., 2009).Concerning the DJ model applied to the NGRIP site, runs with constant ice sheet thickness histories or the one from Vinther et al. (2009) were compared (Buchardt, 2009, Chap. 5, Fig. 5 Johnsen et al., 1995).Tests with a simple DJ model adapted to the GRIP site show that reducing the input surface accumulation rate also reduces the annual layer thicknesses and the total vertical velocity v z , integrated from top to bedrock.In the output timescale the modeled ice age at a certain depth is older.To still get a correct timescale in agreement with GICC05, we would need to deepen the kink height.The same effect would apply to all other Greenland sites.The shape of v z and therefore the kink height is also expected to vary with the ice sheet temperature profile and dust content through changes in ice viscosity (for more details we refer to Cuffey and Paterson, 2010, Chap. 9).During the glacial period, the connection between the Greenland ice sheet and the Ellesmere Island ice could also modify the Greenland ice flow.This effect is expected to affect more the ice flow at NEEM, which is the closest site to Ellesmere Island, than at NGRIP, GRIP and GISP2.In 1969, when creating the DJ model to date the Camp Century ice core, the authors assumed a constant kink height over time due to a lack of information (Dansgaard and Johnsen, 1969).Using a variable kink height will also modify th best guess relationship calculated between input accumulation rate and δ 18 O.
To conclude, there are huge uncertainties on past accumulation rate reconstructions.Comparing three different reconstructions for Summit (Cutler et al., 1995;Cuffey and Clow, 1997;Johnsen et al., 2001), we suggest that the glacial accumulation rate reconstructed by the DJ model has to be taken as a high boundary, the low boundary being at least 50 % lower.This may apply to NGRIP and NEEM.Our firnmodel-based accumulation rates lie in this envelope.We suggest that both thickness and margin location changes should be taken into account in the DJ model.A better agreement between the Cutler et al. (1995) model, the Cuffey and Clow (1997) model and the DJ model may also be found by using a variable kink height in the DJ model.

Conclusions and perspectives
Air and water stable isotopes measurements from four Greenland deep ice cores (GISP2, GRIP, NGRIP and NEEM) have been investigated over a series of Dansgaard-Oeschger events (DO 8, 9 and 10), which are representative of glacial millennial scale variability.We have presented the first δ 15 N data from the NEEM core and combined them with new and previously published δ 15 N data from NGRIP, GRIP and GISP2.Combined with firn modeling, air isotope data allow us to quantify abrupt temperature increases for each ice core site.For DO 8, the reconstructed temperature increase is 8.8 • C for NEEM, 10.4 • C for NGRIP, and 11.1 • C for GISP2.Our data show that for DO 8 and 10, the magnitude of GS-GI increase is up to 2 • C larger in central (GISP2) and North Greenland (NGRIP) than in NW Greenland (NEEM).The reconstructed accumulation increases follow the same spatial pattern.These observations are in agreement with a study of spatial Greenland response to reduced sea-ice extent in the Nordic seas (Li et al., 2010).no spatial gradient is detected for the small DO 9 event.The temporal δ 18 Otemperature relationship varies between 0.3 and 0.6 ‰ • C −1 and is systematically larger at NEEM, possibly due to limited changes in precipitation seasonality compared to GISP2 or NGRIP.The relatively high isotope-temperature relationship for NEEM will have implications for climate reconstructions based on NEEM water isotopes data.Further paleotemperature investigations are needed to assess the stability of this relationship over glacial-interglacial variations.In particular, it would be interesting to compare the presented reconstruction with the temperature-water isotopes relationship over the different climatic context of MIS 5. A better understanding of the causes of the regional isotope and temperature gradients in Greenland requires further investigations of possible source effects (using deuterium excess and O excess), and an improved characterization of atmospheric circulation patterns.We hope that our results will motivate high resolution simulations of DO type changes with climate models equipped with water stable isotopes, in order to test how models capture regional gradients in temperature, accumulation and isotopes, and to understand the causes of these gradients from sensitivity tests (e.g., associated with changes in ice sheet topography, SST patterns, sea ice extent).
The gas age−ice age difference between abrupt warming in water and air isotopes can only be matched with observations when assuming a 26 % (NGRIP) to 40 % (GRIP) lower accumulation rate than derived from the Dansgaard-Johnsen ice flow model.We question the validity of the DJ model to reconstruct past glacial accumulation rate and recommend on the time interval 42 to 36 ka b2k to use our reduced accumulation scenarios.We also suggest that the DJ ice flow model is too simple to reconstruct a correct accumulation rate all along the ice cores and propose to test the incorporation of variable ice sheet margins location and kink height in this model.Our results call for a systematic evaluation of Greenland temperature and accumulation variations during the last glacial-interglacial cycle, combining continuous δ 15 N measurements with firnification modeling.Using a correct accumulation rate is of high importance to reconstruct accurate ice-and gas-age scales and to calculate fluxes based on concentrations of different species in the ice.Moreover, a better estimation of past surface accumulation rates at precise locations in Greenland would help to constrain past changes in ice flow with implications for ice sheet mass balance and dynamics.

Appendix A The Goujon firnification model: method
The firnification model has only one space dimension and calculates the vertical velocity field along the vertical coordinate and the temperature profile across the entire ice sheet for each time step of one year.In the firn, it calculates the density profile from the surface to the close-off depth.The density profile and the accumulation history allow us to obtain the ice age at LID and, assuming gas age equal to zero at LID, the age.The temperature field from surface to bedrock is then used to reconstruct the density profile in the firn, the firn temperature gradient, and from there the δ 15 N at LID.We follow Goujon et al. (2003) where the LID is defined as the depth where the ratio closed to total porosity reaches 0.13.The model is adapted to each ice core site in terms of vertical velocity field, basal melt rate, ice sheet thickness, elevation, surface temperature and accumulation scenarios (Table 1).We assume a convective zone of 2 m at the top of the firn.

A1 Validation of the Goujon firnification model for present day at NEEM
We use the present-day characteristics of the firn at NEEM to validate the Goujon firnification model.During the 2008 summer field season, a shallow core was drilled at the S2 site at NEEM.Firn air was sampled at different depths from the surface to 80 m depth in this borehole (for more details see Buizert et al., 2012).From these air samples, δ 15 N was measured at LSCE (Fig. A1).The increasing δ 15 N with depth reflects the gravitational fractionation.Given the vertical resolution in the data, we do not see a clear convective zone.Below 62 m depth, δ 15 N is constant: the nondiffusive zone is reached.We thus have a LID of 62 m at NEEM for present-day according to these δ 15 N data only.In Buizert et al. (2012), using measurements of different gases in the firn and several diffusion models, the S2 borehole is described as follows: a convective zone of 3 m, a diffusive zone of 59 m down to 63 m depth (LID), and a nondiffusive zone down to 78.8 m depth (total pore closure depth).Following this description and assuming no thermal effect, we calculated the corresponding gravitational fractionation affecting δ 15 N (Fig. A1, blue line).Annual layer counting of the corresponding shallow core and matching with the GICC05 timescale gives an ice age at LID of 190.6 a b2k ±1 a, and 252.5 a at the total pore closure depth.The age of CO 2 is calculated to be 9.6 a at LID and 69.6 a at the total pore closure depth, producing a age of 181 a and 183 a, respectively.The best estimate for the true age is estimated to be 182 +3/−9 a (Buizert et al., 2012).We observe that from the LID, the age becomes constant within uncertainties.Considering the diffusion coefficient to be 1 for CO 2 and using 1.275 for N 2 as in Buizert et al. (2012), the age of N 2 is 7.5 a at the LID, giving a age of 183 a.  (Buizert et al., 2012).Black dots: modeled δ 15 N with the Goujon firnification model.

Clim
We run the Goujon model using the NEEM07S3 shallow core age scale and δ 18 O for the top 60 m (Steen-Larsen et al., 2011) and the NEEM main core below.The δ 18 O record is used to reconstruct the past temperature variations using α = 0.8 (Sjolte et al., 2011); we use β = 9.8 in order to obtain the measured average present-day temperature of −29 • C (Steen-Larsen et al., 2011).Using the density profile measured along the NEEM07S3 core and the corresponding age scale, the past accumulation history was reconstructed (Steen-Larsen et al., 2011) and used here as input for the firnification model.The firnification model estimates the LID at 61.4 m depth.Modeled δ 15 N values agree well with the measured ones in this region (Fig. A1).At 63 m depth, the estimated ice age is 189 a (according to the NEEM07S3 core dating, which agrees very well with the S2 core dating).Ignoring the gas age at LID thus results in an overestimation of the age by less than 10 yr for present-day.

A2 Reconstruction of the past gas age scale
Simulations with the Goujon model show that at the onset of DO events 8, 9 and 10, the heat diffusion in the firn is slow enough so that the peaks of maximum temperature gradient in the firn are synchronous with the δ 18 O peaks.We thus consider that the peaks of δ 15 N occur at the same time as the δ 18 O peaks.However, the Goujon model has no gas diffusion component and this has two consequences: (a) the gas age at LID, due to the time for air to diffuse in the firn, is assumed to be zero; (b) any broadening of the initial δ 15 N peak by gas diffusion in the firn is not taken into account.For present-day, the gas age at LID is 9.6 a for CO 2 (Buizert et al., 2012) and we calculate it to be 7.5 a for N 2 .The Schwander model calculates a N 2 age up to 20 a at the LID over DO 8 to 10 for NGRIP (Huber et al., 2006).The Goujon model thus systematically overestimates the age by 10 to 20 yr in the glacial period, which is within the mean age uncertainty of 60 yr (Table 2).For the NGRIP core, our temperature reconstruction with the Goujon model (without gas diffusion) is in agreement with the temperature reconstruction from Huber et al. (2006) where the Schwander model (with gas diffusion) is used (Fig. D1).We thus consider that the lack of gas diffusion in the Goujon model has an impact that stays within the error estimate (Appendix Sect.B).
For DO 8 to 10 at NEEM, we present the measured and modeled δ 15 N data plotted on an age scale in Fig. A3.The age calculated by the model (Fig. A3, subplot d) is used to synchronize the gas record to the ice record.We have also reported here the age tie-points from Table 2 and we can see that the modeled age reproduces these points, within the error bar.

A3.1 Vertical velocity field
In the firnification model, we used two different parameterizations for the vertical velocity field: the analytical solution from Lliboutry (1979), as in the original model from Goujon et al. (2003), and a Dansgaard-Johnsen type vertical velocity field (Dansgaard and Johnsen, 1969).In the later case, we used the same parametrization as used in the DJ model to calculate past accumulation rates (ice sheet thickness, kink height, fraction of basal sliding, basal melt rate) and then tried different kink heights between 1000 m and 1500 m above bedrock.All these tests produce the same modeled LID and, hence, the same modeled δ 15 N.The different parameterizations actually produce very similar vertical velocity fields in the firn.Because δ 15 N is only sensitive to processes occurring in the firn, huge modification of the vertical velocity field deep in the ice (for example by modifying the kink height) has no impact here.

A3.2 Basal temperature
We also varied the basal temperature between −2.99 • C as measured at present in the borehole (Simon Sheldon, personal communication, 2012) and −1.68 • C, which is the melting temperature as calculated in Ritz (1992) and can be considered as a maximum basal temperature.There is no difference in the modeled LID.Indeed, the relatively high accumulation rate even in the glacial period makes the burial of the snow layers quite fast.As a result, the firn temperature is mostly influenced by the surface temperature but not by the bedrock temperature.

A4 Effect of a convective zone
The Goujon model is written in a way that there is neither heat diffusion nor densification of the firn in the convective zone.No densification of the top firn may be acceptable on the first 2 m but becomes unrealistic when increasing the convective zone.Alternatively, instead of a "true" convective zone, we can prescribe a nongravitational zone, where heat transfert in ice still occurs.This reduces the diffusive column height (DCH) used to calculate the gravitational enrichment of δ 15 N but does not modify its thermal fractionation.
A convective zone deeper than 2 m during the glacial period at NEEM is possible, created by strong katabatic winds due to a steep ice sheet flank, like the 14 m convective zone at YM85 site in Antarctica (Fig. A2 and Kawamura et al., 2006).However, during the glacial period, the Greenland ice sheet may have been connected to Ellesmere Island and the lateral margins were extended compared to present.This would create a flatter surface at the NEEM site, possibly also NGRIP, and would not favor the existence of strong katabatic winds.Marshall and Koutnik (2006) modeled the icebergs delivery from the Northern Hemisphere ice sheets over DO events and showed that the ice sheet margins from Greenland and the Canadian Arctic do not particularly respond to DO events, because these regions remain too cold even during GI.This is not in favor of abrupt change in the convective zone due to ice sheet shape changes at the GS-GI transition.
A convective zone may also be created by a low accumulation rate (as observed at Vostok or Dome F, Fig. A2).The deepest known convective zone (23 m) is reported at the zero-accumulation site Megadunes in Antarctica (Severinghaus et al., 2010).However, note that there is no convective zone at Dome C (Landais et al., 2006) where the present-day annual mean accumulation rate is 2.5 cm w.e.a −1 , slightly higher than at Vostok and Dome F. Our best guess accumulation rate for NEEM, NGRIP and GRIP during MIS3, using a 2 m convective zone, is always higher than at Dome C today (Fig. A2).This is also true for the GISP2 site (Orsi et al., 2013).All these observations are in favor of no deep convective zone at NEEM, NGRIP, GRIP and GISP2 during MIS3.
The existence of a convective zone would affect the average level of δ 15 N, through the reduction of the diffusive zone, but not the modeled depth that is a function of the total firn thickness (LID), itself dependant of surface temperature and accumulation.For NEEM, to reproduce both the measured depth and δ 15 N values, using the original DJ accumulation rate, we need to reduce the temperature scenario d3 by 9 • C everywhere and use a 50 m convective zone (not shown).The obtained system of temperature-accumulationconvective zone is inconsistent with present-day observations in Greenland and Antarctica (accumulation rate much too large compared to the temperature and convective zone).
As a sensitivity test, we have made simulations with the Goujon model using a constant convective zone of 12 m during MIS3 for NEEM (Fig. 3) and NGRIP (not shown), where we still need to reduce the accumulation rate to match the measured δ 15 N profile (by 28 % for NEEM and 19 % for NGRIP).

B1 Temperature increase
The uncertainty associated with temperature reconstruction arises from the contribution of several sources of uncertainties: analytical uncertainties associated with δ 15 N measurements, uncertainty associated with the estimation of the δ 15 N temperature sensitivity ( parameter), uncertainty related to modeling of firn heat diffusion and firnification.In a simple way, based on Eq. (2), we can write the temperature increase T as in the ice, and is the thermal diffusion sensitivity (Grachev and Severinghaus, 2003).
To sum the uncertainties we use the general formula (Press et al., 2007): where x is a function of a, b, and c associated with, respectively, σ a , σ b , and σ c as standard errors.We can thus sum the uncertainties associated to the temperature increase: σ T , δ 15 N therm : this uncertainty results from the analytical uncertainty for δ 15 N measurements and from firnification modeling uncertainty.The pooled standard deviation of our NEEM δ 15 N measurement is 0.006 ‰.We use n GI points to define the GI value (respectively 1, 2, and 3 points for DO 8, 9, and 10) and n GS points for the GS value (respectively 9, 3, and 6 points).For a GS to GI increase, the δ 15 N uncertainty is thus: which gives, respectively, 0.006, 0.006, and 0.004 ‰ for DO 8, 9, and 10.We run the firnification model with a modified temperature scenario in order to exceed the δ 15 N peak value by 0.006 ‰ maximum, for each DO event.The accumulation scenario is kept unchanged.The obtained temperature increase is 0.58 • C larger.If we calculate σ T , δ 15 N as given by Eq. (B3) we obtain 0.52 • C. We conclude that the maximum associated temperature uncertainty is 0.58 • C. Concerning the validity of the firnification modeling, we have already shown in Sect.3.1 that numerous tuning tests performed with the Goujon model do not modify the estimated temperature increase.When using different firnification models (Schwander or Goujon) with similar inputs scenarios, the modeled δ 15 N profiles are similar.Moreover, the duration of temperature increase is well constrained by the GICC05 chronology and the high resolution δ 18 O data.The GICC05 dating and the identification of numerous age tie points (Table 2) between gas-and ice phases gives strong constrains on the accumulation scenario.We are thus quite confident in the validity of our firnification model to reconstruct past surface temperature and accumulation variations.σ T ,D : since the duration of the temperature increase is very well known, the uncertainty on the heat diffusion effect is thus rather small.In our case, it decreases the firn temperature gradient by 1.66 • C with respect to a surface temperature increase of 9.02 • C for DO 8 for NEEM.The major uncertainty in the heat diffusion model is linked to snow/ice conductivity modelisation.For the snow conductivity, we use the formulation from Schwander et al. (1997) where it is a function of the ice conductivity.We have tried different formulations for the ice conductivity (Weller and Schwerdtfeger, 1971;Yen, 1981) and modeled δ 15 N are very similar.
σ T , : we calculate the uncertainty linked to uncertainty (±3 %, Grachev and Severinghaus, 2003) to be 0.22, 0.13, and 0.17 • C for DO 8, 9, and 10 respectively for the NEEM core.This uncertainty increases with the estimated temperature increase and is therefore higher for DO 8.
Summing up all these uncertainties, we estimate the error on the reconstructed NEEM temperature increase for each DO event to be ∼ 0.6 • C (1σ ).Following the same approach, we estimate the uncertainty to be ∼ 1.5 • C (1σ ) for NGRIP and ∼ 0.6 • C for GRIP.For GISP2, we also add an uncertainty of 0.07 • C linked to the amplitude of the accumulation increase (Orsi et al., 2013) and obtain 0.6 • C. The larger uncertainty at NGRIP is mainly caused by the larger analytical uncertainty (0.006 ‰ for NEEM and GRIP δ 15 N data, 0.02 ‰ for NGRIP δ 15 N data).

B2 δ 18 O ice temperature sensitivity
The temperature sensitivity of δ 18 O measured in the ice is defined by the parameter α = δ 18 O/ T .For NEEM, δ 18 O data are measured with an accuracy of 0.07 ‰.We use 2 δ 18 O points to estimate the GI value and 12 (DO 8 and 10) or 9 (DO 9) data points to estimate the GS δ 18 O value.The uncertainty associated with the δ 18 O increase is which gives 0.05 ‰ for DO 8, 9, and 10.Following the same approach, we obtain uncertainties for, respectively, DO 8, 9, and 10 of 0.04, 0.05, and 0.04 ‰ for the NGRIP ice core, and 0.04, 0.08, and 0.06 ‰ for the GRIP ice core.We obtain 0.02 ‰ for the GISP2 core.

B3 Confidence intervals
We propose here two different approaches to determine whether two increases (δ 18 O or temperature), i and j , are different from each other.
First approach: for each calculated i we have calculated an associated error estimate σ i (see previous section).Assuming a Gaussian distribution, each i can be defined by a Gaussian distribution function centered in i: The probability to get a value x ± a i for this i is the corresponding integral: The integral from −∞ to +∞ equals one.For each pair of increases i and j , we calculate the respective maximum confidence intervals [−a i ; a i ] and [−a j ; a j ] in a way that these two intervals do not overlap and that −a j g(x).The integrals p i = p j correspond to the probability to have the value i ± a i and j ± a j .We consider that if p i = p j ≥ 0.9, the two increases are significantly different from each other.
Second approach: we take a couple of temperatures increases T i and T j , assume both to have Gaussian distribution, and calculate the probability of the difference T i -T j to have a value X: Equation (B11) can then be used to find the probability that X = 0, hence that T i = T j .The probability that T i is exactly equal to T j is nil by definition and we have to define the depth interval over which to integrate Eq. (B11).We will thus calculate the probability that T i = T j ± a by assigning the integration range for Eq.(B11) to [−a; a].There is a necessary subjectivity in the choice of this interval.We have chosen to base this estimate of a on the uncertainty associated with our data so that we take a = σ 2 i + σ 2 j : Finally, we consider T i and T j to be significantly different when p(− √ σ i + σ j < x < √ σ i + σ j ) is equal or less than 0.1.We now apply this probability calculation for DO 8 for the following comparison: -NEEM vs. NGRIP: a = 1.6 • C, T NEEM -T NGRIP = 1.6 • C so that the probability that T NEEM = T NGRIP ± 1.6 is calculated with Eq. (B14) as p = 0.48.
These calculations give the same conclusion as obtained from approach 1: for DO 8, T NEEM is significantly different from T GISP2 with 94 % confidence (first approach) or, according to approach 2, these two temperature increases are significantly equal with 5 % confidence.
The results are reported in Table 3 in the following way: for each DO event, if two values are written in two different colors, they are significantly different from each other.A value written in italic is not significantly different from all the others.

Appendix C NEEM firn modelisation with the Herron Langway model
The Herron Langway model (hereafter HL) is an empirical firnification model where the density profile and the ice age in the firn are calculated based on surface temperature and accumulation (Herron and Langway, 1980).We use a surface snow density of 0.350 g cm −3 , as in the Goujon firnification model.Based on the HL density profile in the firn, we calculate the ratio closed to total porosity along the firn column (Goujon et al., 2003).To allow comparison with the Goujon model, we use the same definition for the LID: the depth where the ratio close to total porosity reaches 0.13.At this depth, the HL model gives us the ice age, that we use as a age estimate, and we calculate δ 15 N grav assuming a convective zone of 2 m.This model has no heat diffusion component and we thus use it on periods where the Goujon model shows negligible thermal fractionation for δ 15 N (within the δ 15 N measurement uncertainty), meaning where the surface temperature is stable, without temperature gradient in the firn.We thus can use δ 15 N grav as an estimate for δ 15 N tot .
Here, we apply this model on the stadial periods at NEEM to investigate the surface temperature and accumulation scenarios that match the right δ 15 N level and age.We use the δ 15 N and age values just preceding the DO events (see Table 2) as target values and tune the surface temperature and accumulation.
For DO 8, we use a target δ 15 N value of 0.382 ± 0.006 ‰ and age value of 1198 ± 79 a (Table 2, NEEM tie point n. 2).The HL model can reproduce these values using a surface temperature of −46.76 ± 0.3 • C and an accumulation rate of 0.043 ± 0.004 m.i.e.a −1 (58 % reduction of the one determined by the DJ ice flow model).The LID is at 76 real meters depth (or 52.8 m.i.e.meters ice equivalent).
Modeled surface temperature and accumulation for the onset of DO 8 and 10 are plotted in Fig. A3 with green dots.For DO 8, the HL and Goujon models produce very similar surface temperature scenarios but the HL accumulation rate is lower.For DO 10, the HL and Goujon accumulation rates are similar but the HL temperature is much higher.It is very likely that the differences are due to the strong assumption of no thermal gradient in the firn for the HL model.In order to fit the measured GS level of δ 15 N and the age at the onset of DO 8 and 10 with the HL model, we need to use a reduced accumulation rate by 58 and 51 %, respectively.We here confirm the finding from the Goujon model: decreasing significantly the accumulation rate estimated by the DJ ice flow model is necessary to match both δ 15 N and age data.

Appendix D Reprocessing NGRIP δ 15 N data
To allow the comparison between NEEM and NGRIP, we reconstruct here past temperature and accumulation at NGRIP following the same approach as for the NEEM site.We use the Goujon firnification model adapted to the NGRIP site.
To constrain the model, we minimize the distance between the measured δ 15 N and the modeled one (Fig. D1, c).The corresponding temperature and accumulation scenarios are reported in Fig. D1.For comparison, we also report here the temperature reconstruction from Huber et al. (2006), using the firnification model from Schwander et al. (1997) and the ss09sea06bm timescale.Direct comparison is possible because over DO 8, 9, and 10 (2020 to 2140 m depth), durations proposed by the GICC05 and the ss09sea06bm timescales agree with each other with 5 % difference.Note that the two reconstructions agree well with each other, both for absolute temperature level and temperature variations with time.We use an accumulation rate reduced by 26 % (20 % for Huber et al., 2006) and thus need to reduce the mean temperature level slightly more than Huber et al. (2006) to still match the δ 15 N data.
by Copernicus Publications on behalf of the European Geosciences Union.source: https://doi.org/10.7892/boris.51768| downloaded: 21.9.2017M. Guillevic et al.: Past spatial gradients of temperature, accumulation and δ 18 O in Greenland reduced North Atlantic Deep Water formation

1032M.
Guillevic et al.: Past spatial gradients of temperature, accumulation and δ 18 O in Greenland 2.1.2δ 18 O water isotope data
Fig. A1.Present-day firn at NEEM.Red dots: δ 15 N data measured at LSCE.Blue line: gravitational fractionation for δ 15 N, assuming a convective zone of 3 m and the LID at 63 m depth(Buizert et al., 2012).Black dots: modeled δ 15 N with the Goujon firnification model.
Fig. A3.δ 15 N reconstruction for NEEM (left) and GRIP (right).(a) δ 18 O profile used to reconstruct the surface temperature profile.(b) Surface temperature scenario used in the Goujon model (red line) and in the Herron Langway model (NEEM, 2 green dots).(c) Accumulation scenario for the Goujon model (pink line) and the HL model (NEEM, 2 green dots).(d) δ 15 N data measured at LSCE, this study (red (NEEM) and green (GRIP) line and markers) and modeled δ 15 N by the Goujon model (orange line), using temperature and accumulation scenario shown in (b) and (c), plotted on the ice age scale using the age produced by the Goujon model (e, black line).age tie points are numbered as in Table 2. (e) Markers: tie points between δ 15 N and δ 18 O (blue markers) or accumulation (pink markers), used to constrain the firnification models.Black line: age modeled by the Goujon model.Subplots (a), (b) and (c), in black: GS and GI mean states calculated by the rampfit method.
Fig. D1 δ 15 N reconstruction for NGRIP.(a), δ 18 O profile used to reconstruct the surface temperature profile, NGRIP members (2004).(b) Surface temperature scenario used in the Goujon model (red line, this study) and in Huber et al. (2006), black line.Note that the temperature reconstruction covering DO 8 was missing at that time.(c) Accumulation scenario for the Goujon model.(d) Measured δ 15 N data at the University of Bern, Switzerland, with error bar shown on the last point to the right (dark blue dots Huber et al. (2006), new points in light blue dots).Modeled δ 15 N by the Goujon model (orange line), using temperature and accumulation scenario shown in (b) and (c).Measured and modeled δ 15 N are plotted on the ice age scale using the age produced by the Goujon model (e, black line).age tie points are numbered as in Table 2. Black line: modeled δ 15 N from Huber et al. (2006) using the temperature scenario in (b), black line.(e) Blue markers: tie points between δ 15 N and δ 18 O, used to constrain the firnification models.Black line: age modeled by the Goujon model.Subplots (a), (b) and (c), in black: stadial and interstatial mean states calculated by the rampfit method.

Table 3 .
DO increase (GI-GS) in δ 18 O, temperature and accumulation (mm.i.e.a −1 ) for NEEM, NGRIP, GRIP and GISP2.For δ 18 O, T and α, for each given DO event, two numbers written with different colors are significantly different from each other.A number in italic is not significantly different from all other numbers.See text for details and Sect.B for uncertainties quantification.
, all the other parameters being kept constant (basal sliding, basal melt rate, kink height).The best guess input accumulation rate is a tuned exponential function of δ 18 O (e.g., c.The DJ model assumes that the vertical velocity field (v z ) changes only with surface accumulation rate varia-tions