More efficient North Atlantic carbon pump during the Last Glacial Maximum.

During the Last Glacial Maximum (LGM; ~20,000 years ago), the global ocean sequestered a large amount of carbon lost from the atmosphere and terrestrial biosphere. Suppressed CO2 outgassing from the Southern Ocean is the prevailing explanation for this carbon sequestration. By contrast, the North Atlantic Ocean-a major conduit for atmospheric CO2 transport to the ocean interior via the overturning circulation-has received much less attention. Here we demonstrate that North Atlantic carbon pump efficiency during the LGM was almost doubled relative to the Holocene. This is based on a novel proxy approach to estimate air-sea CO2 exchange signals using combined carbonate ion and nutrient reconstructions for multiple sediment cores from the North Atlantic. Our data indicate that in tandem with Southern Ocean processes, enhanced North Atlantic CO2 absorption contributed to lowering ice-age atmospheric CO2.

Introduction

The North Atlantic Ocean (>~35°N, including the ass="Disease">Nordic Seas and Arctic Ocean) is a major atmospheric ass="Chemical">pan class="Chemical">CO2 sink, which has been mitigating anthropogenic atmospheric CO2 increases[1]. Preindustrial North Atlantic surface water partial pressure of CO2 (pCO2) was up to ~100 μatm lower than the contemporary atmospheric pCO2 of ~280 μatm, which caused substantial atmospheric CO2 invasion[2,3]. Despite its modest area, the North Atlantic Ocean accounts for at least ~30% of the global ocean CO2 uptake today and during preindustrial times[1,4]. Over longer timescales, large-scale oceanic carbon sequestration also occurred during Plio-Pleistocene glaciations[5-7]. This is commonly attributed to reduced glacial Southern Ocean CO2 outgassing[6,8,9], while even the sign of past North Atlantic CO2 uptake efficiency changes remains unconstrained. Here, we present a novel proxy approach to trace atmospheric CO2 invasion in the North Atlantic and thereby evaluate its role in carbon sequestration in ice-age oceans. We find that the last glacial North Atlantic carbon absorption became more efficient, highlighting a critical role of the North Atlantic Ocean in regulating glacial–interglacial atmospheric CO2 changes.

Results

Air–sea CO2 exchange tracers

Any effect of ocean processes on atmospheric ass="Chemical">pCO2 must occur via air–sea ass="Chemical">pan class="Chemical">CO2 exchange. In the North Atlantic, high-nutrient utilization decreases surface-water dissolved inorganic carbon (DIC) and causes surface-water pCO2 to be lower than atmospheric pCO2 (Supplementary Fig. 1). This leads to net air-to-sea CO2 transfer, creating an air–sea exchange signature of DIC (DICas). DICas signals can be distinguished by accounting for within-ocean DIC redistributions that are heavily mediated by biology (Fig. 1). Biological cycling of organic matter depletes DIC and nutrients such as phosphate (PO4) in surface waters and enriches them at depth. Seawater mixing also affects DIC and PO4 concentrations in the ocean. Nevertheless, PO4 variations are ultimately determined by biological processes: without biology, PO4 should be the same everywhere in the ocean regardless of ocean circulation (ignoring the small effect from salinity change). Because marine biology incorporates and releases PO4 and DIC in a relatively fixed proportion following Redfield stoichiometry[3,10] and because PO4 is not affected by air–sea exchange, PO4 can be used to estimate biology-driven within-ocean DIC redistributions (Fig. 1). Any within-ocean DIC redistribution associated with CaCO3 cycling can be accounted for using alkalinity (ALK) and nitrate.

Fig. 1

Concepts to distinguish DICas. For simplicity, only CO2 invasion associated with organic matter cycling is considered. In the ocean box, vertical solid and dashed lines (a–d) represent mean PO4 (blue) and DIC (red) in an abiotic ocean (a). Biology redistributes DIC and PO4 following Redfield stoichiometry (curves; b). This decreases surface-ocean DIC and pCO2, and hence causes air-to-sea CO2 transfer (c). Through mixing and ocean circulation, CO2 invasion raises water-column DIC, i.e., shifting dashed curve (equals the red-solid curve in b) to red-solid curve (c). The shaded region in c represents air–sea exchange DICas signatures. After removing carbon redistribution by biology based on PO4-related curvature of the profiles (b), DICas can be revealed by the shaded region in d

Concepts to distinguish ass="Chemical">DICas. For simplicity, only ass="Chemical">pan class="Chemical">CO2 invasion associated with organic matter cycling is considered. In the ocean box, vertical solid and dashed lines (a–d) represent mean PO4 (blue) and DIC (red) in an abiotic ocean (a). Biology redistributes DIC and PO4 following Redfield stoichiometry (curves; b). This decreases surface-ocean DIC and pCO2, and hence causes air-to-sea CO2 transfer (c). Through mixing and ocean circulation, CO2 invasion raises water-column DIC, i.e., shifting dashed curve (equals the red-solid curve in b) to red-solid curve (c). The shaded region in c represents air–sea exchange DICas signatures. After removing carbon redistribution by biology based on PO4-related curvature of the profiles (b), DICas can be revealed by the shaded region in d

Following the established method[3] to account for within-ocean DIC redistributions by soft-tissue and ass="Chemical">CaCO3 cycling, we calculate preindustrial Atlantic ass="Chemical">pan class="Chemical">DICas using the GLODAP dataset[2] (Fig. 2a). See Methods for details to calculate DICas. More positive DICas values indicate a greater degree of atmospheric CO2 invasion. At basin-scale, the preindustrial DICas of North Atlantic deep water (NADW) is ~50–80 µmol/kg higher than for Antarctic bottom water (AABW) and Antarctic intermediate water (AAIW). This difference reflects North Atlantic CO2 uptake and Southern Ocean release[3,11]. North Atlantic CO2 absorption is driven by (i) an efficient solubility pump due to strong cooling of northward-flowing Gulf Stream waters and (ii) a strong biological pump associated with high nutrient utilization[12-14]. NADW thus represents an efficient pathway for atmospheric CO2 sequestration[6,15]. Through global deep ocean circulation, CO2 absorbed in the North Atlantic is transported throughout the world ocean[1,3], with profound implications for the global carbon cycle.

Fig. 2

Preindustrial Atlantic air–sea exchange tracers. a DICas. b [CO32−]as. Circles represent studied sediment cores. Inset: GLODAP hydrographic data[2] used to generate the sections[96]. NADW North Atlantic deep water, AABW Antarctic bottom water, AAIW Antarctic intermediate water. See Methods for calculation details

Preindustrial Atlantic air–sea exchange tracers. a ass="Chemical">DICas. b [ass="Chemical">pan class="Chemical">CO32−]as. Circles represent studied sediment cores. Inset: GLODAP hydrographic data[2] used to generate the sections[96]. NADW North Atlantic deep water, AABW Antarctic bottom water, AAIW Antarctic intermediate water. See Methods for calculation details

