Introduction

In a seminal paper published in 1976 Hays, Imbrie and Shackleton1 provided the first convincing geologic evidence in support of Milankovitch theory2, which posits a direct link between Earth’s climate and changes in the geometry of our orbit around the Sun. One of the most remarkable observations made by Hays et al. was that the dominant period of climate variability over the last half million years or so (the glacial cycles of the Late Pleistocene) is close to 100 kyr (Fig. 1), even though the forcing at this period (changes in the circularity of Earth’s orbit, eccentricity) is negligible when compared to that related to obliquity (tilt of the rotational axis) or precession (timing of solstice relative to orbital position; Fig. 1). Whether glacial cycles are driven ultimately by changes in eccentricity or more likely by some combination of obliquity and precession3,4,5 (which itself is dependent on eccentricity) the question remains as to why glacial terminations should occur every ~100 kyr and not more frequently as they did during the early Pleistocene6. A clue to this question lies in the asymmetric nature of records of benthic foraminiferal δ18O (a crude proxy for continental ice volume; Fig. 1), which suggest a stepwise build-up of ice sheets over a period of ~90 kyr, followed by a shorter interval (~10 kyr) of relatively fast ice sheet decay during termination7.

Fig. 1: The 100 kyr problem.
figure 1

The geologic record of benthic foraminiferal δ18O65 contains significant spectral power in the eccentricity band (as well as those of precession and obliquity) while the external insolation forcing at this period (~100 kyr) is negligible. Late Pleistocene glacial cycles are ‘terminated’ by relatively short deglacial periods, giving rise to a characteristic saw-tooth character. T1–T5 represent Terminations 1 to 5.

Accordingly, most models of glacial-interglacial (G-IG) variability have called on the long-term build-up of continental ice sheets, and their propensity to collapse beyond some critical threshold, in order to explain the long-lived nature of late Pleistocene glacial cycles e.g. refs. 3,4,8,9,10,11,12. By this reasoning, only once ice sheets have exceeded some critical size, do they become susceptible to a potentially modest (external) forcing to which they respond in a non-linear way (Box 1). Indeed, it has long been appreciated that glacial terminations must involve non-linearities within the climate system in order to amplify the rather modest forcing from insolation1,7,13.

The specific suggestion that abrupt reorganisations of the ocean-atmosphere system might play an active role in the mechanism of termination was introduced by Mix et al.14 (although Ewing and Donn15 had predicted the importance of ocean circulation on the waxing and waning of continental ice sheets several decades earlier) and augmented to include a link to atmospheric CO2 by Broecker and Denton16. The most recent incarnation of this idea was summarised by Denton et al.17 as a sequence of events initiated by the inflow of freshwater to the surface North Atlantic (generated by summer melting of continental ice sheets) suppressing ocean circulation in the Atlantic and ultimately causing the release of CO2 to the atmosphere. At some critical point, the rise in CO2 would pass a threshold beyond which interglacial conditions could be stabilized and maintained17,18,19.

A number of recent studies provide evidence to support Broecker and Denton’s16 original proposition, with several independent lines of proxy evidence emphasising the ubiquitous association between glacial termination and millennial-scale oscillations in ocean circulation, in particular the Atlantic Meridional Overturning Circulation, AMOC20,21,22,23,24,25,26,27. Moreover, ice-core measurements and quantitative carbon cycle models support an active role for abrupt shifts in ocean circulation within the mechanism of glacial termination. For example atmospheric CO2 measurements in Antarctic ice cores combined with proxy reconstructions in marine sediments suggest that CO2 rises on a millennial timescale whenever Atlantic Ocean circulation is in a substantially weakened state25,28,29,30 with additional century-scale increases (at least occasionally) associated with abrupt transitions between weak and strong modes of ocean circulation23,31,32. Modelling studies further suggest that such changes in CO2 may be associated mechanistically with the corresponding variations in ocean circulation, either directly or indirectly through e.g. biological activity33,34,35,36,37,38,39.

However, while the simple model described by Denton et al.17 provides an appealing explanation for the magnification of insolation forcing related to deglaciation it leaves open several questions related to the interplay between millennial and orbital-timescale changes and the anatomy of terminations themselves, for example:

  • Is the temporal correspondence between abrupt shifts in ocean circulation and glacial termination merely coincidence or are there systematic interactions between variations on these timescales that necessarily come into play during deglaciation?

  • How do we explain the occurrence of terminal interruptions like the Bølling-Allerød/Younger Dryas (B/A-YD) oscillation during Termination 1 (T1) while other terminations proceed apparently without such interruption?

  • Why do some terminations end with such high levels of atmospheric CO2 while others (notably the most recent) do not?

