Hindcasting the continuum of Dansgaard–Oeschger variability: mechanisms, patterns and timing

. Millennial-scale variability associated with Dansgaard–Oeschger events is arguably one of the most puzzling climate phenomena ever discovered in paleoclimate archives. Here, we set out to elucidate the underlying dynamics by conducting a transient global hindcast simulation with a 3-D intermediate complexity earth system model covering the period 50 to 30 ka BP. The model is forced by time-varying external boundary conditions (greenhouse gases, orbital forcing, and ice-sheet orography and albedo) and anomalous North Atlantic freshwater ﬂuxes, which mimic the effects of changing northern hemispheric ice volume on millennial timescales. Together these forcings generate a realistic global climate trajectory, as demonstrated by an extensive model/paleo data comparison. Our results are consistent with the idea that variations in ice-sheet calving and subsequent changes of the Atlantic Meridional Overturning Circulation were the main drivers for the continuum of glacial millennial-scale variability seen in paleorecords across the globe.

An important element of DO events is the corresponding variability of ice-rafted debris (an indicator for iceberg surges) (Bond and Lotti, 1995;Sarnthein et al., 2001), illustrated here in a high-resolution ice-rafted debris (IRD) composite record from the Irminger and Iceland seas cores SO82-5 ( van Kreveld et al., 2000) and PS2644 (Voelker et al., 2000) (Fig. 1, upper panel, orange line). We observe that within the age uncertainties, all DO stadials between 30-50 ka BP were accompanied by iceberg surges, which originated from the adjacent northern hemispheric ice sheets. According to these marine records and other data sets from the northern North Atlantic (e.g., Elliot et al., 1998;Grousset et al., 2001;Elliot et al., 2002), the iceberg calving increased during interstadial periods and peaked at the end of the stadials, after which it decayed rapidly. This sequence of events is consistent with the notion of a freshwater-driven throttling of oceanic convection (Sarnthein et al., 2001), meridional mass and heat transport and subsequent sea-ice expansion, which caused the gradual cooling during interstadials. Moreover, possible feedbacks between the Atlantic meridional overturning circulation (AMOC) and ice sheets may have further modulated the evolution of iceberg and freshwater discharges into the  van Kreveld et al., 2000) and PS2644  from the northern North Atlantic (see map inlay). Integral of freshwater flux forcing used in the LOVECLIM MIS3 hindcast simulation (blue) (see Fig. 2, upper panel). Lower panel: GISP2 reconstructed central Greenland temperatures (Alley, 2000). Greenland interstadials (GI) are highlighted by red labels. The main Heinrich events are represented by gray bars and the light green bar marks are the Greenland stadial C7 just prior to GI7. The age shift between the GISP2 record (Alley, 2000) and the NGRIP record on the GICC05 timescale Svensson et al., 2006) is determined. To project the paleorecords onto the common GICC05 age model, we apply this shift to the GISP2 record and to IRD records of cores SO82-5 and PS2644, whose original age models were partly based on correlations with GISP2. Subsequently the IRD composite (orange) was calculated. Timmermann et al., 2003;Shaffer et al., 2004;Alvarez-Solas et al., 2010;Marcott et al., 2011;Alvarez-Solas et al., 2013).
Based on their very distinctive sedimentological characteristics, Heinrich stadials (HS) have often been regarded as dynamically different from other DO stadial-interstadial transitions. However, given the fact that both Heinrich and DO stadials were accompanied by (i) large-scale oceanic changes ( Fig. 2), (ii) IRD layers ( Fig. 1) and (iii) similar global teleconnections, we hypothesize here that Heinrich and DO stadials are part of a continuum of variability that is generated through ice-sheet-driven AMOC changes.
In this paper we set out to simulate the time evolution of DO and Heinrich variability for the period 50-30 ka BP using an intermediate complexity global climate model. We will compare the simulated variability with high-resolution paleoclimate reconstructions. The model simulation is based on the underlying assumption that the continuum of MIS3 climate variability on centennial to millennial timescales can be generated by a suitable North Atlantic freshwater forcing and the associated AMOC response.
The paper is organized as follows: in Sect. 2 the model and experimental setup are described. In Sect. 3 we discuss the patterns of variability associated with DO cycles, the abruptness of stadial-interstadial transitions as well as the timing of Heinrich stadials. We also derive a common age scale that allows for a better comparison between paleoproxy records and model simulations. The paper concludes with a synthesis and discussion of the main results.

