The first evidence for the succession of rapid climatic events in the northern hemisphere during the last glacial period has been provided by the Greenland ice cores [9,21]. Records for that rapid climatic variability were then obtained in numerous marine [2,10,27] and continental archives [5,15,45] in the northern hemisphere and associate events are depicted in the Antarctic ice cores [1]. Despite this sum of information at lower latitudes, Greenland ice cores still provide the more detailed and continuous climatic record of the last glacial period. Among other things, (i) water isotopes profiles (Fig. 1) depict surface-temperature changes [21], (ii) ice chemical contents reflect changes in atmospheric circulation and source conditions [13], and (iii) greenhouse gases concentration (CH₄, N₂O) fluctuations parallel temperature changes [12]. The water isotopes enable one to clearly describe the succession, over a DO, of a rapid ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$/temperature increase, followed by a warm period (interstadial) and a slower ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$/temperature decrease toward a cold period (stadial).

The conventional reconstruction of the temperature profile is based on the spatial linear relationship observed today between surface temperature and ice isotopic composition $({\delta }^{18}{\mathrm{O}}_{\mathrm{ice}})$ [8] with a slope ${\alpha }_{\mathrm{s}}$ defined by: ${\alpha }_{\mathrm{s}}=\mathrm{\Delta }{\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}/\mathrm{\Delta }T=0.67$. Applying this relationship to reconstruct the past temperature evolution implicitly assumes that the spatial and temporal relationships are similar and that the condensation temperature is linked to the surface one with the same relationship back in time. However, the evaporative origin and the seasonality of precipitation have likely varied during the succession of glacial and interglacial periods through the influences of obliquity and ice-sheet extent [4,11,24,26,46].

Another palaeothermometry method was developed from the direct temperature measurements in the borehole and it is based on the heat equation inversion. The most important result of this method for the glacial period is the estimate of the last glacial maximum temperature to $23\pm 2\text{°C}$ less than today in central Greenland, while the conventional temperature estimate from the water isotopes gave a 12 °C difference [6,7,22]. This dramatic underestimation was attributed to the seasonality of the precipitation that is displaced toward summer during the last glacial maximum because of the ice-sheet extent reducing the winter moisture advection to Greenland [11,26,44]. The borehole method is independent of the different biases affecting the water isotopes temperature reconstruction. However, because of heat diffusion, the temporal resolution of the reconstructed temperature profile is limited and cannot resolve the temperature scenario of the rapid events during the last glacial period.

More recently, an alternative method has been proposed to infer the amplitude of the rapid temperature changes [39,40]. It takes advantage of the thermal diffusion that fractionates the atmospheric gases in the firn (constituted of porous snow in the upper 50–120 m of the polar ice-sheet before the pore close-off). The precise temperature estimate requires a high analytical precision, an accurate calibration of the thermal effect and the use of a heat-diffusion model to link the thermal gradient in the firn to the surface temperature change [23]. This method has already been applied to different rapid climatic events leading to the same conclusion than the one drawn from the borehole temperature measurements: the conventional use of the water isotopes underestimates by a factor 1.5–2 the amplitude of the rapid temperature changes during the last glacial [32,33,39,40].

We concentrate here on this new method to study the temperature evolution over DO in Greenland. We first describe the method and its potential to determine precisely the amplitude of the rapid temperature changes. We then apply it to the surface temperature reconstruction over sequences of DO at the onset of the glacial cycle at NorthGRIP. Finally, we present a comparison of rapid events on which air isotopic measurements were performed to compare the temporal evolution of surface temperature during a variety of DO at the beginning, middle or end part of the last glacial period.

The repartition of air isotopes before enclosure at the bottom of the firn is imposed by thermal and gravitational fractionations if the gas species have a constant atmospheric isotopic composition. The nitrogen isotopic composition can possibly show temporal variations through biological processes (denitrification). However, since the residence time of N₂ in the atmosphere is more than 10⁷ years, even the large consumption of nitrogenous fertilizer during recent time does not modify the atmospheric ${\delta }^{15}$N [42]. Argon is inert with respect to the biogeochemical cycling and therefore, it has been shown that the atmospheric ${\delta }^{40}$Ar was constant over the last million years [42].