Below we address each of these questions in turn, highlighting new results and outstanding uncertainties requiring future research focus. We begin by investigating the occurrence of millennial-scale climate variability with respect to orbital-timescale changes in ice volume and atmospheric CO2.

Feedbacks between millennial and orbital-timescale variations

McManus et al.40 first described the pervasive nature of millennial-scale climate variability (Box 2) over the past 500 kyr and suggested that such variability was most pronounced when continental ice volume lay somewhere between full interglacial and extreme glacial values (Box 1); millennial activity was somehow dampened when ice sheets were relatively small (similar to modern) or particularly large. Since then, numerous studies have confirmed that millennial-scale climate variability experiences its greatest amplitude within intermediate climate states and transitions between states (notably terminations)25,27,41,42,43,44.

While geological observations suggest that there exist critical boundary conditions between which millennial-scale activity is enhanced, climate models can help us to identify which climatic parameters are most critical. Recent modelling studies have made some headway into this problem; Zhang et al.45 suggested that the large size of continental ice sheets could indeed be responsible for the (mono) stability of AMOC during full glacial periods, while high (interglacial) concentrations of atmospheric CO2 seem to be responsible for an equivalent stability during interglacial periods, a suggestion also supported by more recent studies46,47. When both parameters occupy the space between full glacial and interglacial conditions a ‘window of opportunity’ arises (Box 1) within which the AMOC exists in a bistable regime and may therefore be expected to experience more frequent and or pronounced variability. Abrupt variations can also occur within full glacial or interglacial conditions when we expect the AMOC to be monostable e.g. Heinrich Event 2 during the LGM or the 8.2ka event in the Early Holocene, but these are most likely transient perturbations forced by distinct events such as freshwater floods e.g. ref. 48 (Box 1).

Thus, millennial-scale variability seems to be dependent on orbital-timescale changes. Indeed, previous studies implied and demonstrated that gradual global warming can directly trigger abrupt transitions in the AMOC49,50,51,52, and it has now been shown that both ice sheet height and atmospheric CO2 can influence the strength of AMOC e.g. refs. 46,53,54,55,56 and the occurrence of abrupt changes45,47,57,58,59. Moreover, it is thought that abrupt variations in ocean circulation can have a reciprocal influence on CO2 (through resultant biophysicochemical changes33,34,35,36,37,38) and possibly sea level (although the actual relationship between variations in AMOC and ice volume remains controversial60,61), highlighting the possibility that millennial-scale variability may be considered both slave and master to orbital-timescale changes, an idea we develop below.

Records of changing surface ocean conditions in the North Atlantic can be used to illustrate the interplay between millennial and orbital-timescale variability over the past 800 kyr (Fig. 2). Based on proxy evidence from the last glacial period22,30 we use variations in North Atlantic SST to infer changes in the AMOC. We quantify millennial power using a Hilbert transform (amplitude envelope) of the 0.5–7 kyr bandpass filtered %NPS (percentage of the polar affiliated foraminifera, Neogloboquadrina pachyderma, within the total assemblage) record from ODP Site 98325 (Fig. 2b). As expected, we observe enhanced power in the millennial-band associated with intermediate climate conditions (the window of opportunity, constrained here using benthic δ18O, atmospheric CO2 and relative sea level; Fig. 2e).

Fig. 2: The sweet spot of millennial variability.
figure 2

a Records of benthic δ18O65, atmospheric CO297, sea level93 and Antarctic δD98. Yellow bars are same as for cf. b Records from ODP Site 98325: %NPS (a proxy for sea surface temperature), millennial power (Hilbert transform of the 0.5–7 kyr Taner bandpass filtered %NPS record) and Ice Rafted Debris (IRD) accumulation. c Binned mean absolute rates of change for the records in a. d Distribution histograms of the records in a. e Mean millennial power for the same bins utilized in c, d. f Mean IRD accumulation rates for the same bins utilized in c, d. In each panel fainter colours represent deglacial intervals (terminations) while bolder colours represent all other times. Yellow bars highlight the respective peaks in millennial power in e (the sweet spot as described in the text and Box 1) and grey shaded regions highlight the broader windows of opportunity.

Less expected is the observation that the peak in millennial power (the sweet spot; Box 1) is associated with a local minimum in the distribution histograms of both benthic δ18O and atmospheric CO2 (Fig. 2d; the bimodal nature of CO2 variability was noted previously62). We suggest that these minima are manifestations of the reciprocal influence of millennial activity on the boundary parameters themselves. Specifically we propose that the feedbacks between abrupt changes in the AMOC and the global parameters of benthic δ18O and CO2 are sufficiently strong when millennial variability is at its maximum that the climate system is effectively forced away from the sweet spot and towards a more stable climate regime. If this is case then the occurrence and active role of abrupt climate shifts during glacial termination (and inception for that matter) is no mere coincidence, but rather represents the necessary and ubiquitous interaction between these two distinct modes of climate variability.