No proxy exists to reconstruct past seaass="Chemical">water DIC and ALK at acceptable precision for direct application, so we employ a linked ass="Chemical">pan class="Chemical">carbonate system parameter for palaeoceanographic studies. Everything else being equal, atmospheric CO2 invasion would decrease seawater carbonate ion concentration ([CO32−]), because CO2 reacts with carbonate ion to form bicarbonate[16]. We thus develop a new tracer, [CO32−]as, which essentially reflects seawater [CO32−] contrasts for the same biological (i.e., PO4) and physical (i.e., temperature–salinity–pressure; T–S–P) conditions (Fig. 2b; see Methods for calculation details). To extract air–sea exchange signals, it is necessary to compare [CO32−] at the same PO4–T–S–P conditions because we must first remove influences on [CO32−] from (i) within-ocean DIC and ALK redistributions by biology and (ii) T−S–P variations via their effects on CO2 system dissociation constants[16]. In the preindustrial Atlantic, the strong negative correlation between [CO32−]as and DICas (Fig. 2, Supplementary Fig. 2) indicates that [CO32−]as variations are affected only by DICas, and thus are ultimately linked to air–sea CO2 exchange.

The Gulf Stream is a major ass="Chemical">NADW source[17]; thus, comass="Chemical">paring the [ass="Chemical">pan class="Chemical">CO32−]as gradient between the Gulf Stream and NADW can provide a measure of CO2 sequestration intensity during transformation of Gulf Stream waters into NADW. Because Gulf Stream waters are more or less in equilibrium with atmospheric pCO2 from ~10°N to 35°N1,2, the Gulf Stream−NADW [CO32−]as gradient mainly reflects North Atlantic (>~35°N) air–sea CO2 exchange (Supplementary Fig. 3). Physical oceanographers have shown that the path of Gulf Stream waters, rather than being a direct conveyor to the polar North Atlantic, is instead a “corkscrew”, where Gulf Stream waters are recirculated south in the subtropical gyre and subduct after being made more dense by air–sea heat loss (e.g., refs. [18,19]). However, our interest lies in net CO2 uptake by the North Atlantic region, and variations in spatial pathways from Gulf Stream to NADW formation sites[18,19] should not significantly complicate our conclusion. The greater the [CO32−]as gradient between Gulf Stream and NADW (instead of their absolute [CO32−]as values), the more efficient air–sea CO2 absorption by the North Atlantic. Linked to large-scale overturning circulation, Gulf Stream−NADW [CO32−]as gradient changes regulate long-term CO2 sequestration into the deep ocean.

Downcore reconstructions

Next, we reconstruct past Gulf Stream–ass="Chemical">NADW [ass="Chemical">pan class="Chemical">CO32−]as gradients to investigate North Atlantic carbon pump efficiency during the LGM (18–27 ka). Previous work suggests that most of North Atlantic subtropical gyre water circulates through the Caribbean Sea before being transported to the subpolar North Atlantic via the Gulf Stream[20]. We, therefore, use Caribbean Sea ODP Site 999 (12.8°N, 78.7°W) to constrain past Gulf Stream physicochemical conditions (Fig. 3, Supplementary Figs. 4 and 5). The feasibility of using ODP Site 999 to reflect the first-order Gulf Stream carbonate chemistry changes between the Holocene and LGM is supported by observations that (i) Caribbean surface waters have similar [CO32−]as values to hydrographic sites located within Gulf Stream during the preindustrial (Supplementary Fig. 3), and (ii) cores from the broader western subtropical Atlantic show comparable Holocene and LGM [CO32−]as signatures as those from ODP 999 (Supplementary Fig. 6). Surface-water T and S are estimated from Globigerinoides ruber Mg/Ca and sea level fluctuations, respectively[21,22]. Previously published G. ruber δ11B (ref. [21]) is used to calculate surface-water pH, while ALK is estimated from S using the modern relationship between S and ALK[21,22]. Along with T, S, and ALK estimates, pH is then used to calculate surface-water [CO32−] and DIC. Given the constraint from pH, seawater ALK and DIC must vary systematically within the ocean carbonate system (Supplementary Fig. 5). This allows precise estimation of [CO32−], because even large ALK uncertainties (100 μmol/kg; ± 2σ, used throughout) only have a minor effect on [CO32−] (~14 μmol/kg). Given its oligotrophic setting, past surface-water PO4 at ODP 999 is assumed to be zero[2,21,22].

Fig. 3

Down core reconstructions. a ODP 999. b BOFS 17 K. c BOFS 14 K. d BOFS 11 K. Seawater [CO32−] values are derived from benthic B/Ca (empty circles) and δ11B (solid circles). Light gray envelopes and error bars: 2σ. Note different y-scales for surface- (ODP 999) and deep-water (BOFS cores) reconstructions. See Methods for reconstruction details

Down core reconstructions. a ass="Chemical">ODP 999. b BOFS 17 K. c BOFS 14 K. d BOFS 11 K. Seaass="Chemical">pan class="Chemical">water [CO32−] values are derived from benthic B/Ca (empty circles) and δ11B (solid circles). Light gray envelopes and error bars: 2σ. Note different y-scales for surface- (ODP 999) and deep-water (BOFS cores) reconstructions. See Methods for reconstruction details

Three cores are used to reconstruct ass="Chemical">deep-ass="Chemical">pan class="Chemical">water conditions of northern-sourced waters (Fig. 3). BOFS 17 K (58°N, 16.5°W, 1150 m) and BOFS 14 K (58.6°N, 19.4°W, 1756 m) are located close to the previously surmised center of Glacial North Atlantic intermediate water (GNAIW)[23], while BOFS 11 K (55.2°N, 20.4°W, 2004 m) is thought to be affected by glacial Nordic Sea overflows[24]. We employ benthic foraminiferal δ11B and B/Ca to reconstruct deep-water [CO32−] with an uncertainty of ~10 μmol/kg[25]. δ11B and B/Ca give consistent downcore [CO32−] reconstructions. Benthic Cd/Ca is used to estimate deep-water Cd and PO4 based on an established approach (Supplementary Fig. 7)[26,27]. Past deep-water T and S changes are estimated from foraminiferal δ18O and sea level fluctuations; use of other methods negligibly affects our conclusion. In total, we present 180 new measurements for benthic foraminiferal δ11B, B/Ca, and Cd/Ca. Details of core materials, methods, new and compiled data, and fully propagated uncertainties are given in Methods and Supplementary Data 1–9.

A pragmatic recipe to estimate [CO32−]as change

Surface-ass="Chemical">water [ass="Chemical">pan class="Chemical">CO32−] at ODP 999 is ~150 μmol/kg higher than deep-water values at BOFS cores (Fig. 3), but this [CO32−] contrast includes influences from physical (via dissociation constants) and biological (via within-ocean DIC and ALK redistributions) changes in addition to any air–sea CO2 changes between surface and deep waters. Below, we present a pragmatic recipe to estimate [CO32−]as gradients between water masses. We take advantage of well-defined sensitivities of [CO32−] to T–S–P (Fig. 4) to calculate normalized seawater [CO32−] ([CO32−]Norm) at conditions of T = 3 °C, S = 35‰, and P = 2500 dbar (Methods). Any variation in T–S–P would affect seawater [CO32−] via (i) changing CO2 system dissociation constants, and (ii) altering the solubility pump and thereby air–sea exchange component CO2 concentrations in seawater. Calculation of [CO32−]Norm only corrects for influences from (i), without affecting any air–sea CO2 signal. After normalization to constant T–S–P conditions and assuming no net air–sea exchange, biological activity drives changes in both [CO32−]Norm and PO4 along the biological trend (green curves in Fig. 5; Methods). Note that along a certain biological trend, seawater [CO32−]Norm and PO4 are only affected by within-ocean DIC and ALK redistributions (Fig. 1b). A net air–sea CO2 change would cause changes in [CO32−]Norm and PO4 across biological curves. At the same PO4, [CO32−]Norm contrasts reflect [CO32−]as gradients due to air–sea CO2 exchange between water masses.

