Dansgaard–Oeschger events in climate models: review and baseline Marine Isotope Stage 3 (MIS3) protocol

Abstract

Dansgaard–Oeschger (D–O) events, millennial-scale climate oscillations between stadial and interstadial conditions (of up to 10–15 ${}^{\circ }$C in amplitude at high northern latitudes), occurred throughout the Marine Isotope Stage 3 (MIS3; 27.8–59.4 ka) period. The climate modelling community up to now has not been able to answer the question of whether our climate models are too stable to simulate D–O events. To address this, this paper lays the ground-work for a MIS3 D–O protocol for general circulation models which are used in the International Panel for Climate Change (IPCC) assessments. We review the following: D–O terminology, community progress on simulating D–O events in these IPCC-class models (processes and published examples), and evidence about the boundary conditions under which D–O events occur. We find that no model exhibits D–O-like behaviour under pre-industrial conditions. Some, but not all, models exhibit D–O-like oscillations under MIS3 and/or full glacial conditions. Greenhouse gases and ice sheet configurations are crucial. However most models have not run simulations of long enough duration to be sure which models show D–O-like behaviour, under either MIS3 or full glacial states. We propose a MIS3 baseline protocol at 34 ka, which features low obliquity values, medium to low MIS3 greenhouse gas values, and the intermediate ice sheet configuration, which our review suggests are most conducive to D–O-like behaviour in models. We also provide a protocol for a second freshwater (Heinrich-event-preconditioned) experiment, since previous work suggests that this variant may be helpful in preconditioning a state in models which is conducive to D–O events. This review provides modelling groups investigating MIS3 D–O oscillations with a common framework, which is aimed at (1) maximising the chance of the occurrence of D–O-like events in the simulations, (2) allowing more precise model–data evaluation, and (3) providing an adequate central point for modellers to explore model stability.

1 Introduction

During a Dansgaard–Oeschger (D–O) event, Greenland transitions between cold Greenland Stadial (GS) and warmer Greenland Interstadial (GI) conditions. The warming can occur within a decade , whilst cooling occurs over a much longer period that is typically several centuries in length. During a warming phase, surface air temperatures over Greenland increase by 10–15 ${}^{\circ }$C . D–O events are best documented during Marine Isotope Stage 3 MIS3; between 27.8–59.4 kyr BP, hereafter ka;, including being recorded in several ice cores from Greenland Fig. ;. Whilst the D–O events recorded in these cores are renowned, the events are global in nature , with known climate signatures including imprints in surface temperature and the hydrological cycle at high northern latitudes , in the tropics , in Eurasia , and in North and South America . While there are no Greenland ice core records of the previous glacial (MIS6 around 140–190 ka), speleothems and Antarctic ice cores indicate that it is extremely likely that D–O events also occurred during MIS6 and earlier glacial periods . This observational evidence shows that D–O events do not occur under interglacial or full Last Glacial Maximum conditions .

In 2011, argued that climate models used in the assessments of the Intergovernmental Panel on Climate Change (IPCC) have not proved their ability to simulate D–O events. This has several implications for the delivery of accurate projections of climate change within the context of tipping points and abrupt climate change . Whilst in the intervening years a number of models have captured key features of D–O events through Atlantic Meridional Overturning Circulation (AMOC) hysteresis behaviour and/or produced D–O-type millennial-scale variability under a range of forcings , we still do not know if climate models are too stable because too few models have been run and led to the publication of an appropriate simulation. This deficiency is related to both the computational expense which prevents models from being run for the longer time periods needed for investigating D–O events and to the lack of an agreed appropriate experimental set-up. The limited knowledge of pre-Last Glacial Maximum (LGM) boundary conditions, in particular in the case of the ice sheet height and distribution, makes it challenging to generate an appropriate MIS3 experimental set-up.

An important question is whether model stability is caused by the model parameters and whether MIS3 conditions are such that the models are in a mono-stable state or in an oscillatory state or if the models exhibit bi-modality, where noise-induced transitions are not induced due to too-low model variability . Previous studies have questioned the significance of the periodic occurrence of D–O events in MIS3 ($\sim \mathrm{1470}$) . If the full glacial period is included, the distribution of waiting time between D–O events is consistent with a random process . Durations of stadials vs. interstadials indicate correlations with global ice volume and orbital parameters , thus underpinning the decision to focus on MIS3 boundary configurations.

Whether models can simulate abrupt changes is a crucial research question: if the current IPCC-class models are too stable to simulate D–O events, their ability to predict future abrupt transitions and their use in identifying tipping points are doubtful. For example, a tipping point may have been recently reached in the Arctic’s Barents Sea ; sea ice loss in the area is linked with enhanced heat transport via an intensified throughflow, or “Atlantification” . In addition, future enhanced precipitation, decline in Arctic sea ice, and melting of glaciers and ice sheets could intensify the supply of freshwater to the North Atlantic and Arctic, which could lead to the reorganisation of the Atlantic circulation and tip the energy distribution between south and north in a similar way to that during D–O events . If climate models do not reliably simulate past tipping events, it suggests that simulations of the coming century may be giving us a false sense of security.

The Coupled Model Intercomparison Project (CMIP) coordinates and designs climate model protocols for the past, present, and future climates and has become an indispensable tool to facilitate our understanding of climate change . The Paleoclimate Model Intercomparison Project 4 (PMIP4) is one of the individual model intercomparison projects which took part in CMIP6 . The design of a common MIS3 experimental protocol would allow the modelling community to address the questions posed above.

This paper compiles current information about unforced D–O-like oscillations in CMIP5/CMIP6 models and discusses the boundary conditions and mechanisms responsible for these oscillations. Given the nomenclature on D–O events varies throughout the literature, we first provide a framework for a more consistent terminology for use within this proposed MIS3 D–O protocol (Table  and Fig. ). Secondly, we review the literature to ascertain whether models reproduce D–O-like events under MIS3 or other climate conditions. We then use this information to develop a protocol for the simulations of D–O events. This protocol focuses on Marine Isotope Stage 3 (MIS3) partly because of the excellent records of D–O events during this period but also because, as our synthesis shows, MIS3 conditions are also conducive to promoting D–O-like events in some models. Given that D–O events did not occur under full glacial conditions in the Last Glacial Period, the proposed modelling protocol is an important improvement on the use of an LGM PMIP protocol. It will undoubtedly help to shed light on the mechanism and processes involved in millennial-scale oscillations during MIS3. The common MIS3 climate modelling protocol is aimed at (1) maximising the chance of the occurrence of D–O-like events in the simulations, (2) improving model–data evaluation, and (3) providing an adequate central point for modellers to also explore model stability. In addition to the protocol for a baseline simulation, we also outline a protocol for a Heinrich-event-preconditioned (freshwater) experiment. These protocols provide a common framework for model experiments to explore cold-period instabilities using commonly specified greenhouse gas (GHG), ice sheet, insolation, and freshwater-related forcings.

Figure 1

MIS3 ice core records and nomenclature. Stable water isotope and CO${}_{\mathrm{2}}$ measurements from Antarctic and Greenland ice cores . See also Table  for D–O nomenclature. The “cnt” red box indicates the period 38 to 32 ka proposed for the MIS3 baseline experiment.

[Figure omitted. See PDF]

Table 1

Terminology.

