# Dansgaard-Oeschger and Heinrich event temperature anomalies in the North Atlantic set by sea ice, frontal position and thermocline structure

J. B. Pedro, C. Andersson, G. Vettoretti, A. H. L. Voelker, Claire Waelbroeck, T. M. Dokken, M. F. Jensen, S. O. Rasmussen, E. G. Sessford, M. Jochum, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

J. B. Pedro, C. Andersson, G. Vettoretti, A. H. L. Voelker, Claire Waelbroeck, et al. Dansgaard-Oeschger and Heinrich event temperature anomalies in the North Atlantic set by sea ice, frontal position and thermocline structure. Quaternary Science Reviews, 2022, 289, 10.1016/j.quascirev.2022.107599. insu-03746437

HAL Id: insu-03746437

https://insu.hal.science/insu-03746437

Submitted on 5 Aug 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ELSEVIER

Contents lists available at [ScienceDirect](https://www.sciencedirect.com)

Quaternary Science Reviews

journal homepage: www.elsevier.com/locate/quascirev

<figure><img src="image_4.png" /></figure>

# Dansgaard-Oeschger and Heinrich event temperature anomalies in the North Atlantic set by sea ice, frontal position and thermocline structure

<figure><img src="image_5.png" /></figure>

Check for
updates

J.B. Pedro<sup>a, b, c, *</sup>, C. Andersson<sup>d</sup>, G. Vettoretti<sup>c</sup>, A.H.L. Voelker<sup>e, f</sup>, C. Waelbroeck<sup>g</sup>, T.M. Dokken<sup>h, i</sup>, M.F. Jensen<sup>h</sup>, S.O. Rasmussen<sup>c</sup>, E.G. Sessford<sup>h</sup>, M. Jochum<sup>c</sup>, K.H. Nisancioglu<sup>h, i, j</sup>

## ARTICLE INFO

**Article history:**

Received 19 December 2021

Received in revised form

26 May 2022

Accepted 3 June 2022

Available online 23 June 2022

Handling Editor: I Hendy

**Keywords:**

Dansgaard-Oeschger event

Heinrich event

abrupt climate change

marine sediment core

ice core

sea ice

## ABSTRACT

We use eighteen timescale-synchronised near-surface temperature reconstructions spanning 10–50 thousand years before present to clarify the regional expression of Dansgaard-Oeschger (D-O) and Heinrich (H) events in the North Atlantic. The North Atlantic Drift region shows D-O temperature variations of ca. 2–5° with Greenland-like structure. The Western Iberian Margin region also shows Greenland-like structure, but with more pronounced surface cooling between interstadials and Heinrich stadials (ca. 6–9 °C) than between interstadials and non-Heinrich stadials (ca. 2–3 °C). The southern Nordic Seas show smaller D-O temperature anomalies (ca. 1–2 °C) that appear out of phase with Greenland. These spatial patterns are replicated in a new global climate model simulation that features unforced (D-O-like) and freshwater forced (H-like) abrupt climate changes. The model simulations and observations suggest consistently that the spatial expression and amplitude of D-O and H event temperature anomalies are dominated by coupled changes in the Atlantic Meridional Overturning, sea ice extent, polar front position and thermocline structure.

© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Dansgaard-Oeschger (D-O) events are abrupt changes in the glacial climate state that were first identified in Greenland ice cores records (Dansgaard et al., 1984). The events are marked in Greenland by surface temperature changes of up to 16 °C (Kindler et al., 2014), sometimes occurring within decades (Erhardt et al., 2019). A large body of work shows abrupt sea surface

temperature, ocean circulation and sea ice changes contemporaneous with the D-O events in the high latitude North Atlantic and Nordic Seas (e.g. Bond et al., 1999; Shackleton et al., 2000; Voelker, 2002; Dokken et al., 2013). Beyond the North Atlantic, the D-O events manifest globally as coupled changes in the ocean, atmosphere and hydrological cycle extending to the tropics and southern high latitudes (e.g. Blunier and Brook, 2001; Wang et al., 2001; Markle et al., 2017; Brovkin et al., 2021).

In Greenland ice cores, the period of warmth that comes with the D-O events is termed a Greenland Interstadial (GI); temperatures gradually cool during the GI, often terminating in an abrupt drop to a cold Greenland Stadial (GS) state (Johnsen et al., 1992;

<footnote><sup>a</sup> Australian Antarctic Division, Kingston, Tasmania, Australia</footnote>

<footnote><sup>b</sup> Australian Antarctic Program Partnership, Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, Tasmania, Australia</footnote>

<footnote><sup>c</sup> Physics of Ice, Climate and Earth, University of Copenhagen, Copenhagen, Denmark</footnote>

<footnote><sup>d</sup> NORCE Norwegian Research Centre AS, Bjerknes Centre for Climate Research, Bergen, Norway</footnote>

<footnote><sup>e</sup> Divisao Geologia e Georecursos Marinhos, Instituto Portugues Do Mar e da Atmosfera, Alges, Portugal</footnote>

<footnote><sup>f</sup> Centre of Marine Sciences (CCMAR), University of the Algarve, Faro, Portugal</footnote>

<footnote><sup>g</sup> Sorbonne Université, CNRS-IRD-MNHN (LOCEAN UMR 7159), Paris, France</footnote>

<footnote><sup>h</sup> Department of Earth Science, University of Bergen, Bergen, Norway</footnote>

<footnote><sup>i</sup> Bjerknes Centre for Climate Research, Bergen, Norway</footnote>

<footnote><sup>j</sup> Center for Earth Evolution and Dynamics, University of Oslo, Oslo, Norway</footnote>

<footnote>* Corresponding author. Australian Antarctic Division, Kingston, Tasmania, Australia.</footnote>

<footnote>E-mail address: joel.pedro@aad.gov.au (J.B. Pedro).</footnote>
Kindler et al., 2014). During Marine Isotope Stage 3 (MIS 3; ca. 27–60 ka) a sequence of eighteen prominent D-O events are observed (North Greenland Ice Core Project members, 2004). Our focus here is on constraining, within a consistent framework, the temperature expression of these events in the North Atlantic.

There is general agreement that GI are associated with an active (comparable to modern) state of the Atlantic Meridional Overturning Circulation (AMOC; e.g. Ganopolski and Rahmstorf, 2001; Böhm et al., 2015), with increased ocean heat transport into the high-latitude North Atlantic and Nordic Seas (e.g. Sessford et al., 2019; Li and Born, 2019) and reduced sea ice extent (e.g. Hoff et al., 2016; Sadatzki et al., 2019). In contrast, GS are thought to be associated with a weak AMOC in which the sea ice margin and sites of deep-water formation are displaced further south (Ganopolski and Rahmstorf, 2001; Hoff et al., 2016; Sadatzki et al., 2019). Some stadials contain Heinrich (H) events, during which layers of ice-rafted debris in North Atlantic marine sediment cores attest to massive iceberg discharges, of 3–15 m sea level equivalent, from the Laurentide Ice Sheet, likely via Hudson Strait (see e.g. Hemming, 2004). The icebergs drifted as far south as the Iberian Margin (e.g. Eynaud et al., 2009), cooling and freshening surface waters on a basin wide scale (e.g. Hillaire-Marcel and de Vernal, 2008).

Many challenges have complicated developing a coherent picture of D-O and H event temperature anomalies in the North Atlantic from the existing body of research. Among these are the range of chronologies in use for North Atlantic marine records, challenges in developing chronologies consistent with the Greenland ice cores, inconsistencies in temperature reconstructions between proxy types and between calibration equations, insufficient temporal resolution at many sites, and the brute task of compiling and consistently dating the disparate data sets.

The observational component of the current study builds on Jensen et al. (2018), which compiled fourteen North Atlantic sea surface temperature reconstructions across the interval 30–40 kyr (GS/GI cycles 5–8). The key differences compared to the Jensen et al. (2018) are that we here extend the time interval to 10–50 ka, transfer all records to a GICC05-consistent timescale (Section 2.3), and use a single method (the maximum likelihood (ML) technique) for all temperature reconstructions.

Our paper has three main objectives: 1), to compile a set of MIS-3 near-surface temperature reconstructions from the North Atlantic and marginal seas, that are internally consistent in terms of proxy type, chronology and transfer function/calibration technique; 2), to use this compilation to constrain the spatial extent and amplitude of D-O and H event temperature anomalies; and 3), to compare these temperature anomaly patterns to a new global climate model simulation, branched from Vettoretti and Peltier (2018), that features unforced (D-O-like) and freshwater forced (H-like) abrupt climate changes. Following this comparison, we draw on the data and model results to identify key factors controlling the observed temperature spatial expression of these events. Previous proxy and model studies have not to our knowledge sought to compare D-O and H event temperature anomalies within the same framework.

## 2. Materials and methods

We use planktonic foraminifera faunas from eighteen marine sediment cores (Table 1) to generate our temperature reconstructions for the North Atlantic and marginal seas. All records are synchronised to the Greenland Ice Core Chronology 2005 (GICC05). The resulting temperature time series are provided in the Appendix Supplement (Figs A1–A3). To identify the primary regional temperature signals in these records, we constructed

stacked time series (Section 2.5) using data respectively from all eighteen sites and from three sub-regions: the Nordic Seas gateways ($n = 2$); North Atlantic Drift and branches ($n = 5$); and the Western Iberian Margin ($n = 8$). Hereafter, we refer to these as the ‘All’, ‘Gateways’, ‘Drift’ and ‘Margin’ stacks (Fig. 1).

The subsections below describe the core selection criteria, the temperature reconstruction technique, the core dating and stacking methodologies, and the definitions of the three sub-regions. The section finishes with a description of the general circulation model experiments.

### 2.1. A note on nomenclature: GS, GI, stadial, interstadial, HS, non-HS

Our nomenclature follows the recommendations of the INTIMATE (Integration of Ice core, Marine and Terrestrial records of the North Atlantic) working group (Rasmussen et al., 2014). The GS and GI terms are defined from the Greenland ice core record and we reserve their use for that archive. We use the more general terms stadial and interstadial, with the same numbering system, for their counterparts in the marine records and model runs. We classify stadials which contain H events as Heinrich stadials (HS) and those that do not contain a known H event as non-Heinrich stadials (non-HS). This classification is not intended to imply that ice rafted debris and meltwater discharge spans the entire stadial duration, it likely does not (Barker et al., 2015).

### 2.2. Core selection criteria