Model and experimental setup
One of the key goals of our study is to simulate the sequence of millennial-scale events during the period 50-30 ka BP and to determine the corresponding global teleconnections. For this task we have chosen the intermediate complexity earth system model LOVECLIM (Timm and Timmermann, 2007;Menviel et al., 2008;Timmermann et al., 2009b;Goosse et al., 2010). The ocean component of LOVE-CLIM (CLIO) consists of a free-surface primitive equation model with a horizontal resolution of 3 • longitude, 3 • latitude, and 20 depth layers. The 3-D atmospheric component (ECBilt) is a spectral T21, three-level model based on quasigeostrophic equations of motion and ageostrophic corrections. LOVECLIM also includes a dynamic-thermodynamic sea-ice model, a land surface scheme, a dynamic global vegetation model (VECODE, Brovkin et al., 1997) and a marine carbon cycle model (LOCH, Menviel et al., 2008;Mouchet, 2011).
Initial conditions for the transient run were obtained by conducting an equilibrium spin-up simulation using an atmospheric CO 2 content of 207.5 ppmv, orbital forcing for the time 50 ka BP and an estimate of the 50 ka BP ice-sheet orography and albedo which were obtained from a 130 ka off-line ice-sheet model simulation (Abe-Ouchi et al., 2007). In the subsequent transient run greenhouse gases, orbital and ice-sheet forcing were updated continuously following the methodology of Timm et al. (2008). Note that our coupled model does not include an interactive ice sheet. Therefore, freshwater withholding from the ocean during phases of icesheet growth and freshwater release into the ocean as a result of ice-sheet calving and ablation are not explicitly captured. To mimic the time evolution of these terms and their effect on the ocean circulation, we apply an anomalous North Atlantic freshwater forcing F (t) to the North Atlantic region 55-10 • W, 50-65 • N. Negative forcing anomalies can be interpreted as periods of ice-sheet growth and excess evaporation over precipitation, whereas positive freshwater anomalies represent times of negative net mass balance of the northern hemispheric ice sheet, associated for instance with massive iceberg calving events or surface ablation.
The freshwater forcing time series F (t) is obtained through an iterative procedure, that optimizes the anomalous freshwater flux such that the simulated temperature anomalies T s (t) in the eastern subtropical North Atlantic best match the target alkenone-based SST anomalies T r (t) reconstructed from the Iberian margin core MD01-2443 (Martrat et al., 2007) (Fig. 2, lower panel, orange line). Starting from an initial guess of the freshwater forcing F i (t) = −αT r (t), a series of about j < 5 experiments was conducted every 1000 yr using additional freshwater flux perturbations δF (j ) (t). In each of these 1000 yr long chunks, the freshwater forcing scenario (j ) was selected with F (j ) (t) = F i (t) + δF (j ) (t) that minimized the cost function within this window [t, t + τ ] with τ = 1000 yr. The simulated temperature evolution for T s (t ) in the integral is a function of the applied freshwater forcing F (j ) (t ). The resulting concatenated freshwater forcing time series F (t) is shown in Fig. 2 (upper panel). Similar to data-assimilation methods that adjust parameters and/or dynamical variables to reduce the mismatch between observations and models, F (t) has the sole purpose to force LOVECLIM into a realistic trajectory with respect to millennial-scale subtropical North Atlantic SST anomalies during MIS3. We did not choose the Greenland temperature reconstruction as an optimization target, because it shows only very weak differences between DO and Heinrich stadials, in contrast to the North Atlantic SST reconstructions.