Term	Description
Abrupt change	We follow the IPCC Assessment Report 4 (IPCC AR4) definition of abrupt event/change . This term refers to a large-scale change, which is much faster than the change in the pertinent forcing (e.g. rising atmospheric CO${}_{\mathrm{2}}$ concentrations).
Tipping point	This term refers to a critical threshold at which a small perturbation can qualitatively modify the development or state of a system .
Tipping element	This term describes large-scale components of the Earth system that could pass a tipping point . Earth system components are the ocean, atmosphere, cryosphere, anthroposphere, and biosphere, which have further important sub-components, e.g. the meridional ocean circulation, the monsoon systems, sea ice, and various ecosystems .
D–O event	During the Last Glacial Period, a series of dramatic climatic fluctuations occurred in the North Atlantic. These are known as D–O events, during which atmospheric and oceanic conditions alternated between relatively mild (interstadial) and full glacial (stadial) conditions . Around 25 abrupt transitions (each completed within a decade) from stadial to interstadial conditions occurred during the Last Glacial Period, and their amplitudes vary from 5 to 16 ${}^{\circ }$C . The duration of interstadials varies from approximately a century to many millennia .
D–O-type oscillations	For the purpose of this MIS3 D–O protocol, the term of D–O-type oscillation refers to D–O-scale climate variability reproduced by climate models, comparable to the D–O events observed in the Greenland ice core record.
Greenland Stadial and Interstadial	We follow the INTIMATE (Integration of Ice core, Marine and Terrestrial records of the North Atlantic) definition of stadial/interstadial terms . The Greenland Interstadial (GI) and Greenland Stadial (GS) period terms are the Greenland expressions of the D–O events and represent warm and cold phases of the North Atlantic (NA) area, respectively.
Heinrich events	These are defined by the presence of layers of ice-rafted debris (IRD) of primarily (but not exclusively) Laurentide origin in North Atlantic sediment cores e.g.. Heinrich events have been observed during some of the longer stadials but likely do not cover the entire period of these longer stadials .
Heinrich Stadial (HS)	This term refers to a stadial containing a specific Heinrich event. indicate that the term HS can refer to the complete stadial period or to part of a stadial only, characterised by changes shown in proxies for IRD, AMOC, or sea surface temperature (SST) .
Bond cycle	D–O events tend to follow a pattern of diminishing amplitude (or a general cooling trend of the GSs) following each Heinrich event (HE) . These cycles of HE-grouped D–O events were named Bond cycles by and . The average gap between HEs is around 7 kyr, so this is the average length of a Bond cycle .

2 Review of spontaneous D–O-type quasi-oscillations in coupled climate models

We compile published evidence of long unforced quasi-oscillations (in the Atlantic Meridional Circulation; AMOC) in IPCC-class models under all climate states in Table , alongside glacial boundary condition simulations which do not show D–O-type oscillations. This permits us to explore the questions of what proportion of models exhibit D–O-like behaviour, which boundary conditions are most conducive to this, and what mechanisms are common to the modelled D–O behaviours. A number of pre-industrial (PI)/present-day model simulations exhibit spontaneous centennial-length cold events (Table ); however, they do not appear to be D–O-like events. We deal with these first.

Under pre-industrial greenhouse gas (GHG) forcing and present-day ice sheets, spontaneous centennial-length cold events that last around 100–200 years occur in four IPCC-class models (Table ). EC-Earth and Community Climate System Model version 4 (CCSM4) show high-atmospheric blocking over the eastern subpolar gyre that causes a cold event under pre-industrial boundary conditions (; Table ). ECHAM6-FESOM also produces cooling events under pre-industrial conditions due to sudden reductions in deep-water convection and increase in sea ice cover in the Labrador Sea . Changes in convection also occur in the Kiel Climate Model (KCM; ); however here centennial-scale variability in the AMOC is linked to variability in Southern Ocean convection. Unlike the CCSM4 and the EC-Earth models, the KCM and ECHAM6-FESOM studies do not indicate an active role of the atmosphere. Although these four models all show abrupt spontaneous cooling events under pre-industrial boundary conditions, these events do not have the typical sawtooth characteristics or longer timescales of D–O-type events.

Regular cycles of D–O-type quasi-oscillations are found in UofT-CCSM4 under LGM boundary conditions . The initiation of the abrupt D–O-type warming events is associated with the opening of a large polynya over the Irminger Sea (Table ). The AMOC spontaneously exhibits D–O-like quasi-oscillations . The salt oscillator is maintained by the salinity gradient between the subtropical gyre and the northern North Atlantic. Although UofT-CCSM4 is the only model to show long unforced quasi-oscillations in the AMOC under full glacial conditions, most of the other PMIP4 LGM simulations have not been run long enough to be sure that such oscillations would not arise if they were run for longer (see Table ). Having said that, ideally models should not show oscillatory D–O-type behaviour when configured under a full glacial climate state, given that in reality D–O events do not occur under full glacial conditions .

D–O-type quasi-oscillations are also found in MIROC4m under mid-glacial conditions . Some aspects of the D–O warming mechanism observed in the UofT-CCSM4, in particular the spatial location of the opening of a big polynya in the Irminger Sea, determining the stadial-interstadial transition, is also identified in MIROC4m (Table ).

Under late-glacial conditions, at 30 ka, a quasi-oscillating AMOC is produced by the HadCM3 model and results from a North Atlantic salt oscillator mechanism similar to that in UofT-CCSM4 . The HadCM3 model also shows millennial-scale climate oscillations triggered by deglacial meltwater discharge in LGM simulations . Under intermediate glacial conditions (MIS3: 40–32 ka), the COSMOS model shows spontaneous millennial-scale climate oscillations triggered solely by orbitally driven insolation changes . Variations in either obliquity or eccentricity-modulated precession lead to climate variations over the tropical and subpolar North Atlantic which exert opposite effects on AMOC strength and hence result in an oscillatory climate regime . The CM2Mc model also produces somewhat smoothed quasi-oscillating AMOC under intermediate MIS3-like boundary conditions, with a present-day ice sheet distribution in combination with a CO${}_{\mathrm{2}}$ concentration of 180 ppm and low obliquity (22${}^{\circ }$) (Table ). The MPI-ESM model exhibits more abrupt D–O-like quasi-oscillations with a present-day ice sheet distribution in combination with CO${}_{\mathrm{2}}$ concentrations ranging between 190–217 ppm Table ;.

In contrast to the above, neither NorESM nor CCSM3 produces D–O-type events or quasi-oscillations under MIS3 conditions (38 ka) (). The NorESM – with a reasonably simulated AMOC and Arctic sea ice distribution in the PI and historical simulations (as documented by ) – simulates a MIS3 climate that is in a stable regime with relatively strong convections in the Norwegian and Labrador seas. Indeed, NorESM sensitivity experiments including large reductions in atmospheric CO${}_{\mathrm{2}}$ levels and Laurentide Ice Sheet heights, aimed at perturbing the system into a cold-stadial-like climate, indicate that the model state appears to be far from a possible threshold . use the LGM ICE-5G ice sheet configuration , with a high Laurentide Ice Sheet (at just over 4000 m), which may have contributed to a strong AMOC in the CCSM3 simulation, alongside its particular background climate.

In summary, IPCC-class models set up with pre-industrial or present-day conditions do not exhibit D–O-type warming events but can feature shorter centennial-length cooling and warming events. This model behaviour is consistent with observations, since millennial-timescale D–O events do not occur under interglacial conditions, but periods of centennial-scale AMOC variability are present throughout several interglacials . Some models which are set up with more MIS3-like conditions exhibit D–O-type warming events, but some do not. Under full LGM conditions only 1 model (UofT-CCSM4) out of 10 (PMIP4 LGM simulations: ) shows spontaneous D–O-type oscillations Tables  and ;.

Since it can take some time for D–O-type oscillations to evolve, it is unclear if some models would develop such oscillations if they were run for longer (at least for 2000 model years). Of the 40 LGM-/MIS3-like simulations Table ;, 16 simulations have been run for less than 2000 years (Table ), which makes it difficult to tell whether any of these simulations are capable of, or likely to, exhibit D–O-like behaviour under specific boundary conditions. In addition the duration of LGM/MIS3 simulations is currently inadequate; we note that the majority of CMIP6 models appear not to have performed any form of glacial period simulation (Table ). Thus, it is difficult to ascertain what proportion of or indeed which models are capable of capturing D–O-like behaviour under any form of glacial period state (Table ).

2.1 The role of ocean and sea ice feedbacks