For inclusion in our temperature database, all cores were required to satisfy the following criteria: 1) temporal resolution of 500 years or better through at least part of the 10–50 ka BP period; 2) complete planktic foraminifera counts; 3) ability to estimate age and age uncertainty relative to the Greenland Ice Core Chronology 2005; and 4) public availability of the required data for 2) and 3) or provision of the data to the authors upon request. The pool of cores satisfying these criteria (Table 1) include nine of those used in Jensen et al. (2018). There is a trade off between spatial coverage and temporal resolution that must be balanced in the core selection criteria. The 500 year cut-off for temporal resolution results in a median resolution across all records of 150 years. The observed signals must therefore be seen as smoothed representations of the actual temperature signals at the respective sites. Nevertheless, as we see further below, signals associated with the majority of GI and GS are resolved in the stacked proxy data.

### 2.3. Temperature reconstruction

Summer season (JAS) temperature reconstructions were generated for this work by applying the maximum likelihood technique (ter Braak and van Dame, 1989) to foraminifera abundances from each core site. We used an Atlantic core-top calibration (Kucera et al., 2005) for all except six sediment cores from along the Western Iberian Margin (MD99-2339, MD95-2042, MD01-2444, MD95-2040, MD95-2041, MD99-2331). For the latter cores, we calculated temperatures by extending the Kucera et al. (2005) Atlantic core-top calibration with additional census counts from core tops from the Iberian Margin (those of Salgueiro et al., 2008; Salgueiro et al., 2014). The motivation for this extended calibration is that due to seasonal upwelling, the western Iberian margin experiences mesoscale variability that is not resolved in the Kucera et al. (2005) data set. The expanded calibration data used for the Iberian margin records is consistent with the calibration of the other records (with regard to species identification and WOA 1998 temperatures for the modern sea surface temperature), while
<caption>**Table 1**</caption>

Core names and basic metadata for the maximum likelihood (ML) near-surface temperature reconstructions. "Res." refers to the median core resolution over the 10–50 ka BP interval. "W19" refers to Waelbroeck et al. (2019).

<table><thead><tr><th>Core name</th><th>°N, °E</th><th>Res.</th><th>Dating</th><th>Foram. counts (ref)</th><th>Foram. counts (data)</th></tr></thead><tbody><tr><td>MD95-2010</td><td>66.68, 4.57</td><td>105</td><td>W19</td><td>Jensen et al. (2018)</td><td>T. Dokken pers. comm. July 2018</td></tr><tr><td>SU90-24</td><td>62.67, -37.38</td><td>192</td><td>W19</td><td>Elliot et al. (1998)</td><td>https://doi.org/10.1594/PANGAEA.881,871</td></tr><tr><td>MD99-2284</td><td>62.37, -0.98</td><td>6</td><td>W19</td><td>Dokken et al. (2013)</td><td>https://doi.org/10.1594/PANGAEA.882,431</td></tr><tr><td>MD99-2281</td><td>60.34–9.46</td><td>147</td><td>W19</td><td>Wary et al. (2015)</td><td>https://www.clim-past.net/11/1507/2015/</td></tr><tr><td>SO82-5</td><td>59.19, -30.90</td><td>86</td><td>W19</td><td>van Kreveld et al. (2000)</td><td>https://doi.org/10.1594/PANGAEA.261,313</td></tr><tr><td>MD04-2829CQ</td><td>58.95, -9.57</td><td>70</td><td>this work</td><td>Hall et al. (2011)</td><td>https://doi.pangaea.de/10.1594/PANGAEA.881,437</td></tr><tr><td>GIK23415-9</td><td>53.18, -19.15</td><td>242</td><td>W19</td><td>Vogelsang et al. (2001)</td><td>https://doi.org/10.1594/PANGAEA.186,156</td></tr><tr><td>BOFS-5k</td><td>50.69, -21.68</td><td>319</td><td>this work</td><td>Maslin (1993)</td><td>M. Maslin pers. comm. July 2018;</td></tr><tr><td>MD04-2845</td><td>45.35, -5.22</td><td>132</td><td>W19</td><td>Sánchez Goni et al. (2008)</td><td>https://doi.org/10.1594/PANGAEA.881,433</td></tr><tr><td>SU92-03</td><td>43.20, -10.11</td><td>690</td><td>W19</td><td>Salgueiro et al. (2010)</td><td>https://doi.org/10.1594/PANGAEA.743,086</td></tr><tr><td>MD99-2331</td><td>42.15, -9.68</td><td>316</td><td>W19</td><td>Sánchez Goni et al. (2008); Naughton et al. (2009)</td><td>M. Sánchez Goni, pers. comm. Jan 2018</td></tr><tr><td>CH69-K09</td><td>41.76, -47.35</td><td>185</td><td>W19</td><td>Labeyrie et al. (1999); Waelbroeck et al. (2001)</td><td>https://doi.org/10.1594/PANGAEA.881,450</td></tr><tr><td>MD95-2040</td><td>40.58, -9.86</td><td>154</td><td>W19</td><td>de Abreu et al. (2003)</td><td>https://doi.org/10.1594/PANGAEA.66,811</td></tr><tr><td>SU90-03</td><td>40.51, -32.05</td><td>368</td><td>W19</td><td>Chapman et al. (2000); Chapman and Shackleton (1998)</td><td>A. Voelker pers. comm., Jul 2018</td></tr><tr><td>MD95-2041</td><td>37.83, -9.51</td><td>111</td><td>W19</td><td>Voelker and de Abreu (2011)</td><td>https://doi.org/10.1594/PANGAEA.737,130</td></tr><tr><td>MD95-2042</td><td>37.80, -10.17</td><td>232</td><td>W19</td><td>Penaud et al. (2011)</td><td>A. Voelker pers. comm. Aug 2018</td></tr><tr><td>MD01-2444</td><td>37.57, -10.13</td><td>111</td><td>W19</td><td>Vautravers and Shackleton (2006)</td><td>F. Bassinot pers. comm. Jan 2018</td></tr><tr><td>MD99-2339</td><td>35.89, -7.53</td><td>78</td><td>this work</td><td>Voelker et al. (2006)</td><td>A. Voelker pers. comm. Oct 2018</td></tr></tbody></table>

adding more modern analogs that are better suited to resolve upwelling related conditions along the Iberian margin (see Salgueiro et al., 2008; Salgueiro et al., 2010, for discussions).

The choice of a transfer function for paleoenvironmental reconstructions is most often based on the root mean square error (RMSE), used to assess the predictive power of the transfer function between the measured and inferred variables. Choosing a transfer function with low RMSE will, in a statistical sense, provide the best transfer function model. Estimating the predictive power of transfer functions assumes that the test sites are independent of the modelling sites. However, autocorrelation is common in ecological data. Using an independent dataset, Telford and Birks (2005) explored and quantified how auto-correlation affects the statistical performance of different transfer functions used to reconstruct Atlantic SSTs from planktic foraminifera. They concluded that the true RMSE is approximately twice the value of previously published estimates for some methods, e.g., the modern analogue technique. In this study, we have chosen the ML transfer function technique, which is based on a unimodal species response. We base this decision on the results by Telford and Birks (2005) that suggest that this method, and other methods based on the assumption of a unimodal species-environment response model, are more robust with regards to spatial structure in the data.

The RMSE of the transfer function model for our ML temperature reconstructions is 1.8 °C and 1.9 °C, based on cross-validation using bootstrapping, for the Atlantic and Iberian Margin cores, respectively.

Time series of the resulting temperature reconstructions from each of the eighteen sediment cores, overlain with the δ¹⁵N-based surface air temperature reconstruction from NGRIP Greenland, are provided in the Supplement (Figs. A1–A3). Although nominally calibrated to 10-m water depth, at some sites the temperature reconstructions may be influenced by conditions at deeper levels (see discussion of Nordic Seas sites in the main text); to acknowledge this caveat we describe the data sets as ‘near-surface temperature reconstructions’ rather than SST reconstructions.

## 2.4. Dating

The cores selected by our criteria all fall north of 35°N. Reliable radiocarbon dating during MIS-3 is generally precluded for these cores by poorly constrained and climate dependent changes in the mixing of surface waters with deeper (radiocarbon depleted)

waters. Specifically, the combined carbon reservoir age, calibration curve and ¹⁴C measurement uncertainties result in age uncertainties on the order of 1000 years. With this level of age uncertainty, radiocarbon dating is not suitable to align stadiad and interstadial intervals between the marine and ice core records. In place of radiocarbon dating at these mid-to-high latitude sites, we use the strategy recommended by the INTIMATE group (Austin and Hibbert, 2012), namely alignment of stratigraphic markers with the ice-core-based event stratigraphy (Rasmussen et al., 2014).

For fifteen of the eighteen cores, we adopt the GICC05-synchronised age models published in Waelbroeck et al. (2019). For the remaining three cores (MD04-2829CQ, BOFS-5k, MD99-2339), we produced GICC05-synchronised age models using the same methods as Waelbroeck et al. (2019). Three types of tie points to GICC05 are used: 1) south of the Greenland Scotland Ridge, rapid increases in proxy sea surface temperature (SST) in the marine cores are aligned with the stadiad-interstadial transitions in NGRIP δ¹⁸O; 2) along and to the north of the Greenland-Scotland Ridge (where SST signals cannot be confidently aligned with NGRIP), magnetic susceptibility signals are aligned instead; 3) dated tephra layers, where available, provide additional tie points to NGRIP (specifically, Faroe Marine Ash Zone (FMAZ) and the North Atlantic Ash Zone (NAAZ) II; Abbott and Davies, 2012).

Inherent in this age modelling approach, is the assumption that the respective SST and/or magnetic susceptibility records vary in phase with the NGRIP δ¹⁸O air temperature proxy. Several lines of evidence exist in support of this assumption. First, the expectation of thermal equilibrium between ocean surface water and overlying air. Second, the observation of synchronous air and sea surface temperature changes in the North Atlantic during the millennial-scale climate variations of the last deglaciation (e.g. Waelbroeck et al., 2019, and references therein). Third, when tephra layers are aligned between ice and marine cores, the marine core magnetic susceptibility signals also become aligned to NGRIP δ¹⁸O (and vice-versa Dokken et al., 2013). Fourth, general circulation models (e.g. Liu et al., 2009) and the current work exhibit synchronous (within years) changes in North Atlantic SST and Greenland air temperature following AMOC perturbations. This rationale supports that internally consistent chronologies are obtained, despite that all three types of tie point are rarely available in a given core (see discussion in Waelbroeck et al., 2019).

