Millennial and sub-millennial scale climatic variations recorded in polar ice cores over the last glacial period

. Since its discovery in Greenland ice cores, the millennial scale climatic variability of the last glacial period has been increasingly documented at all latitudes with studies focusing mainly on Marine Isotopic Stage 3 (MIS 3; 28–60 thousand of years before present, hereafter ka) and characterized by short Dansgaard-Oeschger (DO) events. Recent and new results obtained on the EPICA and NorthGRIP ice cores now precisely describe the rapid variations of Antarctic and Greenland temperature during MIS 5 (73.5–123 ka), a time period corresponding to relatively high sea level. The results display a succession of abrupt events associated with long Greenland InterStadial phases (GIS) enabling us to highlight a sub-millennial scale climatic variability depicted by (i) short-lived and abrupt warming events preceding some GIS (precursor-type events) and (ii) abrupt warming events at the end of some GIS (rebound-type The long accompanying Antarctic seesaw questioning the Greenland ice core as a for


346
E. Capron et al.: Millennial scale climatic variations and bipolar seesaw pattern The "iconic" DO event structure is depicted as a GIS, beginning with an abrupt warming of 8 to 16 • C in mean annual surface temperature within a few decades (Severinghaus et al., 1998;Lang et al., 1999;Landais et al., 2004a;Landais et al., 2006; see also Wolff et al., 2009a for a review). The GIS is then usually characterized by a gradual cooling phase lasting several centuries and its end is marked by a rapid cooling towards a relatively stable cold phase (GS) that persists over several centuries to a thousand of years. This description originates mainly from the DO events occurring over Marine Isotopic Stage 3 (MIS 3, 28-60 ka; See Voelker, 2002 for a review) that benefit from a robust chronology ( Fig. 1; e.g. Wang et al., 2001;Shackleton et al., 2003a;Fairbanks et al., 2005;Svensson et al., 2008).
A dynamical combination between ocean, cryosphere (continental ice sheets and sea ice cover), vegetation and atmosphere is at play during this millennial scale variability (Hendy and Kennett, 1999;Peterson et al., 2000;Kiefer et al., 2001;Wang et al., 2001;Broecker, 2003; but the triggering processes of such a variability are still under discussion (Wunsch, 2006;Friedrich et al., 2009). Current theories point to external forcing mechanisms such as periodic changes in solar activity (Bond et al., 1992;Braun et al., 2008), periodic calving of ice sheets (van Kreveld et al., 2000) and to internal oscillations of the ice sheetocean-atmosphere system through freshwater perturbations (Broecker, 1990;MacAyeal, 1993;Ganopolski and Rahmstorf, 2001;Schulz et al., 2002).
In Antarctic ice cores, millennial scale temperature changes are gradual and out of phase with their abrupt northern counterpart ( Fig. 1; Bender et al., 1994;Blunier et al., 1998;Blunier and Brook, 2001;EPICA-communitymembers, 2006;Capron et al., 2010). Such a "bipolar seesaw" pattern is understood as reflecting changes in the strength of the Atlantic Meridional Oceanic Circulation (AMOC; Broecker, 1998). The physical mechanism for this bipolar seesaw pattern has been explored through a large number of conceptual and numerical models of various complexities (e.g. Stocker et al., 1992;Rind et al., 2001;Velinga and Wood, 2002;Knutti et al., 2004;Kageyama et al., 2009;Liu et al., 2009;Swingedouw et al., 2009). Using the simplest possible model, Stocker and Johnsen (2003) successfully described the Antarctic millennial variability in response to the abrupt temperature changes in the north by involving a southern heat reservoir associated with AMOC variations. Such an important role of the Southern Ocean for the bipolar seesaw mechanism is supported by marine records (Barker et al., 2009).
MIS 5 (∼73.5-130 ka; Shackleton, 1987) includes the last glacial inception and the early glacial which are a time period of great interest since they represent an intermediate stage between full interglacial conditions (defined as MIS 5e, Shackleton et al., 2003b;Fig. 1) and glacial conditions encountered during MIS 2-3. At that time, continental ice sheets are extending, corresponding to sea level variations from 20 to 60 m below present-day sea level (Waelbroeck et al., 2002) compared to 120 m below present-day sea level during MIS 2-3. MIS 5 is also marked by a different orbital configuration with stronger eccentricity and therefore larger seasonal insolation changes compared to MIS 3 (Fig. 1).
The NorthGRIP ice core (Greenland, 75 • 10 N, 42 • 32 W; 2917 m a.s.l.; present accumulation rate of 17.5 cm water equivalent per year (cm w.e. yr −1 )) expands the Summit records back to the last interglacial period (∼123 ka) and offers high resolution (1 cm yr −1 ) due to the basal melt limiting thinning processes (NorthGRIP c.m., 2004). The δ 18 O ice profile unveiled GIS 23, 24, and 25 in the early time of the glacial period (∼101-113 ka; Fig. 1; NorthGRIP c.m., 2004). The discovery of only six abrupt climatic events during MIS 5 (GIS 20 to 25) reveals a longer pacing than the ∼1.5 thousand year (hereafter kyr) approximate DO event frequency suggested by Grootes and Stuiver (1997) between 12 and 50 ka which has been strongly debated (e.g. Grootes and Stuiver, 1997;Schulz, 2002;Rahmstorf, 2003;Ditlevsen et al., 2005Ditlevsen et al., , 2007. We also use here an ice core drilled within the European Project for Ice Coring in Antarctica (EPICA) in the interior of Dronning Maud Land (hereafter, denoted EDML, 75 • S, 0 • E, 2892 m a.s.l., present accumulation rate of 6.4 cm w.e. yr −1 ). It represents a South Atlantic counterpart to the Greenland records (EPICA c.m., 2006) and provides a resolution of ∼30 yrs during MIS 3 and of ∼60 yrs during MIS 5 that makes the EDML core particularly suitable for studying millennial scale climatic variations in Antarctica. During MIS 3, it has been shown that the amplitudes of EDML Antarctic Isotopic Maxima (AIM) are linearly related to the duration of the concurrent GS (EPICA c.m., 2006). Moreover, taking advantage of a new common timescale between EDML and NorthGRIP ice cores based on the global signals  (Meyer et al., 2000). (c) MIS 3 (dotted light blue curve) and MIS 5 (dotted orange curve) sea level variations (Bintanja et al., 2005) reconstructed from the LR04 δ 18 O benthic stack (MIS 3, solid dark blue line; MIS 5, solid pink line; Lisiecki and Raymo, 2005). Both the sea level curve and the LR04 δ 18 O benthic stack are displayed on EDC3 timescale. The timescale synchronisation is done in Parrenin et al. (2007a). (d) MIS 3 (dark blue curve) and MIS 5 (pink curve) CO 2 concentration. Composite CO 2 from EDC and Vostok ice cores (Petit et al., 1999;Monnin et al., 2001;Lüthi et al., 2008). (e) MIS 3 and MIS 5 orbital contexts: 65 • N summer insolation (full line) and obliquity (dotted line) (Laskar et al., 2004   In this paper, we combine the full climatic information available from the NorthGRIP and EPICA ice cores in order to provide a complete description of the abrupt climatic oscillations recorded in the polar regions over MIS 5 in comparison to MIS 3. In Sect. 2, we present the common timescales used for the NorthGRIP and EPICA ice cores over MIS 3 and MIS 5. Section 3 deals with new high resolution measurements of EDML δ 18 O ice and of NorthGRIP δ 15 N and δ 40 Ar which are then used to characterize MIS 5 GS and GIS in terms of structure, temperature amplitude and their relationships with their Antarctic counterparts. In particular, we bring new evidence of sub-millennial scale variability at the onset and the end of the MIS 5 long interstadials. We then discuss the robust features and peculiarities of GIS events of MIS 3 and MIS 5 in relationship with their different climate background (i.e. ice volume, orbital contexts) in Sect. 4. Fi-nally, in Sect. 5, we test the general applicability of the thermal bipolar seesaw concept for the entire glacial period by comparing our results with north-south time-series generated through the conceptual model of Stocker and Johnsen (2003).

Synchronising NorthGRIP and EDML ice cores
An accurate timescale is necessary to characterize the duration and pacing of climatic events and the sequence of events between the Northern and the Southern Hemispheres.
For MIS 3, we use the most recent GICC05 age scale (Greenland Ice Core Chronology 05) extended back to 60 ka  for the NorthGRIP ice core. GICC05 is a timescale based on annual layer counting and has been shown to be compatible within the given counting error with absolutely dated reference horizons in the 0-60 ka period Fleitmann et al., 2009) with GICC05 tending to generally underestimate the age. Synchronization (e.g. Bender et al., 1994;Blunier et al., 1998)   . The uncertainty determination is based on the comparison of DO duration inferred from the EDML-NorthGRIP timescale with their duration on other timescales (marine sediment cores: MD95-2042, Shackleton et al. (2004) andNEAP18K, Chapman andShackelton (2002); Sanbao Cave speleothem record, ; lake record from Monticchio, Brauer et al. (2007)). a DO duration on each chronology: (1) EDML-NorthGRIP timescale (2) MD95-2042 timescale (3) NEAP18K timescale, (4) Sambao Cave timescale, (5) Lago di Monticchio timescale. b For each rapid event, we calculate the mean DO duration deduced from DO durations estimated on each chronology (a). c The uncertainty of each event duration is given as a percentage of error calculated as the ratio between the standard deviation of DO durations on each timescale (a) and the mean event duration (b).

Mean DO b
Event c DO duration on each chronology (yrs) a duration (yrs) duration uncertainty (%) (1) EDML-NGRIP (2) MD95-2042 (3) NEAP-18 K (4) Sambao Cave (5)  performed using the isotopic composition of atmospheric oxygen (δ 18 O atm ) and CH 4 records from air entrapped in EDML and NorthGRIP ice (EPICA c.m., 2006; allowed us to place the two ice cores on the same timescale. From 25 to 50 ka, the maximum uncertainty between the two records is estimated to reach 500 yrs at the onset of GIS 12. To study MIS 5, we transfer the NorthGRIP record onto the EDML1 timescale (Parrenin et al., 2007a;Ruth et al., 2007;Severi et al., 2007) between 75 and 123 ka using also a δ 18 O atm and CH 4 synchronisation . The new EDML-NorthGRIP timescale enables us to quantify the exact phasing between the onsets of AIM and GIS with an accuracy of a few centuries except for the onset of GIS 25, where the uncertainty is higher than 1 kyr due to the lack of high-resoluted methane records both on EDML and North-GRIP that would enable one to determine a precise gas age marker between the two records. Unlike for MIS 3, we are not using an absolute timescale for MIS 5 however focussing on the duration and sequence of events only requires a relative timescale. NorthGRIP basal melting induces a timescale almost linearly proportional to depth by reducing the ice thinning (NorthGRIP c.m., 2004;Dahl-Jensen et al., 2003). For the new EDML-NorthGRIP age-scale, we obtain a smooth evolution of age as a function of depth with less than 10% deviation from the slope deduced from the age/depth relationship of the NorthGRIP glaciological timescale (NorthGRIP c.m., 2004). We thus consider our age markers as being consistent with ice flow conditions at the NorthGRIP site.
We compare this new timescale with independent chronologies from other paleoclimatic archives (Table 1) and this enables us to derive uncertainties associated with the duration of each GIS/GS succession over MIS 5. In the North Atlantic region, marine cores show rapid cooling events (C events) (McManus et al., 1994) that were associated with the GS (i.e. event C 24 is associated with GS 25, McManus et al., 1994;NorthGRIP c.m., 2004;Rousseau et al., 2006). Using such associations, NorthGRIP DO event duration is compared to the one deduced from two marine sediment cores: (i) MD95-2042, providing an age scale with two absolute age markers derived from the Hulu cave between 115 and 81 ka (Shackleton et al., 2003a;Wang et al., 2001) and (ii) NEAP18K, whose age model was constructed by correlation of the benthic δ 18 O records with an orbitally tuned δ 18 O stratigraphy (Shackleton and Pisias, 1985).
Then, the same exercise is carried out by comparing our age-scale with the chronology from a lacustrine sediment core from Lago grande di Monticchio (Brauer et al., 2007) whose chronology is based on lamination counting. Finally, assuming synchronous climatic shifts at low and high latitudes enables us to compare rapid event succession recorded in ice cores through independent dating from speleothem records (i.e. Wang et al., 2001. Such a comparison is difficult because few speleothem records display a clear sequence of rapid events over MIS 5 except the record from Sanbao cave on which we can identify the onset of each GIS event . Uncertainties on DO event durations (i.e. durations of GIS plus GS) obtained from the comparison of the five records are summarized in Table 1. In the following we limit our study to the sequence of events 24, 23, 22, and 21 since the EDML-NorthGRIP synchronisation lacks of robust chronological constraints around GIS 25 . The uncertainties associated with the durations of GIS/GS 24, 23, 22 and 21 represent less than 19% of the duration of each DO event. This general agreement makes the Table 2. North-South rapid events recorded in NorthGRIP and EPICA ice cores. a NorthGRIP δ 18 O ice amplitude ( δ 18 O ice ) at the onset of abrupt events (NorthGRIP c.m., 2004). b Accompanying warming amplitude ( T ) estimated from δ 15 N data with associated uncertainty. c Huber et al. (2006) and Landais et al. (2006) provide a quantification of abrupt temperature change through air isotopes measurements of most of the rapid events over the last glacial period and our study provides new results of temperature estimates at the onset of GIS 21 and GIS 22. d Spatial slope deduced at the onset of each rapid event as α= δ 18 O ice / T . Note that a ±2.5 • C uncertainty on temperature change is translated to an error of ∼0.2 in the calculation of α. e GS durations are given (1) on the GICC05 age scale for GS 2 to GS 12 Svensson et al., 2008), (2) on ss09sea age scale for GS 18 to GS 20 (NorthGRIP c.m., 2004) and (3)on the EDML-NorthGRIP synchronised timescale for GS 21 to GS 24 . GS duration is defined by the interval between the midpoint of the stepwise temperature change at the start and end of a stadial. The errors associated with stadial duration are estimated by using different splines through the data that affect the width of the DO transitions and are linked to visual determination of maxima and minima during transitions. No estimate of GS is given where the beginning of the GS is hard to pinpoint due to the particular structure of events and the corresponding events are labeled with #. f AIM warming amplitudes are given for the EDML ice core based on the temperature reconstruction of Stenni et al. (2010). AIM 22 is damped after d-excess corrections and labeled as * (See Stenni et al. (2010) for details). The amplitude is determined following EPICA c.m. (2006) from the Antarctic δ 18 O ice maximum to the preceding minimum of each event. Uncertainties on MIS 3 AIM amplitudes are determined in EPICA c.m. (2006). For AIM events during MIS 2, MIS 4 and MIS 5, we consider that an error bar of ±0.4 • C encompasses the uncertainty on the determination of the warming amplitude by using different splines through the data.
This study 3 5. possibility of a large (greater than 1.6 kyr) error in interstadial duration unlikely and provides a firm basis to confidently analyse the interstadial structure and pacing of these events over MIS 5.

Past temperature reconstructions
The stable isotopic composition of precipitation at mid and high latitudes is related to local air temperature through the so-called spatial slope (Dansgaard, 1964) with an average of 0.67‰-δ 18 O per • C in Greenland (Johnsen et al., 1989) and 0.75‰-δ 18 O per • C in Antarctica Merlivat, 1977, Masson-Delmotte et al., 2008). This spatial slope, which results from air mass distillation processes, can be used as a surrogate for the temporal slope and may vary over time due to past changes in evaporation conditions or atmospheric transport. However, the ice isotopic composition as a tool for past temperature reconstructions (Dansgaard, 1964) has to be interpreted carefully. Indeed, the hypothesis of similar spatial and temporal slopes is only valid if certain assumptions are satisfied, in particular those concerning the origin and seasonality of the precipitation . In fact, past Greenland temperature reconstructions based on the spatial isotope-temperature slope have been challenged by alternative paleothermometry methods such as (i) the inversion of the borehole temperature profile (e.g. Dahl-Jensen et al., 1998;Cuffey and Clow, 1995;Johnsen et al., 1995) and (ii) the thermal diffusion of air in the firn arising during abrupt climate changes (Severinghaus et al., 1998). The quantitative interpretation of δ 18 O ice in term of temperature variations using this spatial relationship leads to a systematic underestimation of temperature changes (Jouzel, 1999). The amplitudes and shapes of temperature changes can be biased by past changes in precipitation seasonality (Fawcett et al., 1996;Krinner et al., 1997), changes in moisture sources during rapid events (Charles et al., 1994;Boyle, 1997;Masson-Delmotte et al., 2005a) and surface elevation (i.e. Vinther et al., 2009).
In contrast, available information suggests that reconstruction of local surface temperature in Antarctic ice cores using the present-day spatial relationship between isotopic composition of the snow and surface temperature ("isotopic thermometer") is correct, with a maximum associated uncertainty of 20% at the glacial-interglacial scale . As a consequence, the isotopic thermometer is commonly used to quantify past changes in temperature based on the stable isotopic composition measured in deep Antarctic ice cores (e.g. EPICA c.m., Jouzel et al., 2007). More recently, a modelling study dedicated to the stability of the temporal isotope-temperature slope suggests that the classical interpretation of the ice core stable isotopes on EDC may lead to an underestimation of past temperatures for periods warmer than present conditions (Sime et al., 2009).

Greenland temperature reconstruction
Measurements of the isotopic composition of nitrogen (δ 15 N) and/or argon (δ 40 Ar) from air trapped in ice cores provide independent quantitative reconstructions of abrupt local temperature changes (Severinghaus et al., 1998;Lang et al., 1999;Severinghaus and Brook, 1999;Landais et al., 2004a, b, c;Huber et al., 2006). Trapped-air δ 15 N variations are only produced by firn isotopic fractionation due to gravitational settling and thermal diffusion since atmospheric δ 15 N is constant over the past million years (Mariotti, 1983). The combined use of δ 15 N and δ 40 Ar data together with the modelling of physical processes (densification, temperature and gas diffusion) enable the estimation of abrupt surface temperature magnitudes with an accuracy of 2.5 • C (Landais et al., 2004a). Table 2  show a warming of 10 • C ±2.5 • C and 16 • C ±2.5 • C, respectively . In the present paper, we complete the quantification of MIS 5 abrupt warming events based on published δ 15 N measurements  with new δ 40 Ar data for the onset of GIS 21 and new δ 15 N measurements over GIS 22 (Fig. 3). Following the approach of Landais et al. (2004a), we find that the GIS 21 onset is marked by a warming of 11±2.5 • C while new high resolution gas measurements over GIS 22 reveal a weak δ 15 N variation (0.063‰) corresponding to a maximum warming amplitude of 5 • C (Fig. 3). Compared to the classical isotope-temperature relationship, the temporal slope between changes in δ 18 O and in δ 15 N-derived temperature at the onset of GIS events is systematically lower and varies from 0.29 to 0.55‰/ • C. Previous studies suggest an effect of obliquity and ice-sheet on the temporal δ 18 O ice /temperature slope mainly via the seasonality of the precipitation and/or moisture source (Denton et al., 2005;Masson-Delmotte et al., 2005a;Flückiger et al., 2008). Here, we do not find any systematic relationship between the evolution of the temporal slope and the long term evolution of components such as ice sheet volume or orbital parameters. We suggest instead that the temporal slope and thus seasonality and/or moisture source change at the GIS scale.
Note that the amplitude of temperature change at the onset of the different GIS in NorthGRIP should not be used as a quantitative reference for Greenland climate. As an example, NorthGRIP and GISP2, only 325 km apart, present different temperature changes at the onset of GIS 19: 16 • C and 14 • C, respectively (Landais et al., 2004a;Landais, 2004). These regional differences are probably due to a more continental climate at NorthGRIP as has been suggested by the comparison of GRIP and NorthGRIP δ 18 O ice curves (North-GRIP c.m., 2004). Indeed, it has been observed that the difference in δ 18 O ice between NorthGRIP and GISP2 is highly correlated with past continental ice volume reconstructions. This suggests that larger ice sheets enhance the remoteness of NorthGRIP from low latitudes air masses while Summit (GRIP/GISP2) is less affected by such continentality effect. Regional differences in moisture origin during the current interglacial period have also been identified in deuterium excess profiles (Masson-Delmotte et al., 2005b). Finally changes in ice sheet topography can also generate regional elevation changes impacting regional temperature at ice core sites (Vinther et al., 2009). However, these differences modulate the regional expressions of climate variability and do not prevent us to use the NorthGRIP δ 18 O ice profile to qualitatively characterise the relative amplitudes and shapes NorthGRIP c.m., 2004).

Antarctica temperature reconstruction
Few studies have used so far the gas fractionation paleothermometry method on Antarctic records (Caillon et al., 2001;Taylor et al., 2004). Based on δ 15 N and δ 40 Ar data performed on the Vostok ice core, Caillon et al. (2001) estimate a temperature change at the MIS 5d/5c transition consistent within 20% of the classical interpretation of water stable isotope fluctuations. In general, the smoother shape of millennial scale temperature changes prevents the use of the isotopic composition of the air to infer local temperature change during AIM events using thermal diffusion. For Antarctic temperature reconstruction, we thus use the δ 18 O ice -based method improved by using the deuterium excess data (d-excess=δD−8×δ 18 O ice ; Dansgaard, 1964) which allows us to correct for the changes in moisture source conditions (e.g. Stenni et al., 2001Stenni et al., , 2010Vimeux et al., 2002).
Here we use two different temperature reconstructions to characterize AIM warming amplitudes in East Antarctica. First, new high resolution and already published EDML δ 18 O ice and EDC δD profiles (Fig. 3) are corrected for global seawater isotopic composition (Bintanja et al., 2005) following Jouzel et al. (2003) and then converted to past temperatures using the observed spatial slope of 0.82×δ 18 O ice per • C for EDML (Oerter et al., 2004) and 6.04×δD per • C for EDC (Lorius and Merlivat, 1977;Masson Delmotte et al., 2008). These temperatures are corrected for elevation variations and changes in ice origin (upstream and elevation correction; for EDML: Huybrechts et al., 2007;for EDC: Parrenin et al., 2007b).
Secondly, we use the EDML and EDC temperature profiles derived using d-excess data recently published by Stenni et al. (2010). Based on isotopic profiles corrected for sea water isotopic composition and upstream effects, they estimate changes in source conditions and site temperature at both sites (Stenni et al., 2001;Masson-Delmotte et al., 2004). These site temperature reconstructions are expected to be more robust because they account for changes in evaporation conditions. However, for some rapid events, this "inverted" temperature is noisier and makes the beginning and the end of the warming more difficult to identify ( Fig. 3; Stenni et al., 2010). Thus, both "classical" and "site" temperature reconstructions are used for EDC and EDML to identify AIM and associated warming amplitudes.
The AIM temperature change amplitudes are presented in Table 2. We estimate a maximum uncertainty of 0.4 • C based on (i) the warming amplitude difference between the two reconstructions and (ii) the determination of the minima and maxima for Antarctic temperature changes. The moisture correction does not have any major impact on AIM amplitudes nor on their shapes.
Note that Antarctic events with amplitudes of less than 1‰ (equivalent to 0.5 • C) remain delicate to interpret since corrections based on deuterium excess profiles add noise to the temperature reconstruction. We choose hereafter to discuss only small events that have a Greenland counterpart and to ignore the other small δ 18 O ice fluctuations.

Structure of MIS 5 abrupt climate variability
The exceptionally long duration of the GIS during MIS 5 reveals an additional variability within the classical GIS/GS succession. Three types of sub-millennial scale events are identified: (i) short-lived and sharp warming preceding GIS 21 and GIS 23, (ii) abrupt warming during the cooling phase of GIS 21 and (iii) abrupt cooling phases during GIS 24. 46 4) CH 4  and δ 15 N (this study).

Precursor-type peak events
The NorthGRIP isotopic profile contains dramatic reversals in δ 18 O ice before the two longest GIS of the entire glacial period: GIS 21 and GIS 23 (Fig. 4). The first occurs within 200 yrs with a 2.2‰ variation in δ 18 O ice . After a short (100 yrs) return to cold conditions, GIS 21 onset occurs with a 4.2‰ increase in δ 18 O ice . Such a precursor-type struc-ture is also visible before GIS 23 onset: the δ 18 O ice rises by 3.8‰ in 125 yrs at ∼102.5 ka and then drops by 3.6‰ in 100 yrs (hereafter, denoted as GIS 23b). The return to stadial conditions lasts ∼300 yrs before δ 18 O ice increases by 3‰ at the onset of GIS 23. The occurrence of precursor events is confirmed by parallel δ 18 O ice variations measured in GRIP and GISP2 cores Grootes and Stuiver, 1997; and by their detection in high-resolution records of δ 15 N data and CH 4 . The precursor-type peak event leading GIS 23 exhibits a 200 ppbv variation in CH 4 and a 0.050‰ rapid increase in δ 15 N. The reversal prior to GIS 21 is weaker in the NorthGRIP δ 18 O ice profile (2.2‰) than GISP2 (3.4‰) whereas NorthGRIP δ 15 N data over this reversal indicate a 0.08‰ variation in NorthGRIP comparable to the δ 15 N variation of 0.09‰ measured in GISP2 .
Note that these very short δ 18 O ice variations are not only visible during MIS 5. Indeed, the sequence of DO events 13 to 17 during MIS 3 is also extremely unstable with short temperature peaks of 200-400 yrs accompanied by fast shifts in CH 4 concentration (Blunier and Brook, 2001;Huber et al., 2006). This highlights that abrupt climatic variability over the glacial period is more complex than the millennial scale variations expressed by a GIS/GS succession.

Rebound-type events
At the end of the regular cooling phase of GIS 21, δ 18 O ice increases abruptly (∼2‰ in less than 100 yrs) 1.2 kyrs before the sharp return to stadial conditions (Fig. 5). The large scale imprint of the GIS 21 sub-event is detected through GISP2 high resolution CH 4 data (showing a 71 ppbv increase in 140 yrs; Grachev et al., 2009). This "rebound" pattern is identical in δ 18 O ice magnitude, duration and structure to GIS 22. GIS 23 ends with a smooth cooling making it impossible to clearly identify a GS 23 phase. Finally, the GIS 23-22 sequence of events shows exactly the same type of structure as the one observed over GIS 21.
Rebound-type events are not only restricted to MIS 5 as they are also occurring at the end of GIS 11, 12 and 16 (Fig. 5). GIS 13 also appears as a rebound event after the long GIS 14 without a clear GS 14. These rebound-type features are therefore recurrent over the glacial period and except for GIS 11 and GIS 12, are associated with the precursor-type events before GIS.
We also observe that rebound-type events are occurring at the end of a particularly long cooling phase during the GIS. Figure 5 highlights a linear link between the duration of GIS gradual cooling and the rebound event duration from multi-millennial (e.g. sequence of GIS 22-23 events and GIS 21) to few century timescales (e.g. GIS 11 and 16).
www.clim-past.net/6/345/2010/ rebound with associated uncertainties (R 2 =0.95). 1223  5. (a) Sub-millennial scale climatic variability characterised by GIS preceded by precursor-type peak events, a rebound structure after GIS regular cooling phase and rapid cooling phases during GIS 24 (NorthGRIP c.m., 2004). 65 • N insolation is superimposed to NorthGRIP isotopic records (Laskar et al., 2004). Dotted grey arrows indicate the gradual GIS cooling phase followed by the abrupt warming depicted as a "rebound" event. (b) Linear relationship between GIS cooling phase duration and associated duration of the rebound with associated uncertainties (R 2 =0.95).

Abrupt coolings during GIS 24
GIS 24 presents a square wave structure beginning with an abrupt temperature warming of 16 • C  and ending 3.2 kyrs later by a sudden return to stadial conditions (Fig. 6). The warm phase is punctuated by rapid cold events i.e. the slow δ 18 O ice decrease is interrupted by a first drop of δ 18 O ice by 3‰ lasting 200 yrs before a return to interstadial δ 18 O ice level. A second cooling phase occurs 500 yrs later with a 2.5‰-δ 18 O ice decrease. Finally, a stable phase is observed with a duration of 500 yrs followed by the final return to stadial conditions in less than 200 yrs. The abrupt changes in δ 18 O ice are due to changes in surface temperature as confirmed by the associated two 0.04‰ drops in δ 15 N coincident with δ 18 O ice abrupt variations (Fig. 6).
The first rapid cooling over GIS 24 appears to be accompanied by low latitudes counterparts as documented by a simultaneous drop in CH 4 concentration over 150-200 yrs. In addition, sub-millennial scale variations in δ 18 O of O 2 have been identified during GIS 24  reflecting significant changes of biosphere and hydrological cycles at these short timescales .
The very rapid climatic variability observed during the sequence GIS 23-24 with rapid events occurring in addition to the classical succession of GIS and GS, shares some similarities with the sequence of GIS/GS 15-17 that include 6 to 8 individual warming events depending on what one counts as a distinct warming event (Fig. 2).
The particularity of GIS 24 is that the first short cold spell occurs only ∼1380 yrs after the beginning of the GIS. The general picture of sub-millennial variability for this period is thus one of a cold event interrupting a long warm phase (GIS). By contrast, the later sub-millennial variability is better described in terms of brief warm events (GIS or precursor events) interrupting a long glacial phase (GS). With this view of a sharp cold spell interrupting a rather long warm phase, the sub-millennial variability of GIS 24 can only be compared with the 8.2 ka-event that occurred at the beginning of the Holocene    , 1999;Kobashi et al., 2007;Thomas et al., 2007). These two cold events occur during two different periods of transition (glacial inception for the cold events of GIS 24, end of deglaciation for the 8.2 ka-event), but both at a time when ice sheets are relatively small. The AMOC during transitional periods is expected to be subject to rapid instabilities leading to sub-millennial variability because of strong modifications of the freshwater input linked to (i) freshwater discharge (von Grafenstein et al., 1998;Clarke et al., 2004) and/or (ii) enhanced precipitation (Khodri et al., 2001) and favoured by small ice sheets (Eisenman et al., 2009).

