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

Dansgaard–Oeschger (D–O) events, millennial-scale climate oscillations between stadial and interstadial conditions (of up to 10–15 ^∘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 (Kindler et al., 2014; Huber et al., 2006), 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 ^∘C (Andersen et al., 2006; Kindler et al., 2014; Huber et al., 2006). D–O events are best documented during Marine Isotope Stage 3 (MIS3; between 27.8–59.4 kyr BP, hereafter ka; Goni and Harrison, 2010), including being recorded in several ice cores from Greenland (Fig. 1; Johnsen et al., 2001). Whilst the D–O events recorded in these cores are renowned, the events are global in nature (Voelker et al., 2002; Sanchez Goñi and Harrison, 2010; Sánchez Goñi et al., 2017), with known climate signatures including imprints in surface temperature and the hydrological cycle at high northern latitudes (Andersen et al., 2004; Thomas et al., 2009; Seierstad et al., 2014), in the tropics (Deplazes et al., 2013; Baumgartner et al., 2014; Adolphi et al., 2018), in Eurasia (Genty et al., 2003; Wang et al., 2008; Jacobel et al., 2017; Rousseau et al., 2017), and in North and South America (Wang et al., 2004; Wagner et al., 2010; Asmerom et al., 2010; Deplazes et al., 2013; Vanneste et al., 2015). 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 (Lang et al., 1999; Uriarte, 2019; Landais et al., 2004; Turner and Marshall, 2011; Barker et al., 2011; Lambert et al., 2012). This observational evidence shows that D–O events do not occur under interglacial or full Last Glacial Maximum conditions (Galaasen et al., 2014; Tzedakis et al., 2018; Galaasen et al., 2020).

In 2011, Valdes (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 (Brovkin et al., 2021). 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 (Brown and Galbraith, 2016; Galbraith and de Lavergne, 2019; Klockmann et al., 2018; Peltier et al., 2020; Armstrong et al., 2021; Zhang et al., 2021; Vettoretti et al., 2022), 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 (Ditlevsen and Johnsen, 2010). Previous studies have questioned the significance of the periodic occurrence of D–O events in MIS3 (∼1470) (Ditlevsen et al., 2007). If the full glacial period is included, the distribution of waiting time between D–O events is consistent with a random process (Ditlevsen et al., 2005). Durations of stadials vs. interstadials indicate correlations with global ice volume and orbital parameters (Lohmann and Ditlevsen, 2018; Mitsui and Crucifix, 2017), 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 (Barton et al., 2018; Tesi et al., 2021); sea ice loss in the area is linked with enhanced heat transport via an intensified throughflow, or “Atlantification” (Årthun et al., 2012; Polyakov et al., 2017). 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 (Lenton et al., 2008). 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 (IPCC, 2013; Eyring et al., 2016). The Paleoclimate Model Intercomparison Project 4 (PMIP4) is one of the individual model intercomparison projects which took part in CMIP6 (Kageyama et al., 2018). 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 1 and Fig. 1). 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 (Schulz et al., 1999) 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₂ measurements from Antarctic and Greenland ice cores (Bauska et al., 2021; NGRIP Project Members, 2004; Kindler et al., 2014). See also Table 1 for D–O nomenclature. The “cnt” red box indicates the period 38 to 32 ka proposed for the MIS3 baseline experiment.

(Meehl et al., 2007; Brovkin et al., 2021)(Lenton et al., 2008)(Lenton et al., 2008)(Brovkin et al., 2021)(Dansgaard et al., 1982; Johnsen et al., 1992)(Landais et al., 2004; Huber et al., 2006; Kindler et al., 2014)(Rasmussen et al., 2014)(Rasmussen et al., 2014)(e.g. Heinrich, 1988; Hemming, 2004)(Roche et al., 2004; Marcott et al., 2011)Rasmussen et al. (2014)(Barker et al., 2009)(Bond et al., 1992; Alley, 1998; Alley et al., 1999; Clark et al., 2007; Rousseau et al., 2022)Broecker (1994)Alley (1998)(Clark et al., 2007)

Table 1 Terminology.

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 A1, 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 A1); 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 A1). 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 (Drijfhout et al., 2013; Kleppin et al., 2015; Table A1). 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 (Sidorenko et al., 2015). Changes in convection also occur in the Kiel Climate Model (KCM; Martin et al., 2015); 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 (Peltier and Vettoretti, 2014). The initiation of the abrupt D–O-type warming events is associated with the opening of a large polynya over the Irminger Sea (Vettoretti and Peltier, 2016) (Table A1). The AMOC spontaneously exhibits D–O-like quasi-oscillations (Peltier et al., 2020). The Peltier et al. (2020) 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 (Kageyama et al., 2021a) have not been run long enough to be sure that such oscillations would not arise if they were run for longer (see Table B1). 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 (Huber et al., 2006; Galaasen et al., 2014; Kindler et al., 2014; Tzedakis et al., 2018).