Freshwater forcing
The applied North Atlantic freshwater forcing F (t) captures the dominant meltwater pulses associated with Heinrich stadials (Fig. 2). It compares well with a recent freshwater forcing estimate (Jackson et al., 2010) obtained with a North Atlantic box model through Bayesian inversion methods 1 . In both cases stadial-interstadial transitions are triggered by negative forcing anomalies, which increase North Atlantic surface densities and subsequently strengthen the AMOC (Fig. 2, middle panel). Negative freshwater forcing can be interpreted to represent a positive ice-sheet mass balance, which in our modeling framework mimics a reduction of the continental runoff as well as excess evaporation over precipitation over the North Atlantic region. As an independent validation of our freshwater forcing, we compare the time-integral of F (t), which represents the corresponding global sea level changes, with the composite IRD records from the Nordic Seas cores PS2644 and SO82-5 ( Fig. 1) on the GICC05 timescale (see Sect. 3.6 for more details on the synchronization of the cores and the model results). The rationale of this comparison is that high values of IRD correspond to additional freshwater discharge and sea level rise. Furthermore, rising sea level can amplify iceberg calving through ice-shelf instabilities. Except for the simulated sea level rise associated with Heinrich event 5, which is not captured in these eastern North Atlantic paleorecords, we find a relatively good match between model simulation and reconstruction, thus supporting the realism of the applied freshwater forcing. It should be noted here that the simulated DOrelated sea level changes are a factor 2-3 smaller than those reconstructed from the Red Sea during this time (Siddall et al., 2003).

AMOC response
As a result of the applied anomalous North Atlantic freshwater fluxes, the AMOC weakens and strengthens on millennial timescales. Heinrich stadials correspond to a complete shutdown of the AMOC, whereas DO stadials are associated with a 50 % weakening of the AMOC, relative to the interstadial periods (Fig. 2). The resulting AMOC time series compares reasonably well with a reconstruction of Atlantic bottom currents obtained from mass-normalized anhysteretic remanent magnetization (ARM) data (Kissel et al., 2008) (Fig. 2, middle panel), even though the model and paleoceanography time series are based on different underlying age models (a more detailed discussion of age-scale uncertainties is provided is Sect. 3.6).

Temperature response
The excellent agreement between simulated and reconstructed SST anomalies in the Iberian margin area ( Fig. 2 lower panel) is to be expected, because the latter has been used as the target for the optimization of the freshwater fluxes. One important finding is that the temperature drop in the northeast Atlantic around 36 ka BP (referred to as C7, adopting the Chapman and Shackleton, 1999 terminology) can be obtained in our simulation only by a complete shutdown of the AMOC, which is induced by a prolonged freshwater flux of ∼ 0.1 Sv. This is consistent with the presence of an IRD pulse in the Greenland and Irminger seas (see Fig. 1), a drop in sea-surface salinity in SO82-5 ( van Kreveld et al., 2000) and changes in benthic δ 18 O (Margari et al., 2010). Whereas the IRD pulse is well pronounced in the composite IRD time series (Fig. 1), as well as in records from the Irminger basin (SU 90-24) (Elliot et al., 2002), the southern Gardar Drift (JPC-13) (Hodell et al., 2010) and from the Iberian margin (MD95-2040) (Schönfeld et al., 2003); it appears to be absent in other southwestern Atlantic IRD records (e.g., Grousset et al., 1993;Rashid et al., 2003;Nave et al., 2007).
We move on to a more detailed comparison with other temperature reconstructions from the North Atlantic and Mediterranean realm. Figure 3 shows the comparison between simulated surface temperatures in Greenland and NGRIP temperature reconstructions (Huber et al., 2006) on the SS09 timescale. In accordance with paleodata, the simulated Heinrich stadials in Greenland attain very similar minimum temperatures than the DO stadials. This behavior in Greenland is quite distinct from SST reconstructions ( d. c. tracks the underlying AMOC variability much more accurately. Moving further into the eastern Mediterranean region, we find a reasonable qualitative match between the simulated temperature anomalies in Turkey and the δ 18 O record from independently dated speleothems from Sofular Cave (Fleitmann et al., 2009), which capture a combined temperature/hydroclimate signal.
A more comprehensive spatial view of the simulated DO/Heinrich dynamics is obtained from an EOF analysis of global surface air temperatures (Fig. 4a). The dominant EOF mode is characterized by a meridional temperature seesaw in accordance with numerous other modeling studies Timmermann et al., 2009a;Kageyama et al., 2013) and paleoclimate data sets (Blunier et al., 1998;Barker et al., 2009;Stenni et al., 2011). Interstadial conditions are characterized by northern hemispheric warming with strongest amplitudes over the Greenland, Iceland, and Norway seas and the Arctic Ocean. The warming extends into North Africa, Asia and the western North Pacific. This warming pattern is in general agreement with pollenderived temperature reconstructions for GI8 (38 ka BP) and GI6 (33 ka BP) . The corresponding principal component time series (Fig. 4c) clearly features the enhanced cooling during massive Heinrich stadials (HS5 and HS4) and the C7 stadial, in contrast to the weaker cooling associated with DO stadials. Simulated southern hemispheric cooling during interstadials is consistent with the presence of a bipolar temperature seesaw Stocker and Johnsen, 2003).

