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

J. Yu, R F Anderson, Z D Jin, X. Ji, D J R Thornalley, L. Wu, Nicolas

Thouveny, Y. Cai, L. Tan, F. Zhang, et al.

► To cite this version:

J. Yu, R F Anderson, Z D Jin, X. Ji, D J R Thornalley, et al. Millennial atmospheric CO2 changes linked to ocean ventilation modes over past 150,000 years. Nature Geoscience, 2023, 16 (12), pp.1166-1173. 10.1038/s41561-023-01297-x. hal-04562591

HAL Id: hal-04562591

https://hal.science/hal-04562591

Submitted on 29 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Millennial atmospheric CO2 changes linked to ocean ventilation modes over past 150,000 years

J. Yu, R F Anderson, Z D Jin, X. Ji, D J R Thornalley, L. Wu, Nicolas

Thouveny, Y. Cai, L. Tan, F. Zhang, et al.

► To cite this version:

J. Yu, R F Anderson, Z D Jin, X. Ji, D J R Thornalley, et al. Millennial atmospheric CO2 changes linked to ocean ventilation modes over past 150,000 years. Nature Geoscience, 2023, 16 (12), pp.1166-1173. 10.1038/s41561-023-01297-x. hal-04562591

HAL Id: hal-04562591

https://hal.science/hal-04562591v1

Submitted on 29 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Millennial atmospheric CO₂ changes linked to ocean ventilation modes over past 150,000 years

Received: 23 January 2023

Accepted: 20 September 2023

Published online: 26 October 2023

<figure><img src="image_3.png" /></figure>

Check for updates

J. Yu Ⓓ1,2,3 Ⓒ, R. F. Anderson Ⓓ4, Z. D. Jin Ⓓ3,5,6, X. Ji Ⓓ2, D. J. R. Thornalley Ⓓ7, L. Wu Ⓓ1,8, N. Thouveny⁹, Y. Cai Ⓓ6, L. Tan Ⓓ3,6, F. Zhang Ⓓ3, L. Menviel Ⓓ10, J. Tian Ⓓ11, X. Xie¹¹, E. J. Rohling²,¹² & J. F. McManus Ⓓ4

Ice core measurements show diverse atmospheric CO₂ 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 CO₂, 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 CO₂ 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 CO₂ 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 CO₂. Our findings demonstrate that millennial-scale changes in deep-sea carbon storage and atmospheric CO₂ 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.

Ice-core atmospheric CO₂ records show ubiquitous millennial-scale fluctuations superimposed on longer-term glacial–interglacial cycles during the past 800 thousand years (ka)<fnref n="1" />⁻⁴. Ventilation of the carbon-rich deep ocean can critically affect atmospheric CO₂ on various timescales⁵⁻⁷.

A prevailing view⁸⁻¹⁰ to explain past rises in atmospheric CO₂ invokes strengthened Southern Ocean ventilation, which facilitates carbon release from the deep ocean. On the basis of a process called the bipolar seesaw¹⁰,¹¹, Antarctic ice-core temperature may be used to infer changes

<footnote n="1">¹Laoshan Laboratory, Qingdao, China. ²Research School of Earth Sciences, The Australian National University, Canberra, Australia. ³SKLLQG, Institute of Earth Environment, Chinese Academy of Sciences, Xi’an, China. ⁴Lamont–Doherty Earth Observatory, Columbia University, Palisades, NY, USA. ⁵Open Studio for Oceanic–Continental Climate and Environment Changes, Laoshan Laboratory, Qingdao, China. ⁶Institute of Global Environmental Change, Xi’an Jiaotong University, Xi’an, China. ⁷Department of Geography, University College London, London, UK. ⁸Physical Oceanography Laboratory, Ocean University of China, Qingdao, China. ⁹CEREGE, Aix-Marseille Univ, CNRS, IRD, INRA, Coll de Fr, Aix en Provence, France. ¹⁰Climate Change Research Centre, Earth and Sustainability Science Research Centre, University of New South Wales, Sydney, New South Wales, Australia. ¹¹State Key Laboratory of Marine Geology, Tongji University, Shanghai, China. ¹²Ocean and Earth Science, University of Southampton, National Oceanography Centre, Southampton, UK. ☒e-mail: jiminyuanu@gmail.com</footnote>
<figure><img src="image_2.png" /></figure>

# Millennial atmospheric CO₂ changes linked to ocean ventilation modes over past 150,000 years

Received: 23 January 2023

Accepted: 20 September 2023

Published online: 26 October 2023

<figure><img src="image_3.png" /></figure>

Check for updates

J. Yu Ⓓ1,2,3 Ⓒ, R. F. Anderson Ⓓ4, Z. D. Jin Ⓓ3,5,6, X. Ji Ⓓ2, D. J. R. Thornalley Ⓓ7, L. Wu Ⓓ1,8, N. Thouveny⁹, Y. Cai Ⓓ6, L. Tan Ⓓ3,6, F. Zhang Ⓓ3, L. Menviel Ⓓ10, J. Tian Ⓓ11, X. Xie¹¹, E. J. Rohling²,¹² & J. F. McManus Ⓓ4

Ice core measurements show diverse atmospheric CO₂ 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 CO₂, 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 CO₂ 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 CO₂ 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 CO₂. Our findings demonstrate that millennial-scale changes in deep-sea carbon storage and atmospheric CO₂ 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.

Ice-core atmospheric CO₂ records show ubiquitous millennial-scale fluctuations superimposed on longer-term glacial–interglacial cycles during the past 800 thousand years (ka)<fnref n="1" />⁻⁴. Ventilation of the carbon-rich deep ocean can critically affect atmospheric CO₂ on various timescales⁵⁻⁷.

A prevailing view⁸⁻¹⁰ to explain past rises in atmospheric CO₂ invokes strengthened Southern Ocean ventilation, which facilitates carbon release from the deep ocean. On the basis of a process called the bipolar seesaw¹⁰,¹¹, Antarctic ice-core temperature may be used to infer changes

<footnote n="1">¹Laoshan Laboratory, Qingdao, China. ²Research School of Earth Sciences, The Australian National University, Canberra, Australia. ³SKLLQG, Institute of Earth Environment, Chinese Academy of Sciences, Xi’an, China. ⁴Lamont–Doherty Earth Observatory, Columbia University, Palisades, NY, USA. ⁵Open Studio for Oceanic–Continental Climate and Environment Changes, Laoshan Laboratory, Qingdao, China. ⁶Institute of Global Environmental Change, Xi’an Jiaotong University, Xi’an, China. ⁷Department of Geography, University College London, London, UK. ⁸Physical Oceanography Laboratory, Ocean University of China, Qingdao, China. ⁹CEREGE, Aix-Marseille Univ, CNRS, IRD, INRA, Coll de Fr, Aix en Provence, France. ¹⁰Climate Change Research Centre, Earth and Sustainability Science Research Centre, University of New South Wales, Sydney, New South Wales, Australia. ¹¹State Key Laboratory of Marine Geology, Tongji University, Shanghai, China. ¹²Ocean and Earth Science, University of Southampton, National Oceanography Centre, Southampton, UK. ☒e-mail: jiminyuanu@gmail.com</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Deep-water [CO₃²⁻] at site MD95-2039 from the Iberian Margin. **a**, North Greenland Ice Core Project (NGRIP) ice-core δ¹⁸O (ref. 25). **b**, *G. bulloides* δ¹⁸O (three-point-running mean) at MD95-2039 (orange; this study) and ODP 976 (olive²⁶). **c**, Deep-water [CO₃²⁻] at MD95-2039 (this study). **d**, Compiled North Atlantic sediment Pa/Th, an AMOC strength proxy²⁹,³²,³³. **e**, EDML ice-core δ¹⁸O (ref. 27), reflecting temperature changes in the Atlantic sector of Antarctica.</figcaption></figure>

**f**, Atmospheric CO₂ (refs. 1–4). On the basis of temporal [CO₃²⁻] evolutions alongside direction and magnitude of contemporary atmospheric CO₂ changes, we identify five ventilation modes during stadials (vertical colour bandings). ODP, ocean drilling program; EDML, European Project for Ice Coring in Antarctica (EPICA) ice core from Dronning Maud Land.