D–O-type quasi-oscillations are also found in MIROC4m under mid-glacial conditions (Kuniyoshi et al., 2022). 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 (Kuniyoshi et al., 2022) (Table A1).

Under late-glacial conditions, at 30 ka, a quasi-oscillating AMOC is produced by the HadCM3 model (Armstrong et al., 2021) and results from a North Atlantic salt oscillator mechanism similar to that in UofT-CCSM4 (Peltier and Vettoretti, 2014; Vettoretti and Peltier, 2016; Peltier et al., 2020). The HadCM3 model also shows millennial-scale climate oscillations triggered by deglacial meltwater discharge in LGM simulations (Romé et al., 2022). Under intermediate glacial conditions (MIS3: 40–32 ka), the COSMOS model shows spontaneous millennial-scale climate oscillations triggered solely by orbitally driven insolation changes (Zhang et al., 2021). 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 (Zhang et al., 2021). 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₂ concentration of 180 ppm and low obliquity (22^∘) (Brown and Galbraith, 2016; Galbraith and de Lavergne, 2019) (Table A1). The MPI-ESM model exhibits more abrupt D–O-like quasi-oscillations with a present-day ice sheet distribution in combination with CO₂ concentrations ranging between 190–217 ppm (Table A1; Klockmann et al., 2018, 2020).

In contrast to the above, neither NorESM nor CCSM3 produces D–O-type events or quasi-oscillations under MIS3 conditions (38 ka) (Guo et al., 2019b; Zhang and Prange, 2020). The NorESM – with a reasonably simulated AMOC and Arctic sea ice distribution in the PI and historical simulations (as documented by Guo et al., 2019a) – 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₂ 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 (Guo et al., 2019b). Zhang and Prange (2020) use the LGM ICE-5G ice sheet configuration (Peltier, 2004), 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 (Galaasen et al., 2014; Tzedakis et al., 2018; Galaasen et al., 2020). 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: Kageyama et al., 2021a) shows spontaneous D–O-type oscillations (Tables A1 and B1; Kageyama et al., 2021a; Peltier and Vettoretti, 2014).

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 B1; Kageyama et al., 2021a; Armstrong et al., 2021; Klockmann et al., 2018), 16 simulations have been run for less than 2000 years (Table B1), 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 C1). 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 C1).

2.1 The role of ocean and sea ice feedbacks

Changes in the AMOC are crucial to the correct simulation of D–O events (Broecker and Peteet, 1985). 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. Rahmstorf, 2002; Clark et al., 2002). As a consequence, there is a large transport of heat northward across the hemispheres (e.g. Feulner et al., 2013; Buckley and Marshall, 2016), 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. Klockmann et al., 2018). 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 3 Schematic depicting the transition from GI to GS conditions, i.e. a D–O cooling event.

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 (Li et al., 2005, 2010). 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 (Rasmussen and Thomsen, 2004; Singh et al., 2014; Dokken et al., 2013). 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 (Rahmstorf, 1995; Huisman et al., 2010). The interaction of subpolar and tropical salinity anomalies at the surface and in the sub-surface (Jackson and Vellinga, 2013) and possible roles of the intertropical convergence zone and freshwater export through the Fram Strait are also important in D–O-related salinity feedbacks. Klockmann et al. (2018) 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₂ concentrations (Brown and Galbraith, 2016; Klockmann et al., 2018; Guo et al., 2019b; Muglia and Schmittner, 2015; Sherriff-Tadano et al., 2018), 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 (Lynch-Stieglitz, 2017), with the upward mixing of heat from the sub-surface and importation of salt from the tropical Atlantic via gyre mechanisms (Peltier and Vettoretti, 2014).