The gravitational fractionation is simply imposed by the gravity field on the Earth, g. It drives the heavier isotopes towards the bottom end resulting in the barometric equation: ${\delta }^{15}\mathrm{N}={\delta }^{40}\mathrm{Ar}/4=gz/RT$, where z is the depth, R the gas constant and T the mean firn temperature. More interestingly, during a rapid warming, the top of the firn becomes warm, while the bottom end remains cool creating a thermal gradient. The coldest bottom end is enriched in heavier isotopes according to ${\delta }^{15}\mathrm{N}={\alpha }_{\mathrm{N}}\mathrm{\Delta }T$ and ${\delta }^{40}\mathrm{Ar}={\alpha }_{\mathrm{Ar}}\mathrm{\Delta }T$ (${\alpha }_{\mathrm{N}}$ and ${\alpha }_{\mathrm{Ar}}$ are the thermal diffusion coefficients associated to nitrogen and argon isotopic mixings [17,18]). The positive anomaly is then recorded during the pore close-off.

The ${\delta }^{15}$N and ${\delta }^{40}$Ar measurements in the air enclosed in the ice hence result from two effects. To quantify the thermal effect, two methods were proposed. The first method consists in isolating the thermal fractionation from the gravitational one [40]. Indeed the combination of ${\delta }^{15}$N and ${\delta }^{40}$Ar gives access to the ${\delta }^{15}{\mathrm{N}}_{\mathrm{excess}}=({\alpha }_{\mathrm{N}}-{\alpha }_{\mathrm{Ar}}/4)\mathrm{\Delta }T$, directly proportional to the firn thermal gradient. Because of heat diffusion, the firn thermal gradient is smaller than the surface temperature change. A model of heat diffusion in the firn [16,38] is therefore needed to determine the amplitude of the surface temperature change. The second method [32,33] is based on the development of a powerful firnification and heat-diffusion model [38] predicting the firn depth evolution as well as its instantaneous thermal gradient. It should reproduce both the thermal and gravitational signals. The surface temperature scenario used as input is then adjusted so that the model reproduces the measured ${\delta }^{15}$N signal.

Here, we use the first method owing to a high analytical precision both on ${\delta }^{15}$N and ${\delta }^{40}$Ar measurements (respectively 0.006‰ and 0.025‰). This method offers the advantage of ignoring uncertainty in the gravitational signal influenced by changes in accumulation rate and mean temperature [30–32] still badly constrained (uncertainty up to a factor of two in the accumulation rate [44]). Our approach should be almost accumulation-rate independent [30]. Note that the accumulation rate influence on the close-off depth slightly modifies the heat diffusion length but such influence is negligible (the uncertainty in accumulation rate can create a 10% uncertainty in the heat-diffusion length and consequently a maximum error of 0.3 °C in the inferred temperature increase).

Our ${\delta }^{15}$N and ${\delta }^{40}$Ar measurements were performed on the newly drilled NorthGRIP ice core [37] following the analytical procedure previously described [28,29]. Previous studies have already depicted the rapid changes and their amplitudes during the deglaciation (Bølling Allerød, Younger Dryas [39,40]) and full glacial [30]. We therefore concentrated on the glacial inception captured by the new NorthGRIP ice core (DO 23 and 24, whose rapid temperature increases occurred at $-103.3\text{kyr}$ and $-106.9\text{kyr}$, respectively [37]) and on the Marine Isotopic Stage 4, characterized by a rapid growth of the ice sheet [41] ($\sim -62$ to $-78\text{kyr}$, DO 18, 19 and 20, whose rapid temperature increases occurred at −78.4, −74.5, and $-66\text{kyr}$, respectively [37]).

Fig. 2 depicts the ${\delta }^{15}$N and ${\delta }^{40}$Ar records over the two climatic sequences of interest. As expected, the gas records show a positive/negative anomaly corresponding to each warming/cooling. Because of the close-off depth and of the progressive thinning with increasing depth, the warming is recorded a few meters deeper in the gas record than in the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ one (i.e., at each depth level, the air is younger than the surrounding ice, since it is trapped at the bottom of the firn; note that when the water and air isotopic profiles are drawn on a same timescale, the rapid temperature increase is recorded simultaneously as a rapid ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ increase and sharp ${\delta }^{15}$N and ${\delta }^{40}$Ar peaks). The thermal anomalies are characterized by smallest amplitude for ${\delta }^{40}$Ar/4 than for ${\delta }^{15}$N, in agreement with the different thermal diffusion coefficients for the different isotopic pairs [17,18]. Using this property, the ${\delta }^{15}$N_excess evolution (Fig. 2) can be calculated and a firnification and heat-diffusion model [16] is adjusted to reproduce the amplitude of the main rapid warming events. The NorthGRIP surface temperature is found to have increased by 10 °C over DO 23, 11 °C over DO 18 and 20, and up to 16 °C over DO 19 and 24. The associated uncertainty is of 2.5 °C, taking into account the analytical and the dating errors [30] (we suppose that the timing of the temperature increase can be twice more rapid or twice slower than suggested by the proposed timescale [37] because of the accumulation-rate uncertainty). The choice of a different firnification and heat-diffusion model does not influence our results (not shown here).