For example, atmospheric CO2 experienced a long-term decrease from the last interglacial period (Marine Isotope Stage, MIS 5e) towards the LGM (MIS 2; grey headed arrow in Fig. 3), taking it from high concentrations (during which the AMOC was monostable strong and millennial activity was weak) to within the window of opportunity (during MIS 5d) and eventually the sweet spot (around the beginning of MIS 5b). According to our hypothesis at that point the feedbacks between millennial and orbital-timescale changes became more powerful and CO2 was forced upwards and away from the sweet spot (a negative feedback; black headed arrow in Fig. 3). However, CO2 continued its long-term decrease until MIS 3, during which it varied about an intermediate value of ~220 ppmv (within the sweet spot as indicated by our analysis; Fig. 2e ii). Indeed, paleo-reconstructions confirm (or at least suggest) that large and abrupt fluctuations in ocean circulation occurred repeatedly within MIS 322,30 and these are thought to have strongly influenced atmospheric CO2, causing it to rise when circulation was particularly weak (e.g. during Heinrich events) and vice versa25,36. Thus, the effect of accentuated millennial-scale activity at that time was to continually drive CO2 in one direction or another depending on the mode of ocean circulation (i.e. feedbacks could be both negative and positive). On the other hand, the resulting changes in CO2 may themselves have reversed the oceanic conditions that produced them in the first place (e.g. ref. 45) thus forcing CO2 to remain within range of the sweet spot throughout MIS 3. Note that we are not suggesting that all D-O events were caused by changes in CO2 (or vice versa) and while this may have been the case for some of the larger events, we believe that high frequency D-O variability represents the inherent tendency of the AMOC to oscillate from one mode to another as an inevitable consequence of an intermediate climate background63.

Fig. 3: Feedbacks between millennial and orbital-timescale variability.
figure 3

From top to bottom: Benthic δ18O65 (filled regions indicate monostable AMOC as a function of maximum ice volume); Atmospheric CO297 (filled regions indicate monostable AMOC as a function of high CO2); Greenland ice-core δ18O99 with synthetic record27 beyond 110ka; %NPS (temperature proxy) from ODP Site 983; Millennial activity at Site 983 (see text). Pink and yellow shaded boxes represent the window of opportunity and sweet spot respectively (see Fig. 2). Grey headed arrows represent orbital-timescale trends in atmospheric CO2; Black (white) headed arrows represent examples of negative (positive) feedback of millennial-scale variability on orbital changes. Feedbacks are strongest within the sweet spot. Dashed vertical lines indicate timing of AMOC recovery associated with glacial terminations T1 and T2. Lower numbers are Marine Isotope Stages.

CO2 declined further throughout MIS 3 and eventually ice sheets attained their maximum configuration during the LGM (and equivalently MIS 6), inducing a period of reduced millennial-scale activity associated with a monostable strong AMOC. Subsequently, during glacial termination, when climate again passed through the sweet spot, the effect of accentuated millennial-scale activity was to help push CO2 straight through the window of opportunity (a positive feedback; white headed arrow in Fig. 3), at least in the case of T2. Early recovery of the AMOC associated with the B/A during T1 (Fig. 4) slowed the rise of CO2 (a negative feedback) but the inevitable return to a weak mode of circulation during the YD (see sub-section: ‘solutions in the time domain’) enabled CO2 to increase beyond the window of opportunity with the onset of interglacial conditions, eventually returning to a monostable strong mode of AMOC during the Holocene (and MIS 5e). Notably, the occurrence of high amplitude millennial-scale variability associated with high (interglacial-like) CO2 is a peculiar feature of deglacial periods (Fig. 2e ii), which we ascribe to the influence of abrupt change on CO2 as opposed to the influence of CO2 on the occurrence of abrupt change (see also penultimate section).

Fig. 4: The last two terminations.
figure 4

From top to bottom: Deglacial timeseries of atmospheric CO231,97, benthic δ18O (black curve65) and sea level reconstructions (orange100; green92; blue93), NW Atlantic sedimentary εNd21,23 (sensitive to ocean circulation change), GLT_syn_hi27 (a low order approximation of AMOC strength25) and %NPS (a temperature proxy) from the NE Atlantic68. Blue boxes are cold (weak AMOC) intervals. T1, T2 are Terminations 1 and 2; MIS 1, MIS 5e are Marine Isotope Stages; B/A Bølling/Allerød, YD Younger Dryas, HS1, HS11 are Heinrich Stadials, within which Heinrich events, HE1 and HE11 occurred. SLE is Sea Level Equivalent.