in southern westerly winds and, thus, Southern Ocean ventilation¹² (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₂ as seen in ice-core records⁶,⁸,¹⁰,¹¹. Yet it remains unclear why atmospheric CO₂ did not rise during many stadials despite clear Antarctic warmings²,³. The latest high-resolution atmospheric CO₂ record, whose timing and magnitude can be firmly constrained by paired CH₄ measurements from the same ice core, shows no net CO₂ increase during a considerable number of stadials when the bipolar seesaw persistently operated³,¹⁴ (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¹⁵⁻²⁰ 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₂ changes, is yet to be assessed. Even less is known about the impacts on atmospheric CO₂ 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₂ variations on millennial timescales.

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₂ exchange in the deep-water formation regions and respired carbon accumulation during water-mass transportation (Methods). We employ carbonate ion concentrations ([CO₃²⁻]; 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₂ changes¹⁹,²¹,²². Everything else being equal, enhanced Southern Ocean ventilation would reduce deep-sea DIC, raising deep-water [CO₃²⁻] and atmospheric CO₂ (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₃²⁻] and atmospheric CO₂ (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₂—is hampered by a paucity of high-quality [CO₃²⁻] records. Thus, high-resolution [CO₃²⁻] records with robust chronologies are warranted to better understand abrupt ocean ventilation mechanisms.

## Benchmark [CO₃²⁻] record

In this Article, we present the first high-resolution deep-water [CO₃²⁻] 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₃²⁻] record is reconstructed using B/Ca in benthic
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Deep-water [CO₃²⁻] at site MD95-2039 from the Iberian Margin. **a**, North Greenland Ice Core Project (NGRIP) ice-core δ¹⁸O (ref. 25). **b**, *G. bulloides* δ¹⁸O (three-point-running mean) at MD95-2039 (orange; this study) and ODP 976 (olive²⁶). **c**, Deep-water [CO₃²⁻] at MD95-2039 (this study). **d**, Compiled North Atlantic sediment Pa/Th, an AMOC strength proxy²⁹,³²,³³. **e**, EDML ice-core δ¹⁸O (ref. 27), reflecting temperature changes in the Atlantic sector of Antarctica.</figcaption></figure>

**f**, Atmospheric CO₂ (refs. 1–4). On the basis of temporal [CO₃²⁻] evolutions alongside direction and magnitude of contemporary atmospheric CO₂ changes, we identify five ventilation modes during stadials (vertical colour bandings). ODP, ocean drilling program; EDML, European Project for Ice Coring in Antarctica (EPICA) ice core from Dronning Maud Land.

in southern westerly winds and, thus, Southern Ocean ventilation¹² (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₂ as seen in ice-core records⁶,⁸,¹⁰,¹¹. Yet it remains unclear why atmospheric CO₂ did not rise during many stadials despite clear Antarctic warmings²,³. The latest high-resolution atmospheric CO₂ record, whose timing and magnitude can be firmly constrained by paired CH₄ measurements from the same ice core, shows no net CO₂ increase during a considerable number of stadials when the bipolar seesaw persistently operated³,¹⁴ (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¹⁵⁻²⁰ 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₂ changes, is yet to be assessed. Even less is known about the impacts on atmospheric CO₂ 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₂ variations on millennial timescales.

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₂ exchange in the deep-water formation regions and respired carbon accumulation during water-mass transportation (Methods). We employ carbonate ion concentrations ([CO₃²⁻]; 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₂ changes¹⁹,²¹,²². Everything else being equal, enhanced Southern Ocean ventilation would reduce deep-sea DIC, raising deep-water [CO₃²⁻] and atmospheric CO₂ (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₃²⁻] and atmospheric CO₂ (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₂—is hampered by a paucity of high-quality [CO₃²⁻] records. Thus, high-resolution [CO₃²⁻] records with robust chronologies are warranted to better understand abrupt ocean ventilation mechanisms.

## Benchmark [CO₃²⁻] record

In this Article, we present the first high-resolution deep-water [CO₃²⁻] 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₃²⁻] record is reconstructed using B/Ca in benthic
<figure><img src="image_5.png" /></figure>

**Fig. 2** | Detailed view over 30–65 thousand years ago. **a**, Greenland NGRIP ice-core $\delta^{18}$O (ref. 25). **b**, *G. bulloides* $\delta^{18}$O (three-point-running mean) at MD95-2039 (orange; this study). **c**, Deep-water [CO$_3$$^{2-}$] at MD95-2039 (this study). **d**, Compiled North Atlantic sediment Pa/Th, an AMOC strength proxy$^{29,32,33}$. **e**, EDML ice-core $\delta^{18}$O (ref. 27), reflecting temperature changes in the Atlantic sector of Antarctica. **f**, Atmospheric CO$_2$ (ref. 3). **g**, Deep-water [CO$_3$$^{2-}$] at MD07-3076$^{17}$. The curved arrow indicates that the [CO$_3$$^{2-}$] minimum at -45 thousand years ago could be shifted to GS12 (Extended Data Fig. 6). **h**, $\varepsilon$Nd in cores ODP 1063 (triangles)$^{29}$ and TNO57-21 (circles)$^{45}$ from the North and South Atlantic, respectively. See Extended Data Fig. 2 for core locations. During Mode II stadials (yellow bandings), increases in [CO$_3$$^{2-}$] before stadial terminations (c,g) suggest deep-sea carbon loss driven by strong Southern Ocean ventilation, conducive for large atmospheric CO$_2$ increases (**f**). During Mode IV stadials (grey bandings), sustained low [CO$_3$$^{2-}$] (c,g) indicate increased accumulation of respired carbon in the deep ocean driven by AMOC weakening (**d**), helping to explain the lack of any atmospheric CO$_2$ rise (**f**). See text for details.
<figure><img src="image_5.png" /></figure>

**Fig. 2** | Detailed view over 30–65 thousand years ago. **a**, Greenland NGRIP ice-core $\delta^{18}$O (ref. 25). **b**, *G. bulloides* $\delta^{18}$O (three-point-running mean) at MD95-2039 (orange; this study). **c**, Deep-water [CO$_3$$^{2-}$] at MD95-2039 (this study). **d**, Compiled North Atlantic sediment Pa/Th, an AMOC strength proxy$^{29,32,33}$. **e**, EDML ice-core $\delta^{18}$O (ref. 27), reflecting temperature changes in the Atlantic sector of Antarctica. **f**, Atmospheric CO$_2$ (ref. 3). **g**, Deep-water [CO$_3$$^{2-}$] at MD07-3076$^{17}$. The curved arrow indicates that the [CO$_3$$^{2-}$] minimum at -45 thousand years ago could be shifted to GS12 (Extended Data Fig. 6). **h**, $\epsilon$Nd in cores ODP 1063 (triangles)$^{29}$ and TNO57-21 (circles)$^{45}$ from the North and South Atlantic, respectively. See Extended Data Fig. 2 for core locations. During Mode II stadials (yellow bandings), increases in [CO$_3$$^{2-}$] before stadial terminations (c,g) suggest deep-sea carbon loss driven by strong Southern Ocean ventilation, conducive for large atmospheric CO$_2$ increases (f). During Mode IV stadials (grey bandings), sustained low [CO$_3$$^{2-}$] (c,g) indicate increased accumulation of respired carbon in the deep ocean driven by AMOC weakening (d), helping to explain the lack of any atmospheric CO$_2$ rise (f). See text for details.
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Five modes of stadial ventilation.** a–e, Anomalies (Δ) in NGRIP ice-core δ¹⁸O²⁵ (a), *G. bulloides* δ¹⁸O (b) and deep-water [CO₃²⁻] (c) at MD95-2039, atmospheric CO₂ (d) and Antarctic EDML ice-core δ¹⁸O (e). Bold curves show probability maxima, with envelopes representing 2σ uncertainties. The combinations of [CO₃²⁻] (c) and atmospheric CO₂ (d) evolutions define five ventilation modes (Modes I to V). The delayed atmospheric CO₂ 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*</figcaption></figure>

δ¹⁸O²⁶ and adjusted by a factor of ~3 to facilitate plotting. Note different y-axis scales, especially for atmospheric ΔCO₂, 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. 2–7 for statistical details.

foraminifer *Cibicoides wuellerstorfii*, with a reconstruction uncertainty (2σ) of ~10 µmol kg⁻¹ (ref. 24). Altogether, 588 measurements were obtained, yielding an average temporal resolution of 256 years over the past 150 ka. We also constructed a robust chronology for MD95-2039 based on six radiocarbon dates and aligning of MD95-2039 *Globigerina bulloides* δ¹⁸O (δ¹⁸O = [(¹⁸O/¹⁶O)sample/¹⁸O/¹⁶O]standard – 1] × 1,000‰; *n* = 1,257) to benchmark records²⁵,²⁶ (Extended Data Fig. 3). See Methods for details.

Figure 1 shows our [CO₃²⁻] record alongside ice-core δ¹⁸O and atmospheric CO₂, all placed on the same chronology⁴,²⁵,²⁷, 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₃²⁻] fluctuations accompanied by diverse atmospheric CO₂ 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₃²⁻] and contemporaneous ice-core CO₂ during the progress of stadials (Methods and Supplementary Figs. 2–7). On the basis of their [CO₃²⁻]-CO₂ 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²⁸ or marine²⁹ archives.

## Five modes of stadial ocean ventilation

We first assess Heinrich stadials (HSs) associated with ~40–55 ppm atmospheric CO₂ rises at glacial terminations. During HS1 and HS11, our record shows early [CO₃²⁻] increases followed by decreases (Mode I; Fig. 3). On millennial timescales when the global alkalinity change is relatively minor³⁰, a [CO₃²⁻] increase (decrease) probably reflects a DIC decrease (increase)²². The early [CO₃²⁻] increases and associated DIC decreases cannot be explained by northward expansion of carbon-rich SSW³¹. During these stadials, Pa/Th data²⁹,³²,³³ indicate a slowdown of AMOC (Fig. 1d), which would promote respired carbon accumulation³⁴ with an effect to decrease [CO₃²⁻] 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¹⁰,¹¹,³⁵, which would enhance deep-sea ventilation via the Southern Ocean as suggested by proxy and model results⁸,²⁰,³⁶,³⁷. 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⁹,¹²,³⁸ 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¹⁵,¹⁶,¹⁸,²¹,³⁷,³⁹. 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²⁶,⁴⁰,⁴¹ (Fig. 1d). Indeed, Antarctic ice-core records²⁷ 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₂, marked Southern Ocean ventilation enhancement may have played a critical role in initiating and/or sustaining glacial terminations (Fig. 4a).

During late HS1, deep-water [CO₃²⁻] declined at MD95-2039, with lowest values occurring close to the peak of ice-rafted debris (IRD) deposition at the Iberian Margin⁴² (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⁴³ (Extended Data Fig. 4). If so, our data imply a DIC increase in the Atlantic below ~3 km, which would slow atmospheric CO₂ rise, consistent with ice-core data¹. The late HS11 [CO₃²⁻] decline can also be attributed to further AMOC weakening, supported by a coeval benthic δ¹³C decrease at the Iberian Margin⁴⁴ (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¹⁰,¹¹ (Fig. 1e). Although further investigation is needed, we speculate that a substantial deep-sea volume to the south of 40° N
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Five modes of stadial ventilation.** a–e, Anomalies (Δ) in NGRIP ice-core δ¹⁸O²⁵ (a), *G. bulloides* δ¹⁸O (b) and deep-water [CO₃²⁻] (c) at MD95-2039, atmospheric CO₂ (d) and Antarctic EDML ice-core δ¹⁸O (e). Bold curves show probability maxima, with envelopes representing 2σ uncertainties. The combinations of [CO₃²⁻] (c) and atmospheric CO₂ (d) evolutions define five ventilation modes (Modes I to V). The delayed atmospheric CO₂ 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*</figcaption></figure>

δ¹⁸O²⁶ and adjusted by a factor of ~3 to facilitate plotting. Note different y-axis scales, especially for atmospheric ΔCO₂, 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. 2–7 for statistical details.

foraminifer *Cibicidoides wuellerstorfii*, with a reconstruction uncertainty (2σ) of ~10 µmol kg⁻¹ (ref. 24). Altogether, 588 measurements were obtained, yielding an average temporal resolution of 256 years over the past 150 ka. We also constructed a robust chronology for MD95-2039 based on six radiocarbon dates and aligning of MD95-2039 *Globigerina bulloides* δ¹⁸O (δ¹⁸O = [(¹⁸O/¹⁶O)sample/¹⁸O/¹⁶O]standard – 1] × 1,000‰; *n* = 1,257) to benchmark records²⁵,²⁶ (Extended Data Fig. 3). See Methods for details.

Figure 1 shows our [CO₃²⁻] record alongside ice-core δ¹⁸O and atmospheric CO₂, all placed on the same chronology⁴,²⁵,²⁷, 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₃²⁻] fluctuations accompanied by diverse atmospheric CO₂ 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₃²⁻] and contemporaneous ice-core CO₂ during the progress of stadials (Methods and Supplementary Figs. 2–7). On the basis of their [CO₃²⁻]-CO₂ 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²⁸ or marine²⁹ archives.

## Five modes of stadial ocean ventilation

We first assess Heinrich stadials (HSs) associated with ~40–55 ppm atmospheric CO₂ rises at glacial terminations. During HS1 and HS11, our record shows early [CO₃²⁻] increases followed by decreases (Mode I; Fig. 3). On millennial timescales when the global alkalinity change is relatively minor³⁰, a [CO₃²⁻] increase (decrease) probably reflects a DIC decrease (increase)²². The early [CO₃²⁻] increases and associated DIC decreases cannot be explained by northward expansion of carbon-rich SSW³¹. During these stadials, Pa/Th data²⁹,³²,³³ indicate a slowdown of AMOC (Fig. 1d), which would promote respired carbon accumulation³⁴ with an effect to decrease [CO₃²⁻] 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¹⁰,¹¹,³⁵, which would enhance deep-sea ventilation via the Southern Ocean as suggested by proxy and model results⁸,²⁰,³⁶,³⁷. 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⁹,¹²,³⁸ 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¹⁵,¹⁶,¹⁸,²¹,³⁷,³⁹. 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²⁶,⁴⁰,⁴¹ (Fig. 1d). Indeed, Antarctic ice-core records²⁷ 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₂, marked Southern Ocean ventilation enhancement may have played a critical role in initiating and/or sustaining glacial terminations (Fig. 4a).

During late HS1, deep-water [CO₃²⁻] declined at MD95-2039, with lowest values occurring close to the peak of ice-rafted debris (IRD) deposition at the Iberian Margin⁴² (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⁴³ (Extended Data Fig. 4). If so, our data imply a DIC increase in the Atlantic below ~3 km, which would slow atmospheric CO₂ rise, consistent with ice-core data¹. The late HS11 [CO₃²⁻] decline can also be attributed to further AMOC weakening, supported by a coeval benthic δ¹³C decrease at the Iberian Margin⁴⁴ (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¹⁰,¹¹ (Fig. 1e). Although further investigation is needed, we speculate that a substantial deep-sea volume to the south of 40° N
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**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₃²⁻] 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 water-mass boundaries before and during stadials, respectively. Circle indicates the</figcaption></figure>

location of MD95-2039. Changes in Southern Ocean ventilation may be linked to fluctuations in southern westerly wind and Southern Ocean sea-ice cover (grey rectangles), affecting air-sea CO₂ exchange (red arrows; thicker lines indicate stronger CO₂ exchange). In addition, weakened AMOC would reduce ventilation and promote respired carbon accumulation in the deep ocean when southern ventilation was relatively weak (d,e). Thus, Southern Ocean and North Atlantic processes jointly govern ventilation modes and respired carbon storage of the ocean interior and thereby magnitudes and directions of atmospheric CO₂ changes. SO, Southern Ocean.

(latitude of MD95-2039) was well ventilated, which counteracted any increased carbon sequestration in the deep Atlantic volume north of ~40°N and thereby contributed to maintaining the high atmospheric CO₂ rise rate throughout HS11.

Mode II stadials were characterized by ~5–15 ppm atmospheric CO₂ rises at intermediate climate states (Figs. 1–3). During HS4–6 and 7b, deep-water [CO₃²⁻] at MD95-2039 underwent a brief decline followed by a net increase (Fig. 3). The initial [CO₃²⁻] 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²⁹,³²,³³ (Fig. 2d). The subsequent [CO₃²⁻] 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. wuellerstorfii B/Ca (used for deep-water [CO₃²⁻]) and

G. bulloides δ¹⁸O (used for chronology) in the same core (Extended Data Fig. 5 and Methods). Existing εNd data²⁹,⁴⁵ 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¹⁰,¹¹,²⁷,³⁶ (Fig. 2e), the DIC decrease can be readily explained by carbon loss owing to enhanced Southern Ocean ventilation. This is supported, within uncertainties, by independent deep South Atlantic [CO₃²⁻], radiocarbon and oxygenation reconstructions¹⁵⁻¹⁷,²²,⁴⁶ (Fig. 2g and Extended Data Figs. 6 and 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₂ 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₃²⁻] declines at MD95-2039 (Fig. 4b).
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**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₃²⁻] 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 water-mass boundaries before and during stadials, respectively. Circle indicates the</figcaption></figure>

location of MD95-2039. Changes in Southern Ocean ventilation may be linked to fluctuations in southern westerly wind and Southern Ocean sea-ice cover (grey rectangles), affecting air-sea CO₂ exchange (red arrows; thicker lines indicate stronger CO₂ exchange). In addition, weakened AMOC would reduce ventilation and promote respired carbon accumulation in the deep ocean when southern ventilation was relatively weak (d,e). Thus, Southern Ocean and North Atlantic processes jointly govern ventilation modes and respired carbon storage of the ocean interior and thereby magnitudes and directions of atmospheric CO₂ changes. SO, Southern Ocean.

(latitude of MD95-2039) was well ventilated, which counteracted any increased carbon sequestration in the deep Atlantic volume north of ~40°N and thereby contributed to maintaining the high atmospheric CO₂ rise rate throughout HS11.

Mode II stadials were characterized by ~5–15 ppm atmospheric CO₂ rises at intermediate climate states (Figs. 1–3). During HS4–6 and 7b, deep-water [CO₃²⁻] at MD95-2039 underwent a brief decline followed by a net increase (Fig. 3). The initial [CO₃²⁻] 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²⁹,³²,³³ (Fig. 2d). The subsequent [CO₃²⁻] 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. wuellerstorfii* B/Ca (used for deep-water [CO₃²⁻]) and

*G. bulloides* δ¹⁸O (used for chronology) in the same core (Extended Data Fig. 5 and Methods). Existing εNd data²⁹,⁴⁵ 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¹⁰,¹¹,²⁷,³⁶ (Fig. 2e), the DIC decrease can be readily explained by carbon loss owing to enhanced Southern Ocean ventilation. This is supported, within uncertainties, by independent deep South Atlantic [CO₃²⁻], radiocarbon and oxygenation reconstructions¹⁵⁻¹⁷,²²,⁴⁶ (Fig. 2g and Extended Data Figs. 6 and 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₂ 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₃²⁻] declines at MD95-2039 (Fig. 4b).
Although accompanied by -15–20 ppm atmospheric CO₂ rises, deep-water [CO₃²⁻] at MD95-2039 decreased and remained low during Mode III stadials, including the Younger Dryas (YD) and HS8 (Figs. 1 and 3). The YD occurred under relatively warm Antarctic conditions (like late HS11) during the last deglaciation, while HS8 stands as the longest stadial within Marine Isotope Stage (MIS) 5. During the YD, significant increases in deep-water [CO₃²⁻] 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¹⁵,⁴⁷. Increases in deep-water [CO₃²⁻] and oxygenation are also documented at TNO57-21 during HS8²² (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⁸,¹²,²⁷ (Fig. 1e). The associated carbon releases would have contributed to contemporaneous atmospheric CO₂ rises. Yet sustained low [CO₃²⁻] 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₂ generally displays muted net changes during Greenland stadials (GS) 4, 6–8, 10–12, 16, 17 and 25 and HS2, 3, 7a and 10 under intermediate climate states (Figs. 1–3). During these stadials, deep-water [CO₃²⁻] at MD95-2039 varied similarly to Mode III stadials, but no deep-water [CO₃²⁻] or oxygenation increase is seen at South Atlantic sites TNO57-21 and MD07-3076 (44° S, 14° W, 3,770 m)¹⁷,²² (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)²⁷ 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¹³, Antarctica warming would reduce Antarctic Zone sea-ice extent and SSW solubility pump efficiency. Thus, Southern Ocean changes probably promoted CO₂ outgassing with an effect to increase deep Atlantic [CO₃²⁻] during Mode IV stadials. Although well-dated, high-resolution and paired εNd-[CO₃²⁻] reconstructions are desired, published deep Atlantic εNd records²⁹,⁴⁵ show no consistent variations during transitions into Mode IV stadials (Fig. 2h), arguing against increased mixing of low-[CO₃²⁻] SSW as an exclusive explanation for the observed [CO₃²⁻] declines at MD95-2039. We suggest that much of these [CO₃²⁻] declines and associated DIC rises must have been caused by increased accumulation of respired carbon driven by weakened AMOC²⁹,³²,³³ (Fig. 2d), which is supported by model results³¹,⁴⁸. 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⁹,²³, this would counteract the effect of carbon loss from the shallow Southern Ocean and thus curb atmospheric CO₂ rises during Mode IV stadials (Fig. 4d). Indeed, published model simulations⁴⁸ show increased atmospheric CO₂ 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₃²⁻] decrease at MD95-2039 during GS26 when atmospheric CO₂ 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₃²⁻] 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¹¹, both Antarctica and Greenland cooled during GS26²⁵,²⁷ (Fig. 1a,e). In a wider perspective, GS26 was characterized by a global cooling⁴⁹, sea-ice expansion surrounding Antarctica and plausibly an equatorward shift of southern westerlies, which would suppress Southern Ocean ventilation¹⁰,¹² and thus deep-sea carbon leakage. Moreover, concomitant North Atlantic solubility pump enhancement due to cooling²⁵,²⁷ would have promoted atmospheric CO₂ sequestration in NSW. We therefore suggest that a bipolar synergy operated to lower atmospheric CO₂ during GS26, providing a positive feedback to global cooling at the last glacial inception (Fig. 4e).

## Bipolar control on millennial atmospheric CO₂

Our [CO₃²⁻] reconstructions, alongside published data⁸,¹⁵,¹⁶, suggest Southern Ocean ventilation as a key factor to modulate millennial deep Atlantic carbon storage and atmospheric CO₂ 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₂ 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₂ 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₂. As demonstrated by Mode IV data and supported by models³¹,³⁴,⁵⁰, 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₂ (Supplementary Fig. 1), a mechanism widely invoked to explain atmospheric CO₂ declines on glacial–interglacial timescales⁹,¹⁵,²³,³⁴,⁵¹ but largely overlooked on millennial timescales. In fact, all stadials appear to be accompanied by reduced deep-water [CO₃²⁻] 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₃²⁻] (Fig. 3). Our inferred effect on atmospheric CO₂ is supported by model simulations⁴⁸, which show increased North Atlantic CO₂ absorption during all studied stadials (Extended Data Figs. 9 and 10). Under such conditions, if Southern Ocean CO₂ outgassing is not strong enough, atmospheric CO₂ could decrease. In other words, atmospheric CO₂ 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₂ variations.