Changes in the AMOC are crucial to the correct simulation of D–O events . The AMOC features stabilising positive feedbacks: a strong AMOC transports warm and salty water into the subpolar North Atlantic, thus weakening the stratification and also keeping the sea ice cover reduced e.g.. As a consequence, there is a large transport of heat northward across the hemispheres e.g., strong heat loss in the North Atlantic and Arctic, and active deep convection that sustains the strong AMOC. A weak AMOC, on the other hand, is associated with a weaker northward transport of salt and heat. This increases the stratification in the subpolar North Atlantic and thus favours the expansion of sea ice. The weak northward heat transport and the insulating effect of the sea ice keep the density gain due to heat loss small and the AMOC in a weak state e.g.. This weak AMOC state is stable when Antarctic Bottom Water becomes dense and salty enough to replace North Atlantic Deep Water (NADW) in the deep North Atlantic.

Figure 2

Schematic depicting the transition from GS to GI conditions, i.e. a D–O warming event.

[Figure omitted. See PDF]

Figure 3

Schematic depicting the transition from GI to GS conditions, i.e. a D–O cooling event.

[Figure omitted. See PDF]

Sea ice can act as both a slow and fast positive feedback on AMOC-induced changes in climate. Extensive stadial sea ice cover during a weak AMOC state cools Greenland and suppresses atmosphere–ocean exchange of heat and oceanic convection in the North Atlantic . This also leads to a slow build-up of heat in the North Atlantic sub-surface. Foraminifera from marine sediment cores offer evidence to back up the fact that this sub-surface warming occurred before the onset of fast D–O warming events . This heat build-up sets up the conditions for subsequent fast losses of GS sea ice.

Wind-driven and AMOC- and sea-ice-linked salinity changes also play a crucial role in D–O positive and negative feedbacks. Indeed the net freshwater transport in the Atlantic Basin by the AMOC can be used to assess the stability regime of the AMOC . The interaction of subpolar and tropical salinity anomalies at the surface and in the sub-surface and possible roles of the intertropical convergence zone and freshwater export through the Fram Strait are also important in D–O-related salinity feedbacks. note that if the subtropical gyre shifts northward, and the subpolar gyre contracts, an inflow of salty subtropical water extends over the entire Atlantic Basin east of the Mid-Atlantic Ridge. This inflow can supply salty water to the deep-convection sites in the Iceland Basin and Irminger Sea and help maintain continuous deep convection and a strong AMOC even at low CO${}_{\mathrm{2}}$ concentrations , thus preventing the initiation of GS-like conditions. Where the AMOC does enter a weak state for a prolonged period, and the climate enters a GS, a build-up of heat in sub-surface waters and salt in the tropical Atlantic can enable the very rapid resumption of the AMOC , with the upward mixing of heat from the sub-surface and importation of salt from the tropical Atlantic via gyre mechanisms .

The importance of vertical (diapycnal) mixing in the ocean for these long-timescale, D–O-type instabilities has long been recognised . However, we note that the different ocean and climate models (Table ) parameterise diapycnal mixing in very different ways e.g.. The lack of a single consistent parameterisation and differences in the strength of diapycnal mixing across climate models mean it is to be expected that some models will produce D–O-like oscillations “out of the box” under MIS3 boundary conditions, but others may require changes or tuning to their diapycnal mixing parameters. Even within the same model, a large range of diapycnal diffusivities may yield steady states that satisfy common plausibility constraints such as AMOC transport and sea ice distribution (see e.g. ). This is partly because wind-driven Southern Ocean upwelling plays a complementary role in diapycnal mixing in setting the steady-state overturning and partly because surface buoyancy forcing controls the relative strength of the upper (Atlantic) and lower (Antarctic) overturning cells . A realistic AMOC transport may be obtained due to compensating biases in these processes, which has serious implications for whether AMOC feedbacks (necessary for capturing D–O behaviour) are represented in an adequate manner within these models.

Figures  and show some of the key states, processes, and ocean–sea ice feedbacks that enable D–O events. Following , D–O events can be broken down into four periods: (1) cold-stadial state (Fig. a), (2) rapid-warming phase governed by very-fast-timescale mechanisms (Fig. b), (3) warm-interstadial state (Figs. c and a), and (4) gradual-cooling phase (Fig. b) followed by a faster abrupt transition into a cold-stadial phase (Fig. c). For some of the D–O events, the magnitudes of the warming transitions are on the order of 10 ${}^{\circ }$C in a decade, while the slow cooling in the interstadials is on the order of a few degrees in a millennium (the sawtooth shape) . This picture of rapid retreat of North Atlantic sea ice associated with the resumption of convection and the AMOC, alongside an upwards mixing of salt and heat, followed by a slower-cooling phase back into stadial conditions, matches accumulation, temperature, and water isotopes retrieved from Greenland ice core records of D–O warming events .

2.2 The role of Northern Hemisphere ice sheets

Section  and Table  suggest that large Northern Hemisphere ice sheets and the wind regime associated with these can contribute to a strong AMOC, which stabilises the North Atlantic and prevents D–O events. Thus ice sheets have a critical role to play in setting up the conditions for D–O events . Figures  and show some of the key mechanisms and feedbacks that are behind a state of reduced likelihood for D–O events and a potentially D–O-type oscillating state, respectively.

Figure 4

Schematic showing a state of reduced likelihood for D–O events.

[Figure omitted. See PDF]

Figure 5

Schematic of a potentially D–O-type oscillating state.

[Figure omitted. See PDF]

The Northern Hemisphere Eurasian ice sheet was most probably limited to mountainous areas during mid-MIS3 , and its impact on D–O dynamics was probably relatively small. However, the size and presence (or absence) of the Laurentide Ice Sheet (LIS), which has elevations reaching a maximum of approximately 3000 m at the LGM, does appear to cause important and robust (across multiple models) changes to Northern Hemisphere atmospheric circulation and resultant wind forcing of the ocean. LIS-dependent wind changes influence the subpolar gyre and the stability of the atmosphere–ice–ocean coupled system .

A larger LIS (especially its height) causes stronger Northern Hemisphere winds as well as an amplified stationary wave over North America , and it causes the North Atlantic glacial jet to be more stable due to differences in wave–mean flow feedbacks and alters variability in the large-scale atmospheric circulation, especially in the North Atlantic . In addition, LIS height could control the sea ice coverage and gyre circulation by shifting the westerlies over the North Atlantic region .

LIS altered winds that have wide implications for D–O-relevant tipping elements . note that, first, the presence of a large LGM-type LIS is linked to a strong, more zonal and equatorward-shifted North Atlantic jet, which weakens atmospheric heat transport into the North Atlantic and favours episodes of Greenland blocking . Both could trigger the atmosphere–ice–ocean feedbacks that cause abrupt climate change in this area. Second, a steadier and stronger North Atlantic jet strengthens the wind-driven component of the subpolar gyre . Given that at latitudes north of about 45${}^{\circ }$ N, the subpolar gyre, which is essentially wind-driven, plays a crucial role in the northward transport of heat and salt and is strongly linked to the AMOC e.g.; wind-driven changes in this gyre have a strong impact on the density gain in the North Atlantic.

In many simulations with a large LIS (LGM-like ice sheets), the subtropical gyre can shift northward and cause an inflow of salty subtropical water over deep-convection sites, contributing to continuous deep convection and a strong AMOC even at low CO${}_{\mathrm{2}}$ concentrations . Similarly, note that a higher LIS can promote less South Labrador Sea ice export to the north-eastern North Atlantic (which reduces sea ice concentration) to permit deep convection and shift the core of westerlies northwards, strengthening subtropical gyre for heat and salt transport .

For these reasons, large LGM-type ice sheets, particularly a large LIS, tend to lead to a density gain over the North Atlantic, and the northward salt transport is enhanced with respect to the PI ice sheet case. For many models, but not all, this tends to lead to more active convection in the North Atlantic and a strong AMOC (across a wide range of CO${}_{\mathrm{2}}$ concentrations). That said, the AMOC in many LGM simulations is likely too strong . Thus the AMOC is far away from a tipping point with LGM-size ice sheets for many models .