We clearly confirm that the spatial relationship applied to the abrupt changes of surface ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ during the last glacial period underestimates the amplitude of the rapid temperature changes by a factor, $\alpha /{\alpha }_{\mathrm{S}}$, of 1.3 (DO 20) to 2.4 (DO 24) ($\alpha =\mathrm{\Delta }{\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}/\mathrm{\Delta }T$ over the rapid warming, ${\alpha }_{\mathrm{S}}=0.67$). It moreover suggests strong atmospheric and/or sea-ice feedbacks around Greenland to explain a 16 °C temperature increase in 100 years, triggered by a change in the thermohaline circulation mode as classically suggested for the occurrence of DO [14]. Simulations with intermediate complexity models can only produce a maximum 8 °C change [14].

The use of air isotopic measurements can moreover be extended to the reconstruction of the complete temperature history over a sequence of DO if a detailed record of ${\delta }^{15}$N and ${\delta }^{40}$Ar is available. Indeed, even during the glacial period, the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ and the surface temperature records do not show the same apparent temporal slope over the different DO ($\alpha =\mathrm{\Delta }{\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}/\mathrm{\Delta }T=0.5$ for DO 20, 0.4 for DO 19 and around 0.3 for DO 18, 23, 24). As a consequence, the proportionality is probably neither respected inside a DO and the surface temperature scenario should not parallel the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ record. To favour this argument, deuterium excess measurements have been performed over different Greenland ice cores [20,25,31]. They suggest that the source for precipitation largely varied between the cold and the warm phases of a DO, such influencing as well the final ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$.

Some small but significant ${\delta }^{15}$N and ${\delta }^{15}{\mathrm{N}}_{\mathrm{excess}}$ anomalies have no clear counterpart in the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ record. We suggest however that they should be the result of rapid temperature changes. Sensitivity studies driven with a firnification and heat-diffusion model suggest that the mean firn temperature has a stronger influence on the gravitational signal of ${\delta }^{15}$N than the accumulation rate. We have used the firnification and heat-diffusion model [16] in a slightly different way than previously described to best reproduce the ${\delta }^{15}$N and ${\delta }^{15}{\mathrm{N}}_{\mathrm{excess}}$ measurements not only over the rapid temperature increase but over the entire DO [31]. This method has been applied to DO 19, where the highly detailed measurements provide enough constraints for this temperature reconstruction (Fig. 2). The reconstructed surface temperature has a significant different shape than the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ profile. The stadial preceding DO 19 (V) is not simply a cold plateau, but a short temperature minimum, followed by a slow warming (IV) before a last cold dip (III). After the main 16 °C increase, our ${\delta }^{15}$N profile imposes that the temperature should have decreased relatively fast (II). Such a scenario was independently confirmed by the combined use of the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ and δD measurements that enables to account for the source temperature influence (Fig. 2, [31]). The stability of ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ during cold phases appears to be an artefact linked to a site-temperature/source-temperature compensation. Note however that this last method does not take into account the seasonality of precipitation and could therefore not reproduce the full amplitude of the temperature change [35].

Finally, the palaeothermometry method based on air isotopic measurements requires a precise analytical setup and the development of a precise firn modelling tool. However, this method provides unique results on the temperature scenario for rapid events that neither the water isotopic measurements nor the borehole temperature measurements can resolve. Using this method enables us to propose a new shape for the high-latitude-temperature scenario over a DO that should be considered for further modelling studies.

Using the ${\delta }^{15}$N profile over a DO provides a temperature evolution without the different biases affecting the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ profile. Since 1998, measurements of ${\delta }^{15}$N and ${\delta }^{40}$Ar were performed over the different DO in central Greenland to first infer the amplitude of the major temperature increases. We present here a synthesis of the published ${\delta }^{15}$N and ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ profiles over DO as well as the inferred temperature changes to describe the variety of DO in Greenland (Fig. 3).

