Millennial atmospheric CO2 changes linked to ocean ventilation modes over past 150,000 years

Ice core measurements show diverse atmospheric CO2 variations—increasing, decreasing or remaining stable—during millennial-scale North Atlantic cold periods called stadials. The reasons for these contrasting trends remain elusive. Ventilation of carbon-rich deep oceans can profoundly affect atmospheric CO2, but its millennial-scale history is poorly constrained. Here we present a well-dated high-resolution deep Atlantic acidity record over the past 150,000 years, which reveals five hitherto undetected modes of stadial ocean ventilation with different consequences for deep-sea carbon storage and associated atmospheric CO2 changes. Our data provide observational evidence to show that strong and often volumetrically extensive Southern Ocean ventilation released substantial amounts of deep-sea carbon during stadials when atmospheric CO2 rose prominently. By contrast, other stadials were characterized by weak ventilation via both Southern Ocean and North Atlantic, which promoted respired carbon accumulation and thus curtailed or reversed deep-sea carbon losses, resulting in diminished rises or even declines in atmospheric CO2. Our findings demonstrate that millennial-scale changes in deep-sea carbon storage and atmospheric CO2 are modulated by multiple ocean ventilation modes through the interplay of the two polar regions, rather than by the Southern Ocean alone, which is critical for comprehensive understanding of past and future carbon cycle adjustments to climate change. The variable intensity of Southern Ocean as well as North Atlantic deep-water ventilation explains differences in atmospheric CO2 trends and magnitudes during cold stadials over the past 150,000 years, according to a record of deep-ocean acidity.


Article
https://doi.org/10.1038/s41561-023-01297-x In this study, we use the term ventilation to describe effects on ocean interior dissolved inorganic carbon (DIC) driven by changes in air-sea CO 2 exchange in the deep-water formation regions and respired carbon accumulation during water-mass transportation (Methods).We employ carbonate ion concentrations ([CO 3   2−   ]; a parameter tightly linked to seawater pH and acidity) to constrain DIC variations with which to infer deep-sea ventilation states and thereby atmospheric CO 2 changes 19,21,22 .Everything else being equal, enhanced Southern Ocean ventilation would reduce deep-sea DIC, raising deep-water [CO 3 2− ] and atmospheric CO 2 (Supplementary Fig. 1).Conversely, increased respired carbon accumulation in the ocean interior, which could be driven by suppressed ventilation from either polar region, would decrease deep-water [CO 3   2− ] and atmospheric CO 2 (refs.9,15,23).Unfortunately, investigation of millennial-scale variations in deep-sea carbon storage-including their sign, magnitude and timing in relation to atmospheric CO 2 -is hampered by a paucity of high-quality [CO 3 ] records with robust chronologies are warranted to better understand abrupt ocean ventilation mechanisms.