Antarctic sub-millennial scale variability
The new detailed δ 18 O ice measurements on the EDML ice core allow the identification of an Antarctic counterpart to the stadial phase between the precursor and GIS 23, as a 1‰ δ 18 O ice variation within a few decades (Fig. 3). This AIM shows a ∼1 • C temperature increase simultaneous to the cold Greenland phase lasting ∼400 yrs. As for the rapid variability during GIS 24, Antarctic δ 18 O ice and T site reconstructions also exhibit sub-millennial counterparts. After reaching a relative temperature maximum corresponding to AIM 24, the general trend shows a regular decrease interrupted by a 1 kyr plateau that may correspond to the short cold spell occurring during GIS 24 (Fig. 3). Note that we do not identify an Antarctic counterpart to the cold phase between the precursor and GIS 21 (Fig. 3). Two hypotheses could explain such a result: (i) a lack of resolution in the EDML δ 18 O ice profile (ii) the damping of Greenland temperature signals when transferred to Antarctica through the Southern Ocean.

Millennial to sub-millennial scale GIS variability
The detailed analysis of the long GIS of MIS 5 provides evidence for sub-millennial scale variations during these phases. During GIS 21 and GIS 23, we depict a specific structure composed of a precursor-type warming event leading the GIS and a "rebound-type" abrupt event before the GIS abruptly ends. Such a structure is recurrent during MIS 3 at shorter timescales and Fig. 5 displays a linear relationship between the durations of the "rebound-type event" and of the preceding GIS regular cooling.
Inspired by the factors previously proposed for explaining the classical DO variability, we present here some of the possible mechanisms for favouring these additional submillennial scale features: (i) ice sheet size controlling iceberg discharges (MacAyeal, 1993) and the North Atlantic hydrological cycle (Eisenman et al, 2009) and (ii) 65 • N insolation affecting temperature, seasonality, hydrological cycle and ice sheet growth in the high latitudes (e.g. Gallée et al., 1992;Crucifix and Loutre, 2002;Khodri et al., 2003;. Note that these influences may also be enhanced through feedbacks. In particular, sea ice extent variations are often given as trigger (Wang and Mysak, 2006) or amplifiers (Li et al., 2005) of abrupt warming events.
We first discuss the link between the occurrence of the submillennial variability and the ice sheet volume. The length of the GIS displayed on Fig. 5 appears to be related to the mean sea level with the long GIS 23 and 21 being associated with the highest sea level while GIS 11, 12, 14 and 16 are associated with lower sea level during MIS 3 (Fig. 1). Such a link between the GIS length and sea level is expected from a simple Binge-Purge mechanism (MacAyeal, 1993): largest icesheets are expected to be easier to destabilize. However, such a Binge-Purge mechanism is unlikely to explain the existence of sub-millennial scale climatic events during sequences of events 21-24 and 15-17 since they occurs during relative ice sheet volume minima (Bintanja et al., 2005). A more plausible mechanism for these precursor events would be that the smaller ice sheets as observed during MIS 5 (equivalent to sea level of about 20 to 60 m above present sea level; Bintanja et al., 2005) are more vulnerable than large ice sheets observed during MIS 2-3-4 (sea level between 60 and 120 m above present sea level; Bintanja et al., 2005) to local radiative perturbations. If so, a strong 65 • N summer insolation would lead to intermittent freshwater outputs and trigger fast changes in the AMOC intensity. The influence of the Milankovitch insolation forcing on the sub-millennial variability can also be explored (Fig. 5). During MIS 5, the GIS 21 precursor-type event and GIS 24 are both in phase with two relative maxima in summertime insolation at 65 • N. GIS 23 precursor-type event occurs during a relatively strong 65 • N insolation and lags the preceding insolation maximum by only ∼2.5 kyrs (Fig. 5). During MIS 3, we again observe that precursor-type events GIS 14 and 16 are associated with secondary insolation maxima. On the contrary, GIS 11 and 12 are not preceded by a precursor and occur at a time without a marked anomaly in 65 • N summer insolation. Our data therefore suggest a link between high 65 • N insolation and the presence of a sub-millennial scale climatic variability in addition to the GS-GIS succession. This hypothesis also applies to the last deglaciation. Indeed, centennial-scale variations in the NorthGRIP δ 18 O ice profile are superimposed to the Bølling-Alleröd warm phase followed by the Younger-Dryas cooling (Björck et al., 1998) while the 65 • N insolation during those events is equivalent to the one observed during the sequence of events 15-17.
Finally, rebound-type events tend to be associated with long GIS intervals characterized by a slow cooling. We speculate that the rebound at the end of the GIS could be explained by an enhancement of the AMOC. Indeed, a progressive cooling could increase sea ice formation and reduce precipitation amount/runoff, increasing salinity in the North Atlantic region.

The bipolar seesaw pattern
In the above discussion, we described rapid climatic variations over Greenland. Here, we use our common dating of Antarctic and Greenland ice cores to study the north-south millennial scale variability over the whole glacial period and test the general applicability of the conceptual thermal bipolar seesaw of Stocker and Johnsen (2003) especially over the new types of rapid events identified over MIS 5.

Millennial scale variations
Synchronised EDML and NorthGRIP isotopic records emphasized the close link between the amplitude of MIS 3 AIM warming and their concurrent stadial duration in Greenland (EPICA c.m., 2006). To complete this description we have added on Fig. 7 DO/AIM 2, 21, 23 and 24 using our MIS 5 timescale (EPICA c.m., 2006;Capron et al., 2010). Finally, DO/AIM events 18, 19 and 20 have also been added despite a lack of precise north-south common age scale over this period (Fig. 8).
EPICA c.m. (2006) reveal a linear dependency between the amplitude of the AIM warming and the duration of the concurrent stadial in the north for the shorter GIS and GS events during MIS 3. We observe that the linear fit established over MIS 3 also captures the characteristics of DO/AIM events 19, 20, 23 and 24  but 49 -Evolutions of the relationship between Greenland stadial durations and AIM warming 1237 amplitudes inferred from the conceptual model for a thermal bipolar seesaw 1238 Johnsen, 2003;Equation 1) depending on (i) different characteristic timescales (500 yrs, thin 1239 curve; 1000 yrs, dotted thick curve; 1500 yrs, thick curve) and (ii) different values for TN (-1240 1/+1 amplitude, black curves; -2/+2 amplitude; yellow curves). -Evolutions of the relationship between Greenland stadial durations and AIM warming amplitudes inferred from the conceptual model for a thermal bipolar seesaw (Stocker and Johnsen, 2003;Eq. (1)) depending on (i) different characteristic timescales (500 yrs, thin curve; 1000 yrs, dotted thick curve; 1500 yrs, thick curve) and (ii) different values for T N (−1/+1 amplitude, black curves; −2/+2 amplitude; yellow curves).
does not apply for DO/AIM events 2, 18 and 21. In fact, these DO exceptions are all preceded by exceptionally long cold periods in the NorthGRIP record. Exceptionally high temperature amplitudes would be expected from the linear regression as a GS duration of 4 kyr would correspond to an AIM warming of ∼5 • C, much stronger than the observed warming amplitude of the AIM 2, 18 and 21. This shows that for extraordinarily long stadial durations the linear relationship between the stadial duration and the accompanying Antarctic warming amplitude is not longer valid. This feature is indeed expected from the bipolar seesaw concept (Stocker and Johnsen, 2003;EPICA c.m., 2006). Stocker and Johnsen (2003) predict that for long period of reduced AMOC (equivalent to GS duration in their model) a new equilibrium is reached and the Antarctic warming would eventually end. This type of situation could be relevant for the long DO/AIM 21, while DO/AIM events during MIS 3 may be too short for an equilibrium to be reached.  (Stenni et al., 2010) over AIM 2 and the sequence of events from AIM 18 to AIM 20. All are presented on the EDML1 timescale . (b) NorthGRIP δ 18 O ice over 20-28 ka: GICC05 timescale ; over 63-79 ka: ss09sea glaciological timescale (NorthGRIP c.m., 2004). Red Arrows represent warming durations of AIM 2 and AIM 18 and blue arrows represent GS 3 and GS 19 durations.
Here, we make a sensitivity test for the seesaw model in our case using the equation developed in Stocker and Johnsen (2003): Where T S (t) represents the time-dependent temperature variation in the Southern Hemisphere, τ is the characteristic timescale of the heat reservoir in the Southern Hemisphere, T N denotes the time-dependent temperature anomaly of the northern end of the bipolar seesaw. This equation predicts the southern temperature in response to climate signals in the North Atlantic region. The integral form associated with a characteristic time τ for the southern heat reservoir permits to describe the dampened temperature changes in the Southern Ocean in response to abrupt temperature changes in the North Atlantic. Following Stocker and Johnsen (2003), a value of −1 for T N stands for a GS associated with an "off" mode of the AMOC. To model the abrupt GS/GIS transition associated with resumption of the AMOC, T N changes from −1 to +1. A characteristic timescale τ of about 1000-1500 yrs has been determined to fit the Byrd temperature curve using the GRIP data as input. On Fig. 7, we display T S simulations obtained with (i) changes in T N of −1/+1 and −2/+2 and (ii) τ varying between 500 and 1500 yrs. The different results clearly illustrate a saturation level reached in the south when Greenland stadials are particularly long (more than 2000 yrs). The simulations with the largest T N amplitude (−2/+2) permit to fit the AIM amplitude/NorthGRIP stadial duration for MIS 5 events, DO/AIM 8, 12 and 19. However, it is impossible to simulate the behaviours of all events of the glacial period with a fixed amplitude for T N even with very large modifications of τ .
Our analysis suggests that larger amplitudes for T N are needed to explain the Antarctic behaviour when ice sheets are smaller. However, one should be cautious with such interpretation since it is based on the hypothesis that the Antarctic temperature reflects the change in the Southern Ocean. This may not be systematically true. Indeed, AIM events can be linked to millennial scale temperature variations in the subantarctic surface waters (Pahnke et al., 2003) and a recent study based on a marine sediment core from the Southern Ocean shows that, while the amplitude of AIM 21 is clearly larger than the amplitude of AIM 23 in EPICA ice cores, the two respective Sea Surface temperature (SST) increases have the same magnitude (Govin et al., 2009). As a consequence, this change in Antarctic behaviour in regard to rapid variability of SST can be explained by variations in the heat transmission from Southern Ocean SST signals to the interior of Antarctic from one rapid event to the other. Such variations involve many further processes e.g. ocean-atmosphere heat fluxes, polar vortex position, sea-ice formation, ice sheet altitude that are in part related to ice sheets volume (Rind et al., 2001;Velinga and Wood, 2002).
The specific behaviour observed for AIM 2 and AIM 18 is not consistent with the same thermal bipolar seesaw pattern (Fig. 8). In fact, AIM 2 and AIM 18 warming periods are shorter than the corresponding northern stadial phases, ∼700 yrs for each instead of GS durations of 4 kyrs and 5 kyrs, respectively. This highlights that Antarctic warming does not systematically start with the beginning of a GS. Climate conditions of MIS 2 and MIS 4 were particularly cold as recorded in both marine Chapman and Shackleton, 1998;de Abreu et al., 2003) and terrestrial records (e.g. Genty et al., 2003Genty et al., , 2006 and associated with vast ice sheets (Waelbroeck et al., 2002). Numerous studies have already shown that millennial scale climatic variability was reduced during MIS 2 and MIS 4 in relation with ice sheet volume (e.g. McManus et al., 1999;Schulz et al., 2002;Wang and Mysak, 2006;NorthGRIP c.m., 2004;Margari et al., 2010). Our study suggests that the bipolar seesaw was also affected during these cold periods.
Several explanations can be proposed for this particular see-saw pattern: a first possibility could be that the expansion of the Antarctic ice sheet and sea-ice during these two particular periods would increase the isolation of Antarctica and therefore decrease the heat received by the continent from the Southern Ocean (Levermann et al., 2007). A second possibility is linked to the AMOC activity. Marine records have revealed that the AMOC structure and dynamic was different over MIS 4 and end of MIS 2 compared to MIS 3 and MIS 5 in both hemispheres (e.g. Gherardi et al., 2009;Govin et al., 2009;Guihou, 2009). This particular configuration may have lead to an AMOC not strictly in an "off" mode during the whole GS. The AMOC might have been significantly reduced for the entire cold period in the north during GS 3 and 19 but could have collapsed just a few hundred years before the end of the cold phase.

Sub-millennial scale variations
In Sect. 3.3, we have shown that an Antarctic counterpart exists for the sub-millennial variability recorded in Greenland. This is especially obvious for GIS 23b. When displaying the amplitude of the Antarctic warming against the duration in the Greenland cold phase (Fig. 7), we find that it is consistent with the curve representing MIS 5 events. This result highlights that even at sub-millennial scale, the bipolar seesaw model of Stocker and Johnsen (2003) is still valid.
Using an amplitude of ±2 for T N and a characteristic timescale of 1000 yrs for the heat reservoir turned out to be the best way to describe MIS 5 rapid events. We thus apply this tuning for generating T S curves corresponding to the submillennial scale structures highlighted during MIS 5 (GIS 24 and 23).
When we use a stadial duration of ∼1150 yr with a ∼300 yr cold phase between the precursor type peak event and the main abrupt warming, the conceptual model reproduces the same singular structure in the Antarctic counterpart as observed in the data (Fig. 9).
We then construct a time-series of T N corresponding to GIS 24 characterised by an abrupt cooling phase lasting 200 yr (Fig. 10). We observe a plateau interrupting the regular cooling phase after AIM 24 as depicted by the  (Stocker and Johnsen, 2003) for GIS 23 associated with the precursor event. Different configurations of the northern perturbation (T N , superimposed to NorthGRIP δ 18 O ice data) are used to simulate the response of the Southern Hemisphere (T S ). One configuration (dark grey curve) represents the evolution of T S in response to T N that corresponds to an "onoff" signal with 1150 years off (amplitude −2) and 7340 years on (amplitude +2). Two additional configurations (blue and red curves) are superimposed to illustrate the thermal bipolar seesaw pattern at a sub-millennial scale. T N is structured as an "on-off" signal with 1150 years off (amplitude −2), 240 years on (red curve, amplitude +2; blue curve, amplitude +1), 300 years off (amplitude

Fig. 10. (a)
North-south time-series generated through the conceptual thermal bipolar seesaw model (Stocker and Johnsen, 2003) for GIS 24. Different configurations of the northern perturbation (T N , superimposed to NorthGRIP δ 18 O ice data) are used to simulate the response of the Southern Hemisphere (T S ). One configuration (dark grey curve) represents the response of the Southern Hemisphere (T S ) to T N that corresponds to an "on-off" signal with 1400 years off (amplitude −2), 3330 years on (amplitude +2) and 700 years off (amplitude −2). Two other configurations (blue and red curves) are superimposed to illustrate the thermal bipolar seesaw pattern at a sub-millennial scale. T N is structured as an "on-off" signal with 1400 years off (amplitude −2), 1590 years on (amplitude +2), 200 years off (red curve, amplitude −1; blue curve amplitude 0), 1530 years on (amplitude +2) and 700 years off (amplitude −2). (b) EDC δ 18 O ice  and T site reconstruction (Stenni et al., 2010). (c) EDML δ 18 O ice (EPICA c.m., 2006) and T site reconstruction (Stenni et al., 2010). temperature reconstruction of both EPICA cores. Results obtained on both events 23 and 24 emphasize the ability of the model tuned on MIS 5 to explain the sub-millennial scale variability depicted in Antarctic isotopic records.

Summary and perspectives
In this paper, we present the most recent and accurate Greenland-Antarctica common dating over the last 123 ka using the NorthGRIP and EPICA ice cores. We used new and published measurements of air isotopic composition in the NorthGRIP ice core to compare the local amplitudes of temperature changes for GIS of MIS 5 and MIS 3. A study of the δ 18 O ice /temperature slope at the onset of each rapid event shows a strong variability from one GIS event to another but no systematic difference between MIS 3 and MIS 5 events. For Antarctica, we have combined new and published water isotope records to present detailed temperature reconstructions of Antarctic temperature based on EPICA isotopic records.
NorthGRIP records enable us to depict the sub-millennial scale variability during the GIS of MIS 5 and thus, to highlight new type of features (GIS 21, 23) observed also during MIS 3 (GIS 11,12,(13)(14)16). These new patterns appear as (i) precursor-type events prior to the onset of GIS (ii) rebound events at the end of GIS and (iii) centennial-scale cooling during the long and warm GIS 24. In addition to the internal forcing of ice-sheets on the climatic evolution during these events, we have proposed the external influence of the summertime insolation at 65 • N. Disentangling the main processes leading to these sub-millennial scale structures (ice-sheet, insolation, sea-ice, and hydrological cycle forcing) will require dedicated modelling studies. Through our results, we assume that orbital-scale variations play a role in rapid climate change but, also, the millennial-scale variability may hold clues to the long term climatic changes (i.e. Weirauch et al., 2008;Wolff et al., 2009b).
Comparing Antarctic and Greenland behaviour over the succession of AIM/DO back to MIS 5 provides a more complete description of the bipolar seesaw pattern. As expect from the bipolar seesaw concept, a linear relationship between AIM amplitude and preceding GS duration only holds for shorter events, while for extraordinary long GS a new heat flux equilibrium between the Northern and Southern Hemisphere is obtained (EPICA c.m., 2006, Stocker andJohnsen, 2003) and the Southern Ocean warming ceases. The conceptual model of Stocker and Johnsen (2003) for a thermal bipolar seesaw is able to represent most of the variability of the north-south relationship depicted in Greenland and Antarctic isotopic records, even at sub-millennial timescale. However, it is not able to depict the delay of Antarctic warming after the beginning of the GS during the periods associated with large ice sheets (i.e. during MIS 2 and the end of MIS 4). It shows that Greenland ice core temperature proxy records cannot be taken as direct proxy for AMOC changes as suggested from the conceptual model.
To go beyond our description and the conceptual model of Stocker and Johnsen (2003), the new types of DO events identified during MIS 5 should be studied with more complex models (e.g. Ganopolski and Rahmstorf, 2001;Knutti et al., 2004). This would allow quantification of the influence of insolation, ice-sheet volume, sea-ice and hydrological cycle on sub-millennial-scale variability (precursor and rebound events). This should provide also a better understanding of the response of Antarctica to these types of events.