All age models are generated using the *Undatable* age-depth modelling routine (Lougheed and Obrochta, 2019). The
<figure><img src="image_6.png" /><figcaption>**Fig. 1.** Regional near-surface summer (JAS) temperature stacks (blue curves) spanning 10–50 ka BP, showing spatial variability in the structure and amplitude of the D-O events: a) All sites; b) the Nordic Seas gateways; c) the North Atlantic Drift and branches; d) the Western Iberian Margin. Maps on the right show the locations (blue points) and names of cores used in the respective stacks (for the complete list of cores used in the stack of all records see Table 1). The ice core δ¹⁵N-based temperature anomaly at NGRIP Greenland (Kindler et al., 2014) is overlain as a black line on each stack for comparison. The blue shading around each stack estimates its 1-σ uncertainty (defined as the RMSE of the calibration uncertainty and stack standard deviation at each time step). The HS are marked by vertical dark grey shading and numbered at the top of panel a). The GS are marked by vertical</figcaption></figure>
components of the age model uncertainty are the depth uncer-
tainty of each tie point (which we define as half the width of the
transition in SST or magnetic susceptibility and which we estimate
for the tephra layers based on the primary literature) and the 1-σ
uncertainty in the GICC05 timescale (defined as half the GICC05
Maximum Counting Error (MCE)). In addition to these errors,
*Undatable* uses a bootstrap approach to recompute the age models
after withholding individual tie points and sampling over a range of
sedimentation rate uncertainties, giving us quantified age un-
certainties at and between tie points.

Given the above-mentioned dating assumptions, we do not
address in this work the phasing of climate change between cores
or regions; rather we focus on the amplitude and spatial pattern of
the temperature anomalies.

2.5. Regional temperature stacks

Stacking regional temperature proxy data is demonstrated to
reduce local noise and more robustly represent regional climate
trends (e.g. White et al., 1997). Principle component analysis was
not used for our study as there are too many data gaps for it to be
robustly applied.

The All stack is comprised of all records. The regional stacks are
defined as follows: the Nordic Seas gateways (bounded by 62–70°N
and 10°W–5°E; n = 2) chosen to sample the entry to the Nordic
Seas around the Greenland–Scotland ridge; the North Atlantic Drift
and branches (bounded by 48–62°N and 40°W–5°E; n = 5); and
the Western Iberian Margin (bounded by 35–46°N and
15°W–5°W; n = 8). Before stacking, the individual sediment core
temperature reconstructions are placed on a common time grid by
binning into 200-year bins. This width is marginally coarser than
the 150-year median resolution of the compiled proxy data and was
chosen to preserve signal and while minimising oversampling. The
re-gridded records are normalised by subtracting their mean over
the 30–40 ka interval. The subsets of records from within each
region are then averaged (with equal weighting) to generate the
stack. We define the (1-σ) temperature uncertainty of the stack (at
each 200-yr time step) as the standard deviation of the binned data
combined in quadrature with the 1-σ uncertainty of the tempera-
ture calibration. These uncertainties are depicted by the shading
around the respective stacks in Fig. 1, for all four stacks the
respective median uncertainties are less than 2.0 °C. The largest
component of the uncertainty is the calibration component
(1.8–1.9 °C), which is predominantly systematic. Therefore, tem-
peratures anomalies similar to the quoted stack uncertainty should
not be dismissed.

The Iberian Margin stack does not include SU90-03 (32.1°W, 40.5°N) located west of the Azores on the northern margin of the subtropical gyre nor CH69–K09 (47.4°W, 41.8°N), located on the present day Gulf Stream path. These sites are expected to be sensitive to changes in the Gulf Stream and Azores Current dynamics respectively (Labeyrie et al., 1999; Chapman et al., 2000). Both show cold anomalies of 2–6 °C during Heinrich Stadial (HS) intervals (lower in amplitude than most sites on Iberian Margin) and little to no temperature response during non-HS intervals (see Figs A2 and A3 for the respective temperature reconstructions). Inclusion of these sites in the Margin stack was tested and does not significantly change its structure. Core SU90-24 (37.4° W, 62.7° N) falls west of the boundary of our Gateways stack, on the Greenland margin of the Irminger Basin. We did not extend the stack boundary to include it on the basis that frequent IRD pulses from the Iceland ice

sheet (Elliot et al., 1998) complicate its interpretation and dating
relative to the cores in the southern Norwegian seas. Short-lived
(centennial) cold anomalies (of 2–4 °C) likely linked to these IRD
pulses are seen throughout this record, more commonly during HS
than during non-HS (see Fig. A1).

2.6. Earth system model simulations

The new glacial climate simulation used in this study is based
upon a slightly modified (approximately 1° × 1° resolution) version
of the Community Climate System Model Version 4 (CCSM4; Gent
et al., 2011). The model and set up is the same as presented in
Peltier and Vettoretti (2014) and (Vettoretti and Peltier, 2018). The
land (CLM4) and atmosphere (CAM4) share the same Cartesian
latitude and longitude grid, while the ocean (POP2) and sea ice
(CICE4) share the same curvilinear grid with poles centred on
Greenland and Antarctica. The glacial climate simulation uses ice-
age boundary conditions provided by the ICE-6G (VM5a) model
for palaeotopography, bathymetry and land-ice cover (Peltier et al.,
2015). The model is not isotope enabled. The boundary and initial
conditions for this simulation have been discussed in detail previ-
ously in (Vettoretti and Peltier, 2013). Continental land ice, atmo-
spheric trace gases and insolation are representative of the
conditions at Last Glacial Maximum (LGM). A vertically varying
profile of diapycnal diffusivity is used in POP2 in the simulation:
diapycnal diffusivity increases below the main thermocline (to
κ = 10⁻⁴ m²s⁻¹) and tapers toward the surface (to κ = 10⁻⁵ m²s⁻¹),
for a detailed discussion see Vettoretti and Peltier (2018).

This is the first publication of a 1° × 1° CESM4 simulation that
features the unforced D-O like oscillations present Peltier and
Vettoretti (2014) and Vettoretti and Peltier (2018), but with the
addition of Heinrich-like events. The simulation is branched,
starting at year 2200, from that of Vettoretti and Peltier (2018) and
run for approximately 3000 model years. The branch point was
chosen during the first regular unforced stadial. To simulate a H-
like event in the unforced D-O simulation, we apply freshwater
(0.05 Sv for 500 years) to the North Atlantic, between 50° N and 70°
N, during the initial model stadial. The volume and duration of
freshwater is informed by previous studies, which estimate that
MIS 3 Heinrich events produce a freshwater flux over the North
Atlantic of approximately 0.05 Sv for 500 years (e.g. Bassis et al.,
2017). This is equivalent to approximately 2 m of relative sea
level rise. Lithic detrital carbonate deposits in marine cores found
throughout the North Atlantic from ice-rafted debris (IRD) indicate
that the MIS-3 Heinrich events emanate from the Laurentide
Hudson Strait Ice Stream (Hemming, 2004). The freshwater forcing
in the model results in a 5 Sv weaker AMOC than seen in the regular
unforced model stadial and a stronger warming transition into the
interstadial after the H event terminates (Fig. 2). Following a longer
and warmer interstadial period, the simulation returns to the reg-
ular unforced D-O oscillation simulated previously with the model
(Peltier and Vettoretti, 2014; Vettoretti and Peltier, 2018).

Previous modelling studies have not investigated the climate
response to H-events and D-O events within a single simulation.
Interactive ice sheet models have typically been used to simulate
the climate response to H-events (e.g. Roberts et al., 2014; Ziemen
et al., 2019), whereas models that have produced unforced DO-like
oscillations typically have prescribed ice sheets (Brown and
Galbraith, 2016; Klockmann et al., 2020).

In comparing near-surface temperature reconstructions be-
tween model and proxy data we use the standard summer season
<figure><img src="image_7.png" /><figcaption>**Fig. 2.** Modelled freshwater-forced (H-like) and unforced (D-O-like) variations in North Atlantic climate, (simulation branched at year 2000 from Vettoretti and Peltier, 2018). a) Modelled AMOC strength (defined as the maximum stream function in the North Atlantic below 500 m); the grey bar marks the timing of freshwater addition to the model in the northern North Atlantic. The subsequent AMOC reductions occur spontaneously, without freshwater addition. b) Modelled mean annual sea ice area in the North Atlantic and marginal seas. c) The surface temperature at the grid cell corresponding to the NGRIP ice core site. d) The δ¹⁵N-based NGRIP air temperature proxy from Kindler et al. (2014) (black line) and the modelled NGRIP temperature data from panel c averaged to 50 year resolution (grey line). The model data in d) is displayed on an adjusted time axis (see Section 2.7), aligned at the black triangles with the NGRIP temperature data.</figcaption></figure>

from the model (JJA) and the standard foraminifera growth season (JAS) for the proxy calibration. The one month offset negligibly affects the spatial pattern and amplitude of the model temperature anomalies described.

## 2.7. Earth system model timescale treatment

The most obvious difference between the model results and the observed GS-9 to GI-6 sequence is that the modelled stadial and interstadial durations are shorter (by ~50%) than their analogs in the data (Fig. 2c&d). The relatively short interstadial duration may in part be the result of the LGM CO₂ level used in this run (185 ppm)

and its influence on feedbacks that reduce the stability of the interstadial state in the model (see Vettoretti et al., 2022, and Sec 4). To facilitate the model to data comparison, we adjust the time axis of the model to align the modelled stadial-interstadial transitions with those in the observed GS-9 to GI-6 sequence in NGRIP. See Fig. 2d for the modelled and observed NGRIP temperatures displayed on this common timescale and the tie points. The timescale adjustment to NGRIP is made once; all subsequent time-series comparisons between model and marine core data use the same timescale. Our rationale for this adjustment is that the spatial patterns and structure of anomalies provide a useful test of the model physics and the timescale adjustment facilitates the comparison. Ice core support for this assumption comes from the fact that stadial and interstadial durations range between one hundred and many thousand years (Rasmussen et al., 2014), yet the amplitudes of temperature anomalies remain relatively consistent (Kindler et al., 2014). In contrast to our approach, previous model-data comparisons have directly forced the timing of D-O events, generally using North Atlantic freshwater perturbations (e.g. Menviel et al., 2014). The latter approach, while allowing the timing of transitions to match observations, has the drawback of potentially distorting the physics, as it is not demonstrated that freshwater fluxes (or any single perturbation) control D-O events (e.g. Alvarez-Solas et al., 2010; Barker et al., 2015; Kleppin et al., 2015; Li and Born, 2019).

# 3. Results

## 3.1. Regional temperature stacks