Our proposed bipolar control helps resolve a long-standing puzzle² regarding mechanisms for disparate atmospheric CO₂ responses among MIS 3 stadials. The persistent anti-phased temperature relationship between polar ice cores²⁷ suggests sustained operation of the bipolar seesaw during MIS 3 (Fig. 2). One would expect¹⁰⁻¹² enhanced Southern Ocean CO₂ outgassing and thus atmospheric CO₂ rises during all MIS 3 stadials, which is at odds with observations²,³ (Fig. 2). The lack of atmospheric CO₂ rises during many stadials must indicate concomitant processes to offset CO₂ outgassing from the Southern Ocean. Atmospheric CO₂ is affected by carbon exchange with both
Although accompanied by -15–20 ppm atmospheric CO₂ rises, deep-water [CO₃²⁻] at MD95-2039 decreased and remained low during Mode III stadials, including the Younger Dryas (YD) and HS8 (Figs. 1 and 3). The YD occurred under relatively warm Antarctic conditions (like late HS11) during the last deglaciation, while HS8 stands as the longest stadial within Marine Isotope Stage (MIS) 5. During the YD, significant increases in deep-water [CO₃²⁻] 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¹⁵,⁴⁷. Increases in deep-water [CO₃²⁻] and oxygenation are also documented at TNO57-21 during HS8²² (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⁸,¹²,²⁷ (Fig. 1e). The associated carbon releases would have contributed to contemporaneous atmospheric CO₂ rises. Yet sustained low [CO₃²⁻] 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₂ generally displays muted net changes during Greenland stadials (GS) 4, 6–8, 10–12, 16, 17 and 25 and HS2, 3, 7a and 10 under intermediate climate states (Figs. 1–3). During these stadials, deep-water [CO₃²⁻] at MD95-2039 varied similarly to Mode III stadials, but no deep-water [CO₃²⁻] or oxygenation increase is seen at South Atlantic sites TNO57-21 and MD07-3076 (44° S, 14° W, 3,770 m)¹⁷,²² (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)²⁷ 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¹³, Antarctica warming would reduce Antarctic Zone sea-ice extent and SSW solubility pump efficiency. Thus, Southern Ocean changes probably promoted CO₂ outgassing with an effect to increase deep Atlantic [CO₃²⁻] during Mode IV stadials. Although well-dated, high-resolution and paired εNd-[CO₃²⁻] reconstructions are desired, published deep Atlantic εNd records²⁹,⁴⁵ show no consistent variations during transitions into Mode IV stadials (Fig. 2h), arguing against increased mixing of low-[CO₃²⁻] SSW as an exclusive explanation for the observed [CO₃²⁻] declines at MD95-2039. We suggest that much of these [CO₃²⁻] declines and associated DIC rises must have been caused by increased accumulation of respired carbon driven by weakened AMOC²⁹,³²,³³ (Fig. 2d), which is supported by model results³¹,⁴⁸. 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⁹,²³, this would counteract the effect of carbon loss from the shallow Southern Ocean and thus curb atmospheric CO₂ rises during Mode IV stadials (Fig. 4d). Indeed, published model simulations⁴⁸ show increased atmospheric CO₂ 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₃²⁻] decrease at MD95-2039 during GS26 when atmospheric CO₂ 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₃²⁻] 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¹¹, both Antarctica and Greenland cooled during GS26²⁵,²⁷ (Fig. 1a,e). In a wider perspective, GS26 was characterized by a global cooling⁴⁹, sea-ice expansion surrounding Antarctica and plausibly an equatorward shift of southern westerlies, which would suppress Southern Ocean ventilation¹⁰,¹² and thus deep-sea carbon leakage. Moreover, concomitant North Atlantic solubility pump enhancement due to cooling²⁵,²⁷ would have promoted atmospheric CO₂ sequestration in NSW. We therefore suggest that a bipolar synergy operated to lower atmospheric CO₂ during GS26, providing a positive feedback to global cooling at the last glacial inception (Fig. 4e).

## Bipolar control on millennial atmospheric CO₂

Our [CO₃²⁻] reconstructions, alongside published data⁸,¹⁵,¹⁶, suggest Southern Ocean ventilation as a key factor to modulate millennial deep Atlantic carbon storage and atmospheric CO₂ 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₂ 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₂ 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₂. As demonstrated by Mode IV data and supported by models³¹,³⁴,⁵⁰, 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₂ (Supplementary Fig. 1), a mechanism widely invoked to explain atmospheric CO₂ declines on glacial–interglacial timescales⁹,¹⁵,²³,³⁴,⁵¹ but largely overlooked on millennial timescales. In fact, all stadials appear to be accompanied by reduced deep-water [CO₃²⁻] 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₃²⁻] (Fig. 3). Our inferred effect on atmospheric CO₂ is supported by model simulations⁴⁸, which show increased North Atlantic CO₂ absorption during all studied stadials (Extended Data Figs. 9 and 10). Under such conditions, if Southern Ocean CO₂ outgassing is not strong enough, atmospheric CO₂ could decrease. In other words, atmospheric CO₂ 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₂ variations.

Our proposed bipolar control helps resolve a long-standing puzzle² regarding mechanisms for disparate atmospheric CO₂ responses among MIS 3 stadials. The persistent anti-phased temperature relationship between polar ice cores²⁷ suggests sustained operation of the bipolar seesaw during MIS 3 (Fig. 2). One would expect¹⁰⁻¹² enhanced Southern Ocean CO₂ outgassing and thus atmospheric CO₂ rises during all MIS 3 stadials, which is at odds with observations²,³ (Fig. 2). The lack of atmospheric CO₂ rises during many stadials must indicate concomitant processes to offset CO₂ outgassing from the Southern Ocean. Atmospheric CO₂ is affected by carbon exchange with both
ocean and land biosphere reservoirs, but how the terrestrial biosphere carbon changed during stadials remains poorly constrained²,⁴¹. Here we provide critical data evidence (Fig. 2) to demonstrate that respired carbon accumulation, driven by reduced North Atlantic ventilation³¹,⁴⁸,⁵⁰, at times cancelled or surpassed the effect of Southern Ocean outgassing, explaining the observations of little net change or even declines in atmospheric CO₂ during some MIS 3 stadials.

In summary, our reconstructions reveal hitherto unrecognized bipolar ocean ventilation modes behind various types of rapid atmospheric CO₂ changes during the past 150,000 years. Since the Industrial Revolution, the ocean has absorbed a substantial amount of atmospheric CO₂ via high-latitude oceans in both hemispheres, helping to slow down global warming⁵². Given ongoing and possible future AMOC weakening⁵³, it is imperative to fully understand how processes in polar oceans, including both the Southern Ocean and the North Atlantic, affect atmospheric CO₂ and associated ocean acidification. Models used to quantify ocean carbon storage must simulate these processes and their interplay accurately. Our [CO₃²⁻] 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.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-023-01297-x.

## References