In some simulations with reduced ice sheets, the jet stream shifts northwards, leading to regional cooling and a rise in seasonal sea ice concentration over the subpolar gyre region . This freshens the area and lowers deep-water formation, which weakens the subpolar gyre, and as a result the simulations are more prone to enter a weak-convection, weak-AMOC mode which is conducive to D–O-type oscillations . Thus, with intermediate MIS3 LIS, i.e. reduced in its height compared to the LGM, multiple AMOC states are more likely .

3 Contours of a baseline MIS3 experiment protocol

Although the choice of a time within MIS3 for a D–O baseline experiment should be unimportant, given that in reality D–O events occurred during the whole of the MIS3, our analysis of existing simulations, boundary conditions, and mechanisms above suggests that there are periods which may be particularly conducive to D–O events occurring in models. Oscillatory D–O-type behaviour appears to be more likely but is not guaranteed when models are run with intermediate or low MIS3 CO${}_{\mathrm{2}}$ values and ice sheets, i.e. reduced in size compared to the LGM , and particularly without a high LIS.

The impact of orbital parameters has been investigated in less detail than the role of GHGs and ice sheets. Using the model COSMOS, demonstrated that under intermediate glacial conditions, obliquity appears to play a significant role in the occurrence of D–O-type behaviour. In particular, the orbital parameters at 40 ka do not produce D–O-type behaviour, whilst at 34 ka lower obliquity ($\sim \mathrm{22.6}$${}^{\circ }$) leads to D–O-type behaviour (see Fig. 2 from ). Additionally, the MIROC4m model produces D–O-type oscillations (under mid-glacial conditions) and low obliquity (22.9${}^{\circ }$) . From these COSMOS and MIROC4m results, we deduce that low obliquity seems to be conducive to D–O behaviour in models.

These considerations suggest that the interval starting at 38 to 32 ka is a good choice for the proposed baseline experiment: it is characterised by (1) a rather regular sequence of D–O events (Fig. ) and (2) has the ideal intermediate MIS3 ice sheet configuration conducive to generating D–O-type quasi-oscillations (Sect. ).

A baseline simulation needs to be run for a sufficient duration to allow the strong positive feedbacks, together with long-timescale negative feedbacks, that enable D–O-type oscillations. The analysis of existing simulations (Sect. ) suggests that this should be a minimum of 5000 years . However, given computational constraints, a minimum duration of around 2000 years, with a spin-up period of 1000 years, may be a more practical minimum requirement for most modelling groups. It would, however, be important to examine and document key metrics for model drift (such as top-of-atmosphere radiation imbalance and deep-ocean or global-mean ocean temperature) during the initial spin-up. The exact length of spin-up is thus subject to discretion of each modelling group based on these key metrics.

There are two obvious possibilities for spinning up the MIS3 control experiment (MIS3-cnt). The baseline experiment could be initialised from the end of either a well-spun-up LGM or PI experiment. Other possibilities could be to spin up from a linear combination of LGM and PI states (as done in ) or spinning up from present-day observations (as done in ). Modelling groups are encouraged to choose whichever option is more feasible/convenient for them. In the event that several spin-up options are available, short spin-ups with diagnosed top-of-atmosphere (TOA) imbalance or global-mean ocean temperature could help distinguish the faster spin-up option. It is worth noting that initial ocean state (i.e. Atlantic salinity stratification) does play a role in abrupt AMOC change and associated feedbacks , the impacts of which shall be considered and evaluated in the future.

We suggest performing a MIS3-cnt experiment centred at 34 ka using GHG and orbital conditions for 34 ka (Fig. ) and an ice sheet configuration as outlined below (Sect.  and ).

We acknowledge that some models might not oscillate under the proposed 34 ka baseline scenario. Indeed, this is expected for NorESM, which under 38 ka conditions is in a stable regime, and the model state seems to be far from a possible tipping point. In spite of that, standardised MIS3 simulations which do not show D–O-like behaviour are still highly valuable for exactly the same reasons that LGM simulations are relevant to the wider modelling community. These standardised MIS3 simulations could contribute to progress in the overarching CMIP6 questions 1 and 2 : how the earth system responds to forcing and what the origins and consequences of systematic model biases are. With a larger number of standardised MIS3 simulations, we would be able to progress in the following research areas/questions:

• Are state-of-the-art climate models capable of representing D–O events under more realistic MIS3 conditions? Benchmarking these simulations will deliver a measure of how well models simulate abrupt changes and tipping events.

• Standardised MIS3 simulations can help explore the existence of a theoretical sweet spot for millennial activity in current climate models . Being close to or within the sweet spot, the AMOC is characterised by high sensitivity to transient and/or noisy climatic forcing or by self-oscillating behaviours .

• If models are too stable to simulate abrupt transitions, what are the processes that contribute to relative levels of model stability?

• In addition, a larger number of standardised MIS3 simulations could encourage the creation of new data sets, improving model–data evaluation.

Figure 6

Monthly zonal-mean MIS3 (34 ka)–PI anomalies of the top-of-atmosphere short-wave incoming radiation (W m${}^{\mathrm{-}\mathrm{2}}$).

[Figure omitted. See PDF]

3.1 Atmospheric trace gases

MIS3 atmospheric CO${}_{\mathrm{2}}$ values varied between a maximum of $\sim$ 233 ppm and a minimum of $\sim \mathrm{187.5}$ ppm Table , Fig. ;. Interestingly, increases of around 5 ppm happened during the abrupt warming of most D–O events, and increases of up to 10 ppm happened within some Heinrich stadials . GHG forcing is critical to model stability. Low (LGM-like) to intermediate (MIS3) CO${}_{\mathrm{2}}$ concentrations tend to be associated with abrupt D–O-type AMOC transitions in models (Sect. ; ). We thus suggest performing the MIS3-cnt experiment using the GHGs values specified in Table  and keeping these values fixed for the whole duration of the simulation, including the spin-up.

3.2 Northern Hemisphere ice sheets

Constraining MIS3 ice sheet boundary conditions is a challenge. Scarcity and fragmentation of evidence are an issue. In particular, it is difficult to determine the size and shape of the ice sheets during MIS3 because subsequent larger LGM configurations have overridden and destroyed evidence of the position of the margins of these smaller ice sheets.

Global sea level fluctuations during the mid-MIS3 were driven nearly exclusively by the LIS . Global-average sea level remained above $-$ 55 m for the period between 30–55 ka. From glacial isostatic modelling and geological constraints, a global-mean sea level between $\mathrm{-}$ 30 and $\mathrm{-}$ 50 m is inferred . For much of MIS3, since the Eurasian ice sheets and the Cordilleran Ice Sheet were likely restricted to mountain-based caps , the primary control on ice volume is assumed to be from the LIS. Recent work in the area of the Hudson Bay suggests that ice-free conditions may have occurred during mid-MIS3. This implies climatic conditions in this region similar to present and a LIS margin removed from the southern Hudson Bay. Similarly show a considerably lower and less extensive LIS compared to ICE-5G and ICE-6G LGM ice sheet reconstructions. sea level curves are consistent with the estimated MIS3 ice sheet volumes from . Using glacial isostatic adjustment modelling, also show that a small LIS can explain high MIS3 sea level estimates alongside the eastern coast of the United States. The synthesis of of numerical modelling results and empirical data provides additional support for a considerable reduction in the MIS3 LIS extent and very minimal European ice sheet.

Figure 7

MIS3 ice sheet reconstruction from (a, b) . Also shown for comparison is the (c, d) LGM ICE-6G ice sheet reconstruction from . The third row shows the differences between the MIS3 ice sheet reconstruction and the LGM ICE-6G reconstruction.

[Figure omitted. See PDF]