We show a wide variety of DO characterized by rapid warming events ranging from 10 to 16 °C both in Summit (GRIP and GISP2) and in NorthGRIP. The temporal slope between water isotopes and the reconstructed surface-temperature variations confirms that the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ profile cannot be directly used to reconstruct the surface temperature scenario over Greenland: the temporal slope varies from 0.3 to 0.5 over the DO presented here. Finally, comparing the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ and ${\delta }^{15}$N profiles over the different DO leads to the conclusion that their shapes and then the transient surface-temperature evolutions differ from event to event (unfortunately, because of the lack of data, we cannot propose a definite temperature scenario from air isotopic measurements). As a conclusion, describing a ‘typical’ DO through a general scheme is not an easy task.

The last point can be deepened by comparing on the one hand DO 12 (GRIP) and 19 (NorthGRIP) and, on the other hand, DO 18. DO 12 and 19 depict similar ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ and ${\delta }^{15}$N profiles. As the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ is biased by hydrological cycle processes, we concentrate on the ${\delta }^{15}$N one, which suggests a short and small cooling before the main rapid increase (also recorded in the ${\delta }^{15}{\mathrm{N}}_{\mathrm{excess}}$ profile ([30,31], III in Fig. 2). After the main peak, the ${\delta }^{15}$N decreases smoothly and a small peak, corresponding to a rapid warming event, occurs 2000 yr after the main temperature increase (I in Fig. 2). On the contrary, DO 18 shows a rapid temperature increase followed only 200 years later by a rapid temperature decrease according to the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ and ${\delta }^{15}$N profiles. Even if it is highly tempting to suggest, on the basis of common Greenland temperature patterns for DO 12 and 19, a common climatic mechanism, things are more complicated. Indeed, the marine records [3] place the onset of DO 12 after a large iceberg discharge from the Laurentide, the Heinrich event 5. Such event is recorded as a layer rich in ice rafted debris (IRD) and poor in foraminifera between 40°N and 55°N in the Atlantic [2]. According to the same marine records, DO 19 is not preceded by any Heinrich event, ruling out the association of a Heinrich event with a DO depicting the Greenland characteristics of DO 12 or 19. DO 19 is most probably followed by such an event (Heinrich event 6 occurs during the maximum of ice volume at the end of the Marine Isotopic Stage 4, roughly corresponding to DO 18, Fig. 1, [3]). In a similar way, if we expect DO 12 and 18 to be similar (based on the existence of a Heinrich event preceding the rapid surface-temperature increase), the air isotopic measurements performed over those events show that a simple association between the shape of a DO and the presence of a Heinrich event is not straightforward.

To end up with the variety of rapid events, the specificity of DO 23 and 24 should be highlighted. Their ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ and ${\delta }^{15}$N profiles depict high-frequency surface temperature variations (i.e., rapid cooling during interstadials) absent from the other profiles presented here. Such decadal to centennial climatic variability does not seem to be recorded in any marine records (no IRD signal [34]), suggesting that no strong icebergs discharges occurred at that time. We therefore suggest that this rapid climatic variability, occurring during interstadials, is most probably linked to hydrological-cycle processes at that particular climatic period, the glacial inception. It once again confirms that no common scheme can be deduced for the Greenland temperature evolution over a DO and that a direct link with the icebergs discharges as the only modulator of the thermohaline circulation activity is not obvious.

We presented here a recently developed palaeothermometry method that applies to the rapid climatic variability recorded in the Greenland ice cores. Whereas (i) the conventional method based on the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ profile is biased by hydrological processes and (ii) the borehole temperature inversion is unable to capture the rapid climatic variability, the measurements of air isotopes (nitrogen and argon) permit to infer the amplitude of rapid surface-temperature changes as well as the shape of the temperature scenario over DO. We illustrated this method on the recently drilled NorthGRIP ice core and found warming intensities of 10 °C for DO 23, 11 °C for DO 18 and 20 and up to 16 °C for DO 19 and 24. The temperature scenario over DO 19 could also be reconstructed and clearly confirmed that the ${\delta }^{18}{\mathrm{O}}_{\mathrm{ice}}$ profiles in the Greenland ice cores do not fully capture the shape of rapid temperature evolutions. Finally, on the basis of a complete comparison of the existing ${\delta }^{15}$N records of DO in Greenland, we conclude to a wide variety of DO with different shapes and amplitudes as well as different links with the North Atlantic climatic patterns. We especially show that no clear association can be depicted between the existence of a Heinrich event or of another iceberg discharge event and the following amplitude and shape of a DO in Greenland. The climate over Greenland during the last glacial was then subjected to many different influences, such as the polar atmosphere and the extension of northern ice sheets that modifies the hydrological cycle. More definite conclusion should be drawn when a complete ${\delta }^{15}$N/${\delta }^{40}$Ar record is performed over the full glacial period in the NorthGRIP ice core.