1. Bereiter, B. et al. Revision of the EPICA Dome C CO₂ record from 800 to 600 kyr before present. *Geophys. Res. Lett.* **42**, 542–549 (2015).

2. Ahn, J. & Brook, E. J. Siple Dome ice reveals two modes of millennial CO₂ change during the last ice age. *Nat. Commun.* https://doi.org/10.1038/Ncomms4723 (2014).

3. Bauska, T. K., Marcott, S. A. & Brook, E. J. Abrupt changes in the global carbon cycle during the last glacial period. *Nat. Geosci.* **14**, 91–96 (2021).

4. Bereiter, B. et al. Mode change of millennial CO₂ variability during the last glacial cycle associated with a bipolar marine carbon seesaw. *Proc. Natl Acad. Sci. USA* **109**, 9755–9760 (2012).

5. Broecker, W. Glacial to interglacial changes in ocean chemistry. *Progr. Oceanogr.* **11**, 151–197 (1982).

6. Toggweiler, J. R. Variation of atmospheric CO₂ by ventilation of the ocean's deepest water. *Paleoceanogr. Paleoclimatol.* **14**, 571–588 (1999).

7. Fischer, H. et al. The role of Southern Ocean processes in orbital and millennial CO₂ variations—a synthesis. *Quat. Sci. Rev.* **29**, 193–205 (2010).

8. Anderson, R. F. et al. Wind-driven upwelling in the Southern Ocean and the deglacial rise in atmospheric CO₂. *Science* **323**, 1443–1448 (2009).

9. Sigman, D. M., Hain, M. P. & Haug, G. H. The polar ocean and glacial cycles in atmospheric CO₂ concentration. *Nature* **466**, 47–55 (2010).

10. Toggweiler, J. R., Russell, J. L. & Carson, S. R. Midlatitude westerlies, atmospheric CO₂, and climate change during the ice ages. *Paleoceanogr. Paleoclimatol.* **21**, PA2005 (2006).

11. Stocker, T. F. & Johnsen, S. J. A minimum thermodynamic model for the bipolar seesaw. *Paleoceanogr. Paleoclimatol.* **18**, 1087 (2003).

12. Ai, X. E. et al. Southern Ocean upwelling, Earth's obliquity, and glacial-interglacial atmospheric CO₂ change. *Science* **370**, 1348–1352 (2020).

13. Pedro, J. B. et al. Beyond the bipolar seesaw: toward a process understanding of interhemispheric coupling. *Quat. Sci. Rev.* **192**, 27–46 (2018).

14. Rhodes, R. H. et al. Enhanced tropical methane production in response to iceberg discharge in the North Atlantic. *Science* **348**, 1016–1019 (2015).

15. Jaccard, S. L., Galbraith, E. D., Martinez-Garcia, A. & Anderson, R. F. Covariation of deep Southern Ocean oxygenation and atmospheric CO₂ through the last ice age. *Nature* https://doi.org/10.1038/nature16514 (2016).

16. Gottschalk, J. et al. Biological and physical controls in the Southern Ocean on past millennial-scale atmospheric CO₂ changes. *Nat. Commun.* **7**, 11539 (2016).

17. Gottschalk, J. et al. Abrupt changes in the southern extent of North Atlantic Deep Water during Dansgaard-Oeschger events. *Nat. Geosci.* **8**, 950–968 (2015).

18. Burke, A. & Robinson, L. F. The Southern Ocean's role in carbon exchange during the last deglaciation. *Science* **335**, 557–561 (2012).

19. Rae, J. W. B. et al. CO₂ storage and release in the deep Southern Ocean on millennial to centennial timescales. *Nature* **562**, 569–573 (2018).

20. Yu, J. M. et al. Millennial and centennial CO₂ release from the Southern Ocean during the last deglaciation. *Nat. Geosci.* **15**, 293–299 (2022).

21. Yu, J. et al. Loss of carbon from the deep sea since the Last Glacial Maximum. *Science* **330**, 1084–1087 (2010).

22. Yu, J. et al. Sequestration of carbon in the deep Atlantic during the last glaciation. *Nat. Geosci.* **9**, 319–324 (2016).

23. Ito, T. & Follows, M. J. Preformed phosphate, soft tissue pump and atmospheric CO₂. *J. Mar. Res.* **63**, 813–CO839 (2005).

24. Yu, J. M. & Elderfield, H. Benthic foraminiferal B/Ca ratios reflect deep water carbonate saturation state. *Earth Planet. Sci. Lett.* **258**, 73–86 (2007).

25. NGRIP members High-resolution record of Northern Hemisphere climate extending into the last interglacial period. *Nature* **431**, 147–151 (2004).

26. Marino, G. et al. Bipolar seesaw control on last interglacial sea level. *Nature* **522**, 197 (2015).

27. Barbante, C. et al. One-to-one coupling of glacial climate variability in Greenland and Antarctica. *Nature* **444**, 195–198 (2006).

28. Ahn, J. & Brook, E. J. Atmospheric CO₂ and climate on millennial time scales during the last glacial period. *Science* **322**, 83–85 (2008).

29. Bohm, E. et al. Strong and deep Atlantic meridional overturning circulation during the last glacial cycle. *Nature* https://doi.org/10.1038/Nature14059 (2015).

30. Broecker, W. S. & Peng, T. H. The role of CaCO₃ compensation in the glacial to interglacial atmospheric CO₂ change. *Glob. Biogeochem. Cycle* **1**, 15–29 (1987).

31. Gu, S. et al. Remineralization dominating the δ¹³C decrease in the mid-depth Atlantic during the last deglaciation. *Earth Planet. Sci. Lett.* **571**, 117106 (2021).

32. McManus, J. F., Francois, R., Gherardi, J. M., Keigwin, L. D. & Brown-Leger, S. Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes. *Nature* **428**, 834–837 (2004).

33. Henry, L. G. et al. North Atlantic ocean circulation and abrupt climate change during the last glaciation. *Science* **353**, 470–474 (2016).

34. Menviel, L., England, M. H., Meissner, K. J., Mouchet, A. & Yu, J. Atlantic-Pacific seesaw and its role in outgassing CO₂ during Heinrich events. *Paleoceanogr. Paleoclimatol.* **29**, 58–70 (2014).
ocean and land biosphere reservoirs, but how the terrestrial biosphere carbon changed during stadials remains poorly constrained²,⁴¹. Here we provide critical data evidence (Fig. 2) to demonstrate that respired carbon accumulation, driven by reduced North Atlantic ventilation³¹,⁴⁸,⁵⁰, at times cancelled or surpassed the effect of Southern Ocean outgassing, explaining the observations of little net change or even declines in atmospheric CO₂ during some MIS 3 stadials.

In summary, our reconstructions reveal hitherto unrecognized bipolar ocean ventilation modes behind various types of rapid atmospheric CO₂ changes during the past 150,000 years. Since the Industrial Revolution, the ocean has absorbed a substantial amount of atmospheric CO₂ via high-latitude oceans in both hemispheres, helping to slow down global warming⁵². Given ongoing and possible future AMOC weakening⁵³, it is imperative to fully understand how processes in polar oceans, including both the Southern Ocean and the North Atlantic, affect atmospheric CO₂ and associated ocean acidification. Models used to quantify ocean carbon storage must simulate these processes and their interplay accurately. Our [CO₃²⁻] 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.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-023-01297-x.

## References

1. Bereiter, B. et al. Revision of the EPICA Dome C CO₂ record from 800 to 600 kyr before present. *Geophys. Res. Lett.* **42**, 542–549 (2015).

2. Ahn, J. & Brook, E. J. Siple Dome ice reveals two modes of millennial CO₂ change during the last ice age. *Nat. Commun.* https://doi.org/10.1038/Ncomms4723 (2014).

3. Bauska, T. K., Marcott, S. A. & Brook, E. J. Abrupt changes in the global carbon cycle during the last glacial period. *Nat. Geosci.* **14**, 91–96 (2021).

4. Bereiter, B. et al. Mode change of millennial CO₂ variability during the last glacial cycle associated with a bipolar marine carbon seesaw. *Proc. Natl Acad. Sci. USA* **109**, 9755–9760 (2012).

5. Broecker, W. Glacial to interglacial changes in ocean chemistry. *Progr. Oceanogr.* **11**, 151–197 (1982).

6. Toggweiler, J. R. Variation of atmospheric CO₂ by ventilation of the ocean's deepest water. *Paleoceanogr. Paleoclimatol.* **14**, 571–588 (1999).

7. Fischer, H. et al. The role of Southern Ocean processes in orbital and millennial CO₂ variations—a synthesis. *Quat. Sci. Rev.* **29**, 193–205 (2010).

8. Anderson, R. F. et al. Wind-driven upwelling in the Southern Ocean and the deglacial rise in atmospheric CO₂. *Science* **323**, 1443–1448 (2009).

9. Sigman, D. M., Hain, M. P. & Haug, G. H. The polar ocean and glacial cycles in atmospheric CO₂ concentration. *Nature* **466**, 47–55 (2010).

10. Toggweiler, J. R., Russell, J. L. & Carson, S. R. Midlatitude westerlies, atmospheric CO₂, and climate change during the ice ages. *Paleoceanogr. Paleoclimatol.* **21**, PA2005 (2006).

11. Stocker, T. F. & Johnsen, S. J. A minimum thermodynamic model for the bipolar seesaw. *Paleoceanogr. Paleoclimatol.* **18**, 1087 (2003).

12. Ai, X. E. et al. Southern Ocean upwelling, Earth's obliquity, and glacial-interglacial atmospheric CO₂ change. *Science* **370**, 1348–1352 (2020).

13. Pedro, J. B. et al. Beyond the bipolar seesaw: toward a process understanding of interhemispheric coupling. *Quat. Sci. Rev.* **192**, 27–46 (2018).

14. Rhodes, R. H. et al. Enhanced tropical methane production in response to iceberg discharge in the North Atlantic. *Science* **348**, 1016–1019 (2015).

15. Jaccard, S. L., Galbraith, E. D., Martinez-Garcia, A. & Anderson, R. F. Covariation of deep Southern Ocean oxygenation and atmospheric CO₂ through the last ice age. *Nature* https://doi.org/10.1038/nature16514 (2016).

16. Gottschalk, J. et al. Biological and physical controls in the Southern Ocean on past millennial-scale atmospheric CO₂ changes. *Nat. Commun.* **7**, 11539 (2016).

17. Gottschalk, J. et al. Abrupt changes in the southern extent of North Atlantic Deep Water during Dansgaard-Oeschger events. *Nat. Geosci.* **8**, 950–968 (2015).

18. Burke, A. & Robinson, L. F. The Southern Ocean's role in carbon exchange during the last deglaciation. *Science* **335**, 557–561 (2012).

19. Rae, J. W. B. et al. CO₂ storage and release in the deep Southern Ocean on millennial to centennial timescales. *Nature* **562**, 569–573 (2018).

20. Yu, J. M. et al. Millennial and centennial CO₂ release from the Southern Ocean during the last deglaciation. *Nat. Geosci.* **15**, 293–299 (2022).

21. Yu, J. et al. Loss of carbon from the deep sea since the Last Glacial Maximum. *Science* **330**, 1084–1087 (2010).

22. Yu, J. et al. Sequestration of carbon in the deep Atlantic during the last glaciation. *Nat. Geosci.* **9**, 319–324 (2016).

23. Ito, T. & Follows, M. J. Preformed phosphate, soft tissue pump and atmospheric CO₂. *J. Mar. Res.* **63**, 813–CO839 (2005).

24. Yu, J. M. & Elderfield, H. Benthic foraminiferal B/Ca ratios reflect deep water carbonate saturation state. *Earth Planet. Sci. Lett.* **258**, 73–86 (2007).

25. NGRIP members High-resolution record of Northern Hemisphere climate extending into the last interglacial period. *Nature* **431**, 147–151 (2004).

26. Marino, G. et al. Bipolar seesaw control on last interglacial sea level. *Nature* **522**, 197 (2015).

27. Barbante, C. et al. One-to-one coupling of glacial climate variability in Greenland and Antarctica. *Nature* **444**, 195–198 (2006).

28. Ahn, J. & Brook, E. J. Atmospheric CO₂ and climate on millennial time scales during the last glacial period. *Science* **322**, 83–85 (2008).

29. Bohm, E. et al. Strong and deep Atlantic meridional overturning circulation during the last glacial cycle. *Nature* https://doi.org/10.1038/Nature14059 (2015).

30. Broecker, W. S. & Peng, T. H. The role of CaCO₃ compensation in the glacial to interglacial atmospheric CO₂ change. *Glob. Biogeochem. Cycle* **1**, 15–29 (1987).

31. Gu, S. et al. Remineralization dominating the δ¹³C decrease in the mid-depth Atlantic during the last deglaciation. *Earth Planet. Sci. Lett.* **571**, 117106 (2021).