The importance of vertical (diapycnal) mixing in the ocean for these long-timescale, D–O-type instabilities has long been recognised (Welander, 1982). However, we note that the different ocean and climate models (Table A1) parameterise diapycnal mixing in very different ways (e.g. Nilsson et al., 2003; de Lavergne et al., 2019). 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. Holden et al., 2010). This is partly because wind-driven Southern Ocean upwelling plays a complementary role in diapycnal mixing in setting the steady-state overturning (Samelson, 2004) and partly because surface buoyancy forcing controls the relative strength of the upper (Atlantic) and lower (Antarctic) overturning cells (Oliver and Edwards, 2008). 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 2 and 3 show some of the key states, processes, and ocean–sea ice feedbacks that enable D–O events. Following Lohmann and Ditlevsen (2019), D–O events can be broken down into four periods: (1) cold-stadial state (Fig. 2a), (2) rapid-warming phase governed by very-fast-timescale mechanisms (Fig. 2b), (3) warm-interstadial state (Figs. 2c and 3a), and (4) gradual-cooling phase (Fig. 3b) followed by a faster abrupt transition into a cold-stadial phase (Fig. 3c). For some of the D–O events, the magnitudes of the warming transitions are on the order of 10 ^∘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) (Lohmann and Ditlevsen, 2019). This picture of rapid retreat of North Atlantic sea ice (Spolaor et al., 2016; Dokken et al., 2013) 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 (Li et al., 2005, 2010; Sime et al., 2019).

2.2 The role of Northern Hemisphere ice sheets

Section 2 and Table A1 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 (Zhang et al., 2014; Klockmann et al., 2018; Brown and Galbraith, 2016; Muglia and Schmittner, 2015; Sherriff-Tadano et al., 2018). Figures 4 and 5 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 5 Schematic of a potentially D–O-type oscillating state.

The Northern Hemisphere Eurasian ice sheet was most probably limited to mountainous areas during mid-MIS3 (Helmens, 2014; Hughes et al., 2016), 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 (Abe-Ouchi et al., 2015) 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 (Li and Born, 2019; Zhang et al., 2014).

A larger LIS (especially its height) causes stronger Northern Hemisphere winds (Li and Battisti, 2008; Pausata et al., 2011; Hofer et al., 2012; Ullman et al., 2014; Löfverström et al., 2014; Merz et al., 2015) as well as an amplified stationary wave over North America (Manabe and Broccoli, 1985; Cook and Held, 1988), and it causes the North Atlantic glacial jet to be more stable due to differences in wave–mean flow feedbacks (Riviere et al., 2010) and alters variability in the large-scale atmospheric circulation, especially in the North Atlantic (Justino and Peltier, 2005; Pausata et al., 2009; Riviere et al., 2010). In addition, LIS height could control the sea ice coverage and gyre circulation by shifting the westerlies over the North Atlantic region (Zhang et al., 2014).

LIS altered winds that have wide implications for D–O-relevant tipping elements (Seager and Battisti, 2007; Wunsch, 2006). Li and Born (2019) 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 (van der Schrier et al., 2010) and favours episodes of Greenland blocking (Madonna et al., 2017). 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 (Li and Born, 2019). Given that at latitudes north of about 45^∘ 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. Jungclaus et al., 2013); 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₂ concentrations (Brown and Galbraith, 2016; Klockmann et al., 2018; Guo et al., 2019b; Muglia and Schmittner, 2015; Sherriff-Tadano et al., 2018; Zhang et al., 2014). Similarly, Zhang et al. (2014) 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 (Zhang et al., 2014).

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₂ concentrations). That said, the AMOC in many LGM simulations is likely too strong (Klockmann et al., 2018; Kageyama et al., 2021b). Thus the AMOC is far away from a tipping point with LGM-size ice sheets for many models (Zhang et al., 2014; Klockmann et al., 2018; Guo et al., 2019b).

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 (Armstrong et al., 2021). 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 (Klockmann et al., 2018; Armstrong et al., 2021). Thus, with intermediate MIS3 LIS, i.e. reduced in its height compared to the LGM, multiple AMOC states are more likely (Zhang et al., 2014; Kawamura et al., 2017; Zhang and Prange, 2020; Armstrong et al., 2021; Klockmann et al., 2018).

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 (Guo et al., 2019b; Zhang and Prange, 2020) when models are run with intermediate or low MIS3 CO₂ values and ice sheets, i.e. reduced in size compared to the LGM (Brown and Galbraith, 2016; Kawamura et al., 2017; Klockmann et al., 2018; Zhang and Prange, 2020; Galbraith and de Lavergne, 2019; Zhang et al., 2014; Vettoretti et al., 2022), 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, Zhang et al. (2021) 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 (∼22.6^∘) leads to D–O-type behaviour (see Fig. 2 from Zhang et al., 2021). Additionally, the MIROC4m model produces D–O-type oscillations (under mid-glacial conditions) and low obliquity (22.9^∘) (Kuniyoshi et al., 2022). 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. 1) and (2) has the ideal intermediate MIS3 ice sheet configuration conducive to generating D–O-type quasi-oscillations (Sect. 2).

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. 2) suggests that this should be a minimum of 5000 years (Peltier and Vettoretti, 2014; Kleppin et al., 2015; Sidorenko et al., 2015; Brown and Galbraith, 2016; Klockmann et al., 2018, 2020). 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 Klockmann et al., 2016, 2018) or spinning up from present-day observations (as done in Guo et al., 2019b). 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 (Zhang et al., 2013; Knorr et al., 2021), 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. 6) and an ice sheet configuration as outlined below (Sect. 3.1 and 3.2).

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 (Eyring et al., 2016): 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 (Barker and Knorr, 2021). Being close to or within the sweet spot, the AMOC is characterised by high sensitivity to transient and/or noisy climatic forcing (Zhang et al., 2014; Lohmann and Ditlevsen, 2018) or by self-oscillating behaviours (Zhang et al., 2021).