Fig. 4

Carbonate system sensitivities to various changes. a Salinity effect. b Temperature effect. c Pressure effect. d Biological effect. e Air–sea CO2 exchange effect. Calculations are based on GLODAP[2] (n = 55,399; blue) and a LGM output from LOVECLIM[58] (n = 71,768; gray). For a–d, calculations assume no net air–sea CO2 change. Best fits of data are shown by red curves. See Methods for calculation details

Fig. 5

[CO32−]Norm vs. PO4. a Preindustrial Atlantic surface (<100 m, north of 10°N) and deep (>1000 m, 65°N–65°S) water data[2]. b Holocene and LGM data. Error bars, 2σ

ass="Chemical">Carbonate system sensitivities to various changes. a Salinity effect. b Temperature effect. c Pressure effect. d Biological effect. e Air–sea ass="Chemical">pan class="Chemical">CO2 exchange effect. Calculations are based on GLODAP[2] (n = 55,399; blue) and a LGM output from LOVECLIM[58] (n = 71,768; gray). For a–d, calculations assume no net air–sea CO2 change. Best fits of data are shown by red curves. See Methods for calculation details

[ass="Chemical">CO32−]Norm vs. ass="Chemical">pan class="Chemical">PO4. a Preindustrial Atlantic surface (<100 m, north of 10°N) and deep (>1000 m, 65°N–65°S) water data[2]. b Holocene and LGM data. Error bars, 2σ

A plot of [ass="Chemical">CO32−]Norm vs. ass="Chemical">pan class="Chemical">PO4 greatly facilitates investigation of air–sea CO2 exchange from combined [CO32−] and PO4 measurements/reconstructions. Compared to the biological trend, preindustrial North Atlantic surface waters have a steeper trend (Fig. 5a), which reflects CO2 absorption during northward transport. Deep-water data lie on a shallower trend, consistent with mixing between low-[CO32−]as (high DICas) NADW and high-[CO32−]as (low DICas) AABW in the deep Atlantic (Fig. 2).

For our downcore reconstructions, benthic Cd/Ca suggests that ass="Chemical">deep-waters at the BOFS sites had ass="Chemical">pan class="Chemical">PO4 values of ~1.2 and ~0.8 μmol/kg during the Holocene and LGM, respectively (Fig. 5b; Supplementary Fig. 8). Assuming no air–sea CO2 exchange, [CO32−]Norm of ODP 999 surface waters at elevated PO4 due to biological processes can be estimated straightforwardly using the H→H′ and G→G′ trajectories in Fig. 5b for the Holocene and LGM, respectively. For the Holocene, ODP 999 [CO32−]Norm is ~56 ± 8 μmol/kg higher than [CO32−]Norm of BOFS cores at PO4 = 1.2 μmol/kg. For the LGM, ODP 999 [CO32−]Norm is ~114 ± 9 μmol/kg higher than [CO32−]Norm of BOFS cores at PO4 = 0.8 μmol/kg. This suggests a Holocene-to-LGM increase of ~58 ± 12 μmol/kg in the ODP 999−BOFS [CO32−]as gradient.

We also present a second approach to calculate [ass="Chemical">CO32−]as gradients, which involves frequent use of the ass="Chemical">pan class="Chemical">CO2sys program[28] and intermediate-step ALK and DIC parameters (Supplementary Note 1; Supplementary Figs. 9 and 10). The approach gives similar results as the above pragmatic recipe, because both methods are essentially based on the same principle, which is to compare [CO32−] of water masses at the same physical and biological conditions.

Enhanced CO2 uptake in the glacial North Atlantic

What caused the greater ass="Chemical">ODP 999−BOFS [ass="Chemical">pan class="Chemical">CO32−]as gradient during the LGM? We consider influences from biogenic matter composition variations, surface-water ALK and PO4 changes, ocean circulation changes, and North Atlantic air–sea exchange. In Fig. 5, we have used a soft-tissue Redfield C/PO4 of 127 and a rain ratio (R, = Corganic:) of 4 (refs. [3,10,29,30]) to predict the biological trend. Raising LGM C/PO4 to 140 (the high end value in today’s North Atlantic[30]) and R to 8 (doubling of the modern value) could lower the LGM [CO32−]as gradient by ~16 μmol/kg (Supplementary Fig. 15), still leaving ~42 μmol/kg [CO32−]as gradient increase to be explained by other processes. Evidence for such large biological changes is lacking. Importantly, any increase in C/PO4 and R would implicitly sequester more atmospheric CO2 via an enhanced soft-tissue pump and weakened carbonate pump[15]. Inclusion of a whole ocean ALK inventory change[6] or any increased glacial surface-water PO4 at ODP 999 would raise the LGM [CO32−]as gradient (Supplementary Figs. 16 and 17).

Regarding ocean circulation changes, most AAIW upwells in the tropics and less than ~25% of today’s ass="Chemical">NADW is fed directly by AAIW without surfacing at low latitudes[17]. Northward AAIW transport is thought to have been reduced substantially in the glacial Atlantic[23,31-33] in the face of vigorous GNAIW production[34]. Assuming a constant total ass="Chemical">pan class="Chemical">carbon uptake by the North Atlantic, a complete shutdown of AAIW contribution would only raise the ODP 999−BOFS [CO32−]as gradient by ~30%, which is much smaller than the ~100% increase from the Holocene (~56 μmol/kg) to LGM (~114 μmol/kg) (Fig. 5b). Any increased mixing of glacial AABW at BOFS sites would reduce the ODP 999−BOFS [CO32−]as gradient during the LGM. Given the proximity of our deep-water sites to the core of GNAIW and Nordic Sea overflow waters[23,24,31,35,36], the larger LGM [CO32−]as gradient between ODP 999 and BOFS cores likely reflects a greater DICas increase from Gulf Stream to GNAIW. North Atlantic CO2 invasion was responsible for the preindustrial Gulf Stream-NADW [CO32−]as gradient (Fig. 2). Therefore, we ascribe the increased ODP 999−BOFS [CO32−]as gradient during the LGM to more efficient atmospheric CO2 uptake via air–sea exchange and subsequent transport to at least ~2 km depth (BOFS 11 K core depth) in the glacial North Atlantic.

Quantification of North Atlantic CO2 uptake