The recent MIS3 ice sheet reconstruction, PaleoMIST 1.0 (Paleo Margins, Ice Sheets, and Topography), was developed independently of far-field sea level records and indirect proxy records by . This reconstruction is based on trying to fit the evolution of ice flow indicators, as well as chronological constraints of ice-free conditions.

provide a maximum and minimum MIS3 reconstruction, specifically for the Laurentide Ice Sheet. The maximum scenario is more consistent with recently discovered eastward-oriented, pre-LGM ice flow direction indicators found in south-eastern Manitoba , so we currently consider it to be more likely. However, at 35 ka, the difference between the two scenarios is minor. The difference is primarily with the thickness (and therefore also topography) of the ice sheet, rather than extent, but it amounts to less than 1 m of sea level equivalent. The 35 ka time slice represents conditions after Heinrich Event 4 , and the ice margin in the Hudson Strait retreated by about 350 km from the edge of the continental shelf. The ice margin elsewhere for the Laurentide Ice Sheet is based on chronological constraints, most of which are documented in the compilation by . The Cordillera Ice Sheet extent is based on evidence of relatively restricted ice cover during MIS3 . The Greenland Ice Sheet margin is set to be between the LGM and present-day extent. The Eurasian ice extent at 35 ka includes an advance of ice into the Baltic Sea, which happened after Heinrich Event 4 . For East Antarctica, the margin is set to be the same as present. In West Antarctica, the margin at 35 ka is close to the shelf edge, as the maximum extent may have been achieved by 30 ka .

Given its strong evidence basis, we thus suggest the use of the maximum 35 ka PaleoMIST ice sheet configuration . We note that the LIS is considerably reduced in size compared to the ICE-6G LGM reconstruction in the south-eastern margin (Fig. a, d); the Eurasian ice sheet (EIS) is also significantly smaller (Fig. a, d).

Table 2

Summary of the boundary conditions (BCs) and forcings for the MIS3-cnt experiment.

BC/forcing	Suggested value MIS3-cnt
Atmospheric trace gases	CO${}_{\mathrm{2}}$: 208 ppm
	CH${}_{\mathrm{4}}$: 420 ppb
	N${}_{\mathrm{2}}$O: 204 ppb
Insolation	Eccentricity: 0.01567
	Obliquity: 22.6${}^{\circ }$
	Precession: $-$ 0.016
Solar constant	Same as PI control
Ice sheets	35 ka ice sheet reconstruction ;
	mean global salinity increased by 0.6 PSU to account for ice volume
Global freshwater budget	Closed to avoid drifts; snow should not accumulate
	over ice sheets, and rivers should flow into the ocean.
	Models need to consider lakes when closing
	the global freshwater budget.
Vegetation	Dynamic or fixed as in PI;
	if fixed vegetation: tundra in new land points
Dust	As in PI control
Heinrich-kicked variant	Initial 0.04–1 Sv over 250–500 years
	followed by standard MIS3-cnt simulation

Whilst the implementation of the ice sheet will differ between models, the steps of describe how to implement a glacial-state ice sheet in the IPSL climate model. For consistency, we likewise recommend that the same steps should be followed as far as possible. Since a reduced sea level can modify river courses, recommend that as a minimum, rivers should reach the oceans. Also, the ocean should be initialised with a salinity 0.6 PSU higher than the PI experiment to account for the sea level difference between the MIS3 and PI experiment (freshwater stored as ice on land) .

The single ice sheet reconstruction MIS3 set-up summarised above contrasts with the PMIP4 LGM protocol, which provides three different possible ice sheet configurations (PMIP3, ICE6G-C, and GLAC-1D) for the tier 1 LGM experiment . This partially reflects the more limited knowledge of ice sheets in pre-LGM periods. Exploration of the effect of MIS3 ice sheet reconstruction uncertainties on climate models, particularly on model stability, would be valuable. For this purpose, further additional 34 ky/MIS3 ice sheet reconstructions would be very valuable.

3.3 Heinrich-event-preconditioned option

The term Heinrich-event-like “kick” was initially introduced by to invoke a “kicked” salt oscillator hypothesis (pseudo-Heinrich-type behaviour) to induce D–O-type oscillations in an LGM simulation performed with the UofT-CCSM4. During the first 1000 years of the simulation as the model is spun up, and the ocean cools to reach a state consistent with glacial boundary conditions, there are two thermal thresholds during which the strength of the AMOC rapidly reduces (see Fig. 2 in ). These abrupt transitions in the AMOC coincide with abrupt reductions in surface temperatures in the North Atlantic and abrupt expansions of sea ice coverage. During the second of these events, the AMOC is reduced to approximately 12 Sv, about half its strength in the pre-industrial control . This event may resemble the impact of a Heinrich-event-like kick to the AMOC, though no freshwater perturbation was imposed .

In a more recent study, examine the CCSM4 simulations that show unforced D–O-type oscillations , but with the addition of a (freshwater) Heinrich-like event. . The freshwater flux (0.05 Sv injected into the NA for 500 years) leads to (1) a 5 Sv weaker AMOC compared to the one seen in the unforced model stadial and (2) a stronger D–O warming transition into the interstadial phase . Thus this type of HE preconditioning can trigger abrupt reductions in the AMOC strength and in NA surface temperatures and sea ice coverage, and it may also help induce a stadial state in other models, which is more conducive to unforced (D–O-type) oscillations .

Given the importance of HEs to starting Bond cycles of D–O events, an additional experiment to investigate how HE meltwater preconditioning impacts the simulation of D–O-like oscillations under MIS3 boundary conditions would be valuable. HE freshwater preconditioning may, as in reality, be more conducive to a (Bond-cycle-like) sequence of spontaneous D–O-type oscillations (see Table 1).

The freshwater delivered during Heinrich event iceberg discharge suppresses the AMOC, leading to accumulation of heat in the Southern Hemisphere and in the North Atlantic sub-surface waters e.g.. Estimates of the meltwater input into the North Atlantic during Heinrich events range between 2 and 15 m of sea level equivalent ice volume . It is logical to presume that these freshwater events are important in preconditioning the climate system with respect to D–O behaviour . Freshwater perturbations can trigger changes between AMOC states e.g., and a relatively small freshwater flux applied over convection areas can lead to a shutdown of the AMOC e.g.. Studies have shown similarities between observed global features of abrupt D–O changes and the behaviour seen in freshwater forcing (FWF) experiments . However, the sensitivity of the AMOC to a wide range of freshwater inputs varies according to model, where the meltwater is added, and the background climate state . Given these various uncertainties, we suggest that it would be useful to run an additional experiment to investigate how preconditioning through a Heinrich-like (freshwater) event impacts the simulation of D–O-like oscillations under MIS3 boundary conditions.

The range of Heinrich event volumes calculated using ice sheet models varies from 24.2 to $\mathrm{125}\times {\mathrm{10}}^{\mathrm{4}}$ km${}^{\mathrm{3}}$ , whilst isotope-based estimates and a precipitation balance approach have yielded 86, 649, and $\mathrm{946}\times {\mathrm{10}}^{\mathrm{4}}$ km${}^{\mathrm{3}}$ of ice volume . and using a sediment modelling approach, estimated a discharge of 30 to $\mathrm{120}\times {\mathrm{10}}^{\mathrm{4}}$ km${}^{\mathrm{3}}$ of ice volume. Some of the spread in these estimates could be because the relationship between the oxygen isotope record, sea level, and meltwater volume is not constant when ice is lost from marine basins, such that the use of oxygen isotopes for calculating Heinrich event volumes may produce unrealistically high values . There is also some uncertainty about the duration of the Heinrich events, with some previous studies suggesting that they could be as short as 250 years and others suggesting that a duration of 500 years is more typical . These considerations suggest that it is possible to justify the use of anywhere between 0.02–0.6 Sv freshwater flux over 500 years or 0.04–1.2 Sv over 250 years. More recent estimates of Heinrich event magnitudes tend to favour the lower end of this range. If all forcings are set to MIS3-cnt values, and the Heinrich event freshwater flux is distributed across the North Atlantic, this could yield a range of stadial climates . After 250–500 years this freshwater forcing should be switched off.