32. McManus, J. F., Francois, R., Gherardi, J. M., Keigwin, L. D. & Brown-Leger, S. Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes. *Nature* **428**, 834–837 (2004).

33. Henry, L. G. et al. North Atlantic ocean circulation and abrupt climate change during the last glaciation. *Science* **353**, 470–474 (2016).

34. Menviel, L., England, M. H., Meissner, K. J., Mouchet, A. & Yu, J. Atlantic-Pacific seesaw and its role in outgassing CO₂ during Heinrich events. *Paleoceanogr. Paleoclimatol.* **29**, 58–70 (2014).
35. Broecker, W. S. Paleocean circulation during the last deglaciation: a bipolar seesaw? *Paleoceanogr. Paleoclimatol.* **13**, 119–121 (1998).

36. Menviel, L. et al. Southern Hemisphere westerlies as a driver of the early deglacial atmospheric CO₂ rise. *Nat. Commun.* **9**, 2503 (2018).

37. Studer, A. S. et al. Antarctic Zone nutrient conditions during the last two glacial cycles. *Paleoceanogr. Paleoclimatol.* **30**, 845–862 (2015).

38. Martinez-Garcia, A. et al. Iron fertilization of the subantarctic ocean during the last ice age. *Science* **343**, 1347–1350 (2014).

39. Skinner, L. C., Waelbroeck, C., Scrivner, A. E. & Fallon, S. J. Radiocarbon evidence for alternating northern and southern sources of ventilation of the deep Atlantic carbon pool during the last deglaciation. *Proc. Natl Acad. Sci. USA* **111**, 5480–5484 (2014).

40. Grant, K. M. et al. Rapid coupling between ice volume and polar temperature over the past 150,000 years. *Nature* **491**, 744–747 (2012).

41. Gottschalk, J. et al. Mechanisms of millennial-scale atmospheric CO₂ change in numerical model simulations. *Quat. Sci. Rev.* **220**, 30–74 (2019).

42. Hodell, D. A. et al. Anatomy of Heinrich Layer 1 and its role in the last deglaciation. *Paleoceanogr. Paleoclimatol.* **32**, 284–303 (2017).

43. Gherardi, J. M. et al. Glacial-interglacial circulation changes inferred from $^{231}$Pa/$^{230}$Th sedimentary record in the North Atlantic region. *Paleoceanogr. Paleoclimatol*. https://doi.org/10.1029/2008pa001696 (2009).

44. Tzedakis, P. C. et al. Enhanced climate instability in the North Atlantic and southern Europe during the Last Interglacial. *Nat. Commun.* **9**, 4235 (2018).

45. Piotrowski, A., Goldstein, S. J., Hemming, S. R. & Fairbanks, R. G. Temporal relationships of carbon cycling and ocean circulation at glacial boundaries. *Science* **307**, 1933–1938 (2005).

46. Yu, J. et al. Last glacial atmospheric CO₂ decline due to widespread Pacific deep-water expansion. *Nat. Geosci.* https://doi.org/10.1038/s41561-020-0610-5 (2020).

47. Yu, J. et al. Deep South Atlantic carbonate chemistry and increased interocean deep water exchange during last deglaciation. *Quat. Sci. Rev.* **15**, 80–89 (2014).

48. Menviel, L., Timmermann, A., Friedrich, T. & England, M. H. Hindcasting the continuum of Dansgaard-Oeschger variability: mechanisms, patterns and timing. *Clim. Past* **10**, 63–77 (2014).

49. Shakun, J. D., Lea, D. W., Lisiecki, L. E. & Raymo, M. E. An 800 kyr record of global surface ocean $\delta^{18}$O and implications for ice volume–temperature coupling. *Earth Planet. Sci. Lett.* **426**, 58–68 (2015).

50. Schmittner, A. & Lund, D. C. Early deglacial Atlantic overturning decline and its role in atmospheric CO₂ rise inferred from carbon isotopes ($\delta^{13}$C). *Clim. Past* **11**, 135–152 (2015).

51. Anderson, R. F. et al. Deep-sea oxygen depletion and ocean carbon sequestration during the last ice age. *Glob. Biogeochem. Cycles* https://doi.org/10.1029/2018GB006049 (2019).

52. Gruber, N. et al. The oceanic sink for anthropogenic CO₂ from 1994 to 2007. *Science* **363**, 1193–1199 (2019).

53. Thornalley, D. J. R. et al. Anomalously weak Labrador Sea convection and Atlantic overturning during the past 150 years. *Nature* **556**, 227–230 (2018).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023
35. Broecker, W. S. Paleocean circulation during the last deglaciation: a bipolar seesaw? *Paleoceanogr. Paleoclimatol.* **13**, 119–121 (1998).

36. Menviel, L. et al. Southern Hemisphere westerlies as a driver of the early deglacial atmospheric CO₂ rise. *Nat. Commun.* **9**, 2503 (2018).

37. Studer, A. S. et al. Antarctic Zone nutrient conditions during the last two glacial cycles. *Paleoceanogr. Paleoclimatol.* **30**, 845–862 (2015).

38. Martinez-Garcia, A. et al. Iron fertilization of the subantarctic ocean during the last ice age. *Science* **343**, 1347–1350 (2014).

39. Skinner, L. C., Waelbroeck, C., Scrivner, A. E. & Fallon, S. J. Radiocarbon evidence for alternating northern and southern sources of ventilation of the deep Atlantic carbon pool during the last deglaciation. *Proc. Natl Acad. Sci. USA* **111**, 5480–5484 (2014).

40. Grant, K. M. et al. Rapid coupling between ice volume and polar temperature over the past 150,000 years. *Nature* **491**, 744–747 (2012).

41. Gottschalk, J. et al. Mechanisms of millennial-scale atmospheric CO₂ change in numerical model simulations. *Quat. Sci. Rev.* **220**, 30–74 (2019).

42. Hodell, D. A. et al. Anatomy of Heinrich Layer 1 and its role in the last deglaciation. *Paleoceanogr. Paleoclimatol.* **32**, 284–303 (2017).

43. Gherardi, J. M. et al. Glacial-interglacial circulation changes inferred from $^{231}$Pa/$^{230}$Th sedimentary record in the North Atlantic region. *Paleoceanogr. Paleoclimatol*. https://doi.org/10.1029/2008pa001696 (2009).

44. Tzedakis, P. C. et al. Enhanced climate instability in the North Atlantic and southern Europe during the Last Interglacial. *Nat. Commun.* **9**, 4235 (2018).

45. Piotrowski, A., Goldstein, S. J., Hemming, S. R. & Fairbanks, R. G. Temporal relationships of carbon cycling and ocean circulation at glacial boundaries. *Science* **307**, 1933–1938 (2005).

46. Yu, J. et al. Last glacial atmospheric CO₂ decline due to widespread Pacific deep-water expansion. *Nat. Geosci.* https://doi.org/10.1038/s41561-020-0610-5 (2020).

47. Yu, J. et al. Deep South Atlantic carbonate chemistry and increased interocean deep water exchange during last deglaciation. *Quat. Sci. Rev.* **15**, 80–89 (2014).

48. Menviel, L., Timmermann, A., Friedrich, T. & England, M. H. Hindcasting the continuum of Dansgaard-Oeschger variability: mechanisms, patterns and timing. *Clim. Past* **10**, 63–77 (2014).

49. Shakun, J. D., Lea, D. W., Lisiecki, L. E. & Raymo, M. E. An 800 kyr record of global surface ocean $\delta^{18}$O and implications for ice volume–temperature coupling. *Earth Planet. Sci. Lett.* **426**, 58–68 (2015).

50. Schmittner, A. & Lund, D. C. Early deglacial Atlantic overturning decline and its role in atmospheric CO₂ rise inferred from carbon isotopes ($\delta^{13}$C). *Clim. Past* **11**, 135–152 (2015).

51. Anderson, R. F. et al. Deep-sea oxygen depletion and ocean carbon sequestration during the last ice age. *Glob. Biogeochem. Cycles* https://doi.org/10.1029/2018GB006049 (2019).

52. Gruber, N. et al. The oceanic sink for anthropogenic CO₂ from 1994 to 2007. *Science* **363**, 1193–1199 (2019).

53. Thornalley, D. J. R. et al. Anomalously weak Labrador Sea convection and Atlantic overturning during the past 150 years. *Nature* **556**, 227–230 (2018).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023
# Methods

## 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₂ 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₂ 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<sup>15,16,18,51,54</sup>.

## 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<sup>55</sup> 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₂ variations<sup>1-4</sup> 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<sup>56,57</sup>. Fourth, the latitude of MD95-2039 makes it useful to constrain the geographic extent of the relative influence duringstadials of ventilation from the Southern Ocean, which tends to increase deep-water [CO₃²⁻] due to outgassing, and via the North Atlantic, which tends to decrease [CO₃²⁻] due to enhanced biogenic matter respiration (Supplementary Fig. 1).

## Samples and analyses

To ensure enough shells for analyses, about 30 cm³ 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 δ¹⁸O measurements. The average analytical error is -0.08‰ in δ¹⁸O. We compared data from overlapping depths and observed negligible analytical offset between

new and published<sup>58</sup> 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₃) were used for each measurement at Research School of Earth Sciences at the Australian National University (Supplementary Table 1).

To strengthen the robustness of our age model, we have generated new abundances of IRD (n = 152; >150 µm size fraction) and *Neogloboquadrina pachyderma* (*N. pachyderma*/*all* planktonics; n = 66; >150 µm size fraction) (Extended Data Figs. 3 and 5), following approaches described in refs. 58,59.

For B/Ca analyses, the epifaunal benthic foraminiferal species *C. wuellerstorfii* 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<sup>60</sup>. 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<sup>62,63</sup> (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* δ¹⁸O in MD95-2039 to the NGRIP ice-core δ¹⁸O record<sup>25</sup> on the AICC2012 age model<sup>64</sup>. Due to the lack of significant structure in NGRIP δ¹⁸O during HS1, we added two tie points to align MD95-2039 *G. bulloides* δ¹⁸O to the Hulu speleothem δ¹⁸O record (on the U–Th age)<sup>65</sup> 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 δ¹⁸O using a 1.0063 scaling factor (WDC2014 age = 1.0063 × AICC2012 age)<sup>66</sup> to facilitate comparisons of our deep-water [CO₃²⁻] with the latest high-resolution atmospheric CO₂ record from the WDC ice core³ (Fig. 2). Because WDC record covers only the past -65 ka, we used the AICC2012 age model for NGRIP δ¹⁸O for -65–120 ka. Between -120 ka and 140 ka, we have tuned MD95-2039 *G. bulloides* δ¹⁸O to a nearby core ODP Site 976 *G. bulloides* δ¹⁸O, which has been mapped onto speleothem U–Th age (equivalent to the AICC2012 age model)<sup>26</sup>. During 140–150 ka, MD95-2039 *G. bulloides* δ¹⁸O is tuned to the synthetic NGRIP δ¹⁸O record<sup>67</sup>. Other ice-core records used in this study, including EDML δ¹⁸O<sup>27</sup> (WDC δ¹⁸O record is too short and thus not used) and compiled atmospheric CO₂ records<sup>1-4</sup>, 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<sup>58,59</sup> 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.

## Deep-water [CO₃²⁻] reconstructions

Deep-water [CO₃²⁻] values at MD95-2039 are reconstructed using *C. wuellerstorfii* B/Ca<sup>22,24</sup> from [CO₃²⁻]downcore = [CO₃²⁻]PI + ΔB/Ca<sub>downcore-coretop</sub>/Sen, where [CO₃²⁻]PI is the preindustrial (PI) deep-water [CO₃²⁻] value (105 µmol kg⁻¹) estimated from the GLODAP dataset<sup>68</sup>, ΔB/Ca<sub>downcore-coretop</sub> represents the deviation of B/Ca of downcore samples from the core-top value, and the term Sen is the B/Ca-[CO₃²⁻] sensitivity of *C. wuellerstorfii* (1.14 µmol mol⁻¹ per µmol kg⁻¹) (ref. 24).
# Methods

## 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₂ 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₂ 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<sup>15,16,18,51,54</sup>.

## 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<sup>55</sup> 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₂ variations<sup>1-4</sup> 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<sup>56,57</sup>. Fourth, the latitude of MD95-2039 makes it useful to constrain the geographic extent of the relative influence duringstadials of ventilation from the Southern Ocean, which tends to increase deep-water [CO₃²⁻] due to outgassing, and via the North Atlantic, which tends to decrease [CO₃²⁻] due to enhanced biogenic matter respiration (Supplementary Fig. 1).

## Samples and analyses

To ensure enough shells for analyses, about 30 cm³ 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 δ¹⁸O measurements. The average analytical error is -0.08‰ in δ¹⁸O. We compared data from overlapping depths and observed negligible analytical offset between

new and published<sup>58</sup> 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₃) were used for each measurement at Research School of Earth Sciences at the Australian National University (Supplementary Table 1).

To strengthen the robustness of our age model, we have generated new abundances of IRD (n = 152; >150 µm size fraction) and *Neogloboquadrina pachyderma* (*N. pachyderma*/*all* planktonics; n = 66; >150 µm size fraction) (Extended Data Figs. 3 and 5), following approaches described in refs. 58,59.