The All stack, covering the mid-to-high latitude North Atlantic (Fig. 1a), shows strong imprints of the D-O and H events on mean near-surface temperatures. A clear difference from the NGRIP temperature reconstruction (Kindler et al., 2014) is that the amplitude of HS to interstadial warming in the All stack (of ca. 4–5 °C) is clearly larger than the amplitude of non-HS to interstadial warming (of ca. 1–3 °C). See Section 2.1 for our definitions of the HS, non-HS, stadial and interstadial terms. In contrast, there is no clear difference between HS and non-HS at NGRIP (e.g. Kindler et al., 2014, and see black curves in Fig. 1). All GI and GS in the 10–50 ka BP period are represented in the All stack, with the exceptions that conditions remain relatively warm in the stack (i.e. comparable to earlier interstadial levels) throughout the GS-2.1 b–c interval and that no clear signals are resolved in the stack for GS-2.2 and GI-2.2, whose durations are similar to the stack's 150-year median resolution.

The Gateways stack shows a lower amplitude and more variable pattern of temperature changes compared to the other regions. The low amplitude is likely explained by the near-surface temperatures in this region being close to (and unable to cool below) the limit set by the freezing point of seawater (ca. −1.8 °C; see Fig. A1). In the interval 0–40 ka, during which the cores comprising the stack are well-constrained by tephra layers (in MD99-2284) and magnetic susceptibility variations (MD95-2010 and MD99-2284) (Waelbroeck et al., 2019), there is tendency for warming (of ca. 1–2 °C) during stadials and cooling during interstadials. This out-of-phase relationship with Greenland has been documented previously (Dokken et al., 2013; Rasmussen et al., 2016) and we discuss the processes involved further below in light of our GCM results.

The Drift stack shows the most similar structure to Greenland temperature across the D-O events: all interstadials are resolved and the HS are similar to, or only marginally colder than the non-HS. Interstadials in the Drift stack are typically 2–5 °C warmer than stadials.

The Iberian Margin stack is dominated by pronounced cold
anomalies during HS. The near-surface temperatures during HS drop by ca. 6–9 °C below the preceding interstadials. A comparatively weaker cooling, of ca. 2–3 °C, is seen between interstadials and non-HS. On multi-millennial timescales, we note that the LGM interval (generally recognised as ca. 18–24 ka BP; Mix et al., 2001) is warmer than MIS 3 in the Iberian Margin stack. A warm LGM on the Iberian Margin reconciles with proxy (Böhm et al., 2015) and model-based (Brady et al., 2013) evidence for an active LGM AMOC, with NGRIP remaining relatively cold due to its proximity to the enlarged northern ice sheets.

## 3.2. Coupled glacial climate simulation with unforced D-O oscillations

Fig. 2 shows three key North Atlantic climate parameters from the glacial climate simulation: AMOC strength, mean annual sea ice area and the surface temperature at the grid cell corresponding to the NGRIP ice core site. The structure and amplitude of the modelled NGRIP temperature variations in the simulation (Fig. 2c) are remarkably similar to those reconstructed (Kindler et al., 2014) across the GS-9 to GI-6 sequence (Fig. 2d). The model freshwater forcing is similar to the proxy record: in the model, freshwater at a rate based on the observational data (0.05 Sv for 500 years), is applied to the North Atlantic during the first ‘H-like’ stadial, thereafter the modelled ‘D-O-like’ events occur spontaneously; in the proxy record, stadial 9 contains the H4 iceberg discharge, whilst the subsequent stadials 8 and 7 do not contain known H events (Hemming, 2004).

We construct modelled summer (JJA) temperature stacks using the 10-m temperature data from the model grid cells corresponding to the sediment core sites. For consistency with the treatment of the proxy data, the model timeseries are binned to 200-year resolution prior to stacking. The model stacks are plotted together with our respective proxy stacks and the GS-9 to GI-6 NGRIP temperature reconstructions in Fig. 3. The modelled All, Drift and Margin stacks fall within the standard error bounds (Section 2.5) of their counterpart observed stacks (blue shading), except at the very start of interstadials, where modelled warming exceeds observed warming by several degrees and during GS-7 where the modelled cooling exceeds the envelope of the observed cooling by up to 1°. The up to 500 yr resolution of the observed data likely contributes to this more muted temperature signal in the observations.

The modelled Margin stack, also in agreement with the observations, shows stronger cooling during the HS-like stadial than during the non-HS. The modelled Drift stack shows Greenland-like structure, and as in the observations, there is no significant difference in the modelled temperature anomalies between the H-like and unforced stadials. However, a discrepancy is seen between the modelled and observed Gateways stack. Here, the observations indicate warming during the stadials, peak warmth at the onset of the interstadials and then cooling as the interstadials progress (see the MD99-2284 panel of Fig. A1 for a high-resolution view of this sequence). In contrast, the model shows 10-m temperature anomalies that vary in phase with Greenland. We propose an explanation for this model–data difference below.

Previous research (Dokken et al., 2013) proposed that the stadial warming suggested by foraminifera-abundance-based temperature proxies in the gateways region could be explained by deepening of the depth habitat of foraminifera to sample a warming Atlantic inflow layer that is located beneath a cold, low saline, sea ice caped stadial surface layer. This hypothesis is supported by research in the Arctic, where the depth of the preferred habitat of forami-nifera used in temperature reconstructions (N. pachyderma and

<figure><img src="image_8.png" /></figure>

Fig. 3. Observed and modelled regional temperature stacks across GS-9 to GI-6. The observed stacks (blue) are the same as in Fig. 1; the modelled stacks (red) are the mean 10-m summer temperatures in the corresponding model grid cells. For the Gateways stack (b) we show an additional model curve; the mean summer temperature over the top 150 m (magenta). The ice core $\delta^{15}$N-based temperatures anomaly at NGRIP Greenland (black line Kindler et al., 2014) is overlain for comparison. The model time axis is adjusted to align the modelled and observed interstadial onsets (see Section 2.7); facilitating model–data comparison of amplitudes and spatial patterns but prohibiting interpretation of the relative timing between data and model. GS and GI labels from Rasmussen et al. (2014) at the top. Vertical shading marks HS-4 (dark grey) and the non-HS (light grey). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

T. quinqueloba) can increase to ca. 50–200 m in the presence of a low saline surface layer (e.g. 30 PSU, a salinity level reached near the surface in the simulated stadials) (Volkmann, 2000).

Hovmöller plots of the modelled temperature and salinity profiles at site MD99-2284 (Fig. 4) lend further support to this
<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Modelled temperature, sea ice and salinity evolution over the D-O events at the grid point closest to Nordic Seas gateways core site MD99-2284. a) Model 10-m summer temperature (red) and model mean temperature over the top 150 m (magenta). The 10-m temperatures vary in phase with NGRIP, while the 150-m averages capture modelled subsurface warming during stadials. b) Winter (DJF, green) and summer (JJA, amber) sea ice area, showing near perennial modelled sea ice during stadials, ice free conditions at the onset of interstadials and return of winter sea ice cover as the interstadials progress. c) Hovmöller plot of model salinity over the top 250 m, showing halocline deepening during stadials and halocline breakdown at the onset of interstadials. d) Hovmöller plot of model temperature over the top 250 m, showing warming confined below the halocline during stadials before reaching the surface when the halocline breaks down at the interstadial onset. The grey shading in panels a and b represent the modelled Heinrich-like stadial in which freshwater is applied; the following stadal-interstadial oscillations are unforced. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

hypothesis. A cold, low saline layer develops beneath perennial sea ice in the top 50–100 m of the water column during the stadials, a steep halocline confines the warm Atlantic inflow beneath this layer. With this stratification preventing vertical exchange, the subsurface warms as the stadial progresses, even as the surface temperatures approach the limit set by the seawater freezing point (Fig. 4a, c&d). The halocline breaks down at the onset of the interstadial, with convection releasing subsurface heat, contributing to the initial overshoot of surface temperature and the complete removal of sea ice. The subsurface heat is depleted, winter sea ice extent recovers and surface temperatures trend back toward stadial conditions as the interstadial progresses. If we calculate

average temperature over the top 150 m of the model at the gateways sites, on the assumption that this is a better representation of the depth interval sampled by the foraminifera, then we capture a temperature trend in line with the proxy data (magenta line, Fig. 3b). Reducing the depth range over which model temperatures are averaged to the upper 100 m or extending it to the upper 200 m results in similar trends.

### 3.3. Observed and modelled spatial maps of the D-O and H events

To further illustrate the spatial structure of near-surface temperature changes between HS, non-HS and interstadials, we map
the respective data and model anomalies (Fig. 5).

Fig. 5a shows the warming from a non-HS to an interstadial; specifically, interstadial 7 minus stadial 8. For comparison, the observed NGRIP record warms by 8.2 °C between these intervals. The Drift region shows significant warming in the reconstructions from all five core sites (range: 1.7–4.3 °C, mean: 3.2 ± 0.9 °C). South of the Drift, in the Western Iberian Margin region, significant warming is seen at four of the seven core sites with data coverage (range: 2.4–11.6 °C, mean: 6.1 ± 1.1 °C). The modelled and observed temperature anomalies in these regions are broadly consistent (see Fig. 5 shading). At the entrance to the Nordic Seas, MD99-2284 records significantly cooler conditions in the interstadial than in the non-HS. This cooling contrasts with weak surface warming in the model, but is consistent with the modelled temperature changes beneath the near-surface halocline (see Fig. 4). Modelled 10-m temperature anomalies further north in the Nordic Seas and the Arctic Ocean are small (< 0.5°C). The modelled warming from the non-HS to the interstadial is correlated with the displacement of the perennial sea ice margin; the Drift region is sea ice covered in the non-HS and becomes sea ice free in the interstadial (Fig. 5a contours).

Fig. 5b shows the warming from a HS to an interstadial; specifically, interstadial 8 minus stadial 9 (which includes HS-4). The proxy and modelled anomalies are generally intensified versions of those seen in Fig. 5a. Eight of the ten sites around the Iberian Margin show significant warming anomalies (range: 3.8–12.1 °C, mean 6.7 ± 0.7 °C). The amplitudes of these warm anomalies approach, or exceed, the temperature change at NGRIP (7.9 °C) at many sites. Warming is several degrees larger in the data at the southern end of the Iberian Margin than it is in the model. In the Drift region, four sites show significant warming anomalies (range: 3.1–6.5 °C, mean: 4.9 ± 1.0 °C), broadly consistent with the model. Both sites in the Nordic Seas gateways cool significantly (range: –0.7 to – 1.7 °C, mean: – 1.2 ± 1.4). As above, this cooling contrasts with surface warming in the model but agrees with the modelled cooling beneath the halocline (Fig. 3d). The modelled perennial sea ice limit retreats further north toward the Nordic Seas in this example of an interstadial following an HS (Fig. 5b contours) than was the case for the interstadial following the non-HS (Fig. 5a contours).