Hydroclimate response
As a result of the very strong North Atlantic cooling during Heinrich stadials (HS5, HS4) and during the C7 stadial (Fig. 3), northern hemispheric trade winds intensify by up to 60 %, which leads to a southward shift of the Intertropical Convergence Zones, extending from South America, into the tropical Atlantic, equatorial Africa and the Indian Ocean. This is illustrated by the EOF analysis of simulated precipitation in Fig. 4b and d.
The lower amplitude cooling during DO stadials weakens the trade winds by only 30 %. The corresponding southward shift of the tropical rainbands is less pronounced than for Heinrich stadials as shown by the leading principal component of the rainfall EOF analysis (Fig. 4d). In spite of a high correlation (0.92) between the principal components of temperature and rainfall, there are some notable differences. The precipitation mode exhibits a more pronounced two-step structure for the interstadial DO12 (around 47-46 ka) and a stronger difference between Heinrich stadials and DO stadials (Fig. 4d) than the temperature EOF mode (Fig. 4c). Qualitatively the patterns of simulated temperature and rainfall changes agree with those obtained from coupled general circulation models subjected to North Atlantic freshwater perturbations (Broccoli et al., 2006;Timmermann et al., 2007;Kageyama et al., 2013).
Comparing the simulated northern hemispheric rainfall changes on a regional scale with hydroclimate reconstructions for the Mediterranean region (Sánchez-Goñi et al., 2002), the Cariaco Basin (Deplazes et al., 2013), the Arabian Sea (Deplazes et al., 2013), eastern China (Wang et al., 2001) (Fig. 5) and Central America (Hodell et al., 2008) (Fig. 6 and data with stadial (interstadial) conditions corresponding to increased aridity (pluvials). While the pollen-derived precipitation estimate (Sánchez-Goñi et al., 2002) and the model output suggest much drier conditions over the Iberian region during HS5, the chronology of the proxy record is based on a graphic correlation with the GISP2 ice core, which places the HS5 much later than the new GICC05 chronology (Obrochta et al., 2014). It should be noted that reflectance and magnetic susceptibility are indirect hydroclimate proxies and thus cannot give quantitative estimates. In addition, speleothem δ 18 O can be potentially affected by other processes such as changes in temperature, soil evaporation and the water vapor sources. Age model uncertainties associated with the Arabian Sea record (Deplazes et al., 2013) could preclude any conclusions with respect to synchronicity with North Atlantic stadials. However, the high level of correspondence between our simulated precipitation changes and the Arabian Sea reflectance record indicates that North Atlantic stadials lead to drier synchronous conditions over the Arabian Sea. The reverse pattern can be found for southern hemispheric hydroclimate proxies in Brazil (Wang et al., 2007), Peru (Kanner et al., 2013) and Ecuador (Mosblech et al., 2012) as well as for simulated rainfall changes (Fig. 6, lower three panels). In Sect. 3.6 we will try to reconcile some age-model discrepancies by projecting model and proxy data onto the common NGRIP GICC05 timescales Svensson et al., 2006).