As stated earlier, we are using variations in North Atlantic SST recorded at a single core site to infer large-scale changes in the AMOC and we admit that use of alternative records might give different results. For example, the coincidence of our observed sweet spot with minima in the histograms of CO2 and benthic δ18O could be mere coincidence and needs to be confirmed using alternative proxies and locations. Furthermore, the minima themselves are not particularly prominent, and could be a result of controlling factors other than millennial-scale AMOC variability. On the other hand, we observe a maximum in the absolute rate of change of Antarctic temperature (also proposed to reflect changes in the AMOC27,64) aligned with the maximum in millennial power at Site 983 (Fig. 2c iv), which also implies its coincidence with the local minima in CO2 and δ18O. This supports our contention that millennial-scale variability at Site 983 reflects variations in the wider climate system. We also note that the absolute rate of change of CO2 itself reaches a maximum (Fig. 2c ii) in line with the minimum in its histogram (Fig. 2d ii) and the sweet spot of millennial activity i.e. CO2 changes fastest when millennial activity is greatest, which we contend results in the local minimum in the histogram of CO2.

We also observe a peak in millennial power associated with intermediate ice volumes (−40 m Sea Level Equivalent, SLE; Fig. 2e iii). However, we do not find a corresponding minimum in the histogram of sea level. Because benthic foraminiferal δ18O (in this case the stacked record of ref. 65) represents a convolution of ice volume and deep ocean temperature this could mean that the local minimum in the histogram of benthic δ18O reflects the existence of two distinct modes in deep ocean temperature, which could correspond to the two equivalent modes in atmospheric CO2 (but note that such bimodality is not observed in Antarctic temperature; Fig. 2d iv), rather than ice volume alone. Alternatively, it could reflect the very large envelope of uncertainty associated with current sea-level reconstructions i.e. perhaps a local minimum in ice volume is obscured by inaccurate sea-level reconstructions. We return to this issue in the final section.

Finally, our analysis can inform us about the interactions between calving ice sheets around the NE Atlantic and millennial-scale climate variability. There has been some debate over whether or not iceberg calving events cause or amplify abrupt changes in ocean circulation (which can be inferred from changes in high northern latitude sea surface temperature)66,67,68. In a previous study we suggested that icebergs arrived too late to trigger cooling across the NE Atlantic. Here we use an extension of that dataset to show that the peak in millennial variability in SST at this site is misaligned with the peak in ice rafted debris accumulation, which we assume provides a first order indication of iceberg calving (e.g. compare Fig. 2e, f). In fact, with respect to all three parameters (benthic δ18O, CO2 and sea level) we note that ice rafting increases away from the sweet spot of millennial activity, reaching its maximum within the glacial mode of CO2 and δ18O and when ice volume exceeds −50 m SLE. While this observation does not negate the potentially significant influence of icebergs on SST (and ocean circulation), the fact that millennial variability in SST reaches its maximum power (maximum amplitude in the 0.5–7 kyr band) when ice rafting is still relatively subdued suggests that iceberg calving is not the main influence on millennial-timescale climate variability.

The Bølling-Allerød/Younger Dryas oscillation

The abrupt climatic oscillation represented by the HS1-B/A-YD triptych during Termination 1 (Box 2; Fig. 4) has led to much debate as to whether or not such features might characterise glacial terminations in general and further, whether their occurrence reflects spatiotemporal variations in meltwater input to the North Atlantic or a fundamental change in the underlying dynamics of AMOC stability during deglaciation69,70,71,72,73. Based on the notion that abrupt oscillations in the AMOC are more pronounced when ice volume is intermediate between full glacial and interglacial conditions (their so-called ‘D-O window’, which is similar to our ‘window of opportunity’), Sima et al.71 argued that such oscillations represent a change in the underlying dynamical behaviour of the AMOC and are an intrinsic feature of deglacial transitions which by definition must pass through the D-O window. Indeed, based on an array of evidence24,25,26,27,69 we are now confident that abrupt mode switches in the AMOC really are a characteristic feature of glacial terminations. However, such a strong feature as the B/A-YD transition (i.e. a return to near glacial conditions following major deglacial warming) seems to be a rare occurrence25, such that it is still not clear whether it really reflects a systematic change in AMOC stability or e.g. a chance rerouting of meltwater runoff specific to T172,74. Here we investigate the changing stability of the AMOC with respect to key boundary parameters during deglaciation, and how this might interact with freshwater perturbations to give rise to the various permutations represented by previous terminations.

A moveable window of AMOC bistability

Many previous studies have pointed to the existence of AMOC bistability with respect to various boundary conditions45,49,57,75,76. Moreover, it has been shown that the region of bistability with respect to one parameter can vary, depending on the value of another (Box 1). For example, Zhang et al.45 showed that one effect of increasing ice sheet height is to push the window of AMOC bistability with respect to CO2 towards lower concentrations of CO2. If ice sheets get too large (i.e. similar to glacial maximum conditions) then the region of bistability no longer exists for the observed range of G-IG CO2 variability and the AMOC enters a region of monostability. Likewise, for increasing values of CO2, the window of bistability with respect to ice sheet height shifts towards smaller and smaller ice sheets until it disappears (monostability) as CO2 approaches interglacial values (Fig. 5).