4 Conclusions

D–O events are abrupt, large climate changes that punctuated the Last Glacial Period. There is uncertainty whether current IPCC-class models can effectively represent the processes that cause D–O events. We have shown that reduced ice sheets relative to LGM, low obliquity values, and low to medium MIS3 CO${}_{\mathrm{2}}$ values are more likely to lead to unforced quasi-oscillatory D–O-type behaviour. However, the simulations need to be run long enough to allow the strong positive AMOC feedbacks, along with negative feedbacks on long timescales, which can then lead to D–O-type oscillations. Around 40 % of the simulations set up with full LGM or more MIS3-like conditions have a run length of less than 2000 model years, which makes it difficult to tell whether any of these simulations are capable of or likely to exhibit D–O-like behaviour. In addition, the vast majority of PMIP4/CMIP6 models have not run LGM or MIS3-like simulations long enough to be sure which models have the capability to oscillate.

We have provided boundary conditions for a baseline MIS3-cnt simulation and a Heinrich-event-preconditioned variant (freshwater-forced experiment). The MIS3-cnt experiment is centred at 34 ka because it yields the ideal combination of intermediate ice sheets (smaller in size compared to LGM), low obliquity values, and medium to low GHG values conducive to oscillatory D–O-type behaviour in models. Ideally, the MIS3 baseline experiment should be run for 5000 years; however, given computational constraints a minimum duration of 2000 years together with a spin-up of at least 1000 years is a more practical minimum requirement. This baseline MIS3-cnt protocol provides a common framework to explore cold-period instabilities using particular GHG-, insolation-, freshwater-, and NH-ice-sheet-related forcings, together with diapycnal mixing. More model simulations run under the MIS3 D–O protocol proposed here together with analyses across models could provide better insights along the lines of atmospheric–ice–ocean feedbacks behind D–O events. These simulations will allow us to answer questions such as whether current climate models are able to reproduce D–O-type behaviour under more realistic MIS3 conditions, how well models simulate tipping events and abrupt changes, and what the mechanisms are that lead to relative levels of model stability. Moreover, a large number of standardised MIS3 simulations could encourage the creation of new data sets, improving model–data evaluation.

Appendix A Table A1

Models that exhibit spontaneous oscillations.

Study	Model	Period	GHG	Icesheet	Insolation	FWF	Runlength	Main findings	Mechanisms for D–O cooling event	Mechanisms for D–O warming event	Similar mechanism to schematic in Fig.  or
	EC-Earth model	PI	PI	PI	PI	None	1125years	Spontaneous cold event that lasts around 100 years	An anomalous high-atmospheric blocking over the eastern subpolar gyre causes the cold event. Ocean currents transport sea ice southwards, and there is a shutdown of deep-water convection in the Labrador Sea.	No warming event is reported.	No. Atmospheric blocking–sea ice–ocean feedback is identified as a main cause behind the cold event.
	CCSM4	PI	PI	PI	PI	None	1000years	Spontaneous cold event that lasts around 200 years	The cooling event has a duration of 200 years and is linked to a weakened state of the subpolar gyre (SPG) and deep-water convection in the Labrador Sea.	The warming event is triggered by a stronger Icelandic low, and therefore deep-water convection recovers, and SPG circulation resumes.	No. Stochastic atmospheric forcing is identified as a potential cause for sea ice variations.
	ECHAM6-FESOM	Present day (PD)	PD –1990	PD	PD –1990	None	350years	Events of sudden reduction in deep-water convection and increase in sea ice cover in the Labrador Sea	An anomalous inflow of warm and saline water into the deep Labrador Sea causes a weakening of the subpolar gyre and modifies the upper freshwater budget over the Labrador Sea.	No warming event is reported.	No. The strong surface winds over the subtropical NA alter the Gibraltar Strait outflow path into the NA and are identified as the main cause behind the saline and warm anomalies in the deep NA.
	KCM	PD	PD	PD	PD	None	1300years	Centennial-scale variability in the AMOC as well as variations in the NA heat content and subpolar gyre strength	As the NA current accelerates, deep convection in the Weddell Sea enables a positive heat content anomaly to propagate northwards in the upper Atlantic Ocean. Eventually, the heat anomaly reaches the northern NA, resulting in a reduction in deep-water formation there.	The retreat of Antarctic Bottom Water (AABW) leads to an enhanced meridional density gradient that results in an increased North Atlantic Deep Water (NADW) cell.	No. Interhemispheric teleconnection and variability in the Southern Ocean deep-water convection are identified as the main causes behind AMOC oscillations.
	UofT-CCSM4	LGM	LGM –21 ka	LGM –ICE-6G (VM5a)	LGM –21 ka	None	5000years	Spontaneous millennial-scale D–O-type oscillations	The continuous flux of sea ice from the Arctic basin into the NA subpolar gyre area across the East Greenland Current favours melting of sea ice as it moves over the warm ocean surface. This freshwater input restratifies the high-latitude NA and results in a considerable decrease in the rate of NA deep-water formation.	The initiation of the abrupt warming events is associated with the opening of a large polynya over the Irminger Sea. The stability of the water column is key and depends on transport of salt to the subpolar gyre along the Irminger Current and Denmark Strait in the decades preceding the warming event.	Yes
	CCSM4	Glacial simulations	CO${}_{\mathrm{2}}$ levels from 190 to 225 ppm	LGM –ICE-6G (VM5a)	LGM	None	8000–10 000years	Spontaneous millennial-scale D–O-type oscillations within a window of CO${}_{\mathrm{2}}$ levels from 190 to 225 ppm	Old sea ice from the Arctic is exported to the NA, sea ice growth is favoured through ice–albedo feedback, high-latitude convection is reduced through sea ice melt, and consequently Antarctica and Greenland cool. The interstadial-to-stadial transition happens with fast NA sea ice expansion, and NADW production collapses.	During a stadial, sea ice thins in the Southern Ocean and Antarctica warms. There are increases in salt convergence in the NA; NADW fluctuations are amplified via salt advection feedback; and the volume of NADW increases, allowed by late-stadial decreases in AABW formation. Late-stadial vertical stratification experiences thermohaline instability, and the Nordic and Irminger seas are destabilised, triggering rapid sea ice loss in the NA and the transition from stadial to interstadial states.	Yes

Table A1

Continued.