For B/Ca analyses, the epifaunal benthic foraminiferal species *C. wuellerstorfii* 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<sup>60</sup>. 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<sup>62,63</sup> (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* δ¹⁸O in MD95-2039 to the NGRIP ice-core δ¹⁸O record<sup>25</sup> on the AICC2012 age model<sup>64</sup>. Due to the lack of significant structure in NGRIP δ¹⁸O during HS1, we added two tie points to align MD95-2039 *G. bulloides* δ¹⁸O to the Hulu speleothem δ¹⁸O record (on the U–Th age)<sup>65</sup> 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 δ¹⁸O using a 1.0063 scaling factor (WDC2014 age = 1.0063 × AICC2012 age)<sup>66</sup> to facilitate comparisons of our deep-water [CO₃²⁻] with the latest high-resolution atmospheric CO₂ record from the WDC ice core³ (Fig. 2). Because WDC record covers only the past -65 ka, we used the AICC2012 age model for NGRIP δ¹⁸O for -65–120 ka. Between -120 ka and 140 ka, we have tuned MD95-2039 *G. bulloides* δ¹⁸O to a nearby core ODP Site 976 *G. bulloides* δ¹⁸O, which has been mapped onto speleothem U–Th age (equivalent to the AICC2012 age model)<sup>26</sup>. During 140–150 ka, MD95-2039 *G. bulloides* δ¹⁸O is tuned to the synthetic NGRIP δ¹⁸O record<sup>67</sup>. Other ice-core records used in this study, including EDML δ¹⁸O<sup>27</sup> (WDC δ¹⁸O record is too short and thus not used) and compiled atmospheric CO₂ records<sup>1–4</sup>, 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<sup>58,59</sup> 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.

## Deep-water [CO₃²⁻] reconstructions

Deep-water [CO₃²⁻] values at MD95-2039 are reconstructed using *C. wuellerstorfii* B/Ca<sup>22,24</sup> from [CO₃²⁻]downcore = [CO₃²⁻]PI + ΔB/Ca<sub>downcore-coretop</sub>/Sen, where [CO₃²⁻]PI is the preindustrial (PI) deep-water [CO₃²⁻] value (105 µmol kg⁻¹) estimated from the GLODAP dataset<sup>68</sup>, ΔB/Ca<sub>downcore-coretop</sub> represents the deviation of B/Ca of downcore samples from the core-top value, and the term Sen is the B/Ca-[CO₃²⁻] sensitivity of *C. wuellerstorfii* (1.14 µmol mol⁻¹ per µmol kg⁻¹) (ref. 24).
The reconstruction uncertainty ($2\sigma$) is 10 µmol kg⁻¹ in $[CO_3^{2-}]$, based on global core-top calibration samples²⁴,⁶⁹.

## Statistical analyses

Uncertainties associated with time series, including our new deep-water $[CO_3^{2-}]$, *G. bulloides* $\delta^{18}$O and ice-core parameters, were evaluated using a Monte Carlo-style approach⁴⁰,⁷⁰, 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⁷¹, while the individual $[CO_3^{2-}]$ error is estimated to be 10 µmol kg⁻¹ ($2\sigma$). 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.5th percentile) probability intervals for the data distributions⁴⁰,⁷⁰. The same data processing approach was used to analyse $\delta^{18}$O and atmospheric CO₂ records using their respective uncertainties.

We constructed stadial evolutions of NGRIP $\delta^{18}$O²⁵, atmospheric $CO_2^{1,3}$, EDML $\delta^{18}$O²⁷ and MD95-2039 *G. bulloides* $\delta^{18}$O and $[CO_3^{2-}]$ 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 $\delta^{18}$O, following ref. 3. Third, anomalies ($\Delta$) of signals and associated $2\sigma$ 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\sigma$ of all individual stadial signals and their associated uncertainties in each mode using the Monte Carlo-based approach⁴⁰,⁷⁰.

## 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⁴⁴,⁵⁵. During Mode II stadials, the robustness of deep-water $[CO_3^{2-}]$ 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⁻¹) 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 $\delta^{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* $\delta^{18}$O from multiple stadials on the basis of thorough statistical analyses to account for both chronology and reconstruction uncertainties (Supplementary Fig. 4).

## Data availability

All data are provided in Supplementary Data 1.

## Code availability

We used published methods for statistical analyses; no new code was used.

## References

54. Chen, P. J., Yu, J. M. & Jin, Z. D. An evaluation of benthic foraminiferal U/Ca and U/Mn proxies for deep ocean carbonate chemistry and redox conditions. *Geochem. Geophys. Geosyst.* **18**, 617–630 (2017).

55. Shackleton, N. J., Hall, M. A. & Vincent, E. Phase relationships between millennial-scale events 64,000–24,000 years ago. *Paleoceanogr. Paleoclimatol.* **15**, 565–569 (2000).

56. Sarnthein, M. et al. in *The Northern North Atlantic: A Changing Environment* (eds Schafer, P. et al.) 365–410 (Springer, 2000).

57. Sarnthein, M. et al. Changes in east Atlantic deepwater circulation over the last 30,000 years: eight time slice reconstructions. *Paleoceanogr. Paleoclimatol.* **9**, 209–268 (1994).

58. Schonfeld, J., Zahn, R. & de Abreu, L. Surface and deep water response to rapid climate changes at the Western Iberian Margin. *Glob. Planet. Change* **36**, 237–264 (2003).

59. Roucoux, K. H., Shackleton, N. J., de Abreu, L., Schonfeld, J. & Tzedakis, P. C. Combined marine proxy and pollen analyses reveal rapid Iberian vegetation response to North Atlantic. *Quat. Res.* **56**, 128–132 (2001).

60. Barker, S., Greaves, M. & Elderfield, H. A study of cleaning procedures used for foraminiferal Mg/Ca paleothermometry. *Geochem. Geophys. Geosyst.* https://doi.org/10.1029/2003GC000559 (2003).

61. Yu, J. M., Day, J., Greaves, M. & Elderfield, H. Determination of multiple element/calcium ratios in foraminiferal calcite by quadrupole ICP-MS. *Geochem. Geophys. Geosyst.* **6**, Q08P01 (2005).

62. Heaton, T. J. et al. Marine20—the marine radiocarbon age calibration curve (0–55,000 cal BP). *Radiocarbon* **14**, 779–820 (2020).

63. Stuiver, M. & Reimer, P. J. CALIB rev. 8. *Radiocarbon* **35**, 215–230 (1993).

64. Bazin, L. et al. An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120–800 ka. *Clim. Past* **9**, 1715–1731 (2013).

65. Cheng, H. et al. Ice age terminations. *Science* **326**, 248–252 (2009).

66. Buizert, C. et al. The WAIS Divide deep ice core WD2014 chronology—part 1: methane synchronization (68–31 ka BP) and the gas age–ice age difference. *Clim. Past* **11**, 153–173 (2015).

67. Barker, S. et al. 800,000 years of abrupt climate variability. *Science* **334**, 347–351 (2011).

68. Key, R. M. et al. A global ocean carbon climatology: results from Global Data Analysis Project (GLODAP). *Glob. Biogeochem. Cycles* https://doi.org/10.1029/2004GB002247 (2004).

69. Yu, J. et al. Responses of the deep ocean carbonate system to carbon reorganization during the Last Glacial–interglacial cycle. *Quat. Sci. Rev.* **76**, 39–52 (2013).

70. Rohling, E. J. et al. Sea-level and deep-sea-temperature variability over the past 5.3 million years. *Nature* **508**, 477–482 (2014).

71. Ramsey, C. B. Deposition models for chronological records. *Quat. Sci. Rev.* **27**, 42–60 (2008).
The reconstruction uncertainty ($2\sigma$) is 10 µmol kg⁻¹ in $[CO_3^{2-}]$, based on global core-top calibration samples²⁴,⁶⁹.

## Statistical analyses

Uncertainties associated with time series, including our new deep-water $[CO_3^{2-}]$, *G. bulloides* $\delta^{18}$O and ice-core parameters, were evaluated using a Monte Carlo-style approach⁴⁰,⁷⁰, 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⁷¹, while the individual $[CO_3^{2-}]$ error is estimated to be 10 µmol kg⁻¹ ($2\sigma$). 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.5th percentile) probability intervals for the data distributions⁴⁰,⁷⁰. The same data processing approach was used to analyse $\delta^{18}$O and atmospheric CO₂ records using their respective uncertainties.

We constructed stadial evolutions of NGRIP $\delta^{18}$O²⁵, atmospheric $CO_2^{1,3}$, EDML $\delta^{18}$O²⁷ and MD95-2039 G. *bulloides* $\delta^{18}$O and $[CO_3^{2-}]$ 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 $\delta^{18}$O, following ref. 3. Third, anomalies ($\Delta$) of signals and associated $2\sigma$ 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\sigma$ of all individual stadial signals and their associated uncertainties in each mode using the Monte Carlo-based approach⁴⁰,⁷⁰.

## 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⁴⁴,⁵⁵. During Mode II stadials, the robustness of deep-water $[CO_3^{2-}]$ 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⁻¹) 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 $\delta^{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* $\delta^{18}$O from multiple stadials on the basis of thorough statistical analyses to account for both chronology and reconstruction uncertainties (Supplementary Fig. 4).

## Data availability

All data are provided in Supplementary Data 1.

## Code availability

We used published methods for statistical analyses; no new code was used.

## References

54. Chen, P. J., Yu, J. M. & Jin, Z. D. An evaluation of benthic foraminiferal U/Ca and U/Mn proxies for deep ocean carbonate chemistry and redox conditions. *Geochem. Geophys. Geosyst.* **18**, 617–630 (2017).

55. Shackleton, N. J., Hall, M. A. & Vincent, E. Phase relationships between millennial-scale events 64,000–24,000 years ago. *Paleoceanogr. Paleoclimatol.* **15**, 565–569 (2000).

56. Sarnthein, M. et al. in *The Northern North Atlantic: A Changing Environment* (eds Schafer, P. et al.) 365–410 (Springer, 2000).

57. Sarnthein, M. et al. Changes in east Atlantic deepwater circulation over the last 30,000 years: eight time slice reconstructions. *Paleoceanogr. Paleoclimatol.* **9**, 209–268 (1994).

58. Schonfeld, J., Zahn, R. & de Abreu, L. Surface and deep water response to rapid climate changes at the Western Iberian Margin. *Glob. Planet. Change* **36**, 237–264 (2003).

59. Roucoux, K. H., Shackleton, N. J., de Abreu, L., Schonfeld, J. & Tzedakis, P. C. Combined marine proxy and pollen analyses reveal rapid Iberian vegetation response to North Atlantic. *Quat. Res.* **56**, 128–132 (2001).

60. Barker, S., Greaves, M. & Elderfield, H. A study of cleaning procedures used for foraminiferal Mg/Ca paleothermometry. *Geochem. Geophys. Geosyst.* https://doi.org/10.1029/2003GC000559 (2003).

61. Yu, J. M., Day, J., Greaves, M. & Elderfield, H. Determination of multiple element/calcium ratios in foraminiferal calcite by quadrupole ICP-MS. *Geochem. Geophys. Geosyst.* **6**, Q08P01 (2005).

62. Heaton, T. J. et al. Marine20—the marine radiocarbon age calibration curve (0–55,000 cal BP). *Radiocarbon* **14**, 779–820 (2020).

63. Stuiver, M. & Reimer, P. J. CALIB rev. 8. *Radiocarbon* **35**, 215–230 (1993).

64. Bazin, L. et al. An optimized multi-proxy, multi-site Antarctic ice and gas orbital chronology (AICC2012): 120–800 ka. *Clim. Past* **9**, 1715–1731 (2013).

65. Cheng, H. et al. Ice age terminations. *Science* **326**, 248–252 (2009).

66. Buizert, C. et al. The WAIS Divide deep ice core WD2014 chronology—part 1: methane synchronization (68–31 ka BP) and the gas age–ice age difference. *Clim. Past* **11**, 153–173 (2015).

67. Barker, S. et al. 800,000 years of abrupt climate variability. *Science* **334**, 347–351 (2011).

68. Key, R. M. et al. A global ocean carbon climatology: results from Global Data Analysis Project (GLODAP). *Glob. Biogeochem. Cycles* https://doi.org/10.1029/2004GB002247 (2004).

69. Yu, J. et al. Responses of the deep ocean carbonate system to carbon reorganization during the Last Glacial–interglacial cycle. *Quat. Sci. Rev.* **76**, 39–52 (2013).