Abruptness of stadial-interstadial transitions
To determine the response time of various climate variables and ocean transport indicators to freshwater changes, we calculate a composite (Fig. 7, black line) based on several stadial-interstadial DO transitions by aligning the model data relative to the maximum temperature derivative of the simulated Greenland temperature (47.25,43.9,41.8,38.3,35.3,33.6,32.34 ka BP). According to this analysis we find that www.clim-past.net/10/63/2014/ the averaged DO transition takes place within 150 to 200 yr for all the climate variables in Fig. 7. However, for Greenland temperatures (Fig. 7a) and northern hemispheric sea-ice area (Fig. 7e), we see a considerable acceleration and an associated increase of abruptness 100 yr into the transition. This is in agreement with previous estimates of the abruptness of DO stadial-interstadial transitions in Greenland ice cores (∼ 125 yr) (Capron et al., 2010), although much higher rates were reported in atmospheric circulation proxies (Steffensen et al., 2008).
As already demonstrated in Figs. 3 and 5 rainfall changes in the Cariaco Basin area and the Arabian Sea clearly track millennial-scale SST variations in the North Atlantic region. This is further supported by the composite analysis which reveals a very similar time evolution of these variables for the averaged DO stadial-interstadial transition. Rainfall changes in eastern China are less well pronounced owing to a much larger level of simulated rainfall variability that is unrelated to DO cycles.
According to Fig. 7h, changes in the barotropic transport across the Indonesian archipelago occur almost in unison with the AMOC. This surprisingly fast adjustment can be attained by two processes: (i) wind changes in the Pacific (Timmermann et al., 2005b), and (ii) fast oceanic adjustment processes involving wave propagation from the Atlantic into the Indian and Pacific oceans, as discussed in Timmermann et al. (2005a). The former can modulate the Indonesian Throughflow via the island rule (Godfrey, 1989), whereas the latter would have to change the joint effect of baroclinicity and relief (JEBAR) term in the barotropic transport equation (Sarkisyan and Ivanov, 1971;Cane et al., 1998). Irrespective of the relative magnitudes of these terms, our analysis clearly documents that the DO variability has far-reaching fast oceanic impacts that extend also into the other ocean basins.