Fig. 5: Moveable window of AMOC bistability.
figure 5

Schematic stability diagram for AMOC strength as a function of continental ice volume57 for different levels of atmospheric CO245 (CO2 (Glacial) < CO2 < CO2 < CO2+ < CO2 (Interglacial)). Varying atmospheric CO2 leads to structural changes in AMOC stability with respect to ice volume (and vice versa) e.g. increasing CO2 pushes the window of bistability with respect to ice volume to the right. Yellow diamond represents monostable strong glacial mode of AMOC (a function of large ice sheets) while yellow square represents monostable strong interglacial mode of AMOC (a function of high CO2). Other symbols and corresponding letters are referred to in the text.

In Fig. 5 we demonstrate how instantaneous transitions in the AMOC could be invoked by structural changes in AMOC stability with respect to one parameter, as another parameter is varied. Beginning in a glacial mode in which the AMOC is monostable strong45,57 (solid yellow diamond in Fig. 5) we decrease ice sheet height (while keeping CO2 constant) until we cross the grey square to reach 1b. At this point the system is bistable and a transition to a stable weak mode could be invoked by a transient forcing, for example a one-time input of freshwater. As ice sheet height continues to decrease beyond point 1d we leave the window of bistability as the white square is crossed and the AMOC (if still in a strong mode) will experience an abrupt transition to a weak mode (point 2a). All else being equal, ice sheets would then have to grow again by ~18m (SLE)57 in order to overcome the inherent hysteresis and transition back into a monostable strong mode (i.e. to cross the white triangle back to point 1a). Alternatively, if CO2 were to increase simultaneously (CO2+ in Fig. 5) then the transition back to monostable strong would occur earlier (i.e. point 2d →1b) as the window of bistability shifts to the right. On the other hand, from bistable strong point 1d an abrupt transition to a monostable weak mode (point 2b) would occur if CO2 decreased (CO2) even if ice sheets were kept constant.

Implications for the sweet spot

The considerations above have implications for the occurrence of abrupt shifts in the AMOC within an intermediate climate state and when the AMOC is close to a bifurcation point (which we consider to be more likely within the window of opportunity of intermediate CO2 and ice sheet size). As discussed, the region of AMOC bistability with respect to either CO2 or ice sheet height is sensitive to changes in the other parameter45. This means for example that starting in a weak mode (point 2c in Fig. 5) a weak-strong transition could be triggered by a relatively modest increase in ice sheet height (2c→2d), combined with a modest increase in CO2 (shifting the bistable window to the right from ‘CO2’ to ‘CO2+’), resulting in a shift to monostable strong (2d→1b). Likewise, a strong to weak mode transition could be triggered by e.g. the combination of ice decline (1c→1d) plus a CO2 decrease (1d →2b). Hence a mixture of ice volume and CO2 changes that could come about as a result of millennial-scale changes in the AMOC could themselves increase the likelihood of further abrupt changes in the mode of AMOC and thence further changes in CO2 and or ice volume (cf. MIS 3; see previous section; Fig. 3). Such an interpretation could provide a dynamical key to link the sweet spot in the millennial-power band to the local minimum in the distribution histograms of benthic δ18O and atmospheric CO2 (Fig. 2). In this sense, feedbacks between abrupt AMOC changes and CO2 and ice volume might become so effective within the sweet spot that the climate system is forced away from this region as soon as it is entered (Fig. 3).

Implications for deglaciation

We have already discussed the possibility that exaggerated feedbacks between ocean circulation and atmospheric CO2 (and possibly ice volume) might help to amplify the pace of deglaciation as climate passes through the sweet spot for millennial-scale variability. But the notion of a moveable window of AMOC bistability has more fundamental implications for the structure of glacial terminations in general, which we discuss below.