Finally, Fig. 5c shows the temperature difference between a HS and non-HS state; specifically, stadial 9 (HS-4) minus stadial 8. In the data the HS is significantly colder on the Western Iberian Margin than the non-HS (range: –3.2 to –5.0 °C, mean – 4.1 ± 0.7 °C), in agreement with the modelled stronger cooling in this region in the freshwater-forced H-like stadial. There is no significant temperature difference at NGRIP and no coherent pattern of HS to non-HS temperature differences at the core sites north of the Iberian Margin; only three of eight cores in the latter region show significant trends, one cooling and two warming. Similarly in the model, there is little to no cooling trend at sites north of the Margin.

## 4. Discussion

The close agreement between observed and modelled spatial patterns and amplitudes of the non-HS, HS and interstadial temperature anomalies (Figs. 3 and 5) suggests that the model is capturing key elements of the physics of these climate states and

the transitions between them. We now proceed to use the model to help explain the observed spatial variability. The AMOC is weak in the non-HS, yet weaker in the HS and active in the interstadials, consistent with proxy data (Böhm et al., 2015; Henry et al., 2016). The AMOC strength is tightly coupled in the model to changes in sea ice extent, the position of ocean fronts, thermocline structure and subsurface heat accumulation. We focus here on the role of these factors in setting the temperature anomalies. As discussed in Li and Born (2019), the subpolar gyre is a key coupling region for the atmosphere-ice-ocean system in the North Atlantic and the modelled gyre strength also correlates in the simulation with the changes in sea ice, surface temperatures and AMOC strength during the D-O oscillations.

There is a close correspondence between the temperature anomalies and the regions over which sea ice retreats between respective intervals (Fig. 5). At the onset of interstadials, the modelled perennial sea ice limit (defined here as the 15% JJA contour) rapidly retreats north across the Drift region into the southern Norwegian and Greenland Seas (Fig. 5a&b). Note that the mapped sea ice contours show the mean states across the respective intervals and that more extensive retreat accompanies the initial overshoot stage of the interstadials; indeed the southern Norwegian Seas remain completely ice-free in the model in both summer and winter in the first century of the interstadials (e.g. at the MD99-2284 site, as shown in Fig. 4b). The open-ocean conditions are consistent with biomarker-based sea ice reconstructions from the MD99-2284 site (Sadatzki et al., 2019) and from ~ 150 km further west on the Faroe Islands Margin (Hoff et al., 2016). In the biomarker records and a dinocyst-based reconstruction, also from the Faroe Margin (Wary et al., 2016), winter sea ice returns as the interstadial progresses. Similarly, in the model, winter sea ice expands following the overshoot, exceeding > 50% coverage in the southern Norwegian Seas several centuries before the subsequent stadial onset (Fig. 4b).

With the stadial onset, the modelled summer sea ice limit advances across the Drift region, over the Bay of Biscay to just northwest of the Iberian Margin (Fig. 5a). Sea ice presence in the Bay of Biscay at 45°N during all stadials is supported observationally by the near saturation (90–100%) of polar foraminifera N. pachyderma in marine core MD04-2845 (Sánchez Goñi et al., 2008). By modern analogue, such concentrations are associated with sea ice coverage in the modern East Greenland and Labrador Currents.

During HS, the modelled summer sea ice limit extends yet further south, reaching around 42°N on the Western Iberian Margin (Fig. 5c). There is evidence to support that the sea ice edge extended at least this far and possibly further south during H events. Proxy-based reconstructions of the position of the Arctic polar front, which is generally aligned with the sea ice edge, indicate that it reached close to 40°N during HS and that Arctic and subarctic waters penetrated as far south as 36°N (Eynaud et al., 2009). Closely related to the position of the polar front, the North Atlantic IRD belt and associated iceberg melt extended to 40–39°N (Hemming, 2004; Voelker and de Abreu, 2011), and the C37:4 alkenone is recorded as far south as 37.5°N (e.g. Bard et al., 2000). The IRD presence and C37:4 alkenone both point to low surface salinity and surface stratification that would be conducive to winter sea ice presence. The advance of the polar front and sea ice limit to the Western Iberian Margin during H events is therefore consistent with the high amplitude of the HS temperature anomalies in these
<figure><img src="image_10.png" /></figure>

cores (Figs. 1d, 3d and 5c), in contrast to non-HS when the sea ice edge does not migrate this far south.

The spatial correspondence between the observed and modelled temperature anomalies and the regional sea ice anomalies supports a dominant role of sea ice in setting the pattern of D-O and H event temperature anomalies. During interstadials, reduced sea ice accentuates warming by allowing the ocean to absorb more heat in summer, which it releases back to the atmosphere in winter (Li et al., 2005). During stadials, and even more so during HS, extensive sea ice accentuates surface cooling through greatly reduced solar gain in summer and reduced heat loss from the (sea ice insulated) ocean in winter. The fact that in data and model there is no difference in surface temperatures between the HS and non-HS in Greenland, neither in the Labrador nor Nordic Seas, is due to the high-latitude sea ice extent (north of ~55°N) being the same for both states. Therefore the regional ocean to atmosphere heat fluxes, which are critical to setting Arctic air temperatures (Rhines et al., 2008), are also expected to be the same. On the Western Iberian Margin, HS are significantly colder (by on average - 4.1 ± 0.7 °C) than the non-HS stadials (e.g. Figs. 1d, 3d and 5c) because additional HS cooling is driven by further southward displacement of both the polar front and the sea ice edge.

In the model, the subsurface temperature anomalies (and thus the subsurface heat storage) are greater during the HS than during the non-HS (Fig. 4d), this phenomenon extends across the Nordic, Irminger and Labrador seas. The heat storage is greater during HS than non-HS because the surface freshwater input makes the thermocline inversion more stable and longer lasting (He et al., 2020). Temperature proxies sampling intermediate-depth waters in the southern Norwegian Seas (Dokken et al., 2013; Sessford et al., 2019), Faroe Margin (Ezat et al., 2014), Denmark Strait (Sessford et al., 2018) and Labrador Sea (Marcott et al., 2011) substantiate subsurface warming during stadials, stronger subsurface warming during HS and subsurface cooling during interstadials. The modelled subsurface warming during stadials is consistent with that invoked in previous hypotheses, in which subsurface warming impinging on an isostatically-depressed Laurentide ice sheet margin (Alvarez-Solas et al., 2010; Bassis et al., 2017) is a precondition for H events (Shaffer et al., 2004). Some stadials are spared H events in this scenario because of the millennial-scale timescale required to build enough ice to depress the ice sheet margin to the level where the subsurface warming can drive the destabilising subsurface melt (Bassis et al., 2017).

Recent work using an ensemble of runs with the same model as used here, but configured to run at coarser 3 × 3° resolution, has probed in detail the timescale of the D-O oscillation (Vettoretti et al., 2022). Their results support that the intrinsic slow timescale of the oscillation results from the build up and release of this subsurface heat in the North Atlantic. They further show that the

**Fig. 5.** Mapping the summer 10-m temperature differences between climate states in the D-O events in proxy data (points) and in the model (shading). a) Interstadial warming following a non-Heinrich stadial. b) Interstadial warming following a Heinrich stadial. c) The additional cooling in a Heinrich stadial compared to a non-Heinrich stadial. In each map, violet and cyan lines mark the model mean summer sea ice extent (15% contour) for the warmer and cooler states respectively. Proxy data are only plotted where there are two or more data points per interval at the core site. Data and model anomalies that are not significant (according to a two-sided Students' t-test at the *p* < 0.05 level) are filled/shaded white. In calculating the temperature means at each site over the stadials and interstadials, we leave out the first and last 200 years of each interval to reduce sensitivity to dating uncertainty (changing the buffer between 0 and 500 years has little affect on the key patterns, though it does cause the significance of some anomalies to drop from the 95%-90% level). Although proxy data are nominally calibrated to 10-m water depth, the reconstructed temperatures may be representative of greater depths in some locations (see discussion of Nordic Seas sites in the main text). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
length of the stadials and interstadials is sensitive to the back-
ground concentration of atmospheric CO2, as well as noise associ-
ated with internal climate variability. Lower CO2 drives feedbacks,
including importantly on sea ice growth, which favour the stability
of the stadial state and thus longer stadials. Higher CO2 drives
feedbacks favouring the stability of the interstadial state and thus
longer interstadials.

**5. Conclusions**

Our results do not identify the underlying cause or trigger of D-O
and H events. However, they lend further support to the view that
D-O events are part of an oscillatory climate mode (Dokken et al.,
2013; Vettoretti and Peltier, 2018; Vettoretti et al., 2022), that is
not reliant on a systematic trigger (Capron et al., 2021) and in which
freshwater release from the northern ice sheets is an important
feedback in the oscillation (Alvarez-Solas et al., 2010; Bassis et al.,
2017).

Model studies seeking to represent D-O and H events have traditionally required tuning of each abrupt event by adjusting freshwater fluxes (e.g. Ganopolski and Rahmstorf, 2001; Menviel et al., 2014) or applying some other direct forcing, e.g. adjusting the height of the northern ice sheets (Zhang et al., 2014). Parsimony calls into question this direct forcing approach now that free-running GCMs can generate unforced D-O like oscillations without recourse to such tuning (Brown and Galbraith, 2016; Peltier and Vettoretti, 2014; Klockmann et al., 2020). The advance here compared to the latter studies has been to show in a single framework 1) that the modelled temperature anomalies are consistent with those seen in compiled North Atlantic near surface temperature data, and 2) that adding freshwater to the model North Atlantic during the stadial stage leads to more intense cooling on the Iberian margin (though not further north), consistent with the compiled proxy data from HS intervals. Future work is planned to extend the model-data comparison performed here to evaluate the temperature expression of DO and H events and the driving physical processes across a range of forced and unforced simulations in different models.

The current model lacks an interactive ice sheet in which there
would be the opportunity for subsurface warming during the D-O
event to internally trigger ice sheet collapse and a self-sustained
HS, such coupling remains out of reach for GCMs running on
palaeoclimate timescales and remains one of the outstanding
challenges in modelling and understanding of these abrupt events.
Our data compilation will be released to aid evaluation and testing
of other models seeking to represent the physics of abrupt climate
change and for intercomparison with other proxy-based
reconstructions.

Data availability

This manuscript includes the following supplementary data.

* *agenmodels.zip*: containing the age model files for cores MD04-2829CQ, BOFS-5K and MD99-2339. For each core there are 3 files: the inputs to *Undateable* ( `<corename>_udinput_Sl.txt`), the age-depth model

* ( <corename> _admodel.txt ) and a figure showing the
  age–depth model ( <corename> _admodel.pdf ).

* stacks.zip: containing time series of the observed and modelled regional stacks (i.e as used to produce Figs. 1 and 3 of the Main Text).

* maps.zip: containing the temperature anomalies between climate states (i.e. as used to produce Fig. 5a–c maps in the Main Text).