70. Rohling, E. J. et al. Sea-level and deep-sea-temperature variability over the past 5.3 million years. *Nature* **508**, 477–482 (2014).

71. Ramsey, C. B. Deposition models for chronological records. *Quat. Sci. Rev.* **27**, 42–60 (2008).
72. Buizert, C. et al. Precise interpolar phasing of abrupt climate change during the last ice age. *Nature* **520**, 661–665 (2015).

73. Gottschalk, J. et al. Southern Ocean link between changes in atmospheric CO₂ levels and Northern-Hemisphere climate anomalies during the last two glacial periods. *Quat. Sci. Rev.* **230**, 106067 (2020).

74. Shin, J. et al. Millennial-scale atmospheric CO₂ variations during the Marine Isotope Stage 6 period (190–135 ka). *Clim. Past* **16**, 2203–2219 (2020).

75. Schlitzer, R. Ocean Data View v.3.2.2 https://odv.awi.de/ (ODV, 2006).

76. Zeebe, R. E. & Wolf-Gladrow, D. A. CO₂ in Seawater: Equilibrium, Kinetics, Isotopes Vol. 65 (Elsevier, 2001).

77. Thomson, J. et al. Implications for sedimentation changes on the Iberian margin over the last two glacial/interglacial transitions from ($^{230}$Th$_{excess}$)$_{0}$ systematics. *Earth Planet. Sci. Lett.* **165**, 255–270 (1999).

78. Röthlisberger, R., Crosta, X., Abram, N. J., Armand, L. & Wolff, E. W. Potential and limitations of marine and ice core sea ice proxies: an example from the Indian Ocean sector. *Quat. Sci. Rev.* **29**, 296–302 (2010).

79. Abram, N. J., Wolff, E. W. & Curran, M. A. J. A review of sea ice proxy information from polar ice cores. *Quat. Sci. Rev.* **79**, 168–183 (2013).

80. Gottschalk, J. et al. Past carbonate preservation events in the deep southeast Atlantic Ocean (Cape Basin) and their implications for Atlantic overturning dynamics and marine carbon cycling. *Paleoceanogr. Paleoclimatol.* **33**, 643–663 (2018).

81. Barker, S. & Diz, P. Timing of the descent into the last ice age determined by the bipolar seesaw. *Paleoceanogr. Paleoclimatol.* https://doi.org/10.1002/2014PA002623 (2014).

82. Sachs, J. & Anderson, R. F. Increased productivity in the subantarctic ocean during Heinrich events. *Nature* **434**, 1118–1121 (2005).

83. Ninnemann, U. S. & Charles, C. D. Changes in the mode of Southern Ocean circulation over the last glacial cycle revealed by foraminiferal stable isotopic variability. *Earth Planet. Sci. Lett.* **201**, 383–396 (2002).

84. Goosse, H. et al. Description of the Earth system model of intermediate complexity LOVECLIM version 1.2. *Geosci. Model Dev.* **3**, 603–633 (2010).

85. Berger, A. Long-term variations of daily insolation and Quaternary climatic changes. *J. Atmos. Sci.* **35**, 2362–2367 (1978).

86. Abe-Ouchi, A., Segawa, T. & Saito, F. Climatic conditions for modelling the Northern Hemisphere ice sheets throughout the ice age cycle. *Clim. Past* **3**, 423–438 (2007).

## Acknowledgements

Core MD95-2039 was collected during the international coring campaign IMAGES Cruise (MD 101) of the French vessel *Marion Dufresne* led by L. Labeyrie, F. Bassinot and Y. Lancelot. This project would not be possible without long-time and dedicated laboratory assistance to process numerous samples by P. Wang, Y. Dai, J. Zhang, L. Wang, P. Chen and L. Ma. This study is supported by NSFC (42076056), Laoshan Laboratory Science and Technology Innovation Project (LSKJ202203300), ARC Discovery Projects (DP190100894) and Future Fellowship (FT140100993) to J.Y. We are indebted to S. Barker for his permission to use stable isotope data collected at Cardiff University by D.J.R.T. as part of NERC grant NE/F002734/1 (to Barker). We thank D. Sigman and J. Gottschalk for discussions.

## Author contributions

J.Y. designed the project, interpreted data and wrote the paper. J.Y. and D.J.R.T. collected core samples with help from N.T., J.T. and X.X. Z.D.J. managed shell picking and IRD counting with help from F.Z. J.Y. supervised B/Ca data collection and analysed data with substantial input from X.J. D.J.R.T. did planktic faunal counts and contributed to stable isotope and B/Ca measurements and data analyses. L.W. provided expertise in physical oceanography. Y.C. and L.T. assisted with stable isotope analyses. L.M. assisted with model outputs and ED Figs. 9 and 10. R.F.A., L.W., E.J.R., J.F.M., D.J.R.T., L.M. and F.Z. commented on the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at https://doi.org/10.1038/s41561-023-01297-x.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41561-023-01297-x.

**Correspondence and requests for materials** should be addressed to J. Yu.

**Peer review information** *Nature Geoscience* thanks Fortunat Joos, Antje Voelker and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: James Super, in collaboration with the *Nature Geoscience* team.

**Reprints and permissions information** is available at www.nature.com/reprints.
72. Buizert, C. et al. Precise interpolar phasing of abrupt climate change during the last ice age. *Nature* **520**, 661–665 (2015).

73. Gottschalk, J. et al. Southern Ocean link between changes in atmospheric CO₂ levels and Northern-Hemisphere climate anomalies during the last two glacial periods. *Quat. Sci. Rev.* **230**, 106067 (2020).

74. Shin, J. et al. Millennial-scale atmospheric CO₂ variations during the Marine Isotope Stage 6 period (190–135 ka). *Clim. Past* **16**, 2203–2219 (2020).

75. Schlitzer, R. Ocean Data View v.3.2.2 https://odv.awi.de/ (ODV, 2006).

76. Zeebe, R. E. & Wolf-Gladrow, D. A. CO₂ in Seawater: Equilibrium, Kinetics, Isotopes Vol. 65 (Elsevier, 2001).

77. Thomson, J. et al. Implications for sedimentation changes on the Iberian margin over the last two glacial/interglacial transitions from ($^{230}$Th$_{excess}$)$_{0}$ systematics. *Earth Planet. Sci. Lett.* **165**, 255–270 (1999).

78. Röthlisberger, R., Crosta, X., Abram, N. J., Armand, L. & Wolff, E. W. Potential and limitations of marine and ice core sea ice proxies: an example from the Indian Ocean sector. *Quat. Sci. Rev.* **29**, 296–302 (2010).

79. Abram, N. J., Wolff, E. W. & Curran, M. A. J. A review of sea ice proxy information from polar ice cores. *Quat. Sci. Rev.* **79**, 168–183 (2013).

80. Gottschalk, J. et al. Past carbonate preservation events in the deep southeast Atlantic Ocean (Cape Basin) and their implications for Atlantic overturning dynamics and marine carbon cycling. *Paleoceanogr. Paleoclimatol.* **33**, 643–663 (2018).

81. Barker, S. & Diz, P. Timing of the descent into the last ice age determined by the bipolar seesaw. *Paleoceanogr. Paleoclimatol.* https://doi.org/10.1002/2014PA002623 (2014).

82. Sachs, J. & Anderson, R. F. Increased productivity in the subantarctic ocean during Heinrich events. *Nature* **434**, 1118–1121 (2005).

83. Ninnemann, U. S. & Charles, C. D. Changes in the mode of Southern Ocean circulation over the last glacial cycle revealed by foraminiferal stable isotopic variability. *Earth Planet. Sci. Lett.* **201**, 383–396 (2002).

84. Goosse, H. et al. Description of the Earth system model of intermediate complexity LOVECLIM version 1.2. *Geosci. Model Dev.* **3**, 603–633 (2010).

85. Berger, A. Long-term variations of daily insolation and Quaternary climatic changes. *J. Atmos. Sci.* **35**, 2362–2367 (1978).

86. Abe-Ouchi, A., Segawa, T. & Saito, F. Climatic conditions for modelling the Northern Hemisphere ice sheets throughout the ice age cycle. *Clim. Past* **3**, 423–438 (2007).

## Acknowledgements

Core MD95-2039 was collected during the international coring campaign IMAGES Cruise (MD 101) of the French vessel *Marion Dufresne* led by L. Labeyrie, F. Bassinot and Y. Lancelot. This project would not be possible without long-time and dedicated laboratory assistance to process numerous samples by P. Wang, Y. Dai, J. Zhang, L. Wang, P. Chen and L. Ma. This study is supported by NSFC (42076056), Laoshan Laboratory Science and Technology Innovation Project (LSKJ202203300), ARC Discovery Projects (DP190100894) and Future Fellowship (FT140100993) to J.Y. We are indebted to S. Barker for his permission to use stable isotope data collected at Cardiff University by D.J.R.T. as part of NERC grant NE/F002734/1 (to Barker). We thank D. Sigman and J. Gottschalk for discussions.

## Author contributions

J.Y. designed the project, interpreted data and wrote the paper. J.Y. and D.J.R.T. collected core samples with help from N.T., J.T. and X.X. Z.D.J. managed shell picking and IRD counting with help from F.Z. J.Y. supervised B/Ca data collection and analysed data with substantial input from X.J. D.J.R.T. did planktic faunal counts and contributed to stable isotope and B/Ca measurements and data analyses. L.W. provided expertise in physical oceanography. Y.C. and L.T. assisted with stable isotope analyses. L.M. assisted with model outputs and ED Figs. 9 and 10. R.F.A., L.W., E.J.R., J.F.M., D.J.R.T., L.M. and F.Z. commented on the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at https://doi.org/10.1038/s41561-023-01297-x.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41561-023-01297-x.

**Correspondence and requests for materials** should be addressed to J. Yu.

**Peer review information** *Nature Geoscience* thanks Fortunat Joos, Antje Voelker and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: James Super, in collaboration with the *Nature Geoscience* team.

**Reprints and permissions information** is available at www.nature.com/reprints.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 1** | Stadial atmospheric CO₂ changes based on latest high-resolution ice-core records during 30–65 ka. **a**, NGRIP δ¹⁸O²⁵. **b**, West Antarctic Ice Sheet Divide ice core (WDC) CH₄ (ref. 14). **c**, WDC δ¹⁸O⁷². **d**, WDC CO₂ (ref. 3). In c–d, bold curves show probability maxima with shading envelopes representing 2σ uncertainty. **e**, Changes (Δ) in CO₂ for stadials. **f**, Total net ΔCO₂ against duration of stadials. **g**, As **f**, but from ref. 73 (circles) and ref. 74 (crosses) based on older ice-core data¹. Note different y-axis scales between **f** and **g**. Vertical banding as shown in Fig. 2. We follow the literature³ to define

starting and ending dates of stadials. Given paired, high-resolution CH₄ and CO₂ measurements from the same ice core to, respectively, better constrain timings and magnitudes of stadial CO₂ changes, the latest WDC data³,¹⁴ likely provide more accurate information about past rapid CO₂ changes. Antarctica warmed during all Northern Hemisphere stadials (a–c), but atmospheric CO₂ 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₂ and stadial duration (f). GS: Greenland stadial; HS: Heinrich stadial.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 1** | Stadial atmospheric CO₂ changes based on latest high-resolution ice-core records during 30–65 ka. **a**, NGRIP δ¹⁸O²⁵. **b**, West Antarctic Ice Sheet Divide ice core (WDC) CH₄ (ref. 14). **c**, WDC δ¹⁸O⁷². **d**, WDC CO₂ (ref. 3). In c-d, bold curves show probability maxima with shading envelopes representing 2σ uncertainty. **e**, Changes (Δ) in CO₂ for stadials. **f**, Total net ΔCO₂ against duration of stadials. **g**, As **f**, but from ref. 73 (circles) and ref. 74 (crosses) based on older ice-core data¹. Note different y-axis scales between **f** and **g**. Vertical banding as shown in Fig. 2. We follow the literature³ to define