With reconstructed ass="Chemical">ODP 999−BOFS [ass="Chemical">pan class="Chemical">CO32−]as gradients, we further quantify North Atlantic air–sea CO2 absorption changes between the Holocene and LGM. [CO32−]as/DICas sensitivities can be precisely estimated (Fig. 4e), making [CO32−]as gradients a useful proxy to calculate DICas changes. The 58 ± 12 μmol/kg Holocene-to-LGM [CO32−]as increase (Fig. 5b) indicates a DICas increase of 91 ± 20 μmol/kg due to enhanced North Atlantic air–sea CO2 absorption (Methods). Compared to the preindustrial Gulf Stream-NADW DICas gradient of ~90 μmol/kg (Fig. 2, Supplementary Fig. 3), this suggests a doubling of CO2 uptake efficiency in the LGM North Atlantic.

Beside ass="Chemical">DICas gradient changes, which indicate air–sea ass="Chemical">pan class="Chemical">CO2 uptake efficiency, knowledge of northern-sourced-water volumes in the global deep ocean is required to determine total North Atlantic carbon sequestration. Figure 6 shows the total extra carbon absorbed by the LGM North Atlantic for a range of northern-sourced-water volumes (Methods). Sedimentary Pa/Th, radiocarbon, neodymium isotopes, and paired benthic Cd/Ca–δ13C suggest[32,34,35,37] vigorous glacial northern-sourced intermediate water production and subsequent transport to the remaining world ocean. Based on previous estimates[35,36,38,39], we tentatively assume that NADW- and GNAIW-derived waters occupy ~50% and ~30%, respectively, of the global deep ocean volume (1 × 1018 m3 for >1 km). In this case, our ~91 μmol/kg Holocene-to-LGM DICas increase yields ~100 Petagrams of carbon (PgC; 1 Pg = 1 × 1015 g) greater CO2 sequestration by the LGM North Atlantic (Fig. 6; Methods). To maintain similar total carbon uptake between the Holocene and LGM, GNAIW would need to be less than ~50% of NADW in volume, which we consider unlikely given evidence for intensive GNAIW export to the global ocean[23,32,34,35,37]. We acknowledge uncertainties associated with our calculations, and encourage future work to better constrain volumes and carbonate chemistry changes of various water masses in the past.

Fig. 6

North Atlantic CO2 budget. The LGM–Holocene extra carbon uptake is based on Holocene-to-LGM DICas increase of 91 μmol/kg. The large red square represents our best estimate of ~100 PgC, assuming that NADW and GNAIW occupied ~50% and ~30% of the global deep ocean (>1 km), respectively[35,36,38,39]. See Methods for calculation details

North Atlantic ass="Chemical">CO2 budget. The ass="Chemical">pan class="Chemical">LGM–Holocene extra carbon uptake is based on Holocene-to-LGM DICas increase of 91 μmol/kg. The large red square represents our best estimate of ~100 PgC, assuming that NADW and GNAIW occupied ~50% and ~30% of the global deep ocean (>1 km), respectively[35,36,38,39]. See Methods for calculation details

Discussion

Previous work[40-42] has tried to constrain air–sea ass="Chemical">CO2 exchange by reconstructing surface conditions. This requires reconstructions of the air–sea ass="Chemical">pan class="Chemical">pCO2 difference (influenced by T, S, and nutrient utilization), the gas transfer velocity (a power function of wind speed), solubility of CO2 in seawater (mainly affected by T), and the area and contact time of surface waters available for air–sea exchange[1]. Sea ice cover[43] possibly expanded, reducing glacial North Atlantic CO2 absorption. A larger LGM meridional surface temperature gradient[43,44] would enhance the North Atlantic solubility pump[13]. Existing planktonic δ15N and Cd/Ca data[40,45] show conflicting results regarding the glacial North Atlantic nutrient conditions, perhaps due to complications associated with surface-water proxies and spatial/seasonal nutrient variations in the North Atlantic. A decreased preformed nutrient in the glacial North Atlantic might be inferred from a lower GNAIW PO4 (Fig. 3), but faster ventilation and/or reduced glacial AAIW could also cause a nutrient decline in GNAIW[23,34,46]. Little is known about past wind intensity and air–sea contact time changes. Consequently, potential North Atlantic glacial CO2 invasion remains poorly understood. Bypassing the necessity to reconstruct surface-water conditions for which some proxies are still lacking (e.g., wind), our new approach, to our knowledge, offers the first proxy-based quantitative estimate of air–sea CO2 uptake efficiency in the glacial North Atlantic.

In contrast to previous calculations[47-49] which concern combined biological (i.e., within-ocean DIC redistribution) and air–sea exchange ass="Chemical">carbon changes (Fig. 1c), our total North Atlantic ass="Chemical">pan class="Chemical">carbon uptake estimate only represents the net air–sea CO2 change that is more directly relevant to atmospheric and terrestrial carbon inventory variations. Our estimated ~100 PgC sequestration constitutes ~15% of the Holocene-LGM ~600 PgC change associated with the atmosphere (~200 PgC) and terrestrial biosphere (~400 PgC)[5,6]. Given this global carbon budget context, our work reinforces the role of other polar regions (e.g., Southern Ocean) in controlling the glacial–interglacial carbon cycle. However, if there were no efficiency enhancement for the LGM North Atlantic, a 40% shrinkage of NADW volume would decrease air–sea component CO2 sequestration by ~240 PgC in the deep ocean (Methods). Therefore, by overcoming this opposing “volume effect”, the improved glacial North Atlantic efficiency increased DICas values of northern-sourced deep waters (termed the “endmember effect”) and thereby contributed substantially to air–sea CO2 sequestration in the LGM deep ocean.

Atmospheric ass="Chemical">pCO2 is controlled by both ass="Chemical">pan class="Chemical">CO2 gains (e.g., via Southern Ocean outgassing) and losses (e.g., via North Atlantic absorption)[2,3,11]. Growing evidence indicates that processes outside the Southern Ocean may have affected past atmospheric CO2 variations[50-52]. Our proxy-based results indicate that the North Atlantic CO2 pump efficiency during the LGM was almost doubled relative to the Holocene. This increased efficiency and associated “endmember effect” effectively outcompeted the opposing “volume effect” due to any shrinkage of northern-sourced deep waters in the world ocean. In addition to the well-recognized role of reduced outgassing in the Southern Ocean[6,8,9,47,53,54], we therefore suggest that variations in the uptake and sequestration of atmospheric CO2 via the North Atlantic Ocean were important contributors to glacial/interglacial carbon cycling.

Methods

CO2 system calculations

For both the preindustrial ocean and down-core ass="Chemical">CO2 system calculations, ass="Chemical">pan class="Chemical">seawater carbonate system variables were calculated using the CO2sys.xls program[28] with dissociation constants K1 and K2 according to Mehrbach et al.[55] and according to Dickson[56]. Seawater total boron concentration was calculated from the boron–salinity relationship of Lee et al.[57]. For the GLODAP dataset, the anthropogenic CO2 contribution was subtracted from the measured DIC to obtain preindustrial DIC values[2].

Preindustrial Atlantic DICas and [CO32−]as