Benchmark [CO 3 2− ] record
In this Article, we present the first high-resolution deep-water [CO 3   2− ] record to span the entire last glacial cycle (Figs. 1 and 2).The record was obtained using core MD95-2039 (40.6°N, 10.3° W, 3,381 m) from the Iberian Margin in the North Atlantic, utilizing multifold advantages of sediments from this setting for palaeo-reconstructions (Extended Data Fig. 2 and Methods).Our [CO 3   2− ] record is reconstructed using B/Ca in benthic in southern westerly winds and, thus, Southern Ocean ventilation 12 (see ref. 13 for further discussions).When Atlantic meridional overturning circulation (AMOC) weakens during stadials, Antarctica and surrounding regions warm up, accompanied by strengthened and/or poleward shift of southern westerlies, which would enhance deep-sea ventilation via the Southern Ocean and thereby increase atmospheric CO 2 as seen in ice-core records 6,8,10,11 .Yet it remains unclear why atmospheric CO 2 did not rise during many stadials despite clear Antarctic warmings 2,3 .The latest high-resolution atmospheric CO 2 record, whose timing and magnitude can be firmly constrained by paired CH 4 measurements from the same ice core, shows no net CO 2 increase during a considerable number of stadials when the bipolar seesaw persistently operated 3,14 (Extended Data Fig. 1).To fully understand the ocean-atmosphere carbon interactions, marine proxy records are required because they provide more direct clues about oceanic processes.The ocean interior is ventilated mainly by deep waters, including southern-sourced waters (SSWs) formed in the polar Antarctic Zone and northern-sourced waters (NSWs) formed in the North Atlantic.Published marine reconstructions [15][16][17][18][19][20] yield valuable information about the signs of deep-sea ventilation changes related to SSWs, although how the relative ventilation strength varied between stadials, which is critical to understand magnitudes and rates of atmospheric CO 2 changes, is yet to be assessed.Even less is known about the impacts on atmospheric CO 2 from ventilation changes linked to NSW production rate and volumetric extent.In this Article, we investigate the role of deep Atlantic ventilation via the two polar regions in modulating deep-sea carbon storage and atmospheric CO 2 variations on millennial timescales.] at MD95-2039 (this study).d, Compiled North Atlantic sediment Pa/Th, an AMOC strength proxy 29,32,33 .e, EDML ice-core δ 18  ] at MD95-2039 (this study).d, Compiled North Atlantic sediment Pa/Th, an AMOC strength proxy 29,32,33 .e, EDML ice-core δ 18  ] at MD07-3076 17 .The curved arrow indicates that the [CO 3
Figure 1 shows our [CO 3 2− ] record alongside ice-core δ 18 O and atmospheric CO 2 , all placed on the same chronology 4,25,27 , which is a prerequisite for confident comparisons of marine and ice-core records.This study focuses on stadial changes.Our record provides sufficient data to resolve striking [CO 3 2− ] fluctuations accompanied by diverse atmospheric CO 2 changes during 25 stadials over the past 150,000 years.We use a Monte Carlo-style approach to identify major features associated with our [CO 3 2− ] and contemporaneous ice-core CO 2 during the progress of stadials (Methods and Supplementary Figs.2-7).On the basis of their [CO 3 2− ]-CO 2 similarities and contrasts, the 25 stadials can be classified into 5 groups, leading us to infer 5 ocean ventilation modes and associated hydrographic changes in the two polar regions (Figs.1-4).Note that our mode classification differs from those based on either ice-core 28 or marine 29 archives.

Five modes of stadial ocean ventilation
We first assess Heinrich stadials (HSs) associated with ~40-55 ppm atmospheric CO 2 rises at glacial terminations.During HS1 and HS11, our record shows early [CO 3 2− ] increases followed by decreases (Mode I; Fig. 3).On millennial timescales when the global alkalinity change is relatively minor 30 , a [CO 3 2− ] increase (decrease) probably reflects a DIC decrease (increase) 22 .The early [CO 3 2− ] increases and associated DIC decreases cannot be explained by northward expansion of carbon-rich SSW 31 .During these stadials, Pa/Th data 29,32,33 indicate a slowdown of AMOC (Fig. 1d), which would promote respired carbon accumulation 34 with an effect to decrease [CO 3   2− ] in the deep Atlantic.Nevertheless, reduced AMOC cooled the North Atlantic and, through the bipolar seesaw, warmed the Southern Ocean (Fig. 1a,e).This was probably accompanied by strengthening and/or poleward shifts of southern westerlies 10,11,35 , which would enhance deep-sea ventilation via the Southern Ocean as suggested by proxy and model results 8,20,36,37 .We thus ascribe our reconstructed DIC decreases at MD95-2039 during early stages of HS1 and HS11 to carbon losses driven by enhanced Southern Ocean ventilation (Fig. 4a), although we do not exclude other mechanisms 9,12,38 such as a weakened biological pump that might also assist deep-sea carbon release.Our inferred Southern Ocean ventilation enhancement is supported by opal flux, radiocarbon and oxygenation data 15,16,18,21,37,39 .Moreover, the location of MD95-2039 suggests that our inferred Southern Ocean ventilation enhancement was sufficiently intense and far-reaching to reduce deep Atlantic carbon all the way up to 40° N.This may have resulted from particularly pronounced bipolar seesaw responses to substantial AMOC weakening at glacial terminations 26,40,41 (Fig. 1d).Indeed, Antarctic ice-core records 27 indicate that early HS1 and HS11 were accompanied by strongest warmings over the past 150 ka (Fig. 1e).By efficiently releasing deep-sea carbon and hence substantially increasing atmospheric CO 2 , marked Southern Ocean ventilation enhancement may have played a critical role in initiating and/or sustaining glacial terminations (Fig. 4a).

2−
] declined at MD95-2039, with lowest values occurring close to the peak of ice-rafted debris (IRD) deposition at the Iberian Margin 42 (Fig. 1c and Extended Data Fig. 3).The associated DIC increase was due, at least partly, to increased accumulation of respired carbon due to further AMOC weakening as registered by Pa/Th in a nearby core 43 (Extended Data Fig. 4).If so, our data imply a DIC increase in the Atlantic below ~3 km, which would slow atmospheric CO 2 rise, consistent with ice-core data 1 .The late HS11 [CO 3 2− ] decline can also be attributed to further AMOC weakening, supported by a coeval benthic δ 13 C decrease at the Iberian Margin 44 (Extended Data Fig. 4).Compared with late HS1, late HS11 was characterized by greater Antarctic warming and perhaps stronger deep-sea ventilation via the Southern Ocean 10,11 (Fig. 1e).Although further investigation is needed, we speculate that a substantial deep-sea volume to the south of 40° N ] (c) and atmospheric CO 2 (d) evolutions define five ventilation modes (Modes I to V).The delayed atmospheric CO 2 decline associated with Mode V may be due to age uncertainties (Supplementary Fig. 7).Data for HS11 shown in a (dashed curve) are from ODP 976 G. bulloides δ 18 O 26 and adjusted by a factor of −3 to facilitate plotting.Note different y-axis scales, especially for atmospheric ΔCO 2 , between different modes.To assist comparison, records are normalized by the same approach to be plotted against the progress of stadials (x axis), with vertical dashed lines indicating start (0%) and end (100%) of the stadial.For clarity, only stacked changes are shown for Modes II (HS4-6, 7b) and IV (HS2, 3, 7a and 10; GS4, 6-8, 10-12, 16, 17 and 25).See Methods and Supplementary Figs.  − ] at MD95-2039 underwent a brief decline followed by a net increase (Fig. 3).The initial [CO 3

2−
] decline and associated DIC rise may reflect increased biogenic matter respiration and/ or northward expansion of carbon-rich SSW in the deep Atlantic, linked to AMOC slowdown 29,32,33 (Fig. 2d).The subsequent [CO 3 2− ] increase and thus DIC decrease occurred significantly before stadial terminations.This timing lead cannot be ascribed to age uncertainties or bioturbation because a similar timing relationship is observed between C. wuellerstorfi B/Ca (used for deep-water [CO 3 2− ]) and G. bulloides δ 18 O (used for chronology) in the same core (Extended Data Fig. 5 and Methods).Existing εNd data 29,45 show little sign of increased mixing of low-DIC NSW in the deep Atlantic within Mode II stadials (Fig. 2h).Given large Antarctic warming and sea-ice retreat in surrounding oceans 10,11,27,36 (Fig. 2e), the DIC decrease can be readily explained by carbon loss owing to enhanced Southern Ocean ventilation.7).The location of MD95-2039 suggests carbon losses from a sizeable mass of deep waters from the Southern Ocean to ~40° N, helping to explain the prominent atmospheric CO 2 rises during Mode II stadials.Compared with the Mode I state, however, Southern Ocean ventilation enhancement during Mode II stadials appears to have been insufficiently intense to prevent the initial [CO 3 2− ] declines at MD95-2039 (Fig. 4b).] and oxygen levels are observed, respectively, at sites TNO57-21 (41° S, 8° E, 4,981 m) and TNO57-13PC (53° S, 5° E, 2,848 m) from the South Atlantic 15,47 .Increases in deep-water [CO 3 2− ] and oxygenation are also documented at TNO57-21 during HS8 22 (Extended Data Fig. 8).These observations suggest that Southern Ocean ventilation during Mode III stadials was sufficiently vigorous to release carbon from the abyssal South Atlantic, in agreement with marked Antarctic warming and enhanced upwelling in the Antarctic Zone 8,12,27 (Fig. 1e).The associated carbon releases would have contributed to contemporaneous atmospheric CO 2 rises.Yet sustained low [CO 3

2−
] and hence high DIC in deep waters at MD95-2039 indicates that Southern Ocean ventilation did not strengthen enough to overcome the effects from respired carbon accumulation and/or increased mixing of carbon-rich SSW driven by AMOC reductions during the YD and HS8.Thus, relative to Modes I and II stadials, deep-sea carbon loss driven by enhancement of Southern Ocean ventilation remained geographically more restricted to the south of ~40° N during Mode III stadials (Fig. 4c).
For Mode IV, atmospheric CO 2 generally displays muted net changes during Greenland stadials (GS) ] or oxygenation increase is seen at South Atlantic sites TNO57-21 and MD07-3076 (44° S, 14° W, 3,770 m) 17,22 (Fig. 2g and Extended Data Figs.6 and 7).These observations suggest that carbon loss via the Southern Ocean must have been restricted to relatively shallow Atlantic waters located to the south of 44° S (Fig. 4d).Nevertheless, the concurrent Antarctic warming (Fig. 2e) 27 implies that Southern Ocean ventilation was possibly enhanced by some degree compared with pre-stadial conditions.Even considering potential complexities to link Antarctic temperature and Southern Ocean ventilation changes 13 , Antarctica warming would reduce Antarctic Zone sea-ice extent and SSW solubility pump efficiency.Thus, Southern Ocean changes probably promoted CO 2 outgassing with an effect to increase deep Atlantic [CO 3 ] reconstructions are desired, published deep Atlantic εNd records 29,45 show no consistent variations during transitions into Mode IV stadials (Fig. 2h), arguing against increased mixing of low-[CO 3 2− ] SSW as an exclusive explanation for the observed [CO 3

2−
] declines at MD95-2039.We suggest that much of these [CO 3 2− ] declines and associated DIC rises must have been caused by increased accumulation of respired carbon driven by weakened AMOC 29,32,33 (Fig. 2d), which is supported by model results 31,48 .Alongside DIC rises during many Mode IV stadials at South Atlantic site MD07-3076 (Fig. 2g), there appeared to be widespread sequestration of respiratory carbon, and hence nutrients, in the ocean interior driven by suppressed ventilation via the North Atlantic.By increasing ocean's biological pump efficiency 9,23 , this would counteract the effect of carbon loss from the shallow Southern Ocean and thus curb atmospheric CO 2 rises during Mode IV stadials (Fig. 4d).Indeed, published model simulations 48 show increased atmospheric CO 2 absorption by the North Atlantic owing to improved biological and solubility pump efficiencies during interstadial-to-stadial transitions under intermediate climate conditions (Extended Data Figs. 9 and 10).
Mode V is illustrated by a deep-water [CO 3 2− ] decrease at MD95-2039 during GS26 when atmospheric CO 2 declined at the onset of the last glaciation (Figs. 1 and 3).A similar change may have occurred at HS9, but its age in core MD95-2039 is less certain.The observed [CO 3   2− ] decrease and thus DIC increase can be caused by increased respired carbon accumulation and/or northward SSW penetration linked to AMOC weakening (Fig. 1d).Different from Modes I-IV stadials, which were mechanistically related to bipolar seesaw processes 11 , both Antarctica and Greenland cooled during GS26 25,27 (Fig. 1a,e).In a wider perspective, GS26 was characterized by a global cooling 49 , sea-ice expansion surrounding Antarctica and plausibly an equatorward shift of southern westerlies, which would suppress Southern Ocean ventilation 10,12 and thus deep-sea carbon leakage.Moreover, concomitant North Atlantic solubility pump enhancement due to cooling 25,27 would have promoted atmospheric CO 2 sequestration in NSW.We therefore suggest that a bipolar synergy operated to lower atmospheric CO 2 during GS26, providing a positive feedback to global cooling at the last glacial inception (Fig. 4e).