Common age scale
To better compare the paleorecords in Figs. 2, 3, 5 and the simulated climate variables, we make an attempt to bring the time series all onto the same age scale. We have chosen the Greenland Ice Core Chronology 2005 (GICC05) Svensson et al., 2006) as the common age model. To project the model simulation (30-50 ka BP) onto this age model, we compare the simulated Greenland temperature with the NGRIP δ 18 O  and identify in a 1800 yr sliding window the lag at which the lag correlation attains its maximum value. Here we allow maximum lags of ±750 yr. To avoid large local discontinuities or reversals in the age model projection, we subsequently filtered the resulting time series of age adjustments for the model simulation using a 500 yr running mean. The resulting Greenland temperature-based age shift is then applied to all other model variables, thus keeping, at least to first order, the leadlag structure within the model intact. We also project the   Svensson et al., 2006) of Nordic Seas salinity from core S082-5 ( van Kreveld et al., 2000) and simulated surface salinity anomalies at this location; simulated maximum meridional overturning circulation in the North Atlantic (Sv) compared to North Atlantic marine sediment cores ARM data (Kissel et al., 2008); simulated SST anomalies off the Iberian margin (15-8 • W, 37-43 • N) compared to alkenone-based SST anomalies from marine sediment core MD01-2444 (Martrat et al., 2007). Model results are in black and paleoproxy records in orange.
sea-surface salinity data from SO82-5 ( van Kreveld et al., 2000), the ARM data and the Iberian margin SST (Fig. 8) as well as the Cariaco and Arabian seas color records (Fig. 9), onto the GICC05 timescale. Here we assume that at least in a 1800 yr sliding window, the proxy data varies in synchrony with the Greenland temperature record at zero lag. This assumption is well justified by the model results that show maximum correlation of 0.92 between the principal components of the leading EOFs of temperature ( Fig. 4a and c) and precipitation ( Fig. 4b and d) at zero lag. Furthermore, our model-based composite analysis of DO stadial-interstadial transitions supports the notion of near synchronicity (within ±100 yr) of the physical variables under consideration. Having synchronized the model and paleoproxy data with the NGRIP δ 18 O record on GICC05, we find a much better agreement between model and paleoproxy records, particularly for the period 50-40 ka BP. The ARM data now nicely features marked AMOC weakening during HS5, HS4 and C7 as well as during most of the DO stadials (Fig. 8). In addition,  Table showing the timing of stadials HS5, HS4, C7/HS3.2 and HS3 as recorded in high-resolution well-dated paleorecords. We calculate the mean (ka BP) and one standard deviation (ka) of the timing of each stadial as well as the duration (yr) of each stadial and its standard deviation (yr). In the Pacupahuain record, the two peaks around HS5 were treated as a single event. A similar method was used for HS3 in the Hulu Cave record.  the precipitation records from the Cariaco Basin and the Arabian Sea are now in better agreement with the model, particularly for HS5 (Fig. 9). We conclude that if forced with a freshwater forcing that leads to simulated salinity anomalies which closely resemble (within the dating uncertainties) paleosalinity reconstructions from the Nordic Seas (Fig. 8), the LOVECLIM model hindcast captures the dominant modes of Heinrich and DO variability found in paleoreconstructions. The model results further support that freshwater forcing triggered changes of the AMOC and North Atlantic SSTs, which subsequently caused the observed hydroclimate shifts across both Hemispheres. This confirms our initial hypothesis that ice-sheetdriven AMOC variations played a crucial role in generating the continuum of millennial-scale DO/Heinrich variability in the North Atlantic during MIS3 (see also Sarnthein et al., 2001). Potential feedbacks of AMOC variability on the mass balance of the major ice sheets will be discussed in Sect. 4.