• 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⁻²).

3.1 Atmospheric trace gases

MIS3 atmospheric CO₂ values varied between a maximum of ∼ 233 ppm and a minimum of ∼187.5 ppm (Table 2, Fig. 1; Bauska et al., 2021). 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 (Bauska et al., 2021). GHG forcing is critical to model stability. Low (LGM-like) to intermediate (MIS3) CO₂ concentrations tend to be associated with abrupt D–O-type AMOC transitions in models (Sect. 2; Klockmann et al., 2018; Zhang et al., 2017, 2014; Brown and Galbraith, 2016; Vettoretti et al., 2022). We thus suggest performing the MIS3-cnt experiment using the GHGs values specified in Table 2 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 (Kleman et al., 2010; Batchelor et al., 2019) 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 (Gowan et al., 2021). 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 − 30 and − 50 m is inferred (Dalton et al., 2022). For much of MIS3, since the Eurasian ice sheets and the Cordilleran Ice Sheet were likely restricted to mountain-based caps (Helmens, 2014; Hughes et al., 2016; Clague and Ward, 2011), the primary control on ice volume is assumed to be from the LIS. Recent work in the area of the Hudson Bay (Dalton et al., 2016, 2019; McMartin et al., 2019; Dalton et al., 2022) suggests that ice-free conditions may have occurred during mid-MIS3. This implies climatic conditions in this region similar to present (Dalton et al., 2017) and a LIS margin removed from the southern Hudson Bay. Similarly Tarasov et al. (2012) show a considerably lower and less extensive LIS compared to ICE-5G and ICE-6G (Peltier, 2004; Peltier et al., 2015) LGM ice sheet reconstructions. Pico et al. (2017) sea level curves are consistent with the estimated MIS3 ice sheet volumes from Batchelor et al. (2019). Using glacial isostatic adjustment modelling, Pico et al. (2017) also show that a small LIS can explain high MIS3 sea level estimates alongside the eastern coast of the United States. The synthesis of Batchelor et al. (2019) 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) Gowan et al. (2021). Also shown for comparison is the (c, d) LGM ICE-6G ice sheet reconstruction from Peltier et al. (2015). The third row shows the differences between the MIS3 ice sheet reconstruction and the LGM ICE-6G reconstruction.

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 Gowan et al. (2021). This reconstruction is based on trying to fit the evolution of ice flow indicators, as well as chronological constraints of ice-free conditions.

Gowan et al. (2021) 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 (Gauthier and Hodder, 2020), 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 (Andrews and Voelker, 2018), 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 Dalton et al. (2019). The Cordillera Ice Sheet extent is based on evidence of relatively restricted ice cover during MIS3 (Clague and Ward, 2011). 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 (Hughes et al., 2016). 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 (Larter et al., 2014).

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

(Bauska et al., 2021)(Loulergue et al., 2008)(Schilt et al., 2010)(Berger et al., 1998)(Berger et al., 1998)(Berger et al., 1998)(Gowan et al., 2021)

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