Study	Model	Period	GHG	Icesheet	Insolation	FWF	Runlength	Main findings	Mechanisms for D–O cooling event	Mechanisms for D–O warming event	Similar mechanism to schematic in Fig.  or
	HadCM3B	MIS3 – 30 ka	30 ka	30 ka	30 ka	None	6000years	Spontaneous millennial-scale D–O-type oscillations	Ocean forcing initiates the stadial phase. The collapse of the salinity gradient between the northern NA and subtropical gyre (STG) leads to a reduced advection in the Nordic Seas and decreased deep-water formation.	The initiation of the interstadial phase is associated with a wind-driven atmospheric forcing in the Nordic Seas due to increased regional temperatures, reduced sea ice cover, and increased sea level pressure, which enhances wind stress and convection.	Yes, partially. The D–O-type oscillations reflect a salt oscillator mechanism in the NA.
	COSMOS	MIS3: 40–32 ka	40 ka	40 ka	One transient simulation of 40–32 ka and one 40 ka snapshot simulation (with 34 ka orbital parameters)	None	$+\mathrm{5000}$years	Spontaneous millennial-scale D–O-type oscillations, orbitally induced AMOC changes	Transitions from warm interstadial to cold stadial are linked to (1) a precession-controlled rise in low-latitude boreal summer insolation by modifying the NA low-latitude hydroclimate and/or (2) an obliquity-controlled reduction in high-latitude annual insolation by altering high-latitude sea ice–ocean–atmosphere interactions.	While the AMOC is in its weak phase, a gradual increase in sub-surface temperature in the subpolar ocean together with enhanced northward transport of salt in the NA drive the AMOC back to its strong phase.	Some similarities are that unforced AMOC oscillations are triggered by the tropical salt impact (linked to precession-controlled summer insolation) and/or the subpolar thermal impact (linked to obliquity-controlled mean annual insolation).
	CM2Mc	Mixed forcing	CO${}_{\mathrm{2}}$ $=\mathrm{180}$ ppm	PI	Obliquity: 22${}^{\circ }$; precession: 90${}^{\circ }$	None	Morethan8000years	Spontaneous millennial-scale D–O-type oscillations	During a weak AMOC phase, NA deep convection is largely reduced, and there is an expansion of sea ice in the north-eastern Atlantic. Heat accumulates at depth in the NA linked to the weak advection of warm waters from the tropics.	During a strong AMOC phase, NA deep convection is intense, and there is a retreat of sea ice in the northeast Atlantic.	Yes, partially. Salt advection is a key driver of the oscillations, specifically the salt exchange between subpolar and subtropical NA.
	MPI-ESM	Mixed forcing	CO${}_{\mathrm{2}}$ $=$ 195–217 ppm; CH${}_{\mathrm{4}}$ $=$ 396–494 ppb; N${}_{\mathrm{2}}$O $=$ 209–227 ppb	PI	LGM – 21 ka	None	8000–12 350years	Spontaneous millennial-scale D–O-type oscillations	Stadial phases correspond to weak AMOC and strong SPG phases. The extensive SPG results in low northward salt transport, and deep convection only occurs sporadically in the Iceland Basin. The Nordic Seas are entirely ice-covered, which results in a weak Icelandic Low and therefore in a weak wind stress curl. Sub-surface waters in the Nordic Seas are around 3 K warmer than during interstadial phases.	During interstadial phases, the AMOC is strong, and the SPG is contracted and weak. There is a broad inflow of salty subtropical water to the subpolar NA. Changes in the SPG are driven by variations in the cross-gyre density difference. The eastern NA is fully ice-free. Deep convection occurs continuously in the Iceland Basin, the Irminger Sea, and the Nordic Seas.	Yes, partially. The proposed mechanism behind the spontaneous AMOC oscillations compromises three components: (1) oscillations in salinity comparable to , (2) a density-driven feedback loop comparable to , and (3) a wind-driven feedback loop comparable to and .
	MIROC4m	Mid-glacial conditions	CO${}_{\mathrm{2}}$ $=$ 220 ppm	LGM – ICE-5G	Obliquity: 22.949${}^{\circ }$; eccentricity: 0.04; perihelion: 270${}^{\circ }$ and 90${}^{\circ }$	None	6000years	Spontaneous millennial-scale D–O-type oscillations	Changes in sub-surface ocean temperature in the NA play an important role in modifying the stratification of the vertical water column and then reversing the AMOC mode (thermohaline oscillator).	The SPG remains weak (strong) when the AMOC is weak (strong), as well as during the transitions between the two AMOC modes. There is a positive feedback between AMOC and SPG, in agreement with .	Yes, partially. The opening of a big polynya determines the stadial–interstadial transition. Abrupt changes in AMOC lead to changes in salt advection with the NA subpolar gyre and work as a positive feedback.

Table A2

List of simulations run under MIS3/mid-glacial conditions.

Study	Model	Period	GHG	Ice sheet	Insolation	FWF	Run length	Main findings
	NorESM1-F	MIS3 –38 ka	CO${}_{\mathrm{2}}$ $=$ 215 ppm, CH${}_{\mathrm{4}}$ $=$ 550 ppb, N${}_{\mathrm{2}}$O $=$ 260 ppb	Data-constrained 38 ka	38 ka	None	2500 years (recently extended to 6000 model years)	The equilibrium MIS3 simulation does not show spontaneous D–O-type oscillations. Attempts at perturbing the system into a cold-stadial state by modifying the height of the LIS and atmospheric CO${}_{\mathrm{2}}$ levels show that the modelled MIS3 interstadial state is rather stable, thus questioning the occurrence of spontaneous D–O-type oscillations in the lack of interactive ice sheet–meltwater dynamics.
	CCSM3	MIS3 –38 ka	CO${}_{\mathrm{2}}$ $=$ 215 ppm, CH${}_{\mathrm{4}}$ $=$ 501 ppb, N${}_{\mathrm{2}}$O $=$ 234 ppb	ICE-5G ice sheet configuration .	38 ka	12 hosing/extraction experiments with freshwater fluxes from $\pm \mathrm{0.005}$ to $\pm \mathrm{0.2}$ Sv, injected in the Nordic Seas for 500 years	2170 years	AMOC is more sensitive to meltwater fluxes under MIS3 conditions than under LGM conditions. The lower AMOC stability under MIS3 conditions proposes that D–O-type oscillations could have been triggered by small perturbations in the ocean surface meltwater forcing; e.g. they linked to ice sheet processes.
	MIROC4m	Mid-glacial state	CO${}_{\mathrm{2}}$ $=$ 215 ppm, CH${}_{\mathrm{4}}$ $=$ 350 ppb, N${}_{\mathrm{2}}$O $=$ 200 ppb	Intermediate-size ice sheet configuration (at 15 ka)	15 ka	Hosing experiments with freshwater fluxes of 0.05 and 0.1 Sv, injected in the North Atlantic Ocean (50 to 70${}^{\circ }$ N) for 500 years	More than 2000 years	The climate response to freshwater perturbations is much lower under LGM conditions than under MIS3 conditions. The unperturbed LGM AMOC is unusually weak (around 6 Sv) and thus could barely be further lessened, such that meltwater hosing does not largely affect the large-scale climate.
	CCSM4	Glacial run	CO${}_{\mathrm{2}}$ $=$ 210 ppm	LGM – ICE-6G (VM5a)	LGM	Hosing experiment (Heinrich-event-like pulse) with freshwater fluxes of 0.05 Sv for 500 years in two separate stadial periods in a glacial simulation run with CO${}_{\mathrm{2}}$ of 210 ppm. The freshwater flux is injected in the North Atlantic (50 to 70${}^{\circ }$ N).	8000 years	The Heinrich simulation has a large Northern Hemisphere temperature and AMOC overshoot after the Heinrich stadial ends. Nevertheless, this fast AMOC rise above regular interstadial levels is in agreement with observations only for a few Heinrich-stadial periods (H4 and H5).

Table A2

Continued.

Study	Model	Period	GHG	Ice sheet	Insolation	FWF	Run length	Main findings
	COSMOS (ECHAM5–JSBACH–MPIOM)	Mixed forcing	LGM	Sensitivity experiments applying different heights of the NH ice sheets	LGM	Freshwater fluxes from 0.00 to $\pm \mathrm{0.02}$ Sv, injected in the NA for 100–300 years	Snapshot and transient simulations of 250–700 years	An AMOC bi-stability regime is found under intermediate CO${}_{\mathrm{2}}$ and ice sheet conditions roughly resembling those of the MIS3 climate. In the bi-stable MIS3 regime, transitions from weak to strong AMOC state and vice versa could be initiated by not only gradual variations in LIS height and atmospheric CO${}_{\mathrm{2}}$ but also freshwater perturbations. Changes in the LIS height can initiate a positive atmosphere–ocean–sea ice feedback, leading to D–O-type climate shifts. A gradual increase in the ice sheet height results in a northward shift in the winds and favours a more saline Labrador Sea by both reducing the sea ice/freshwater import from the Arctic and increasing the advection of salt into the area.
	COSMOS (ECHAM5–JSBACH–MPIOM)	Mixed forcing	CO${}_{\mathrm{2}}$ levels increased from 185 to 205 ppm during 500 years.	40 % of LGM ice sheet configuration	LGM	None	500 years	For a moderate ice volume, 0.25–0.45 times the LGM, two stable AMOC modes are identified. A variation of 15 ppm in atmospheric CO${}_{\mathrm{2}}$ concentration – equivalent to changes during D–O cycles containing HEs – is enough to trigger oscillations between a weak stadial state to a strong interstadial circulation state.
	COSMOS (ECHAM5–JSBACH–MPIOM)	Mixed forcing	CO${}_{\mathrm{2}}$ levels increased from 185 to 239 ppm at a rate of 0.02 ppm per year.	Intermediate ice sheet configuration – ice volume equivalent to a sea level drop of approximately 40 m	LGM	None	600 years
	COSMOS (ECHAM5–JSBACH–MPIOM)	Mixed forcing	CO${}_{\mathrm{2}}$ levels increased from 185 and 245 ppm at a rate of 0.05 ppm per year.	LGM ice volume	LGM	Persistent freshwater flux of 0.15 Sv	600 years