The GLODAP dataset[2] is used to calculate preindustrial ocean ass="Chemical">CO2 system variables. Following the established method of Broecker and Peng[3], we account for ass="Chemical">pan class="Disease">DIC anomalies created by (1) freshwater addition or removal based on S, (2) soft-tissue carbon creation and respiration based on PO4, and (3) CaCO3 formation and dissolution based on ALK and nitrate (NO3). See Fig. 1 for the simplified concept. We adopt the term DICas to represent net air–sea exchange component DIC signatures from:

where the subscript “s” represents values normalized to S of 35 (e.g., DICs = DIC × 35/S); the superscript “mo” denotes mean ocean values at S = 35 (ass="Chemical">PO4mo = 2.2 μmol/kg, ass="Chemical">pan class="Chemical">ALKmo = 2383 μmol/kg, DICmo = 2267 μmol/kg, and NO3mo = 31 μmol/kg)[29]; C/PO4 represents the soft-tissue stoichiometric Redfield ratio; and the arbitrary DICconstant ( = 2285 μmol/kg) is designed to bring zero DICas close to the NADW–AABW boundary (Fig. 2). The term (PO4s − PO4mo) × C/PO4 corrects for DIC changes due to photosynthesis and soft-tissue degradation, and the term ½ × (ALKs − ALKmo + NO3s − NO3mo) accounts for DIC changes caused by CaCO3 formation and dissolution. To be consistent with previous work[3,30], we used C/PO4 = 127 to calculate DICas and [CO32−]as in Fig. 2. Using other C/PO4 values[10] does not significantly affect spatial DICas and [CO32−]as patterns (Supplementary Figs. 11 and 12). Neither are their patterns affected by using other PO4–ALK–NO3 values to replace global mean values in Eq. (1) (Supplementary Figs. 13 and 14). Ideally, DICconstant would be the mean DIC value of an abiotic ocean (Fig. 1), but this value cannot be simply determined from modern observations. Because our interest lies in spatial DICas contrasts instead of absolute values, the choice of DICconstant has no effect on our interpretation.

To obtain [ass="Chemical">CO32−]as, we first calculate [ass="Chemical">pan class="Chemical">CO32−]PO4–T–S–P using (DICas + DICconstant), ALKmo, and PO4mo at T = 3 °C, S = 35, and P = 2500 dbar. [CO32−]as is then calculated by [CO32−]as = [CO32−]PO4–T–S–P − [CO32−]constant, where [CO32−]constant ( = 78 μmol/kg, calculated using DICconstant and ALKmo) is designed to bring zero [CO32−]as close to the NADW–AABW boundary. In essence, the [CO32−]as distribution reflects the variation of [CO32−] when normalized to the same PO4−T−S−P conditions.

CO2 system sensitivities and calculation of [CO32−]Norm

Because the seaass="Chemical">water ass="Chemical">pan class="Chemical">CO2 system is nonlinear, there is currently no simple way to derive these sensitivities based on CO2 system equations[16]. We use GLODAP preindustrial data[2] to calculate numerically [CO32−] sensitivities to various physiochemical parameters. Use of LGM outputs from the LOVECLIM model[58] yields comparable sensitivities. We first use hydrographic data, including T, S, P, DIC, ALK, PO4, and SiO3 to calculate [CO32−]. We then change S to 35‰ and other chemical concentrations proportionally. For example, ALK and DIC will change as follows:We use S = 35‰, ALKS=35, DICS=35, [PO4]S=35, and [SiO3]S=35 along with hydrographic T and P to calculate [CO32−]S=35. The [CO32−] to S sensitivity (Sen_S) is calculated by:To estimate temperature effects, we calculate [CO32−]S=35, T=3 °C using S = 35‰, ALKS=35, DICS=35, [PO4]S=35, [SiO3]S=35, T = 3 °C, and hydrographic P. The sensitivity of [CO32−]S=35 to temperature (Sen_T) is defined by:

Regarding pressure effects, we calculate [ass="Chemical">CO32−]S=35, T=3 °C, P=2500 ass="Chemical">pan class="Chemical">dbar using S = 35‰, ALKS=35, DICS=35, [PO4]S=35, [SiO3]S=35, T = 3°C, and P = 2500 dbar. The sensitivity of [CO32−]S=35, T=3°C to P (Sen_P) is defined by:

To estimate the influence on [ass="Chemical">CO32−] from within-ocean ALK–DIC redistributions by biological processes, we assume a 0.1 μmol/kg increase in ass="Chemical">pan class="Chemical">PO4 (i.e., ΔPO4 = 0.1 μmol/kg) due to biological respiration (photosynthesis has an opposite effect). The resultant ALK (ALKS=35+respiration) and DIC (DICS=35+respiration) can then be calculated from:Resultant [CO32−] ([CO32−]Norm+respiration) values are calculated using DICS=35+respiration, ALKS=35+respiration, and ([PO4]S=35 + ΔPO4) at constant physical conditions of T = 3 °C, S = 35, and P = 2500 dbar. The sensitivity of [CO32−]Norm to PO4 is defined by:We consider four Redfield stoichiometric scenarios: C/PO4 = 127, R = 4 (the reference composition; Fig. 4d); C/PO4 = 140, R = 4; C/PO4 = 127, R = 8; and C/PO4 = 140, R = 8 (Supplementary Fig. 15). In all cases, strong exponential correlations exist between [CO32−]Norm/PO4 sensitivities and [CO32−]Norm. The correlations may reflect the buffering effect of the seawater CO2 system: for seawater with high DIC (low [CO32−] and high buffering capability), [CO32−] would be relatively less sensitive to biological DIC and ALK disturbances. All of the above sensitivity calculations assume no net air–sea CO2 change.

To calculate air–sea exchange sensitivities, we assume a 10 μmol/kg increase in DICS=35 due to atmospheric ass="Chemical">CO2 invasion (i.e., Δass="Chemical">pan class="Chemical">DICas = 10 μmol/kg). We calculate [CO32−]Norm+as using S = 35‰, ALKS=35, DICS=35+as ( = DICS=35 + ΔDICas), [PO4]S=35, [SiO3]S=35, T = 3 °C, and P = 2500 dbar. The sensitivity of [CO32−]as to DICas is defined by:

Using sensitivities shown in Fig. 4, [pan class="Chemical">CO32−]Norm can be calculated by:Excel spreadsheets are provided in Supplementary Data 7–8 to calculate [ass="Chemical">pan class="Chemical">CO32−]Norm and the biological curves shown in Fig. 5.

LGM–Holocene North Atlantic carbon budget