Age models for cores on the Waelbroeck et al. (2019) age models
are at https://doi.pangaea.de/10.1594/PANGAEA.9,000,734/4/
59,554/. All proxy temperature reconstructions and new age
models will be archived on PANGAEA and model data used in this
paper will be archived at https://sid.erda.dk/cgi-sid/ls.py?share_
id=CjgXyMzOlV.

Author contributions

JBP, CA, KHM, TMD conceptualised the work. CA performed the temperature calibrations. JBP compiled the marine sediment core records with assistance from AHLV, MFJ, CW, CA, EGS and TMD. JBP, CW, AHL and TMD performed the dating. JBP carried out the stacking, data model comparisons and made the figures with assistance from GV, SOR, MJ and KHM. GV performed the simulations. JBP wrote the original draft and all authors contributed to review and editing. KHM, TMD and SOR acquired the funding.

Declaration of competing interest

The authors declare that they have no known competing
financial interests or personal relationships that could have
appeared to influence the work reported in this paper.

Acknowledgements

This research was supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement 610055 (the *ice2ice* project) and by the *ChronoClimate* project funded by the Carlsberg Foundation. A.V. acknowledges funding support from FCT - Foundation for Science and Technology through projects UIDB/04326/2020 and IF/01500/2014. J.B.P. acknowledges support from the Australian Government Department of Industry Science Energy and Resources, grant ASCI000002. C.W. acknowledges support from the European Research Council grant ACCLIMATE/339108. S.O.R. acknowledges support from the IceFlow project funded by the VIL-LUM Foundation. The authors thank *ice2ice* project colleagues and the many palaeoceanographers who made their data available. We are grateful to the anonymous reviewers for constructive comments that improved the manuscript.

**Appendix A. Supplementary data**

Supplementary data to this article can be found online at
https://doi.org/10.1016/j.quascirev.2022.107599.
**Appendix A.** Fig. A1 to A3

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Near-surface summer (JAS) temperature reconstructions from sites between 67 and 58°N. The NGRIP $\delta^{15}$N-based surface air temperature reconstruction (Kindler et al., 2014) is overlain for comparison. Greenland stadials containing Heinrich events (Heinrich Stadials: 'HS') are shaded grey. All y-axes for the sediment and ice core temperature reconstructions span a range of 16 °C. The right panels give the core names and locations. The temperature anomalies are all calculated using a foraminifera-abundance based maximum likelihood technique (Section 2.3). Although the temperature reconstructions are nominally calibrated to 10-m water depth, at the sites along and north of the Greenland-Scotland ridge foraminifera are expected to descend beneath the fresh stadial surface layer into the comparatively warm Atlantic layer, so sampling temperature variations over a larger depth range (see the main text for discussion). The combined (1-$\sigma$ RMSE) temperature proxy calibration and measurement uncertainties at all sites range from 1.8 to 1.9 °C. All cores are synchronised to the Greenland Ice Core Chronology 2005 (GICC05). For complete references to the foraminifera data sources and the age models for individual cores refer to Table 1.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. A.2. Near-surface temperature reconstructions (as for Fig. A1) but for sites between 55 and 41° N.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. A.3.** Near-surface temperature reconstructions (as for Fig. A1) but for sites between 41 and 35° N.</figcaption></figure>

## Appendix B. Supplementary Data

This manuscript includes the following supplementary data.

* *agemodels.zip*: containing the age model files for cores MD04-2829CQ, BOFS-5K and MD99-2339. For each core there are 3 files: the inputs to *Undateable* (corename\_udinput_S1.txt), the age-depth model

* (corename\_admodel.txt) and a figure showing the age-depth model (corename\_admodel.pdf).

* *stacks.zip*: containing time series of the observed and modeled regional stacks (i.e as used to produce Figs 1 and 3 of the Main Text).

* *maps.zip*: containing the temperature anomalies between climate states (i.e. as used to produce the Fig 5a-c maps in the Main Text).

## References

Abbott, P.M., Davies, S.M., 2012. Volcanism and the Greenland ice-cores: the tephra record. Earth Sci. Rev. 115, 173–191. https://doi.org/10.1016/j.earscirev.2012.09.001.

Alvarez-Solas, J., Charbit, S., Ritz, C., Paillard, D., Ramstein, G., Dumas, C., 2010. Links between ocean temperature and iceberg discharge during Heinrich events. Nat. Geosci. 3, 122–126. https://doi.org/10.1038/ngeo752.

Austin, W.E.N., Hibbert, F.D., 2012. Tracing time in the ocean: a brief review of chronological constraints (60-8 kyr) on North Atlantic marine event-based stratigraphies. Quat. Sci. Rev. 36, 28–37. https://doi.org/10.1016/j.quascirev.2012.01.015.

Bard, E., Rostek, F., Turon, J.L., Gendreau, S., 2000. Hydrological impact of heinrich events in the subtropical northeast atlantic. Science 289, 1321–1324. https://doi.org/10.1126/science.289.5483.1321.

Barker, S., Chen, J., Gong, X., Jonkers, L., Knorr, G., Thornalley, D., 2015. Icebergs not the trigger for North Atlantic cold events. Nature 520, 333–336. https://doi.org/10.1038/nature14330.

Bassis, J.N., Petersen, S.V., Mac Cathles, L., 2017. Heinrich events triggered by ocean forcing and modulated by isostatic adjustment. Nature 542, 332–334. https://doi.org/10.1038/nature21069.

Blunier, T., Brook, E.J., 2001. Timing of millennial-scale climate change in Antarctica and Greenland during the last glacial period. Science 291, 109–112. https://doi.org/10.1126/science.291.5501.109.
Böhm, E., Lippold, J., Gutjahr, M., Frank, M., Blaser, P., Antz, B., Fohlmeister, J., Frank, N., Andersen, M.B., Deininger, M., 2015. Strong and deep Atlantic meridional overturning circulation during the last glacial cycle. Nature 517, 73–76. https://doi.org/10.1038/nature14059.

Bond, G.C., Showers, W., Elliot, M., Evans, M., Lotti, R., Hajdas, I., Bonani, G., Johnson, S., 1999. The north atlantic's 1-2 kyr climate rhythm: relation to heinrich events, dansgaard/oeschger cycles and the little ice age. Geophys. Monogr.-Am. Geophys. Union 112, 35–58.

Brady, E.C., Otto-Bliesner, B.L., Kay, J.E., Rosenbloom, N., 2013. Sensitivity to glacial forcing in the ccsm4. J. Clim. 26, 1901–1925. https://doi.org/10.1175/JCLI-D-11-00416.1.

Brovkin, V., Brook, E., Williams, J.W., Bathiany, S., Lenton, T.M., Barton, M., DeConto, R.M., Donges, J.F., Ganopolski, A., McManus, J., Praetorius, S., de Vernal, A., Abe-Ouchi, A., Cheng, H., Claussen, M., Crucifix, M., Gallopín, G., Iglesias, V., Kaufman, D.S., Kleinen, T., Lambert, F., van der Leeuw, S., Liddy, H., Loutre, M.F., McGee, D., Rehfeld, K., Rhodes, R., Seddon, A.W.R., Trauth, M.H., Vandeveken, L., Yu, Z., 2021. Past abrupt changes, tipping points and cascading impacts in the Earth system. Nat. Geosci. 14, 550–558. https://doi.org/10.1038/s41561-021-00790-5.

Brown, N., Galbraith, E.D., 2016. Hosed vs. unhosed: interruptions of the Atlantic Meridional Overturning Circulation in a global coupled model, with and without freshwater forcing. Clim. Past 12, 1663–1679. https://doi.org/10.5194/cp-12-1663-2016.

Capron, E., Rasmussen, S.O., Popp, T.J., Erhardt, T., Fischer, H., Landais, A., Pedro, J.B., Vettoretti, G., Grinsted, A., Gkinis, V., Vaughn, B., Svensson, A., Vinther, B.M., White, J.W.C., 2021. The anatomy of past abrupt warmings recorded in Greenland ice. Nat. Commun. 12, 8107–8126. https://doi.org/10.1038/s41467-021-22241-w.

Chapman, M.R., Shackleton, N.J., 1998. Millennial-scale fluctuations in North Atlantic heat flux during the last 150000. Earth Planet Sci. Lett. 159, 57–70.

Chapman, M.R., Shackleton, N.J., Duplessy, J.C., 2000. Sea surface temperature variability during the last glacial-interglacial cycle: assessing the magnitude and pattern of climate change in the North Atlantic. Palaeogeogr. Palaeoclimatol. Palaeoecol. 157, 1–25. https://doi.org/10.1016/S0031-0182(99)00168-6.

Dansgaard, W., Johnsen, S.J., Clausen, H.B., Dahl-Jensen, D., Gundestrup, N., Hammer, C.U., Oeschger, H., 1984. North Atlantic Climatic Oscillations Revealed by Deep Greenland Ice Cores, vol. 29. Washington DC American Geophysical Union Geophysical Monograph Series, pp. 288–298. https://doi.org/10.1029/GM029p0288.

de Abreu, L., Shackleton, N.J., Schönfeld, J., Hall, M., Chapman, M., 2003. Millennial-scale oceanic climate variability off the Western Iberian margin during the last two glacial periods. Mar. Geol. 196, 1–20. https://doi.org/10.1016/S0025-3227(03)00046-X.

Dokken, T.M., Nisancioglu, K.H., Li, C., Battisti, D.S., Kissel, C., 2013. Dansgaard-Oeschger cycles: interactions between ocean and sea ice intrinsically to the Nordic seas. Paleooceanography 28, 491–502. https://doi.org/10.1002/palo.20042.

Elliot, M., Labeyrie, L., Bond, G., Cortijo, E., Turon, J.L., Tisnerat, N., Duplessy, J.C., 1998. Millennial-scale iceberg discharges in the irvingian basin during the last glacial period: relationship with the heinrich events and environmental settings. Paleooceanography 13, 433–446. https://doi.org/10.1029/98PA01792.

Erhardt, T., Capron, E., Rasmussen, S.O., Schüpbach, S., Bigler, M., Adolphi, F., Fischer, H., 2019. Decadal-scale progression of the onset of dansgaard-oeschger forcing events. Clim. Past 15, 811–825. https://doi.org/10.5194/cp-15-811-2019. URL: https://cp.copernicus.org/articles/15/811/2019/.

Eynaud, F., de Abreu, L., Voelker, A., Schönfeld, J., Salgueiro, E., Turon, J.L., Penaud, A., Toucanne, S., Naughton, F., Sánchez Goni, M.F., Malaizé, B., Cacho, I., 2009. Position of the polar front along the western iberian margin during key cold episodes of the last 45 ka. G-cubed 10, Q07U05. https://doi.org/10.1029/2009GC002398.