The last deglaciation (T1; Fig. 4) was characterized by a pronounced weakening of AMOC during its earliest phase (Heinrich Stadial 1, HS1)20. This event is thought to reflect (at least in part) the addition of freshwater to the surface North Atlantic by decaying ice sheets and is temporally related to the ice rafting event known as Heinrich Event 1 (HE1). Equivalent ice sheet wasting events accompanied every deglaciation of the past 1Myr or so24,25, and various studies have attempted to demonstrate that these were also related to a weakening of AMOC25,26,77. A large freshwater event could perturb the AMOC into a weak mode even from its glacial monostable strong state but this would only be transient and the AMOC would recover following the end of freshwater input. However, the decrease in ice sheet height that might accompany such a freshwater event could also push the system into a bistable regime or even toward monostable weak (Fig. 5) and in this case the AMOC would remain in a weak mode even after cessation of freshwater input. From this point there are a few possible ways in which a recovery of the AMOC (which must happen at some point during glacial termination) could be achieved:

  1. (1)

    An increase in ice sheet height pushes the system back toward monostable strong (e.g. Fig. 5 point 2d→ 1a). This is unlikely to occur during a deglaciation, when ice sheets are generally thought to be in recession.

  2. (2)

    Recovery within the bistable window (e.g. 2c→1c). This would require an external forcing such as a negative freshwater (salinification) event, which is also unexpected given the tendency for meltwater generation during the deglacial period.

  3. (3)

    Increasing atmospheric CO2 pushes the bistable window to the right and forces the system toward monostable strong (e.g. 2d→1b for CO2+). This would require that the rate of CO2 increase effectively outpaced the rate of ice sheet decline, which is possible but probably unlikely again given that ice sheets tend to diminish throughout deglaciation. We suggest that this mechanism could explain AMOC recovery during the Bølling-Allerød.

  4. (4)

    CO2 eventually reaches near-interglacial levels and the system is forced into a monostable strong regime. This scenario is inevitable once CO2 reaches the desired threshold. The result would be an AMOC recovery towards the end of deglaciation.

We suggest that the last option (4) is the most likely scenario given the reasoning above for options (1) and (2) and the fact that most terminations are not interrupted by a strong oscillation25 as we might expect from option (3) (see discussion below).

Solutions in the time domain

In Fig. 6 we transfer the stability characteristics from our previous discussion into the time domain for three possible deglacial scenarios (A–C, below). To reflect the decay of ice sheets during deglaciation we make the assumption that a finite flux of freshwater will enter the surface North Atlantic throughout the deglacial period. If the AMOC enters a bistable regime during that time, this flux is sufficient to force it into the weak mode and recovery will not occur unless the system experiences a structural change in stability. It should be stressed that while the bistable regime and deglaciation may overlap in the time domain, they are not synonymous: the bistable regime can exist whenever ice volume and CO2 are at intermediate values (as apparent during MIS 3). Furthermore, if enough freshwater is added, this could induce a weak mode even when the AMOC is in a monostable strong regime (while ice sheets are still very large) and this could modulate the precise timing or occurrence of abrupt transitions during deglaciation beyond those predicted by stability changes alone.

Fig. 6: Deglacial changes in structural stability.
figure 6

Idealised timeseries of CO2, ice volume and North Atlantic freshwater (FW) flux for last two terminations with corresponding variations in AMOC strength according to text and Fig. 4. Solid curves represent T1; Dashed curves represent T2; Dotted curves represent a scenario for T1 in which variations in the AMOC are forced purely by freshwater (not requiring an underlying structural change in stability). Times t0–t6 are referred to in text.

Scenario A: The HS1-B/A-YD triptych (Termination 1; solid curves in Fig. 6). Starting in a glacial mode (t0 in Fig. 6; yellow diamond in Fig. 5), warming northern hemisphere summers induce melting of northern continental ice sheets. The resulting freshwater tends to weaken the AMOC (t1) and decreasing ice sheet height pushes it into a bistable regime so that a weak mode becomes inevitable (HS1): t2 in Fig. 6; 1a→2d in Fig. 5. During this period of weakened circulation increasing atmospheric CO2 pushes the window of bistability to the right fast enough to outpace the decline in ice sheet height, triggering a recovery of AMOC as the system is forced back toward monostable strong (B/A): t3; 2d→1b (CO2+). Note that recovery of the AMOC in this scenario can occur even in the presence of continued freshwater input45. Within Termination 1, the deglacial rise in atmospheric CO2 stagnates during the B/A such that the return to a weak mode of AMOC becomes inevitable as ice sheets continue to decline, pushing the system back towards bistability and possibly to monostable weak (YD): t4; 1b→c2(→2a). Atmospheric CO2 then starts to rise again and eventually reaches a level sufficient to induce an interglacial monostable strong mode of AMOC (yellow square in Fig. 5) and a second recovery is achieved (t5; end of YD).

Scenario B: ‘Standard’ termination (e.g. Termination 2; dashed curves in Fig. 6). Starting in a glacial mode (t0 in Fig. 6; yellow diamond in Fig. 5), warming northern hemisphere summers induce melting of northern continental ice sheets. The resulting freshwater tends to weaken the AMOC (t1) and decreasing ice sheet height pushes it into a bistable regime so that a weak mode becomes inevitable (terminal HS event): t2; 1a→2d. As deglaciation proceeds, increasing CO2 and wasting ice sheets keep pace with one another so that the AMOC remains in a bistable or monostable weak mode until CO2 eventually reaches a level sufficient to induce an interglacial monostable strong mode of AMOC (yellow square) and recovery is achieved (t5). Brief cold events such as the 8.2ka event and the cold event following HS11 recorded at ODP Site 983 (Fig. 4) may be expected as northern ice sheets continue to decay. However, these are distinct from the YD because they occur within an interglacial state (when CO2 is high and the AMOC is monostable strong) and are therefore transient by nature i.e. solely driven by e.g. freshwater fluxes.