Whilst the implementation of the ice sheet will differ between models, the steps of Kageyama et al. (2017) 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, Kageyama et al. (2017) 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) (Guo et al., 2019b).

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 (Kageyama et al., 2017). 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 Peltier and Vettoretti (2014) 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 Peltier et al., 2020). 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 (Peltier et al., 2020). This event may resemble the impact of a Heinrich-event-like kick to the AMOC, though no freshwater perturbation was imposed (Peltier et al., 2020).

In a more recent study, Pedro et al. (2022) examine the CCSM4 simulations that show unforced D–O-type oscillations (Peltier and Vettoretti, 2014; Vettoretti and Peltier, 2018), but with the addition of a (freshwater) Heinrich-like event. (Pedro et al., 2022). 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 (Pedro et al., 2022). 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 (Pedro et al., 2022).

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. Stocker and Johnsen, 2003). Estimates of the meltwater input into the North Atlantic during Heinrich events range between 2 and 15 m of sea level equivalent ice volume (Hemming, 2004; Chappell, 2002; Rohling et al., 2004; Roche et al., 2004; Roberts et al., 2014b; Siddall et al., 2008; Grant et al., 2014). It is logical to presume that these freshwater events are important in preconditioning the climate system with respect to D–O behaviour (Peltier and Vettoretti, 2014; Peltier et al., 2020; Pedro et al., 2022). Freshwater perturbations can trigger changes between AMOC states (e.g. Ganopolski and Rahmstorf, 2001; Timmermann et al., 2003; Stouffer et al., 2006; Hu et al., 2008; Kageyama et al., 2013; Roberts et al., 2014a; Sime et al., 2019), and a relatively small freshwater flux applied over convection areas can lead to a shutdown of the AMOC (e.g. Roche et al., 2010). Studies have shown similarities between observed global features of abrupt D–O changes and the behaviour seen in freshwater forcing (FWF) experiments (Liu et al., 2009; Menviel et al., 2014). 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 (Ganopolski and Rahmstorf, 2001; Timmermann et al., 2003; Stouffer et al., 2006; Hu et al., 2008; Kageyama et al., 2013; Roberts et al., 2014a; Zhang et al., 2014). 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 125×10⁴ km³ (MacAyeal, 1993; Dowdeswell et al., 1995; Marshall and Clarke, 1997; Hulbe, 1997), whilst isotope-based estimates and a precipitation balance approach have yielded 86, 649, and 946×10⁴ km³ of ice volume (Hemming, 2004; Roche et al., 2004; Levine and Bigg, 2008). Roberts et al. (2014b) and using a sediment modelling approach, estimated a discharge of 30 to 120×10⁴ km³ 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 (Gasson et al., 2016; Hemming, 2004; Roberts et al., 2014b). 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 (Hemming, 2004) and others suggesting that a duration of 500 years is more typical (Roberts et al., 2014b). 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 (Sime et al., 2019; Zhang and Prange, 2020). 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₂ 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

Drijfhout et al. (2013)Kleppin et al. (2015)Sidorenko et al. (2015)Martin et al. (2015)Peltier and Vettoretti (2014); Vettoretti and Peltier (2016, 2018)Vettoretti et al. (2022)Armstrong et al. (2021)Zhang et al. (2021)Brown and Galbraith (2016)Klockmann et al. (2018, 2020)Peltier and Vettoretti (2014)Montoya et al. (2011)Drijfhout et al. (2013)Kleppin et al. (2015)Kuniyoshi et al. (2022)Li and Born (2019)

Table A1 Models that exhibit spontaneous oscillations.

Guo et al. (2019b)Zhang and Prange (2020)(Peltier, 2004)Kawamura et al. (2017)Vettoretti et al. (2022)Zhang et al. (2014, 2017)Zhang et al. (2014, 2017)Zhang et al. (2014, 2017)Zhang et al. (2014, 2017)

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

Appendix B

Peltier and Vettoretti (2014)Lohmann et al. (2020)Sidorenko et al. (2019)Tierney et al. (2020)Valdes et al. (2017)Lhardy et al. (2021)Volodin et al. (2018)Sepulchre et al. (2020)Ohgaito et al. (2021)Mauritsen et al. (2019)Armstrong et al. (2021)Zhang et al. (2021)Guo et al. (2019b)Zhang and Prange (2020)Kawamura et al. (2017)Kuniyoshi et al. (2022)Vettoretti et al. (2022)Brown and Galbraith (2016)Klockmann et al. (2018)Zhang et al. (2014)

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

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

Appendix C

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