Bipolar control on millennial atmospheric CO 2
Our [CO 3 2− ] reconstructions, alongside published data 8,15,16 , suggest Southern Ocean ventilation as a key factor to modulate millennial deep Atlantic carbon storage and atmospheric CO 2 variations over the last glacial cycle (Figs.1-4).Although based on data from the Iberian Margin for the sake of high-resolution reconstructions with robust chronologies, our inferred changes in Southern Ocean ventilation perhaps also affected carbon storage in the more voluminous deep Indo-Pacific oceans.Instead of a simple 'on-off switch', our data reveal that Southern Ocean ventilation modes varied in both strength and geographic extent among stadials (Fig. 4).Repeated millennial-scale enhancement of Southern Ocean ventilation appears to have been strong and often far-reaching, and thus helps explain relatively large atmospheric CO 2 rises during many stadials, particularly those at glacial terminations (Modes I-III; Figs. 1 and 4).Conversely, suppressed Southern Ocean ventilation facilitated deep-sea carbon sequestration, contributing to lowering atmospheric CO 2 during GS26 and possibly HS9 (Mode V; Figs. 1 and 4).
Importantly, our data underscore an indispensable but previously underappreciated role of North Atlantic ventilation in governing deep-sea carbon sequestration and atmospheric CO 2 .As demonstrated by Mode IV data and supported by models 31,34,50 , reduced ventilation via the North Atlantic, instead of the Southern Ocean, probably promoted respired carbon accumulation in the ocean interior (Figs. 2 and 4d).By improving the oceanic biological pump efficiency, the associated carbon accumulation would lower atmospheric CO 2 (Supplementary Fig. 1), a mechanism widely invoked to explain atmospheric CO 2 declines on glacial-interglacial timescales 9,15,23,34,51 but largely overlooked on millennial timescales.In fact, all stadials appear to be accompanied by reduced deep-water [CO 3 2− ] at some point, indicative of pervasive influences of changes in North Atlantic ventilation because plausible Southern Ocean ventilation enhancements (except for Mode V) tend to raise deep Atlantic [CO 3 2− ] (Fig. 3).Our inferred effect on atmospheric CO 2 is supported by model simulations 48 , which show increased North Atlantic CO 2 absorption during all studied stadials (Extended Data Figs. 9 and 10).Under such conditions, if Southern Ocean CO 2 outgassing is not strong enough, atmospheric CO 2 could decrease.In other words, atmospheric CO 2 is jointly affected by the relative influence of Southern Ocean outgassing and North Atlantic absorption.We thus propose that the interplay of ventilation via the two polar regions, namely a bipolar control, should be considered when evaluating diverse millennial atmospheric CO 2 variations.
Our proposed bipolar control helps resolve a long-standing puzzle 2 regarding mechanisms for disparate atmospheric CO 2 responses among MIS 3 stadials.The persistent anti-phased temperature relationship between polar ice cores 27 suggests sustained operation of the bipolar seesaw during MIS 3 (Fig. 2).One would expect [10][11][12] enhanced Southern Ocean CO 2 outgassing and thus atmospheric CO 2 rises during all MIS 3 stadials, which is at odds with observations 2,3 (Fig. 2).The lack of atmospheric CO 2 rises during many stadials must indicate concomitant processes to offset CO 2 outgassing from the Southern Ocean.Atmospheric CO 2 is affected by carbon exchange with both Article https://doi.org/10.1038/s41561-023-01297-xocean and land biosphere reservoirs, but how the terrestrial biosphere carbon changed during stadials remains poorly constrained 2,41 .Here we provide critical data evidence (Fig. 2) to demonstrate that respired carbon accumulation, driven by reduced North Atlantic ventilation 31,48,50 , at times cancelled or surpassed the effect of Southern Ocean outgassing, explaining the observations of little net change or even declines in atmospheric CO 2 during some MIS 3 stadials.
In summary, our reconstructions reveal hitherto unrecognized bipolar ocean ventilation modes behind various types of rapid atmospheric CO 2 changes during the past 150,000 years.Since the Industrial Revolution, the ocean has absorbed a substantial amount of atmospheric CO 2 via high-latitude oceans in both hemispheres, helping to slow down global warming 52 .Given ongoing and possible future AMOC weakening 53 , it is imperative to fully understand how processes in polar oceans, including both the Southern Ocean and the North Atlantic, affect atmospheric CO 2 and associated ocean acidification.Models used to quantify ocean carbon storage must simulate these processes and their interplay accurately.Our [CO 3 2− ] record-covering various types of stadials with robust chronologies throughout the last glacial cycle-can be employed by models to assess the relative importance of ventilation changes in the two polar regions under different climate conditions, with critical implications for future carbon cycling and climate changes.