Ezat, M.M., Rasmussen, T.L., Groeneveld, J., 2014. Persistent intermediate water warming during cold stadials in the southeastern Nordic seas during the past 65 k.y. Geology 42, 663–666. https://doi.org/10.1130/G35579.1.

Ganopolski, A., Rahmstorf, S., 2005. Rapid changes of glacial climate simulated in a coupled climate model. Nature 409, 153–158.

Gent, P.R., Danabasoglu, G., Donner, L.J., Holland, M.M., Hunke, E.C., Jayne, S.R., Lawrence, D.M., Neale, R.B., Rasch, P.J., Vertenstein, M., Worley, P.H., Yang, Z.L., Zhang, M., 2011. The community climate system model version 4. J. Clim. 24, 4973–4991. https://doi.org/10.1175/2011JCLI4083.1.

Hall, I.R., Colmenero-Hidalgo, E., Zahn, R., Peck, V.L., Hemming, S.R., 2011. Centennial- to millennial-scale ice-ocean interactions in the subpolar northeast Atlantic 18-41 kyr ago. Paleooceanography 26, PA2224. https://doi.org/10.1029/2010PA002084.

He, C., Liu, Z., Zhu, J., Zhang, J., Gu, S., Otto-Bliesner, B.L., Brady, E., Zhu, C., Jin, Y., Sun, J., 2020. North atlantic subsurface temperature response controlled by effective freshwater input in Heinrich events. Earth Planet Sci. Lett. 539, 116247. https://doi.org/10.1016/j.epsl.2020.116247. URL: https://www.sciencedirect.com/science/article/pii/S0012821X20301904.

Hemming, S.R., 2004. Heinrich events: massive late Pleistocene detritus layers of the North Atlantic and their global climate imprint. Rev. Geophys. 42, RG1005. https://doi.org/10.1029/2003RG000128.

Henry, L.G., McManus, J.F., Curry, W.B., Roberts, N.L., Piotrowski, A.M., Keigwin, L.D., 2016. North Atlantic ocean circulation and abrupt climate change during the last glaciation. Science 353, 470–474. https://doi.org/10.1126/science.aaf5529.

Hillaire-Marcel, C., de Vernal, A., 2008. Stable isotope clue to episodic sea ice

formation in the glacial North Atlantic. Earth Planet Sci. Lett. 268, 143–150. https://doi.org/10.1016/j.epsl.2008.01.012.

Hoff, U., Rasmussen, T.L., Stein, R., Ezat, M.M., Fahl, K., 2016. Sea ice and millennial-scale climate variability in the Nordic seas 90 kyr ago to present. Nat. Commun. 7, 12247. https://doi.org/10.1038/ncomms12247.

Jensen, M.F., Nummelin, A., Nielsen, S.B., Sadatzki, H., Sessford, E., Risebrobakken, B., Andersson, C., Voelker, A., Roberts, W.H.G., Pedro, J., Born, A., 2018. A spatiotemporal reconstruction of sea-surface temperatures in the north atlantic during dansgaard-oeschger events 5-8. Clim. Past 14, 901–922. https://doi.org/10.5194/cp-14-901-2018.

Johnsen, S.J., Clausen, H.B., Dansgaard, W., Fuhrer, K., Gundestrup, N., Hammer, C.U., Iversen, P., Jouzel, J., Stauffer, B., Steffensen, J.P., 1992. Irregular glacial interstadials recorded in a new Greenland ice core. Nature 359, 311–313. https://doi.org/10.1038/359311a0.

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., Leuenberger, M., Spahni, R., Capron, E., Chappellaz, J., 2014. Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core. Clim. Past 10, 887–902. https://doi.org/10.5194/cp-10-887-2014.

Kleppin, H., Jochum, M., Otto-Bliesner, B., Shields, C.A., Yeager, S., 2015. Stochastic atmospheric forcing as a cause of Greenland climate transitions. J. Clim. 28, 7741–7763. https://doi.org/10.1175/JCLI-D-14-00728.1.

Klockmann, M., Mikolajewicz, U., Kleppin, H., Marotzke, J., 2020. Coupling of the subpolar gyre and the overturning circulation during abrupt glacial climate transitions. Geophys. Res. Lett. 47, e2020GL090361. https://doi.org/10.1029/2020GL090361. https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020GL090361. https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020GL090361.e2020GL09036110.1029/2020GL090361.

Kucera, M., Rosell-Melé, A., Schneider, R., Waelbroeck, C., Weinetl, M., 2005. Multipoxy approach for the reconstruction of the glacial ocean surface (MARGO). Quat. Sci. Rev. 24, 813–819. https://doi.org/10.1016/j.quascirev.2004.07.017.

Labeyrie, L., Leclaire, H., Waelbroeck, C., Cortijo, E., Duplessy, J.C., Vidal, L., Elliot, M., Le Coat, B., Auffret, G., 1999. Temporal Variability of the Surface and Deep Waters of the North West Atlantic Ocean at Orbital and Millennial Scales, vol. 112. Washington DC American Geophysical Union Geophysical Monograph Series, pp. 77–98. https://doi.org/10.1029/GM112p0077.

Li, C., Born, A., 2019. Coupled atmosphere-ice-ocean dynamics in Dansgaard-Oeschger events. Quat. Sci. Rev. 203, 1–20. https://doi.org/10.1016/j.quascirev.2018.10.031.

Li, C., Battisti, D.S., Schrag, D.P., Tziperman, E., 2005. Abrupt climate shifts in Greenland due to displacements of the sea ice edge. Geophys. Res. Lett. 32, L19702. https://doi.org/10.1029/2005GL023492.

Liu, Z., Otto-Bliesner, B.L., He, F., Brady, E.C., Tomas, R., Clark, P.U., Carlson, A.E., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., Cheng, J., 2009. Transient simulation of last deglaciation with a new mechanism for Bølling-Allerød warming. Science 325, 310. https://doi.org/10.1126/science.1171041.

Lougheed, B., Obrochta, S., 2019. A rapid, deterministic age-depth modeling routine for geological sequences with inherent depth uncertainty. Paleozoeanogr. Paleoclimatol. 34, 122–133. https://doi.org/10.1029/2018PA003457.

Marcott, S.A., Clark, P.U., Padman, L., Klinkhammer, G.P., Springer, S.R., Liu, Z., Otto-Bliesner, B.L., Carlson, A.E., Ungerer, A., Padman, J., He, F., Cheng, J., Schmittner, A., 2011. Ice-shelf collapse from subsurface warming as a trigger for Heinrich events. Proc. Natl. Acad. Sci. Unit. States Am. 108, 13415–13419. https://doi.org/10.1073/pnas.1104772108.

Markle, B.R., Steig, E.J., Buizert, C., Schoenemann, S.W., Bitz, C.M., Fudge, T.J., Pedro, J.B., Ding, Q., Jones, T.R., White, J.W.C., Sowers, T., 2017. Global atmospheric teleconnections during Dansgaard-Oeschger events. Nat. Geosci. 10, 36–40. https://doi.org/10.1038/ngeo2848.

Maslin, M., 1993. A Study of the Palaeoceaneanography of the N.E. Atlantic during the Pleistocene. Phd Thesis. University of Cambridge, Cambridge, UK.

Menviel, L., Timmermann, A., Friedrich, T., England, M.H., 2014. Hindcasing the continuum of Dansgaard-Oeschger variability: mechanisms, patterns and timing. Clim. Past. 10, 63–77. https://doi.org/10.5194/cp-10-63-2014.

Mix, A.C., Bard, E., Schneider, R., 2001. Environmental processes of the ice age: land, oceans, glaciers (EPILOG). Quat. Sci. Rev. 20, 627–657. https://doi.org/10.1016/S0277-3791(00)00145-1.

Naughton, F., Sánchez Goni, M.F., Kageyama, M., Bard, E., Duprat, J., Cortijo, E., Desprat, S., Malaizé, B., Joly, C., Rostek, F., Turon, J.L., 2009. Wet to dry climatic trend in north-western Iberia within Heinrich events. Earth Planet Sci. Lett. 284, 329–342. https://doi.org/10.1016/j.epsl.2009.05.001.

North Greenland Ice Core Project members, 2004. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature 431, 147–151. https://doi.org/10.1038/nature02805.

Peltier, W.R., Vettoretti, G., 2014. Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: a “kicked” salt oscillator in the Atlantic. Geophys. Res. Lett. 41, 7306–7313. https://doi.org/10.1002/2014GL061413.2014GL061413.

Peltier, W.R., Argus, D.F., Drummond, R., 2015. Space geodesy constrains ice age terminal deglaciation: the global ICE-6C_C (VM5a) model. J. Geophys. Res. 120, 450–487. https://doi.org/10.1002/2014JB011176.

Penaud, A., Eynaud, F., Sánchez-Goni, M., Malaizé, B., Turon, J.L., Rossignol, L., 2011. Contrasting sea-surface responses between the western Mediterranean Sea and eastern subtropical latitudes of the North Atlantic during abrupt climatic events of MIS 3. Mar. Micropaleontol. 80, 1–17. https://doi.org/10.1016/j.marmicro.2011.03.002.
Rasmussen, S.O., Bigler, M., Blockley, S.P., Blunier, T., Buchardt, S.L., Clausen, H.B., Cvikjanovic, I., Dahl-Jensen, D., Johnsen, S.J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W.Z., Lowe, J.J., Pedro, J.B., Popp, T., Seierstad, I.K., Steffensen, J.P., Svensson, A.M., Vallelonga, P., Vinther, B.M., Walker, M.J., Wheatley, J.J., Winstrup, M., 2014. A stratigraphic framework for abrupt climatic changes during the last glacial period based on three synchronized Greenland ice-core records: refining and extending the (INTIMATE) event stratigraphy. Quat. Sci. Rev. 106, 14–28. https://doi.org/10.1016/j.quascirev.2014.09.007.

Rasmussen, T.L., Thomsen, E., Moros, M., 2016. North Atlantic warming during Dansgaard-Oeschger events synchronous with Antarctic warming and out-of-phase with Greenland climate. Sci. Rep. 6, 20535. https://doi.org/10.1038/srep20535.

Rhines, P., Häkkinen, S., Josey, S.A., 2008. Is Oceanic Heat Transport Significant in the Climate System? Springer Netherlands, Dordrecht, pp. 87–109. https://doi.org/10.1007/978-1-4020-6774-7_5.

Roberts, W.H.G., Valdes, P.J., Payne, A.J., 2014. A new constraint on the size of Heinrich Events from an iceberg/sediment model. Earth Planet Sci. Lett. 386, 1–9. https://doi.org/10.1016/j.epsl.2013.10.020.