Appendix B Table B1

Summary of LGM-/MIS3-like simulations discussed in the text. Models that reproduce D–O-type oscillations are highlighted in bold.

Study	Model	Period	N° of simulations	Run length
	UofT-CCSM4	PMIP4 LGM	1	5000
	AWI-ESM1-1-LR	PMIP4 LGM	1	1300
	AWI-ESM-2-1-LR	PMIP4 LGM	1	600
	CESM1.2	PMIP4 LGM	1	1800
	HadCM3B-M2.1aD	PMIP4 LGM	3	400–2900${}^{\mathrm{a}}$
	iLOVECLIM1.1.4	PMIP4 LGM	2	5000
	INM-CM4-8	PMIP4 LGM	1	50
	IPSLCM5A2	PMIP4 LGM	1	1200
	MIROC-ES2L	PMIP4 LGM	1	8960
	MPI-ESM1.2	PMIP4 LGM	1	3850
	HadCM3B-M2.1aD	MIS3 (30 ka)	1	6000
	COSMOS	MIS3 (40–32 ka)	2	5000
	NorESM	MIS3 (38 ka)	1	$+\mathrm{6000}$
	CCSM3	MIS3 (38 ka)	1	2170
	MIROC4m	Mid-glacial conditions	1	$+\mathrm{2000}$
	MIROC4m	Mid-glacial conditions	2	6000
	CCSM4	Glacial conditions	4${}^{\mathrm{b}}$	8000
	CM2Mc	Mixed forcing	1	$+\mathbf{8000}$
	MPI-ESM	Mixed forcing	3	$+\mathbf{8000}$
	COSMOS	Mixed forcing	11${}^{\mathrm{c}}$	300–4000${}^{\mathrm{d}}$

${}^{\mathrm{a}}$ Only one simulation run longer than 2000 model years. ${}^{\mathrm{b}}$ Four simulations run with CO${}_{\mathrm{2}}$ levels: 200, 210, 220, 225 ppm. ${}^{\mathrm{c}}$ We do not consider FWF runs nor transient simulations forced with varying CO${}_{\mathrm{2}}$ and/or NH ice sheet height. ${}^{\mathrm{d}}$ Only two simulations with a duration of $+\mathrm{2000}$ model years.

Appendix C Table C1

Contributing members to PMIP4/CMIP6 that have run simulations under LGM or MIS3 conditions.

Model	Period	Run length
ACCESS-ESM1-5	–	–
AWI-ESM-1-1-LR	LGM	1300
CESM2	–	–
CNRM-CM6-1	–	–
EC-Earth3-LR	–	–
FGOALS-f3-L	–	–
FGOALS-g3	–	–
GISS-E2-1-G	–	–
HadGEM3-GC31-LL	–	–
INM-CM4-8	LGM	50
IPSL-CM6A-LR	–	–
MIROC-ES2L	LGM	$+\mathrm{8000}$
MPI-ESM1-2	LGM	3850
MRI-ESM2-0	–	–
NESM3	–	–
NorESM1-F	MIS3	6000
NorESM2-LM	–	–

Code availability

Access to the Met Office Unified Model source code is available under licence from the Met Office at https://www.metoffice.gov.uk/research/approach/collaboration/unified-model/partnership .

Data availability

The climate model data (HadCM3) are available on request from http://www.bridge.bris.ac.uk/resources/simulations .

PaleoMIST 1.0, which includes the ice sheet margin, paleotopography, and ice thickness datasets and Stokes coefficients used in this study are available on Pangaea (10.1594/PANGAEA.905800, ).

Team list

Ayako Abe-Ouchi (The University of Tokyo, Kashiwa, Japan), Andreas Born (University of Bergen, Bergen, Norway), Nathaelle Bouttes (Laboratoire des sciences du climat et de l’environnement, France), Peter Ditlevsen (Niels Bohr Institute, University of Copenhagen, Denmark), Michael P. Erb (School of Earth and Sustainability, Northern Arizona University, Flagstaff, AZ, USA), Georg Feulner (Potsdam Institute for Climate Impact Research, Germany), Evan J. Gowan (Department of Earth and Environmental Sciences, Kumamoto University, Japan), Lauren Gregoire (University of Leeds, UK), Chuncheng Guo (NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, Bergen, Norway), Sandy P. Harrison (University of Reading, UK), Heather Andres (Memorial University of Newfoundland and Labrador, Canada), Masa Kageyama (Laboratoire des sciences du climat et de l’environnement, France), Marlene Klockmann (Helmholtz-Zentrum Hereon, Germany), Fabrice Lambert (Pontifical Catholic University of Chile, Chile), Allegra N. LeGrande (Columbia University, NASA Goddard Institute for Space Studies), Ute Merkel (MARUM – Center for Marine Environmental Sciences, University of Bremen, 28359 Bremen, Germany), Larissa S. Nazarenko (Columbia University, NASA Goddard Institute for Space Studies), Kerim H. Nisancioglu (University of Bergen, Norway), Kevin Oliver (University of Southampton, UK), Bette Otto-Bliesner (Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, Colorado), William R. Peltier (University of Toronto, Canada), Matthias Prange (MARUM – Center for Marine Environmental Sciences, University of Bremen, 28359 Bremen, Germany), Kira Rehfeld (Department of Geosciences, Department of Physics, University of Tübingen, Germany, and Institute of Environmental Physics, Heidelberg University, Germany), Alexander J. Robinson (Complutense University of Madrid, Madrid, Spain, and Instituto de Geociencias, CSIC-UCM, Madrid, Spain), Lev Tarasov (Memorial University of Newfoundland and Labrador, Canada), Paul J. Valdes (University of Bristol, UK), Guido Vettoretti (Niels Bohr Institute, University of Copenhagen, Denmark), Nils Weitzel (Department of Geosciences, University of Tübingen, Germany, and Institute of Environmental Physics, Heidelberg University, Germany), Qiong Zhang (Stockholm University, Sweden), Xu Zhang (State Key Laboratory of Tibetan Plateau Earth System, Resources and Environment (TPESRE), Chinese Academy of Sciences (CAS), Beijing, China, and Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany).

Author contributions

IMV compiled all tables. LCS and IMV wrote the first draft of this paper. LCS and IMV produced all figures. All authors contributed to the final draft.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

Evan J. Gowan is funded by an international postdoctoral fellowship from the Japan Society for the Promotion of Science. Bette Otto-Bliesner acknowledges funding by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. Xu Zhang acknowledges funding from NSFC (no. 42075047). Matthias Prange and Ute Merkel acknowledge support from the PalMod project (http://www.palmod.de, last access: 14 April 2023; FKZ 01LP1916C), funded by the German Federal Ministry of Education and Research (BMBF). Kira Rehfeld and Nils Weitzel acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), project no. 395588486, and the PalMod project (https://www.palmod.de/, last access: 13 October 2022), subproject no. 01LP1926C. Chuncheng Guo acknowledges funding from the Research Council of Norway under grant no. 325333 (ABRUPT).

Financial support

This project is TiPES contribution no. 123. It has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 820970.

Review statement

This paper was edited by Dominik Fleitmann and reviewed by two anonymous referees.