The use of the term 'ventilation'
To avoid confusion, we clarify the meaning of the term 'ventilation' in the context of the carbon cycle as discussed in the main text.In this study, we use 'ventilation' to describe DIC changes in the ocean interior related mainly to two processes: (1) air-sea CO 2 exchange in the deep-water formation regions and (2) DIC changes related to respired carbon accumulation during the propagation of the water-mass signals (Supplementary Fig. 1).We use 'ventilation' as a shorthand to describe their combined effect on DIC in the ocean interior.Although we do not specifically distinguish the respective roles of the two processes (as under certain circumstances it is difficult to do so), it does not affect our conclusions in the main text.
To exemplify, everything else being equal, enhanced Southern Ocean ventilation would decrease DIC in the deep ocean by strengthened CO 2 release associated with improved air-sea gas exchange in the deep-water formation regions surrounding Antarctica and/or by reduced respired carbon accumulation due to faster transport of SSW within the ocean (Supplementary Fig. 1).Reduced North Atlantic ventilation would increase deep Atlantic DIC via promoting accumulation of respired carbon due to a sluggish AMOC.
It is important to note that in our study an increased mixing proportion of SSW does not mean better ventilation via the Southern Ocean.For example, the deep Atlantic might be occupied by more SSW during the Last Glacial Maximum than during the Holocene, but ventilation via the Southern Ocean was probably reduced during the Last Glacial Maximum as suggested by proxy data 15,16,18,51,54 .