starting and ending dates of stadials. Given paired, high-resolution CH₄ and CO₂ measurements from the same ice core to, respectively, better constrain timings and magnitudes of stadial CO₂ changes, the latest WDC data³,¹⁴ likely provide more accurate information about past rapid CO₂ changes. Antarctica warmed during all Northern Hemisphere stadials (a-c), but atmospheric CO₂ 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₂ and stadial duration (f). GS: Greenland stadial; HS: Heinrich stadial.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 2** | Modern Atlantic meridional distributions of $[CO_3]^{2-}$ (shading) and DIC (contours). Note that around the equator the hydrographic transect changes from the eastern basin in the North Atlantic to the western basin in the South Atlantic. For simplicity, we use NSW and SSW to indicate northern- and 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$^{68}$.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 2** | Modern Atlantic meridional distributions of $[CO_3]^{2-}$ (shading) and DIC (contours). Note that around the equator the hydrographic transect changes from the eastern basin in the North Atlantic to the western basin in the South Atlantic. For simplicity, we use NSW and SSW to indicate northern- and 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$^{68}$.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 3** | Age model and sedimentation rate of core MD95-2039.
a, NGRIP $\delta^{18}$O$^{25}$. Crosses indicate age tie points. b, *G. bulloides* $\delta^{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.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 3** | Age model and sedimentation rate of core MD95-2039.
a, NGRIP $\delta^{18}$O$^{25}$. Crosses indicate age tie points. b, *G. bulloides* $\delta^{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.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 4** | 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**, solid curve shows Pa/Th in cores GGC5 (33.7°N, 57.6°W, 4,550 m) and CDH19 (33.7°N, 57.6°W, 4,541 m) from the Bermuda Rise<sup>29,32,33</sup>, while dashed curve shows data from core SU81-18 (37.8°N, 10.2°W, 3,135 m) on the Iberian Margin<sup>43</sup>. Pa/Th data for HS11 are only available from ODP1063 (33.7°N, 57.6°W, 4,584 m)<sup>29</sup>, 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

<figure><img src="image_17.png" /></figure>

opposite change (**e**), while high-resolution and well-dated records show an IRD peak in MD95-2039 (Extended Data Fig. 3d) and a benthic $\delta^{13}$C decline in Iberian Margin core MD01-2444 (37.66°N, 10.1°W, 2,637 m)<sup>44</sup> (**f**). Compared to the late HS11, the late HS1 is characterized by (i) an overall colder Antarctica and a stalling of Antarctic warming<sup>27</sup> (**d**, **i**), (ii) a larger size but smaller retreat of Southern Ocean sea ice<sup>77</sup> (**b**, **g**), although note different views about using SSNa to reflect sea-ice extent<sup>78,79</sup>, and (iii) a much smaller atmospheric $CO_2$ rise<sup>1</sup> (**c**, **h**). These suggest different Southern Ocean ventilation states during late stages of HS1 and HS11 at the last two glacial terminations.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 4** | 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, solid curve shows Pa/Th in cores GGC5 (33.7°N, 57.6°W, 4,550 m) and CDH19 (33.7°N, 57.6°W, 4,541 m) from the Bermuda Rise<sup>29,32,33</sup>, while dashed curve shows data from core SU81-18 (37.8°N, 10.2°W, 3,135 m) on the Iberian Margin<sup>43</sup>. Pa/Th data for HS11 are only available from ODP1063 (33.7°N, 57.6°W, 4,584 m)<sup>29</sup>, 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

<figure><img src="image_17.png" /></figure>

opposite change (e), while high-resolution and well-dated records show an IRD peak in MD95-2039 (Extended Data Fig. 3d) and a benthic $\delta^{13}$C decline in Iberian Margin core MD01-2444 (37.66°N, 10.1°W, 2,637 m)<sup>44</sup> (f). Compared to the late HS11, the late HS1 is characterized by (i) an overall colder Antarctica and a stalling of Antarctic warming<sup>27</sup> (d, i), (ii) a larger size but smaller retreat of Southern Ocean sea ice<sup>77</sup> (b, g), although note different views about using SSNa to reflect sea-ice extent<sup>78,79</sup>, and (iii) a much smaller atmospheric $CO_2$ rise<sup>1</sup> (c, h). These suggest different Southern Ocean ventilation states during late stages of HS1 and HS11 at the last two glacial terminations.
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 5** | Relative timing of deep-water [CO₃²⁻] and other proxy data during HS4-6 in core MD95-2039. a, *G. bulloides* δ¹⁸O, as an age model control. b, Deep-water [CO₃²⁻]. c, IRD⁵⁹. d, *N. pachyderma*/all planktonics percentage. All data are from the same core and intentionally plotted against

sediment depth. A lead of deep-water [CO₃²⁻] rise (b) relative to *G. bulloides* δ¹⁸O decline (a) is observed during HS4-6. Thick sedimentation during these times helps to minimize any bioturbation effect. Vertical bands show HSs (yellow) and GSs (grey).
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 5** | Relative timing of deep-water [CO₃²⁻] and other proxy data during HS4-6 in core MD95-2039. a, *G. bulloides* δ¹⁸O, as an age model control. b, Deep-water [CO₃²⁻]. c, IRD⁵⁹. d, *N. pachyderma*/all planktonics percentage. All data are from the same core and intentionally plotted against

sediment depth. A lead of deep-water [CO₃²⁻] rise (b) relative to *G. bulloides* δ¹⁸O decline (a) is observed during HS4-6. Thick sedimentation during these times helps to minimize any bioturbation effect. Vertical bands show HSs (yellow) and GSs (grey).
<figure><img src="image_19.png" /></figure>

**Extended Data Fig. 6** | Deep-water proxies and age model in core MD07-3076 from the South Atlantic. Data¹⁶,¹⁷ from MD07-3076 (b-g) are plotted against published age model¹⁷, based on alignment of planktonic faunal census data including *N. pachyderma* and *G. bulloides* abundances to Antarctic temperatures²⁷ (f-h). Curved arrows indicate that the [CO₃²⁻] 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).
<figure><img src="image_19.png" /></figure>

**Extended Data Fig. 6** | Deep-water proxies and age model in core MD07-3076 from the South Atlantic. Data¹⁶,¹⁷ from MD07-3076 (b-g) are plotted against published age model¹⁷, based on alignment of planktonic faunal census data including *N. pachyderma* and *G. bulloides* abundances to Antarctic temperatures²⁷ (f-h). Curved arrows indicate that the [CO₃²⁻] 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).
<figure><img src="image_20.png" /></figure>

**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* δ¹⁸O (f) to Antarctic ice-core temperature proxies (g). Note that the age at -46 ka is possibly too young: *G. bulloides* δ¹⁸O minimum can be shifted older by -2 ka, that is, shifting the dashed rectangle area to HS5 as suggested by the horizontal arrow. This age adjustment would put alkenone changes³² (e) to be consistent with sub-Antarctic Zone export reductions during HSs, as expected from decreased dust deposition associated with a poleward shift of southern

westerlies³⁸. Taking this age uncertainty into account, deep waters at TNO57-21 experienced increases in [CO₃²⁻]²²,⁴⁶ (b) and O₂ (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₃²⁻] changes during many GSs (a, b). Yellow bands: Mode II stadials; Grey bands: Mode IV stadials. Data from refs. 27,80,81.
<figure><img src="image_20.png" /></figure>

**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* δ¹⁸O (f) to Antarctic ice-core temperature proxies (g). Note that the age at -46 ka is possibly too young: *G. bulloides* δ¹⁸O minimum can be shifted older by -2 ka, that is, shifting the dashed rectangle area to HS5 as suggested by the horizontal arrow. This age adjustment would put alkenone changes³² (e) to be consistent with sub-Antarctic Zone export reductions during HSs, as expected from decreased dust deposition associated with a poleward shift of southern

westerlies³⁸. Taking this age uncertainty into account, deep waters at TNO57-21 experienced increases in [CO₃²⁻]²²,⁴⁶ (b) and O₂ (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₃²⁻] changes during many GSs (a, b). Yellow bands: Mode II stadials; Grey bands: Mode IV stadials. Data from refs. 27,80,81.
<figure><img src="image_21.png" /></figure>

**Extended Data Fig. 8 | Deep-water [CO₃²⁻] at TNO57-21 and enhanced Southern Ocean ventilation during HS8.** a, NGRIP δ¹⁸O²⁵. b-h, proxy data in core TNO57-21²²,⁴⁵,⁵⁴,⁸¹,⁸². i, EDML δ¹⁸O²⁵. The age model of TNO57-21 is based on tuning *G. bulloides* δ¹⁸O (*h*) to Antarctic temperature proxy (*i*). The reliability of the age model is supported by benthic δ¹⁸O⁸³ (*h*) and a decline in alkenone fluxes during HS8⁸² (*e*), the latter of which is expected from decreasing dust deposition and hence surface export in the sub-Antarctic Zone³⁸. At TNO57-21,

εNd data⁴⁵ (f) indicate increased mixing of low-[CO₃²⁻] SSW during HS8. Thus, the observed increase in deep-water [CO₃²⁻]²² (b) most likely suggest enhanced ventilation via the Southern Ocean during HS8. This inference is further supported by concomitant benthic δ¹³C increase⁸³ (g) and improvement of deep-water oxygenation based on proxies⁴⁴,⁸² from the same core (c, d). Data from refs. 27,80,81.
<figure><img src="image_21.png" /></figure>

**Extended Data Fig. 8 | Deep-water [CO₃²⁻] at TNO57-21 and enhanced Southern Ocean ventilation during HS8.** a, NGRIP δ¹⁸O²⁵. b-h, proxy data in core TNO57-21²²,⁴⁵,⁵⁴,⁸¹,⁸². i, EDML δ¹⁸O²⁵. The age model of TNO57-21 is based on tuning *G. bulloides* δ¹⁸O (*h*) to Antarctic temperature proxy (*i*). The reliability of the age model is supported by benthic δ¹⁸O⁸³ (*h*) and a decline in alkenone fluxes during HS8⁸² (*e*), the latter of which is expected from decreasing dust deposition and hence surface export in the sub-Antarctic Zone³⁸. At TNO57-21,

εNd data⁴⁵ (f) indicate increased mixing of low-[CO₃²⁻] SSW during HS8. Thus, the observed increase in deep-water [CO₃²⁻]²² (b) most likely suggest enhanced ventilation via the Southern Ocean during HS8. This inference is further supported by concomitant benthic δ¹³C increase⁸³ (g) and improvement of deep-water oxygenation based on proxies⁴⁴,⁸² from the same core (c, d). Data from refs. 27,80,81.
<figure><img src="image_22.png" /></figure>

**Extended Data Fig. 9** | **Changes in northern North Atlantic CO₂ absorption in an Earth system model**⁴⁸. **a**, Timeseries of the simulated Atlantic meridional overturning circulation (AMOC). **b**, Integrated CO₂ 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₂ absorption by the northern North Atlantic. An AMOC weakening leads to greater CO₂ 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⁴⁸. **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⁸⁴. The model is forced from 50 to 30 ka with time-varying boundary conditions including orbital parameters⁸⁵, northern hemispheric ice-sheet orography and albedo⁸⁶, and atmospheric CO₂ concentration²⁸. 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.
<figure><img src="image_22.png" /></figure>

**Extended Data Fig. 9** | **Changes in northern North Atlantic CO₂ absorption in an Earth system model**⁴⁸. **a**, Timeseries of the simulated Atlantic meridional overturning circulation (AMOC). **b**, Integrated CO₂ 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₂ absorption by the northern North Atlantic. An AMOC weakening leads to greater CO₂ 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⁴⁸. **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⁸⁴. The model is forced from 50 to 30 ka with time-varying boundary conditions including orbital parameters⁸⁵, northern hemispheric ice-sheet orography and albedo⁸⁶, and atmospheric CO₂ concentration²⁸. 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.
GS-IS

<figure><img src="image_23.png" /></figure>

SST anomalies

<figure><img src="image_24.png" /></figure>

HS-IS

Surface PO4 anomalies

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

Sea-air CO2 anomalies

<figure><img src="image_27.png" /></figure>

Extended Data Fig. 10 | Simulated climatic and biogeochemical responses to AMOC changes in an Earth system model¹⁸. Anomalies for (a, b) sea surface temperature (SST; unit °C), (c, d) surface [PO₄³⁻] (µmol/kg), and (e, f) sea-air CO₂ (µ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,

<figure><img src="image_28.png" /></figure>

respectively. Enhanced CO₂ absorption during HS and GS are driven by increased solubility pump (related to cooling; a, b) and biological pump (related to [PO₄³⁻] drawdown; c, d). GS = Greenland Stadial at 44.2 ka, IS = Interstadial at 46.6 ka, and HS = Heinrich Stadial at 48.8 ka.
GS-IS

<figure><img src="image_23.png" /></figure>

SST anomalies

<figure><img src="image_24.png" /></figure>

HS-IS

Surface PO4 anomalies

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

Sea-air CO2 anomalies

<figure><img src="image_27.png" /></figure>

Extended Data Fig. 10 | Simulated climatic and biogeochemical responses to AMOC changes in an Earth system model¹⁸. Anomalies for (a, b) sea surface temperature (SST; unit °C), (c, d) surface [PO₄³⁻] (µmol/kg), and (e, f) sea-air CO₂ (µ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,

<figure><img src="image_28.png" /></figure>

respectively. Enhanced CO₂ absorption during HS and GS are driven by increased solubility pump (related to cooling; a, b) and biological pump (related to [PO₄³⁻] drawdown; c, d). GS = Greenland Stadial at 44.2 ka, IS = Interstadial at 46.6 ka, and HS = Heinrich Stadial at 48.8 ka.