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

We use eighteen timescale-synchronised near-surface temperature reconstructions spanning 10 e 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 e 5 (cid:1) 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 e 9 (cid:1) C) than between interstadials and non-Heinrich stadials (ca. 2 e 3 (cid:1) C). The southern Nordic Seas show smaller D-O temperature anomalies (ca. 1 e 2 (cid:1) 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/).


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;Kindler et al., 2014).During Marine Isotope Stage 3 (MIS 3; ca.27e60 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€ ohm 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 3e15 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 30e40 kyr (GS/GI cycles 5e8).The key differences compared to the Jensen et al. (2018) are that we here extend the time interval to 10e50 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.

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 A1eA3).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 INTI-MATE (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).

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 10e50 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.

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 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 d 15 N-based surface air temperature reconstruction from NGRIP Greenland, are provided in the Supplement (Figs.A1eA3).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.

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 14 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 stadial 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 GICC05synchronised 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 stadialeinterstadial transitions in NGRIP d 18 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 d 18 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 millennialscale 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 d 18 O (and viceversa 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 Table 1 Core names and basic metadata for the maximum likelihood (ML) near-surface temperature reconstructions."Res."refers to the median core resolution over the 10e50 ka BP interval."W19"refers to Waelbroeck et al. (2019).1).The ice core d 15 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-s 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 components of the age model uncertainty are the depth uncertainty 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-s 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 uncertainties 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.

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 62e70 N and 10 We5 E; n ¼ 2) chosen to sample the entry to the Nordic Seas around the GreenlandeScotland ridge; the North Atlantic Drift and branches (bounded by 48e62 N and 40 We5 E; n ¼ 5); and the Western Iberian Margin (bounded by 35e46 N and 15 We5 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 30e40 ka interval.The subsets of records from within each region are then averaged (with equal weighting) to generate the stack.We define the (1-s) 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-s uncertainty of the temperature 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.8e1.9C), which is predominantly systematic.Therefore, temperatures 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 CH69eK09 (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 2e6 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 2e4 C) likely linked to these IRD pulses are seen throughout this record, more commonly during HS than during non-HS (see Fig. A1).

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 iceage 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 previously in (Vettoretti and Peltier, 2013).Continental land ice, atmospheric 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 k ¼ 10 À4 m 2 s À 1) and tapers toward the surface (to k ¼ 10 À5 m 2 s À1 ), 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 Hlike 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 regular 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 between model and proxy data we use the standard summer season light shading and numbered at the base of panel a).The very short lived GS-2.2 and the disputed GS-14 are not shaded (see Rasmussen et al., 2014).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)from the model (JJA) and the standard foraminifera growth season (JAS) for the proxy calibration.The one month offset negligibly affect the spatial pattern and amplitude of the model temperature anomalies described.

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 2 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 stadialeinterstadial 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 modeledata 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).

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.4e5 C) is clearly larger than the amplitude of non-HS to interstadial warming (of ca.1e3 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 10e50 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 bec 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 0e40 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.1e2 C) during stadials and cooling during interstadials.This outof-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 2e5 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.6e9 C below the preceding interstadials.A comparatively weaker cooling, of ca.2e3 C, is seen between interstadials and non-HS.On multi-millennial timescales, we note that the LGM interval (generally recognised as ca.18e24 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€ ohm 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.

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 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 modeledata 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 capped stadial surface layer.This hypothesis is supported by research in the Arctic, where the depth of the preferred habitat of foraminifera used in temperature reconstructions (N.pachyderma and T. quinqueloba) can increase to ca. 50e200 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€ oller plots of the modelled temperature and salinity profiles at site MD99-2284 (Fig. 4) lend further support to this hypothesis.A cold, low saline layer develops beneath perennial sea ice in the top 50e100 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.

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.7e4.3C, 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.4e11.6C, 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.8e12.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.1e6.5C, 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.

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€ ohm 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 (90e100%) of polar foraminifera N. pachyderma in marine core MD04-2845 (S anchez 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. Proxybased 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 40e39 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 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%e90% 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 background concentration of atmospheric CO 2 , as well as noise associated with internal climate variability.Lower CO 2 drives feedbacks, including importantly on sea ice growth, which favour the stability of the stadial state and thus longer stadials.Higher CO 2 drives feedbacks favouring the stability of the interstadial state and thus longer interstadials.

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 freerunning 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 modeledata 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.

Fig. 1 .
Fig. 1.Regional near-surface summer (JAS) temperature stacks (blue curves) spanning 10e50 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 Table1).The ice core d 15 N-based temperature anomaly at NGRIP Greenland

Fig. 2 .
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 d 15 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.

Fig. 3 .
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 d 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 modeledata 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.)

Fig. 4 .
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€ oller plot of model salinity over the top 250 m, showing halocline deepening during stadials and halocline breakdown at the onset of interstadials.d) Hovm€ oller 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 stadialeinterstadial 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.)

Fig
Fig. A.2. Near-surface temperature reconstructions (as for Fig. A1) but for sites between 55 and 41 N.