Core selection
We appreciate that the regions such as Southern Ocean and the Pacific Ocean can provide valuable information about deep-ocean ventilation histories due to their proximity to the deep-water formation sites surrounding Antarctica (for example, Antarctic Zone) or large volumes of high-DIC waters (for example, deep Pacific Ocean).However, millennial-scale reconstructions based on sediments from these regions may be complicated by chronological uncertainties, lack of sufficient benthic foraminiferal shells or low sedimentation rates.We chose to work on core MD95-2039 from the Iberian Margin for the following advantageous reasons.First, we can follow the established approach 55 to construct a robust age model for the core, a prerequisite for high-confidence comparisons of marine and ice-core records.The Iberian Margin is arguably the best location in the world ocean to construct reliable chronologies for long time series of seawater property reconstructions.Second, the core location is characterized by high sediment-deposition rates, conducive for detailed reconstructions with minimal influence of bioturbation (Extended Data Fig. 3).Third, MD95-2039 presents as a sensitive location to investigate effects of water-mass changes on millennial atmospheric CO 2 variations [1][2][3][4] during well-documented and extensively studied MIS 3 when the NSW-SSW boundary is thought to be located proximal to the core's water depth 56,57 .Fourth, the latitude of MD95-2039 makes it useful to constrain the geographic extent of the relative influence during stadials of ventilation from the Southern Ocean, which tends to increase deep-water [CO 3 ] due to enhanced biogenic matter respiration (Supplementary Fig. 1).

Samples and analyses
To ensure enough shells for analyses, about 30 cm 3 of sediment for each sample (∼2 cm thickness) from core MD95-2039 was disaggregated in de-ionized water and was subsequently wet sieved through 63 μm sieves.About 10 shells of G. bulloides from each sample were picked from the 250-300 μm size fraction for δ 18 O measurements.The average analytical error is ~0.08‰ in δ 18 O.We compared data from overlapping depths and observed negligible analytical offset between new and published 58 data, so new (n = 846) and published (n = 411) data (Fig. 1a) are combined for consideration.
We have measured six radiocarbon dates for the early Holocene and Bølling/Allerød periods (three dates for each period).About 400 shells of G. bulloides from the 300-355 μm size fraction (equivalent to ~4-6 mg CaCO 3 ) were used for each measurement at Research School of Earth Sciences at the Australian National University (Supplementary Table 1).
For B/Ca analyses, the epifaunal benthic foraminiferal species C. wuellerstorfi was picked from the 250-500 μm size fraction.For each sample, ~10-20 shells were picked and then double checked under a microscope before crushing to ensure that consistent morphologies were used throughout the core.Foraminiferal shells were then crushed and cleaned following the established protocol 60 .Carbonate B/Ca ratios were measured on an inductively coupled plasma mass spectrometer using procedures outlined in ref. 61, with an analytical error better than ~4% (2σ).The Mn/Ca and Al/Ca were also measured, and they showed no correlation with B/Ca, suggesting minimal influences from silicate or diagenetic coatings.