Sadatzki, H., Dokken, T.M., Berben, S.M.P., Muschitiello, F., Stein, R., Fahl, K., Menviel, L., Timmermann, A., Jansen, E., 2019. Sea ice variability in the southern Norwegian sea during glacial dansgaard-oeschger climate cycles. Sci. Adv. 5, eaau6174. https://doi.org/10.1126/sciadv.aau6174.

Salgueiro, E., Voelker, A., Abrantes, F., Meggers, H., Pflaumann, U., Loncarić, N., González-Álvarez, R., Oliveira, P., Bartels-Jónsdóttir, H.B., Moreno, J., Wefer, G., 2008. Planktonic foraminifera from modern sediments reflect upwelling patterns off Iberia: insights from a regional transfer function. Mar. Micropaleontol. 66, 135–164. https://doi.org/10.1016/j.marmicro.2007.09.003.

Salgueiro, E., Voelker, A.H.L., de Abreu, L., Abrantes, F., Meggers, H., Wefer, G., 2010. Temperature and productivity changes off the western Iberian margin during the last 150 ky. Quat. Sci. Rev. 29, 680–695. https://doi.org/10.1016/j.quascirev.2009.11.013.

Salgueiro, E., Naughton, F., Voelker, A.H.L., de Abreu, L., Alberto, A., Rossignol, L., Duprat, J., Magalhães, V.H., Vaqueiro, S., Turon, J.L., Abrantes, F., 2014. Past circulation along the western Iberian margin: a time slice vision from the Last Glacial to the Holocene. Quat. Sci. Rev. 106, 316–329. https://doi.org/10.1016/j.quascirev.2014.09.001.

Sánchez Goni, M.F., Landais, A., Fletcher, W.J., Naughton, F., Desprat, S., Duprat, J., 2008. Contrasting impacts of Dansgaard Oeschger events over a western European latitudinal transect modulated by orbital parameters. Quat. Sci. Rev. 27, 1136–1151. https://doi.org/10.1016/j.quascirev.2008.03.003.

Sessford, E.G., Tisserand, A.A., Risebrobakken, B., Andersson, C., Dokken, T., Jansen, E., 2018. High-resolution benthic mg/ca temperature record of the intermediate water in the Denmark strait across d-o stadia-interstadial cycles. Paleogeogr. Paleoclimatol. 33, 1169–1185. https://doi.org/10.1029/2018PA003370.

Sessford, E.G., Jansen, M.F., Tisserand, A.A., Muschitiello, F., Dokken, T., Nisancioglu, K.H., Jansen, E., 2019. Consistent fluctuations in intermediate water temperature off the coast of Greenland and Norway during Dansgaard-Oeschger events. Quat. Sci. Rev. 223, 105887. https://doi.org/10.1016/j.quascirev.2019.105887.

Shackleton, N.J., Hall, M.A., Vincent, E., 2000. Phase relationships between millennial-scale events 64,000-24,000 years ago. Paleogeography 15, 565–569. https://doi.org/10.1029/2000PA000513.

Shaffer, G., Olsen, S.M., Bjerrum, C.J., 2004. Ocean subsurface warming as a mechanism for coupling Dansgaard-Oeschger climate cycles and ice-rafting events. Geophys. Res. Lett. 31, L24202. https://doi.org/10.1029/2004GL020968.

Telford, R.J., Birks, H.J.B., 2005. The secret assumption of transfer functions: problems with spatial autocorrelation in evaluating model performance. Quat. Sci. Rev. 24, 2173–2179. https://doi.org/10.1016/j.quascirev.2005.05.001.

ter Braak, C.J.F., van Dame, H., 1989. Inferring ph from diatoms: a comparison of old and new calibration methods. Hydrobiologia 178, 209–223.

van Kreveld, S., Sarnthein, M., Erlenkeuser, H., Grootes, P., Jung, S., Nadeau, M.J., Pflaumann, U., Voelker, A., 2000. Potential links between surging ice sheets, circulation changes, and the Dansgaard-Oeschger Cycles in the Irminger Sea, 60-18 Kyr. Paleogeography 15, 425–442. https://doi.org/10.1029/1999PA000464.

Vautravers, M.J., Shackleton, N.J., 2006. Centennial-scale surface hydrology off Portugal during marine isotope stage 3: insights from planktonic foraminiferal fauna variability. Paleogeography 21, PA3004. https://doi.org/10.1029/2005PA001144.

Vettoretti, G., Peltier, W.R., 2013. Last Glacial Maximum sea ice sheet impacts on North Atlantic climate variability: the importance of the sea ice lid. Geophys. Res. Lett. 40, 6378–6383. https://doi.org/10.1002/2013GL058486.

Vettoretti, G., Peltier, W.R., 2018. Fast physics and slow physics in the nonlinear dansgaard-oeschger relaxation oscillation. J. Clim. 31, 3423–3449. https://doi.org/10.1175/JCLI-D-17-0559.1.

Vettoretti, G., Ditlevesen, P., Jochum, M., Rasmussen, S., 2022. Atmospheric CO2 control of spontaneous millennial-scale ice age climate oscillations. Nat. Geosci. 15, 300–306. https://doi.org/10.1038/s41561-022-00920-7.

Voelker, A.H., 2002. Global distribution of centennial-scale records for marine isotope stage (mis) 3: a database. Quat. Sci. Rev. 21, 1185–1212. https://doi.org/10.1016/S0277-3791(01)00139-1. decadal-to-Millennial-Scale Climate Variability. URL: https://www.sciencedirect.com/science/article/pii/S0277379101001391.

Voelker, A.H.L., de Abreu, L., 2011. A Review of Abrupt Climate Change Events in the Northeastern Atlantic Ocean (Iberian Margin): Latitudinal, Longitudinal, and Vertical Gradients, vol. 193. Washington DC American Geophysical Union Geophysical Monograph Series, pp. 15–37. https://doi.org/10.1029/2010GM001021.

Voelker, A.H.L., Lebreiro, S.M., Schönfeld, J., Cacho, I., Erlenkeuser, H., Abrantes, F., 2006. Mediterranean outflow strengthening during northern hemisphere coolings: a salt source for the glacial Atlantic? Earth Planet Sci. Lett. 245, 39–55. https://doi.org/10.1016/j.epsl.2006.03.014.

Vogelsang, E., Sarnthein, M., Pflaumann, U., 2001. d18O Stratigraphy, Chronology, and Sea Surface Temperatures of Atlantic Sediment Records (GLAMAP-2000 Kiel). Institut für Geowissenschaften, Christian-Albrechts-Universität, Kiel, Germany. https://doi.org/10.2312/reports-igf.2001.13. Technical Report 13.

Volkmann, R., 2000. Planktic foraminifers in the outer latVICE sea and the fram strait—modern distribution and ecology. J. Foraminifer. Res. 30, 157–176.

Waelbroeck, C., Duplessy, J.C., Michel, E., Labeyrie, L., Paillard, D., Duprat, J., 2001. The timing of the last deglaciation in North Atlantic climate records. Nature 412, 724–727.

Waelbroeck, C., Lougheed, B.C., Vazquez Riveiros, N., Missiaen, L., Pedro, J., Dokken, T., Hajdas, I., Wacker, L., Abbott, P., Dumoulin, J.P., Thil, F., Eynaud, F., Rossignol, L., Fersi, W., Albuquerque, A.L., Arz, H., Austin, W.E.N., Came, R., Carlson, A.E., Collins, J.A., Dennielou, B., Desprat, S., Dickson, A., Elliot, M., Farmer, C., Giraudeau, J., Gottschalk, J., Henderiks, J., Hughen, K., Jung, S., Knutz, P., Lebreiro, S., Lund, D.C., Lynch-Stieglitz, J., Malaizé, B., Marchitto, T., Martínez-Méndez, G., Mollenhauer, G., Naughton, F., Nave, S., Nürnberg, D., Oppo, D., Peck, V., Peeters, F.J.C., Penaud, A., Portilho-Ramos, R.D.C., Repschläger, J., Roberts, J., Rühlemann, C., Salgueiro, E., Sanchez Goni, M.F., Schönfeld, J., Scussolini, P., Skinner, L.C., Skonieczny, C., Thornalley, D., Toucanne, S., Rooij, D.V., Vidal, L., Voelker, A.H.L., Wary, M., Weldeab, S., Ziegler, M., 2019. Consistently dated Atlantic sediment cores over the last 40 thousand years. Sci. Data 6, 165. https://doi.org/10.1038/s41597-019-0173-8.

Wang, Y.J., Cheng, H., Edwards, R.L., An, Z.S., Wu, J.Y., Shen, C.C., Dorale, J.A., 2001. A high-resolution absolute-dated late pleistocene monsoon record from Hulu Cave, China. Science 294, 2345–2348. https://doi.org/10.1126/science.1064618.

Wary, M., Eynaud, F., Sabine, M., Zaragosí, S., Rossignol, L., Malaizé, B., Palis, E., Zumaque, J., Caille, C., Penaud, A., Michel, E., Charlier, K., 2015. Stratification of surface waters during the last glacial millennial climatic events: a key factor in subsurface and deep-water mass dynamics. Clim. Past 11, 1507–1525. https://doi.org/10.5194/cp-11-1507-2015.

Wary, M., Eynaud, F., Rossignol, L., Lapuyade, J., Gasparotto, M.C., Londeix, L., Malaizé, B., Castéra, M.H., Charlier, K., 2016. Norwegian Sea warm pulses during Dansgaard-Oeschger stadials: zooming in on these anomalies over the 35-41 ka cal BP interval and their impacts on proximal European ice-sheet dynamics. Quat. Sci. Rev. 151, 255–272. https://doi.org/10.1016/j.quascirev.2016.09.011.

White, J.W.C., Barlow, L.K., Fisher, D., Grootes, P., Jouzel, J., Johnsen, S.J., Stuiver, M., Clausen, H., 1997. The climate signal in the stable isotopes of snow from Summit, Greenland: results of comparisons with modern climate observations. J. Geophys. Res. 102, 425–439.

Zhang, X., Lohmann, G., Knorr, G., Purcell, C., 2014. Abrupt glacial climate shifts controlled by ice sheet changes. Nature 512, 290–294. https://doi.org/10.1038/nature13592.

Ziemen, F.A., Kapsch, M.L., Klockmann, M., Mikolajewicz, U., 2019. Heinrich events show two-stage climate response in transient glacial simulations. Clim. Past 15, 153–168. https://doi.org/10.5194/cp-15-153-2019. URL: https://cp.copernicus.org/articles/15/153/2019/.