The total extra ass="Chemical">carbon increase (Δ∑Cass="Chemical">pan class="Chemical">LGM−Holocene) in Fig. 6 is calculated by Δ∑CLGM−Holocene = V × density × %GNAIW × ([CO32−]as_ODP999-BOFSLGM/0.61) × 12 − V × density × %NADW × ([CO32−]as_ODP999-BOFSHolocene/0.59) × 12, where V is the global deep ocean volume (>1 km water depth) at 100.8 × 1016 m3, density = 1027.8 kg/m3 (ref. [29]), %GNAIW and %NADW, respectively, represent their volume fractions in the deep ocean, [CO32−]as_ODP999-BOFSHolocene = 56 μmol/kg, [CO32−]as_ODP999-BOFSLGM = 114 μmol/kg (Fig. 5), terms 0.61 and 0.58, respectively, represent the absolute LGM and Holocene [CO32−]as/DICas sensitivities (Fig. 4e) used to transfer [CO32−]as_ODP999-BOFS into ODP999–BOFS DICas contrasts (LGM: 186 μmol/kg; Holocene: 95 μmol/kg), and the number 12 converts C from moles into weight. Based on previous estimates, %NADW is thought to be ~50% (refs. [38,39]), while %GNAIW remained roughly similar to %NADW or shrank (refs. [35,36]). These estimates are debated and have large uncertainties, and we thus calculate Δ∑CLGM−Holocene for a range of %NADW and %GNAIW values (Fig. 6). Any influence from AAIW is ignored because of its similar [CO32−]as signals to Gulf Stream during the Holocene (Supplementary Fig. 3) and much reduced northward advection during the LGM[23,31-33]. We tentatively treat Δ∑CLGM-Holocene of ~100 PgC using %NADW = 50% and %GNAIW = 30% as our best estimate. Assuming no Holocene–LGM DICas gradient change (i.e., the same CO2 uptake efficiency) and everything else being equal, Δ∑CLGM–Holocene would be −240 PgC at %NADW = 50% and %GNAIW = 30%.

Cores, age models, samples, and analytical methods

We used ass="Chemical">ODP Site 999 for Gulf Stream surface-ass="Chemical">pan class="Chemical">water reconstructions (Fig. 2). The age model is from Schmidt et al.[59]. Planktonic foraminiferal Globigerinoides ruber (sensu stricto, white variety) δ18O, Mg/Ca, and δ11B data are from refs. [21,22,59]. Briefly, about 25 and 55 shells from the 250–350 μm size fraction were used for δ18O and Mg/Ca analyses, respectively. Samples for δ18O analyses were sonicated in methanol for 5–10 s, roasted under vacuum at 375 oC for 30 min, and analyzed on a Fisons Optima IRMS with a precision of <0.06‰. Shells for Mg/Ca were cleaned following the reductive cleaning procedure[60] and measured on an inductively-coupled plasma mass spectrometer (ICP–MS) with a precision of ~1.7%. For δ11B analyses, about 100–120 G. ruber (w) shells from the 300–355 μm size fraction were cleaned following the “Mg-cleaning” procedure[61], to minimize material loss during cleaning[62]. G. ruber (w) δ11B was measured on a Neptune multicollector (MC)–ICP–MS with an analytical error in δ11B of about ±0.25‰ (ref. [21]).

Three cores (BOFS 17, BOFS 11, and BOFS 14 K) from the polar North Atlantic Ocean are used for ass="Chemical">deep-ass="Chemical">pan class="Chemical">water reconstructions (Fig. 3). Their age models are based on published chronologies[24,63-65]. For each sample (~2 cm thickness), ~10–20 cm3 of sediment was disaggregated in de-ionized water and was wet sieved through 63 μm sieves. To facilitate analyses, we picked the most abundant species for measurements. For each B/Ca analysis, ~10–20 monospecific shells of the benthic foraminifera C. mundulus (BOFS 17 K) and C. wuellerstorfi (BOFS 14, 11 K) were obtained from 250 to 500 μm size fraction. The shells were double checked under a microscope before crushing to ensure that consistent morphologies were used throughout the core. On average, following this careful screening the starting material for each sample was ~8–12 shells, which is equivalent to ~300–600 μg of carbonate. For benthic B/Ca analyses, foraminiferal shells were cleaned with either the “Mg-cleaning” method[61] or the “Cd-cleaning” protocol[61], to investigate cleaning effects on trace element/Ca in foraminiferal shells[62,66]. No discernable B/Ca difference is observed between the two cleaning methods[25,62]. Benthic B/Ca ratios were measured on an ICP–MS using procedures outlined in ref. [67], with an analytical error better than ~5%.

For each benthic Cd/Ca analysis, ~10–20 shells of the benthic ass="Species">foraminiferal taxa C. mundulus (BOFS 17 K), ass="Chemical">pan class="Species">C. wuellerstorfi (BOFS 14 K, 11 K), and Uvigerina spp. (BOFS 17 K) were picked from the 250–500 μm size fraction. Previous studies[26,27,68] showed similar Cd/Ca ratios between infaunal Uvigerina spp. and epifaunal Cibicidoides, and we thus combined Cd/Ca data from these taxa to obtain continuous downcore PO4 records. We used the “Cd-cleaning” method[60,69] to clean benthic shells for Cd/Ca measurements. Cd/Ca ratios were measured on an ICP–MS with an analytical error better than ~5% (ref. [67])

For δ11B measurements, about 20 benthic shells from the 250–500 μm size fraction were picked for each sample. Shells used for δ11B analyses were cleaned using the “Mg-cleaning” method, to minimize loss of shell material[61]. After cleaning, shells were dissolved and pure ass="Chemical">boron was extracted using column chemistry as described by Foster[21]. Benthic δ11B was measured on a Neptune multi-collector (MC)–ICP–MS following ref. [21]. The analytical error in δ11B is about ± 0.25‰. Due to the relatively large sample size requirement, shell availability, and lengthy chemical treatments for δ11B, we present low-resolution δ11B for C. mundulus from BOFS 17 K and for ass="Chemical">pan class="Species">C. wuellerstorfi from BOFS 11 K. Note that consistent [CO32−] results from B/Ca and δ11B strengthen the reliability of our reconstructions (Fig. 3).

Published benthic Cd/Ca and B/Ca results are included in Fig. 3. Altogether, we generated 180 new measurements of benthic δ11B, B/Ca, and Cd/Ca. All data are listed in Supplementary Data 1–9.

ODP 999 reconstructions

ass="Chemical">ODP Site 999 was used to constrain ass="Chemical">past physical conditions and ass="Chemical">pan class="Chemical">carbonate chemistry of the Gulf Stream (Supplementary Fig. 4). Following previous approaches[21,22], surface water temperature (Tsurface) and salinity (Ssurface) were estimated based on G. ruber Mg/Ca (ref. [59]) and sea level changes[21,22,59], respectively. We first convert G. ruber δ11B to borate δ11B (δ11Bborate), following the conversion method of ref. [22]. Surface water pH (pHsurface) was calculated from seawater δ11Bborate along with Tsurface and Ssurface. To constrain the CO2 system, two CO2 system variables are necessary[16]. In addition to δ11B-derived pH, literature studies[21,22,41] generally estimate past surface-water ALK (ALKsurface) changes. Following refs. [21,22], we estimate ALKsurface from Ssurface using the modern Ssurface–ALKsurface relationship (ALKsurface = 59.19 × Ssurface + 229.08, R2 = 0.99)[21]. Together with Tsurface and Ssurface, pHsurface, and ALKsurface were used to calculate other CO2 system variables including surface-water [CO32−] ([CO32−]surface) and DIC (DICsurface) using the CO2sys program[28]. Surface-water PO4 concentration at ODP 999 is assumed to be zero over the last 27 ka.