Age model
The age model of MD95-2039 is based on a combination of methods.For samples from <~14 ka, the age model relies on six radiocarbon ages and an assumed core-top age of 0 ka.The radiocarbon dates were converted into calendar ages using the CALIB 8.1 with the Marine20 calibration curve 62,63 (Supplementary Table 1).No radiocarbon date is used during the YD to avoid any complications with past surface reservoir age variations.For samples from ~14.5 to 22 ka and from ~65 to 120 ka, the age model was constructed mainly using the established approach of aligning G. bulloides δ 18 O in MD95-2039 to the NGRIP ice-core δ 18 O record 25 on the AICC2012 age model 64 .Due to the lack of significant structure in NGRIP δ 18 O during HS1, we added two tie points to align MD95-2039 G. bulloides δ 18 O to the Hulu speleothem δ 18 O record (on the U-Th age) 65 at 16 ka and 17.7 ka, following the approach from ref. 39 (Supplementary Table 2).During ~22-65 ka, we adopted the WDC2014 age model for NGRIP δ 18 O using a 1.0063 scaling factor (WDC2014 age = 1.0063 × AICC2012 age) 66 to facilitate comparisons of our deep-water [CO 3

2−
] with the latest high-resolution atmospheric CO 2 record from the WDC ice core 3 (Fig. 2).Because WDC record covers only the past ~65 ka, we used the AICC2012 age model for NGRIP δ 18 O for ~65-120 ka.Between ~120 ka and 140 ka, we have tuned MD95-2039 G. bulloides δ 18 O to a nearby core ODP Site 976 G. bulloides δ 18 O, which has been mapped onto speleothem U-Th age (equivalent to the AICC2012 age model) 26 .During 140-150 ka, MD95-2039 G. bulloides δ 18 O is tuned to the synthetic NGRIP δ 18 O record 67 .Other ice-core records used in this study, including EDML δ 18 O 27 (WDC δ 18 O record is too short and thus not used) and compiled atmospheric CO 2 records [1][2][3][4] , are placed on the same age scales as those used to construct the MD95-2039 chronology.The robustness of our age model is confirmed by new and published 58,59 abundances of IRD and N. pachyderma (Extended Data Figs. 3 and 5).All age tie points for MD95-2039 are given in Supplementary Table 2.

Statistical analyses
Uncertainties associated with time series, including our new deep-water [CO 3 2− ], G. bulloides δ 18 O and ice-core parameters, were evaluated using a Monte Carlo-style approach 40,70 , considering errors with both the chronology (x axis) and reconstructions (y axis) (Fig. 3 and Supplementary Figs.2-7).Age errors are estimated using the OxCal programme 71 , while the individual [CO 3 2− ] error is estimated to be 10 μmol kg -1 (2σ).All data points were sampled separately and randomly 2,000 times within their chronological and [CO 3 2− ] uncertainties, and each iteration was then interpolated linearly.At each time step, the probability maximum and data distribution uncertainties of the 2,000 iterations were assessed.We present probability maxima (bold curves) and ±95% (shading; 2.5th-97.5thpercentile) probability intervals for the data distributions 40,70 .The same data processing approach was used to analyse δ 18 O and atmospheric CO 2 records using their respective uncertainties.
We constructed stadial evolutions of NGRIP δ ] for the five modes discussed in the main text (Fig. 3).First, each record was processed using the Monte Carlo approach described in the preceding.Second, the starting and ending times of stadials were identified on the basis of rapid changes in NGRIP δ 18 O, following ref. 3. Third, anomalies (Δ) of signals and associated 2σ uncertainties were calculated relative to the onset of each stadial using Monte Carlo-derived time series and ±95% probability intervals.Fourth, the duration of stadial was normalized (0-100%) to facilitate comparisons.Fifth, to obtain the overall trends, signals for Modes II and IV (number of stadials ≥ 5) were stacked by averaging the normalized anomalies of all stadials from each mode.To assess the uncertainties of the stacks, we calculated 2σ of all individual stadial signals and their associated uncertainties in each mode using the Monte Carlo-based approach 40,70 .