Scenario C: No change in structural stability of AMOC (dotted curves in Fig. 6). Variations in AMOC are forced purely by changes in freshwater input to the surface North Atlantic.

A test of these predictions requires, at the very least, robust records of atmospheric CO2 and ice volume (or preferably ice sheet height and configuration). However, while CO2 reconstructions are available back to 800ka, reliable sea-level reconstructions are limited to the last deglaciation (note disagreement among SL reconstructions across T2 in Fig. 4). As we discuss in the final section, this is a critical area for future study. Furthermore, Scenarios A and B could be modified by variable addition of meltwater during deglaciation. For example, AMOC recovery could be delayed beyond entering a monostable strong regime until freshwater addition waned sufficiently. The interplay between freshwater and structural changes in stability might also explain some of the more nuanced features observed during previous terminations (e.g. the 2-part division of HS178 or the partial recovery of AMOC during T325). Of relevance here is a newly published study5 in which the authors posit a direct link between insolation forcing and the duration of glacial terminations. Quantifying the duration of a glacial termination using benthic δ18O is not ideal (as alluded to in the previous section) but nevertheless the results are instructive since they suggest that higher levels of insolation might lead to a shorter termination by driving faster ice sheet decay. According to our arguments, in these cases we should probably not expect termination to be interrupted by an early recovery of the AMOC (SL rises too fast with respect to CO2). Furthermore, by their5 calculations Termination 1 is not only relatively long in duration (ranking joint 4th longest of the 11 analysed) but is also longer than might be predicted from its insolation (in this case using integrated summer energy as the relevant quantity). Perhaps its longer duration (i.e. slower SL rise with respect to CO2) made Termination 1 particularly prone to an early recovery of the AMOC.

Early interglacial legacy of deglacial climate instability

So far we have described the interplay between millennial and orbital-timescale variations with respect to glacial terminations but as discussed recently25 these interactions are likely to have lasting effects during the period immediately following termination (i.e. during early interglacial times). This becomes important when making comparisons between interglacials whose preceding terminations may have been very different. Such comparisons are important in light of ongoing changes in e.g. greenhouse gas concentrations, which need to be viewed within the context of natural variability to allow robust differentiation of potential human influences on the climate system79. For example the rise in atmospheric CH4 and CO2 over last ~5 and ~8 kyr (respectively) of the Holocene80,81 has sparked debate over whether or not humans have been influencing climate for many thousands of years82,83,84,85.

At the heart of these discussions lies the contrasting histories of both CH4 and CO2 over the course of the Holocene (which experienced upwards trends in both gases) as compared with previous interglacials (which experienced downward trends in some cases e.g. MIS 7, 9 and 19)86. In a recent study25 we made the case that the earliest portion (first few kyr) of several previous interglacials should not be included when making these comparisons if quasi-equilibrium is a prerequisite for such comparisons to be valid. We defined (quasi-)equilibrium as a (hypothetical) situation in which Earth’s climate is equilibrated with respect to key boundary conditions including orbital configuration, ice volume and atmospheric CO2 concentration and we argued that the earliest parts of many previous interglacials were not at equilibrium thanks to continued adjustment of ocean circulation patterns following deglacial oscillations of the AMOC (note the peculiar coincidence of high amplitude millennial activity and high CO2 associated with deglacial periods as mentioned earlier; Fig. 2e ii). We concluded by questioning whether the downward trends in CO2 observed for MIS 7, 9 and 19 would remain if the initial period of higher and decreasing CO2 associated with AMOC recovery (equivalent to the interval between t5 and t6 in Fig. 6) were omitted from analysis. Below we make the case that the Holocene upward trend in CO2 may be considered anomalous a priori thanks to the early recovery of AMOC midway through Termination 1, in which case the early intervention of humans need not be invoked to explain the apparent difference between the Holocene and previous interglacials.