Following refs. [21,22,59], errors are estimated to be 1 °C, 1‰, 100 μmol/kg, and ~0.43‰ for Tsurface, Ssurface, ALKsurface, and δ11Bborate, respectively. Integrated average uncertainties in [ass="Chemical">CO32−]surface and DICsurface for a single reconstruction are, respectively, ~20 (ass="Chemical">pan class="Chemical">Holocene: ~18, LGM: ~24) and ~90 μmol/kg, based on quadratic addition of all individual errors sourced from Tsurface ([CO32−]surface: 2 μmol/kg, DICsurface: 3 μmol/kg), Ssurface ([CO32−]surface: 2 μmol/kg, DICsurface: 5 μmol/kg), ALKsurface ([CO32−]surface: 14 μmol/kg, DICsurface: 86 μmol/kg), and δ11Bborate ([CO32−]surface: 16 μmol/kg, DICsurface: 24 μmol/kg; note that δ11Bborate leads to an error in [CO32−] via pH). Uncertainties for calculated CO2 system variables at ODP 999 are tabulated in Supplementary Data 1. Use of other methods to estimate ALK would have little impact on our conclusions (Supplementary Figs. 16 and 17).

From pH to [CO32−]

For palaeo-studies, surface-ass="Chemical">water pH is generally obtained from planktonic ass="Chemical">pan class="Species">foraminiferal δ11B. To calculate [CO32−], a second CO2 system variable is needed[16]. Following the previous approach[21,22], past ALKsurface at ODP 999 have been estimated from S using the Ssurface–ALKsurface relationship. Due to limited knowledge about the past Ssurface–ALKsurface relationship, a generous uncertainty has been assigned to ALKsurface at ±100 μmol/kg (ref. [21,22]), which is about half of the entire ALK range in the present global ocean[2]. Using ALKsurface and pHsurface along with Tsurface and Ssurface, [CO32−]surface and DICsurface can be calculated using the CO2sys program[28]. Because of the large uncertainty in ALKsurface, large errors in DICsurface might be expected (Supplementary Fig. 4). However, given the constraint from pHsurface, seawater ALKsurface and DICsurface variations are not random but must vary systematically within ALK–DIC space (Supplementary Fig. 5). Because of the close relationship between pH and [CO32−] (i.e., roughly parallel patterns of pH and [CO32−] within ALK−DIC space; Supplementary Fig. 5), this systematic ALK−DIC variation allows us to confine [CO32−] with acceptable uncertainty. For a given pH at ODP 999, an error of 100 μmol/kg in ALK only leads to an error of about ±14 μmol/kg in [CO32−] (Supplementary Fig. 5).

For clarity, Supplementary Fig. 5a, b only consider the effect of ALK errors on [ass="Chemical">CO32−] estimates assuming constant pH and T–S–P conditions. To fully proass="Chemical">pagate errors from various sources including Tsurface, Ssurface, ALKsurface, and pHsurface, we use a Monte Carlo approach (n = 10,000) to calculate the integrated error in [ass="Chemical">pan class="Chemical">CO32−] (ref. [70]). As can be seen from Supplementary Fig. 5c–f, the final errors (~20–25 μmol/kg) in an individual [CO32−] reconstruction based on the Monte-Carlo are similar to those (~18–24 μmol/kg) based on quadratic addition of individual errors, justifying our major error estimation approach (i.e., quadratic addition).

Subtropical western North Atlantic surface [CO32−]

Because most of North Atlantic subtropical gyre ass="Chemical">waters circulate through the Caribbean Sea before being transported to the subpolar North Atlantic via the Gulf Stream, ass="Chemical">pan class="Chemical">ODP 999 from Caribbean Sea is used to constrain past Gulf Stream carbonate chemistry[20]. To further test the feasibility of using ODP 999 to represent the first-order Gulf Stream [CO32−] changes during the Holocene and LGM, we have estimated surface-water [CO32−] for four sites from the wider subtropical western Atlantic region (latitude: 12–33°N, longitude: 61–91°W). Among these sites, KNR140–51GGC (33°N, 76°W) is located within the Gulf Stream today[71]. Because subtropical surface waters cycle multiple times through the upper ocean gyre circulations, it is possible that surface waters have been close to equilibrium with past atmospheric pCO2 (refs. [21,22]). Therefore, we assume surface-water pCO2 of 270 and 194 ppm for the Holocene and LGM, respectively[72]. We assign a ±15 ppm error to surface-water pCO2 to account for any potential air–sea CO2 disequilibrium. For these sites, we use surface temperature and salinity reconstructions from previous publications[71,73-75]. ALK is calculated based on the same approach for ODP 999. The reconstructed in situ [CO32−] values show some differences between cores, due to local T–S conditions. Since we are interested in air–sea CO2 exchange signals, we convert reconstructed in situ [CO32−] into [CO32−]Norm using Eq. (11). As can be seen from Supplementary Fig. 6 and Supplementary Data 2, these cores show similar [CO32−]Norm values for the Holocene (~260 μmol/kg) and LGM (~300 μmol/kg) as ODP 999. Therefore, we argue that ODP 999 sufficiently records first-order Gulf Stream air–sea exchange carbonate chemistry for the Holocene and LGM. Because we aim to obtain a proxy-based estimates, we use ODP 999 data for calculations in the main text.

Benthic B/Ca and δ11B to deep-water [CO32−]

Most ass="Chemical">deep-ass="Chemical">pan class="Chemical">water [CO32−] values are reconstructed using benthic B/Ca (refs. [25,47]) from [CO32−]downcore = [CO32−]PI + ΔB/Cadowncore–coretop/k, where [CO32−]PI is the preindustrial (PI) deep-water [CO32−] value estimated from the GLODAP dataset[2], ΔB/Cadowncore–coretop represents the deviation of B/Ca of down-core samples from the core-top value, and k is the B/Ca–[CO32−] sensitivity of C. wuellerstorfi (1.14 μmol/mol per μmol/kg) or C. mundulus (0.69 μmol/mol per μmol/kg)[25]. We use a reconstruction uncertainty of ±10 μmol/kg in [CO32−] based on global core-top calibration samples[25,76].

For cores BOFS 17 K and BOFS 11 K, new monospecific epifaunal benthic δ11B values were converted into ass="Chemical">deep-ass="Chemical">pan class="Chemical">water [CO32−] following the approach detailed in ref. [77]. Briefly, benthic δ11B is assumed to directly reflect deep-water borate δ11B, as suggested by previous core-top calibration work[78]. Deep-water pH is calculated using benthic δ11B along with Tdeep and Sdeep, similar to the approach to calculate surface-water pH at ODP 999 (refs. [21,22]). We assume constant ALK at the studied sites (2313 μmol/kg at BOFS 17 K and 2310 μmol/kg at BOFS 11 K) in the past. Following ref. [77], a generous error of 100 μmol/kg is assigned to ALK estimates. We then calculate deep-water [CO32−] from pH and ALK using the CO2sys program[28]. The integrated average uncertainty in deep-water [CO32−] is ~±10 μmol/kg, based on quadratic addition of individual errors of ~±2 μmol/kg sourced from Tdeep (±1 °C), ~±2 μmol/kg from Sdeep (±1‰), ~±5 μmol/kg from ALK (±100 μmol/kg), and ~±8 μmol/kg from δ11Bborate (~±0.25‰). As demonstrated by Supplementary Fig. 5, the large ALK error only contributes a small uncertainty to the final [CO32−] estimate. As shown in Fig. 3, benthic B/Ca and δ11B yield consistent deep-water [CO32−] reconstructions for the Holocene and LGM.