Timing and duration of Heinrich stadials
Here we will take the opportunity to revisit the timing of some important climate events during MIS3. We will focus in particular on Heinrich stadials (HS5, HS4 and HS3) and stadial C7 (see Fig. 1) and use high-resolution paleoclimate reconstructions with independent age control that capture Heinrich and DO variability. The timing from NGRIP on GICC05 , Sofular Cave, Turkey (Fleitmann et al., 2009), Pacupahuain Cave, Peru (Kanner et al., 2013), Santiago Cave, Ecuador (Mosblech et al., 2012) and Hulu Cave (Wang et al., 2001) is summarized in Ta -Solas et al., 2010). Furthermore, sea level changes generated by one ice sheet can trigger ice-shelf instabilities in another ice sheet and subsequent accelerated flow and iceberg calving. Once the ice sheets reach a new mass balance, the freshwater input into the North Atlantic ceases and the AMOC starts its recovery thus initiating a stadial-interstadial DO transition. HS4,.8 ka BP, respectively, which agrees well with our model simulation (Figs. 2 and 3). These ranges also agree well with estimates from North Atlantic marine sediment cores for the timing of HS5 and HS4 (50-47 and 40.2-38.3 ka BP, respectively) (Sánchez-Goñi and . The C7 stadial was accompanied by a considerable IRD pulse in the northeastern North Atlantic, as seen for instance in the sediment cores PS2644 , SO82-5 ( van Kreveld et al., 2000), JPC-13 (Hodell et al., 2010), MD95-2040(Schönfeld et al., 2003), NA87-22 and SU 90-24 (Elliot et al., 2002. Furthermore, we find strongly reduced surface temperatures in the Atlantic and widespread northern hemispheric aridity (Figs. 2, 3 and 4). In our model simulation and in the paleoproxy data, the C7 stadial shares many common characteristics with the typical response for Heinrich stadials 3-5. In fact in the EOF analysis of the model simulation (Fig. 4), this period is basically indistinguishable from the other prominent Heinrich stadials, both in terms of temperature and rainfall. We therefore propose to introduce the term "Heinrich stadial 3.2" (HS3.2), using the same nomenclature introduced for Heinrich stadial 5.2 (Sarnthein et al., 2001).
Paleorecords as well as the model results display little coherency regarding the amplitude and the timing of HS3. It was suggested (Elliot et al., 1998;Scnoeckx et al., 1999;Grousset et al., 2000) that in contrast to other Heinrich events, HS3 may have originated from the Fennoscandian ice sheet. It is thus possible that iceberg discharges during Heinrich event 3 had a different impact on the AMOC than for other Heinrich events.

Conclusions
Here we presented a new transient model simulation that covers the period 30-50 ka BP. This climate model hindcast experiment was designed in such a way that freshwater forcing between 55-10 • W, 50-65 • N generates AMOC changes and subsequently northeastern Atlantic temperature anomalies that are in agreement with alkenone-based temperature reconstructions from the Iberian margin area. With this weak constraint on model/data agreement, we were able to independently evaluate the model performance with numerous other high-resolution climate proxies from both hemispheres. A marked weakening of the AMOC reduces the oceanic and atmospheric poleward heat transport thus leading to a strong cooling over the North Atlantic region (Kageyama et al., 2013). In our model the cooling is centered on Scandinavia, extends over Greenland and northern Europe and is also simulated over southern Europe, North Africa and Asia. The cooling is the strongest at high latitudes due to sea-ice albedo feedbacks. Such temperature changes lead to a stronger surface temperature gradient over the North Atlantic and therefore to a strengthening of the North Easterly trades. The cooler conditions over the North Atlantic and stronger trades induce a southward shift of the ITCZ over the Atlantic region with drier conditions simulated over Europe, the northern part of South America, North Africa and the Middle East.
The resulting high level of agreement between the model and paleoproxy records provides strong support for our initial hypothesis: namely that Heinrich and DO variability during MIS3 were caused by ice-sheet-driven changes in the strength of the AMOC. The relationship between stadials and IRD records from the North Atlantic region further indicates that northern hemispheric ice-sheet calving and freshwater discharges play a major role in disrupting the AMOC during MIS3 (Fig. 1). Figure 10 summarizes a scenario explaining DO/Heinrich variability, which follows some elements of Sarnthein et al. (2001). We suggest that ice-sheet instabilities of low amplitude, originating mostly from the Eurasian, Iceland and Greenland ice sheets (Bond and Lotti, 1995), are the main driver of the fast DO variability. The corresponding low amplitude freshwater flux perturbations triggered a weakening of the AMOC, but not a complete collapse. We further acknowledge the possibility that DO climate variability may have caused ice-sheet mass imbalances and calving from circum-Atlantic ice sheets (Bond and Lotti, 1995;Grousset et al., 2000Grousset et al., , 2001Marshall and Koutnik, 2006), thus contributing to the DO-synchronized delivery of IRD into the North Atlantic. In contrast, instabilities from the Laurentide ice sheet occurred less frequently but were associated with much larger iceberg and freshwater discharges, leading to complete AMOC shutdown and larger SST and hydroclimate changes in the North Atlantic realm and beyond. Changes in sea level during Heinrich events  and subsurface temperatures (Shaffer et al., 2004;Mignot et al., 2007;Alvarez-Solas et al., 2010;Marcott et al., 2011;Alvarez-Solas et al., 2013) (Fig. 10) may have subsequently triggered marine-ice-sheet instabilities, thus increasing the initial freshwater discharge. Such processes may have played a key role in synchronizing ice-sheet dynamics in the Northern Hemisphere and in prolonging ice-sheet instabilities during Heinrich events. Once the ice sheets reach a new mass balance, the freshwater input into the North Atlantic ceases, salinity increases rapidly to the more saline Arctic/North Atlantic glacial background state and the AMOC starts its recovery thus initiating a stadial-interstadial DO transition.
A more detailed view of the underlying mechanisms is provided in Fig. 11, which shows the time evolution of the composite IRD record (Fig. 1) from cores SO82-5 and PS2644 (orange, upper curve), the salinity reconstruction (blue, upper curve) from Irminger Sea core SO82-5 ( van Kreveld et al., 2000), the North Atlantic ARM data (red, middle curve) (Kissel et al., 2008) as a proxy for changes in AMOC strength, summer SST variations from the Irminger Sea (cyan, middle curve) (van Kreveld et al., 2000) and the GISP2 ice-core temperature reconstruction (black, lower curve) (Alley, 2000). All data were interpolated onto the GICC05 timescales (see caption to Fig. 1 and Sect. 3.6 for more details). Here we begin with the high IRD values during HS4, low salinities and cold conditions in the Irminger Sea and a weak AMOC (40-39 ka BP). Around ∼ 39 ka BP the strong freshwater forcing vanishes abruptly. Concomitantly, Fig. 11. Top panel: composite IRD record (orange) obtained from the average of the normalized northern North Atlantic IRD records SO82-5 ( van Kreveld et al., 2000) and PS2644  (see Fig. 1) and Irminger Sea sea-surface salinities (blue) from SO82-5 ( van Kreveld et al., 2000); middle panel: North Atlantic marine sediment cores ARM data (red) (Kissel et al., 2008) and Summer SST anomalies (cyan) from marine sediment core SO82-5 ( van Kreveld et al., 2000); bottom panel: GISP2 reconstructed central Greenland temperatures (black) (Alley, 2000). All data are displayed on the GICC05 timescale (see Fig. 1). sea-surface salinity increases thus initiating the AMOC recovery. The AMOC strengthening leads to Greenland and North Atlantic warming as well northern North Atlantic seaice retreat. Greenland interstadial 8 (GI8) is characterized by a very warm initial period which lasts for about 100-200 yr. IRD is at its minimum, Nordic Seas surface salinity is high and so is the strength of the AMOC (Fig. 8). We consider this period of minimum ice-sheet calving a period of positive northern hemispheric ice-sheet mass balance (i.e. growth). Around 37.5 ka BP calving resumes and increases until 37 ka BP. This evolution is briefly interrupted for about 100 yr between 36.9-36.8 ka BP, before a period of rapid iceberg surging and an associated salinity decrease leads into the stadial cooling phase during HS3.2. The stadial iceberg surging period lasts for 1 ka and comes to an end when the ice-sheet calving has exhausted itself. This is the initiation of GI7. The scenario outlined here for a set of paleoproxy data sets is entirely consistent with the modelingbased evidence from Figs. 8 and 9.
Our paper further highlights the different response characteristics of various climate variables to AMOC changes (Fig. 7). The extraordinary abruptness of Greenland temperature changes during the DO stadial-interstadial transition was identified as a regional phenomenon, which is likely induced by sea-ice feedbacks (Li et al., 2010;Deplazes et al., 2013). Given the fact that the North Atlantic temperature and AMOC composite shown in Fig. 7 have already reached 2/3 of their full DO amplitude at zero lag whereas Greenland temperatures have only attained about 50 %, it may appear as if the Greenland record is lagging the other variables. This is merely a reflection of the nonlinearity of the Greenland temperature response. It should be noted here that this feature may impact the synchronization of high-resolution proxy time series with Greenland climate reconstructions. According to our model and data-based evidence, we conclude that ice-sheet/freshwater-driven AMOC variations and local feedbacks determined the timing and abruptness of DO events as well as their global teleconnections.