Deaney et al.23 compared the evolution of atmospheric CO2 across the last two terminations and concluded that the larger apparent magnitude of CO2 change across T2 was due to the late recovery of AMOC associated with that termination (e.g. Fig. 4). The continuous rise in CO2 associated with Heinrich Stadial 11 (HS11) and subsequent ~20 ppmv jump in CO2 associated with recovery of the AMOC ~129 ka resulted in a larger net change in CO2 across T2 than the corresponding change across T1. Thus, atmospheric CO2 was effectively lower at the onset of the Holocene than it might otherwise have been if the B/A had not occurred (an analogous claim was made recently concerning mean ocean heat content, which was also higher at the onset of MIS 5e than it was at the start of the Holocene87). This logic was borne out in a modelling study by Ganopolski and Brovkin38 in which they found that a pronounced overshoot in CO2 occurs at the beginning of an interglacial when recovery of AMOC happens only at the end of termination. This is followed by a rather constant or even decreasing trend in CO2 during the subsequent interglacial. In contrast (and analogous to T1) if a termination is interrupted by an early recovery of the AMOC, CO2 will tend to be lower at the start of the subsequent interglacial and proceed to increase throughout the next several thousand years.

We note that CO2 also increased throughout MIS 11, which has been used as an analogy for MIS 1 because of the similarity in insolation forcing88. Indeed, insolation forcing across T1 was much weaker than e.g. T2 because of lower eccentricity and this has been suggested as a possible reason why an early AMOC recovery did not occur during T272, However, more recently MIS 19 has been advocated as a better insolation analogue for MIS 186 and it is notable that T9 (into MIS 19) did not experience an early AMOC recovery25 and CO2 decreased throughout MIS 19, following an initial overshoot25.

Earlier we made a case that Termination 1 (T1) was unusual in the sense that it was interrupted by an early recovery of the AMOC (associated with the Bølling/Allerød). Given the tendency for changes in ice volume and CO2 to coincide during deglaciation, together with the ubiquitous input of freshwater, we argue that the most likely case for a typical termination is Scenario B: essentially an uninterrupted period of weakened AMOC and consequently a late recovery of AMOC (as compared with T1). Following this logic and the studies mentioned above23,38 we conclude that the relatively low concentration of CO2 observed at the onset of MIS 1 (compared for example with MIS 5, 7 and 9), and its subsequent increase throughout the Holocene, may also be considered as atypical with respect to other interglacials, regardless of any influence that humans may or may not have had. This should be taken into account when making comparisons between the Holocene and more typical interglacials such as MIS 5.

Summary and outlook

In this perspective we have highlighted progress in our understanding of glacial terminations as the ultimate expression of the interaction between millennial and orbital-timescale variations in Earth’s climate. Specifically we discussed the importance of feedbacks, which appear most accentuated within the sweet spot of millennial activity, and in particular during glacial termination when they help to explain the magnitude and rapidity of these events. We then made the case that structural changes in the stability of the AMOC can explain the occurrence of terminal AMOC oscillations such as the HS1-B/A-YD triptych. Finally we argued that the unusual AMOC oscillation of Termination 1 led to the low initial concentration and subsequent rise of atmospheric CO2 throughout the Holocene, which itself may therefore be considered unusual a priori (without the need to call on human intervention).

However, there are still major gaps in our knowledge, which need to be filled before we can test some of the hypotheses put forward here. For example, sea-level reconstructions beyond the last glacial maximum are uncertain and inconsistent (e.g. throughout MIS 360,61 and across MIS 5a/489,90 and Termination 291,92,93,94; Fig. 4), making it difficult to test our ideas about changing AMOC stability as a driver for deglacial oscillations in ocean circulation or even to assess the boundary controls on the occurrence of millennial-timescale variability on G-IG timescales. More robust sea-level reconstructions are therefore required that will allow us to compare the relative rates of CO2 rise and ice sheet decline during previous terminations and assess our assertions for example about the uniqueness of the B/A-YD oscillation.

An area of gathering momentum is model intercomparison. Many of the exciting theoretical discoveries made over the past few years have been obtained using a single model, leaving their results vulnerable to possible inadequacies in a particular model. Projects such as PMIP94,95 will allow more confidence to be gained (or not) in proposed mechanisms as will efforts to compare unforced (e.g. without freshwater forcing) AMOC oscillations in a variety of different model setups. To date most simulations of abrupt climate variability have been forced with freshwater but a growing number of models display forms of ‘unforced’ oscillations (or at least strongly non-linear responses to gradually changing boundary conditions). Given the apparent importance of factors other than freshwater forcing we advocate a community effort to simulate abrupt transitions through other means. Ultimately, the incorporation of an interactive carbon cycle and ice sheet dynamics (including e.g. solid earth and ice shelf cavity processes) will allow a fuller representation of the mechanisms involved in the coupled earth system. In this way we will be able to vastly improve our understanding of how internal climate system dynamics respond to external forcing across deglacial transitions and finally reach a turning point in our understanding of glacial terminations. This knowledge is not merely of academic interest; uncertainties over the location of stability thresholds with respect to our current climate mean that we cannot rule out the possibility that the AMOC is already within a regime of multiple equilibria with respect to North Atlantic freshwater forcing96.