Timing of deep-water [CO 3 2− ] during Mode II stadials
Due to short durations associated with millennial events, any lead or lag relationship in timing between signals is inherently small, which urges caution about complications from bioturbation when interpreting proxy data.However, timing of geochemical signals can be well constrained by making the best use of high sedimentation and measurements of coexisting surface-and deep-dwelling foraminiferal shells from MD95-2039, as demonstrated previously for other cores from the Iberian Margin 44,55 .During Mode II stadials, the robustness of deep-water [CO 3 ] increases earlier than stadial terminations (Figs. 2  and 3) is supported by multiple lines of evidence.First, the high sedimentation rate (~20-30 cm ka -1 ) and large thicknesses of sediments (~30-70 cm) during these stadials at MD95-2039 not only allow detailed reconstructions but also substantially minimize any bioturbation effect (Extended Data Figs. 3 and 5).The early [CO 3

2−
] rise is unlikely an artefact of downward mixing of shells with high [CO 3 2− ] from the subsequent interstadial; otherwise, similar features would be expected to be more prevalent during shorter GSs with thinner sediment deposition, which would be subject to stronger bioturbation effects, a phenomenon not observed.Second, the [CO 3

2−
] rise before stadial terminations persists when examining changes in B/Ca (used to reconstruct [CO 3 2− ]) and δ 18 O (used to construct age model) against sediment depth (Extended Data Fig. 5).Because these measurements were made on coexisting benthic and planktic shells, bioturbation, if any, would shift depths of both species with minimal effect on their lead-lag relationship.We do acknowledge that any bioturbation effect also depends on the relative abundances of benthic versus planktic shells, which may change between stadials and interstadials, but a reliable means is yet lacking to obtain the 'original' species abundance information before bioturbation.Third, instead of relying on data from any single event, the early [CO 3 2− ] rise is revealed by [CO 3 2− ] and G. bulloides δ 18 O from multiple stadials on the basis of thorough statistical analyses to account for both chronology and reconstruction uncertainties (Supplementary Fig. 4).  1 .Note different y-axis scales between f and g.Vertical banding as shown in Fig. 2. We follow the literature 3 to define starting and ending dates of stadials.Given paired, high-resolution CH 4 and CO 2 measurements from the same ice core to, respectively, better constrain timings and magnitudes of stadial CO 2 changes, the latest WDC data 3,14 likely provide more accurate information about past rapid CO 2 changes.Antarctica warmed during all Northern Hemisphere stadials (a-c), but atmospheric CO 2 declined during many stadials (d-f) despite persistent operation of the bipolar seesaw.For stadials shorter than ~1000 years, no correlation is found between ΔCO 2 and stadial duration (f)., while dashed curve shows data from core SU81-18 (37.8°N, 10.2°W, 3,135 m) on the Iberian Margin 43 .Pa/Th data for HS11 are only available from ODP1063 (33.7°N, 57.6°W, 4,584 m) 29 , but caution the low resolution and potentially significant age errors with the record (e).During the late HS11, AMOC possibly slowed down although one (dubious, possibly related to age uncertainty) Pa/Th data point suggests the opposite change (e), while high-resolution and well-dated records show an IRD peak in MD95-2039 (Extended Data Fig. 3d) and a benthic δ 13 C decline in Iberian Margin core MD01-2444 (37.66°N, 10.1°W, 2,637 m) 44 (f).Compared to the late HS11, the late HS1 is characterized by (i) an overall colder Antarctica and a stalling of Antarctic warming 27 (d, i), (ii) a larger size but smaller retreat of Southern Ocean sea ice 77 (b, g), although note different views about using SSNa to reflect sea-ice extent 78,79 , and (iii) a much smaller atmospheric CO 2 rise 1 (c, h)..Model description: LOVECLIM consists of an OGCM, with a resolution of 3° × 3° and 20 vertical levels, coupled to dynamic-thermodynamic sea-ice model, a quasi-geostrophic atmospheric model (spectral T21 with 3 vertical levels), a dynamic global vegetation model, and a marine carbon cycle model 84 .The model is forced from 50 to 30 ka with time-varying boundary conditions including orbital parameters 85 , northern hemispheric ice-sheet orography and albedo 86 , and atmospheric CO 2 concentration 28 .In order to simulate AMOC changes characteristic of D-O cycles, freshwater is added in the North Atlantic (55°W-10°W, 50°N-65°N).For more details, see ref. 48.