Benthic Cd/Ca to deep-water PO4

We follow the established approach[26,46,79] to convert benthic (ass="Species">C. wuellerstorfi, C. mundulus, and Uvigerina spp.) ass="Chemical">pan class="Species">foraminiferal Cd/Ca into deep-water Cd concentrations. Partition coefficients (DCd) are used to calculate deep water Cd from: Cd (nmol/kg) = [(Cd/Ca)foram/DCd] × 10. Bertram et al.[65] used empirical DCd values of 2.3, 2.2, and 2.7 for BOFS 17, 14, and 11 K, respectively. However, these DCd values would result in Holocene Cd of 0.3–0.4 nmol/kg, higher than the observed value of ~0.25 nmol/kg from modern hydrographic measurements (Supplementary Fig. 7)[80]. This offset may suggest higher DCd values for the North Atlantic Ocean, which has been acknowledged recently[81]. We thus adjust DCd (~25% increase) so that the calculated Holocene deep-water Cd concentrations match modern measurements. This adjustment is supported by consistent Cd reconstructions from this study and previous reconstructions based on Cd/Ca measurements for Hoeglundina elegans. Compared to Cibicidoides, DCd into H. elegans is far less variable[79]. As can be seen from Supplementary Fig. 8, for cores with similar benthic δ13C from similar water depths (i.e., bathed in similar water masses), our Cd reconstructions match favorably with those based on H. elegans measurements[82]. Deep water Cd is converted into PO4 using the relationship based on the latest North Atlantic Ocean measurements (Supplementary Fig. 7)[80]. Using older published Cd–PO4 relationships[26,83] only marginally affects our PO4 estimates.

Uncertainties associated with Cd and ass="Chemical">PO4 reconstructions are estimated as follows. Error for Cd is estimated using 2σCd = , where and (=5%) are errors for ass="Chemical">pan class="Chemical">DCd and Cd/Ca, respectively. Due to poorly defined uncertainty for DCd from the literature, we assume an error of 50%, and then compare our final errors with literature estimates to assess the appropriateness of our calculations. Seawater PO4 is calculated from Cd using: PO4 = , where 2σ and 2σ, respectively, represent 95% confidence errors associated with a and b (Supplementary Fig. 7b). The PO4 uncertainty was calculated from: , where  = ,  = , and  = . Our final errors on individual Cd and PO4 are ~0.12 nmol/kg (~55%) and ~0.5 μmol/kg (~50%), respectively. When compared with previously published uncertainties (~0.08 nmol/kg for Cd and ~0.17 μmol/kg for PO4)[46,68], our error estimates are possibly too generous. Here we use ~50% error to be conservative. We encourage future work to improve uncertainty estimates for the benthic Cd/Ca proxy.

The oceanic residence time of ass="Chemical">PO4 is ~100,000 years[84]. The ass="Chemical">pan class="Chemical">LGM deep ocean was possibly more reducing[85], which might have facilitated sediment organic matter preservation, and, thus, PO4 removal from the ocean. However, this effect might have been compensated by decreased organic burial on continental slopes due to shallower LGM sea levels[86,87]. Considering the short (~10,000 years) last deglacial[84], we assume that global PO4 and Cd reservoirs remained constant between the Holocene and LGM. Our reconstructions (Fig. 3) are consistent with high benthic δ13C and low benthic Cd/Ca at numerous glacial North Atlantic mid-depth sites[23,31,46,65,88,89].

Deep-water temperature and salinity estimates

ass="Chemical">Deep-ass="Chemical">pan class="Chemical">water temperature (Tdeep) is estimated from the ice volume corrected benthic δ18O (δ18OIVC) and the δ18O-temperature equation of Marchitto et al.[90] from Tdeep = 2.5 − (δ18OIVC − 2.8)/0.224, where δ18OIVC = δ18Obenthic − δ18Oglobal_sealevel. δ18Oglobal_sealevel was estimated from sea level curves[86,87] with a global δ18Oseawater−sea level scaling of 0.0085‰/m (ref. [91]). Deep-water salinity (Sdeep) is calculated by: Sdeep = Score_top + 1.11 × δ18Oglobal_sealevel, where Score_top is the modern Sdeep (35.06, 34.926, and 34.893 at BOFS 17 , 11 , and 14 K, respectively[2]) and the term 1.11 is the scaling term for a global S−δ18Oglobal_sealevel relationship[29,91]. We assume ±1 °C and ±1‰ uncertainties in Tdeep and Sdeep, respectively. Use of other methods to estimate Tdeep and Sdeep negligibly affects our conclusions, due to relatively weak sensitivities of [CO32−]Norm to T and S changes (Fig. 4).

Uncertainties and statistical analyses

Uncertainties associated with [ass="Chemical">CO32−] and ass="Chemical">pan class="Chemical">PO4 were evaluated using a Monte-Carlo approach[92,93]. Errors associated with the chronology (x-axis) and [CO32−] and PO4 reconstructions (y-axis) are considered during error propagation. Age errors are assumed to be ±3000 years for the three BOFS cores. Methods to calculate errors associated with individual [CO32−] and PO4 reconstructions (y-axis) are given above. All data points were sampled separately and randomly 5000 times within their chronological and [CO32−] or PO4 uncertainties and each iteration was then interpolated linearly. At each time step, the probability maximum and data distribution uncertainties of the 5000 iterations were assessed. Figure 3 shows probability maxima (bold curves) and ±95% (light gray; 2.5−97.5th percentile) probability intervals for the data distributions, including chronological and proxy uncertainties. For details, see refs. [92,93].

For a time period (e.g., ass="Chemical">Holocene) where multiple analyses are available, uncertainties are calculated following the method from ref. [94] by 2σ = , where n is the number of reconstructions and 2σ is the error associated with individual reconstruction. For [ass="Chemical">pan class="Chemical">CO32−] or [CO32−]Norm offsets between the Holocene and LGM, 2σ = , where and are 2σ of Holocene and LGM values, respectively. Other methods (e.g., weighted mean)[95] would give similar results.

When using Eq. (11) to calculate [ass="Chemical">CO32−]Norm, errors from various sensitivities are <1.5 μmol/kg (see Supplementary Data 8 for crosschecking). Because [ass="Chemical">pan class="Chemical">CO32−] is normalized to a constant condition (i.e., no error with final T–S–P), the error in [CO32−]Norm is largely sourced from [CO32−] reconstruction uncertainties. For surface water [CO32−]Norm calculations, T and S errors are already included in surface [CO32−] reconstructions. For calculations associated with deep waters, [CO32−]Norm errors are ~0.5, ~3.5, and ~0.1 μmol/kg from ± 1 °C in T, ±1‰ in S, and ±50 dbar in P, respectively. Therefore, these uncertainties (already included in error calculations) are relatively less important compared to the reconstruction error of ±10 μmol/kg for deep water [CO32−].