Fig. 4 |
Fig. 4 | Schematic illustrating five stadial ventilation modes involving both Southern Ocean and North Atlantic processes.a-e, Modes I to V. From Mode I to Mode V, the strength of Southern Ocean ventilation changes in sequence: very strong (a), strong (b), strong but geographically restricted (c), weak (d), very weak (e).Shading shows seawater [CO 3 2− ] anomalies (opposite changes expected for DIC) relative to pre-stadial conditions.Blue and grey arrows denote extents of volumes ventilated by SSW and NSW, respectively, with thicker curves indicating stronger ventilation.Dashed and solid curves indicate possible watermass boundaries before and during stadials, respectively.Circle indicates the IV stadials.Although well-dated, high-resolution and paired εNd-[CO 3 2−

2 −]
due to outgassing, and via the North Atlantic, which tends to decrease [CO 3 2−

Fig. 2 |ExtendedFig. 4 |
Modern Atlantic meridional distributions of [and DIC (contours).Note that around the equator the hydrographic transect changes from the basin in the North Atlantic to the western basin in the South Atlantic.For simplicity, use NSW and SSW to indicate northernand southern-sourced deep waters, respectively.As expected from marine CO 2 system 22,76 , DIC and [CO 3 2− ] are inversely correlated.Also shown are sediment cores discussed in the main text.NSW: Northern-Sourced Waters; SSW: Southern-Sourced Waters; DIC: Dissolved Inorganic Carbon.Map is generated using Ocean Data View 75 based on hydrographic data for sites (inset) compiled by the GLODAP dataset .Extended Data Fig. 3 | Age model and sedimentation rate of core MD95-2039.a, NGRIP δ 18 O 25 .Crosses indicate age tie points.b, G. bulloides δ 18 O in MD95-2039 (red; this study) and ODP976 (olive) 26 .c, Sedimentation rate.d, IRD (HS11: this study; other times: ref. 59).e, N. pachyderma/all planktonics percentage ( < 32 ka: ref. 58; >32 ka: this study).Vertical dashed lines indicate some Heinrich events, as illustrated by prominent IRD peaks.f, age vs. depth.Green circles represent age control points.The envelope shows 1σ uncertainties based on OxCal analyses 71 .Data shown in c-f are from MD95-2039.Comparison of changes between HS1 and HS11.a-d, HS1 data.e-i, HS11 data.Vertical dashed lines separate HS1 and HS11 into early and late phases.In a, curve shows Pa/Th in cores GGC5 57.6°W, 4,550 m) CDH19 (33.7°N, 57.6°W, 4,541 m) from the Bermuda Rise

Fig. 6 |
Deep-water proxies and age model in core MD07-3076 from the South Atlantic.Data16,17 from MD07-3076 (b-g) are plotted against published age model 17 , based on alignment of planktonic faunal census data including N. pachyderma and G. bulloides abundances to Antarctic temperatures 27 (f-h).Curved arrows indicate that the [CO 3 2−] minimum at ~46 ka could be shifted to GS12, so that N. pachyderma and G. bulloides abundance peaks can be aligned with an Antarctic warming at ~44 ka (a, b, h).Overall, considering age model and reconstruction uncertainties, the suite of reconstructions at MD07-3076 (b-e) support our conclusions: strong Southern Ocean ventilation during Mode II stadials (HS4-6; yellow bands) and relatively weak ventilation during Mode IV stadials (grey bands).Extended Data Fig.7| Deep-water proxies, surface export, and age model in core TNO57-21 from the South Atlantic.The age model is based on alignment of G. bulloides δ 18 O (f) to Antarctic ice-core temperature proxies (g).Note that the age at ~46 ka is possibly too young: G. bulloides δ18 O minimum can be shifted older by ~2 ka, that the dashed rectangle area to HS5 as suggested by the horizontal arrow.This age adjustment would put alkenone changes82 (e) to be consistent with sub-Antarctic Zone export reductions during HSs, as expected from decreased dust deposition associated with poleward shift of southern westerlies38 .Taking this age uncertainty into account, deep waters at TNO57-21 experienced increases in [CO 3 2− ] 22,46 (b) and O 2 (refs.54,82) (c, d) during HS4-6, supporting enhanced Southern Ocean ventilation during Mode II stadials as proposed in the main text.Due to relatively low sampling resolution (∼10-15 cm per sample), it is yet difficult to resolve [CO 3 2− ] changes during many GSs (a, b).Yellow bands: Mode II stadials; Grey bands: Mode IV stadials.Data from refs.27,80,81.Extended Data Fig. 9 | Changes in northern North Atlantic CO 2 absorption in an Earth system model 48 .a, Timeseries of the simulated Atlantic meridional overturning circulation (AMOC).b, Integrated CO 2 flux in the northern North Atlantic (40-65°N).The horizontal dashed line indicates the mean flux value during 34-50 ka.Negative values indicate atmospheric CO 2 absorption by the northern North Atlantic.An AMOC weakening leads to greater CO 2 uptake in the northern North Atlantic during stadials.Results are from transient simulations of MIS 3 performed with an Earth system model of intermediate complexity, LOVECLIM

Fig. 10 |
Simulated climatic and biogeochemical responses to AMOC changes in an Earth system model 48 .Anomalies for (a, b) sea surface temperature (SST; unit °C), (c, d) surface [PO 4 3− ] (μmol/kg), and (e, f) sea-air CO 2 (μatm) for GS-IS (left column) and HS-IS (right column).Black and red curves show sea ice fronts (defined as 15% of sea ice density) during IS and GS/HS, respectively.Enhanced CO 2 absorption during HS and GS are driven by increased solubility pump (related to cooling; a, b) and biological pump (related to [PO 4 3− ] drawdown; c, d).GS = Greenland Stadial at 44.2 ka, IS = Interstadial at 46.6 ka, and HS = Heinrich Stadial at 48.8 ka.