Our mission is to educate the largest number of citizens possible on the fraudulent science, politics, and business of the deceptive carbon dioxide based anthropogenic global warming hypothesis through rigorous scientific experimentation and analysis, and the false viability of “renewable” energy sources; thereby negating the green energy agenda of America’s elected leaders, leading to the removal of carbon dioxide from the U.S. EPA’s 2015 Endangerment Finding and the elimination of all carbon regulation, taxes, and trading worldwide.
A non-profit 501(c)(3) compliant educational & research institute
More efficient North Atlantic carbon pump during the Last Glacial Maximum
J. Yu, L. Menviel, Z. D. Jin, D. J. R. Thornalley, G. L. Foster, E. J. Rohling, I. N. McCave, J. F. McManus, Y. Dai, H. Ren, F. He, F. Zhang, P. J. Chen & A. P. Roberts
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 Nordic Seas and Arctic Ocean) is a major atmospheric CO2 sink, which has been mitigating anthropogenic atmospheric CO2 increases1. 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 invasion2,3. Despite its modest area, the North Atlantic Ocean accounts for at least ~30% of the global ocean CO2 uptake today and during preindustrial times1,4. Over longer timescales, large-scale oceanic carbon sequestration also occurred during Plio-Pleistocene glaciations5,6,7. This is commonly attributed to reduced glacial Southern Ocean CO2 outgassing6,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 pCO2 must occur via air–sea 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 stoichiometry3,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
Following the established method3 to account for within-ocean DIC redistributions by soft-tissue and CaCO3 cycling, we calculate preindustrial Atlantic DICas using the GLODAP dataset2 (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 release3,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 utilization12,13,14. NADW thus represents an efficient pathway for atmospheric CO2 sequestration6,15. Through global deep ocean circulation, CO2 absorbed in the North Atlantic is transported throughout the world ocean1,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 data2 used to generate the sections96. NADW North Atlantic deep water, AABW Antarctic bottom water, AAIW Antarctic intermediate water. See Methods for calculation details
No proxy exists to reconstruct past seawater DIC and ALK at acceptable precision for direct application, so we employ a linked 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 bicarbonate16. 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 constants16. 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 NADW source17; thus, comparing the [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 sites18,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–NADW [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 Stream20. 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, respectively21,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 ALK21,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 zero2,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
Three cores are used to reconstruct deep-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 overflows24. We employ benthic foraminiferal δ11B and B/Ca to reconstruct deep-water [CO32−] with an uncertainty of ~10 μmol/kg25. δ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-water [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 GLODAP2 (n = 55,399; blue) and a LGM output from LOVECLIM58 (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 data2. b Holocene and LGM data. Error bars, 2σ
A plot of [CO32−]Norm vs. 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 deep-waters at the BOFS sites had 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 [CO32−]as gradients, which involves frequent use of the CO2sys program28 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 ODP 999−BOFS [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:CCaCO3) 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 Atlantic30) 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 pump15. Inclusion of a whole ocean ALK inventory change6 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 NADW is fed directly by AAIW without surfacing at low latitudes17. Northward AAIW transport is thought to have been reduced substantially in the glacial Atlantic23,31,32,33 in the face of vigorous GNAIW production34. Assuming a constant total 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 waters23,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 ODP 999−BOFS [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 DICas gradient changes, which indicate air–sea 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 suggest32,34,35,37 vigorous glacial northern-sourced intermediate water production and subsequent transport to the remaining world ocean. Based on previous estimates35,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 ocean23,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), respectively35,36,38,39. See Methods for calculation details
Discussion
Previous work40,41,42 has tried to constrain air–sea CO2 exchange by reconstructing surface conditions. This requires reconstructions of the air–sea 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 exchange1. Sea ice cover43 possibly expanded, reducing glacial North Atlantic CO2 absorption. A larger LGM meridional surface temperature gradient43,44 would enhance the North Atlantic solubility pump13. Existing planktonic δ15N and Cd/Ca data40,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 GNAIW23,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 calculations47,48,49 which concern combined biological (i.e., within-ocean DIC redistribution) and air–sea exchange carbon changes (Fig. 1c), our total North Atlantic 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 pCO2 is controlled by both 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 variations50,51,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 Ocean6,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 CO2 system calculations, seawater carbonate system variables were calculated using the CO2sys.xls program28 with dissociation constants K1 and K2 according to Mehrbach et al.55 and KSO4
according to Dickson56. 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 values2.
Preindustrial Atlantic DICas and [CO32−]as
The GLODAP dataset2 is used to calculate preindustrial ocean CO2 system variables. Following the established method of Broecker and Peng3, we account for 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:
(1)
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 (PO4mo = 2.2 μmol/kg, 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 work3,30, we used C/PO4 = 127 to calculate DICas and [CO32−]as in Fig. 2. Using other C/PO4 values10 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 [CO32−]as, we first calculate [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 seawater CO2 system is nonlinear, there is currently no simple way to derive these sensitivities based on CO2 system equations16. We use GLODAP preindustrial data2 to calculate numerically [CO32−] sensitivities to various physiochemical parameters. Use of LGM outputs from the LOVECLIM model58 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:
ALKs=35=ALK×35/S,and
(2)DICs=35=DIC×35/S.
(3)
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: Sen_S=([CO32−]−[CO32−]S=35)/(S−35).
(4)
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: Sen_T=([CO32−]S=35,T=3∘C−[CO32−]S=35)/(3−T).
(5)
Regarding pressure effects, we calculate [CO32−]S=35, T=3 °C, P=2500 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: Sen_P=([CO32−]S=35,T=3 ∘C,P=2500dbar−[CO32−]S=35,T=3∘C)/(2500−P)×100.
(6)
To estimate the influence on [CO32−] from within-ocean ALK–DIC redistributions by biological processes, we assume a 0.1 μmol/kg increase in 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: ALKS=35+respiration=ALKs=35+ΔPO4×C/PO4÷R×2−ΔPO4×N/PO4.
(8)
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: [CO32−]Norm/PO4sensitivity=([CO32−]Norm+respiration−[CO32−]Norm)/ΔPO4.
(9)
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 CO2 invasion (i.e., Δ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:
(10)
Using sensitivities shown in Fig. 4, [CO32−]Norm can be calculated by: [CO32−]Norm=[CO32−]+(35−S)×Sen_S+(3−T)×Sen_T+(2500−P)/100×Sen_P.
(11)
Excel spreadsheets are provided in Supplementary Data 7–8 to calculate [CO32−]Norm and the biological curves shown in Fig. 5.
LGM–Holocene North Atlantic carbon budget
The total extra carbon increase (Δ∑CLGM−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 LGM23,31,32,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 ODP Site 999 for Gulf Stream surface-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 procedure60 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” procedure61, to minimize material loss during cleaning62. 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).
Cores, age models, samples, and analytical methods
We used ODP Site 999 for Gulf Stream surface-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 procedure60 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” procedure61, to minimize material loss during cleaning62. 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 deep-water reconstructions (Fig. 3). Their age models are based on published chronologies24,63,64,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” method61 or the “Cd-cleaning” protocol61, to investigate cleaning effects on trace element/Ca in foraminiferal shells62,66. No discernable B/Ca difference is observed between the two cleaning methods25,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 foraminiferal taxa C. mundulus (BOFS 17 K), C. wuellerstorfi (BOFS 14 K, 11 K), and Uvigerina spp. (BOFS 17 K) were picked from the 250–500 μm size fraction. Previous studies26,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” method60,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 material61. After cleaning, shells were dissolved and pure boron was extracted using column chemistry as described by Foster21. 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 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
ODP Site 999 was used to constrain past physical conditions and carbonate chemistry of the Gulf Stream (Supplementary Fig. 4). Following previous approaches21,22, surface water temperature (Tsurface) and salinity (Ssurface) were estimated based on G. ruber Mg/Ca (ref. 59) and sea level changes21,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 necessary16. In addition to δ11B-derived pH, literature studies21,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 program28. 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 [CO32−]surface and DICsurface for a single reconstruction are, respectively, ~20 (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-water pH is generally obtained from planktonic foraminiferal δ11B. To calculate [CO32−], a second CO2 system variable is needed16. Following the previous approach21,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 ocean2. Using ALKsurface and pHsurface along with Tsurface and Ssurface, [CO32−]surface and DICsurface can be calculated using the CO2sys program28. 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 [CO32−] estimates assuming constant pH and T–S–P conditions. To fully propagate 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 [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 waters circulate through the Caribbean Sea before being transported to the subpolar North Atlantic via the Gulf Stream, ODP 999 from Caribbean Sea is used to constrain past Gulf Stream carbonate chemistry20. 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 today71. 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, respectively72. 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 publications71,73,74,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 deep-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 dataset2, Δ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 samples25,76.
For cores BOFS 17 K and BOFS 11 K, new monospecific epifaunal benthic δ11B values were converted into deep-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 work78. 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 program28. 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 approach26,46,79 to convert benthic (C. wuellerstorfi, C. mundulus, and Uvigerina spp.) 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 recently81. 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 variable79. 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 measurements82. 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 relationships26,83 only marginally affects our PO4 estimates.
Uncertainties associated with Cd and PO4 reconstructions are estimated as follows. Error for Cd is estimated using 2σCd = (2σDCd)2+(2σCd/Ca)2−−−−−−−−−−−−−−−−−√, where 2σDCd and 2σCd/Ca (=5%) are errors for 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 = Cd−b±(2σb)a±(2σa), where 2σa and 2σb, respectively, represent 95% confidence errors associated with a and b (Supplementary Fig. 7b). The PO4 uncertainty was calculated from: 2σPO4=(∂PO4/∂a⋅2σa)2+(∂PO4/∂b⋅2σb)2+(∂PO4/∂Cd⋅2σCd)2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√, where ∂PO4/∂a = −(Cd−b)a2, ∂PO4/∂b = −1a, and ∂PO4/∂Cd = 1a.
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 PO4 is ~100,000 years84. The LGM deep ocean was possibly more reducing85, 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 levels86,87. Considering the short (~10,000 years) last deglacial84, 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 sites23,31,46,65,88,89.
Deep-water temperature and salinity estimates
Deep-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 curves86,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, respectively2) and the term 1.11 is the scaling term for a global S−δ18Oglobal_sealevel relationship29,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 [CO32−] and PO4 were evaluated using a Monte-Carlo approach92,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., Holocene) where multiple analyses are available, uncertainties are calculated following the method from ref. 94 by 2σ = [∑ni=1(2σi)2]/n−−−−−−−−−−−−√, where n is the number of reconstructions and 2σi is the error associated with individual reconstruction. For [CO32−] or [CO32−]Norm offsets between the Holocene and LGM, 2σ = (2σHolocene)2+(2σLGM)2−−−−−−−−−−−−−−−−−−−√, where 2σHolocene and 2σLGM 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 [CO32−]Norm, errors from various sensitivities are <1.5 μmol/kg (see Supplementary Data 8 for crosschecking). Because [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−].
Data availability
The data reported in the paper are presented in Supplementary Data.
Additional information
Journal peer review information:Nature Communications thanks Anne Willem, Omta, Shuang Zhang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
References
1.
Takahashi, T. et al. Climatological mean and decadal change in surface ocean pCO2, and net sea-air CO2 flux over the global oceans. Deep Sea Res Part II56, 554–577 (2009).
Gloor, M. et al. A first estimate of present and preindustrial air-sea CO2 flux patterns based on ocean interior carbon measurements and models. Geophys. Res. Lett. 30, 10-11–10-14 (2003).
Keeling, C. D. & Heimann, M. Meridional eddy diffusion model of the transport of atmospheric carbon dioxide: 2. Mean annual carbon cycle. J. Geophys. Res.91(D7), 7782–7796 (1986).
Toggweiler, J. R., Murnane, R., Carson, S., Gnanadesikan, A. & Sarmiento, J. L. Representation of the carbon cycle in box models and GCMs: 2. Organic pump. Glob. Biogeochem. Cycle17, 1027 (2003).
13.
Toggweiler, J. R., Gnanadesikan, A., Carson, S., Murnane, R. & Sarmiento, J. L. Representation of the carbon cycle in box models and GCMs: 1. Solubility pump. Glob. Biogeochem. Cycle17, 1026 (2003).
14.
Palter, J. B., Lozier, M. S. & Barber, R. T. The effect of advection on the nutrient reservoir in the North Atlantic subtropical gyre. Nature437, 687–692 (2005).
Talley, L. D. Closure of the global overturning circulation through the Indian, Pacific, and Southern Oceans: schematics and transports. Oceanography26, 80–97 (2013).
Foukal, N. P. & Lozier, M. S. No inter-gyre pathway for sea-surface temperature anomalies in the North Atlantic. Nat. Commun. 7, 11333 (2016).
20.
Johns, W. E., Townsend, T. L., Fratantoni, D. M. & Wilson, W. D. On the Atlantic inflow to the Caribbean Sea. Deep Sea Res. Part I: Oceanogr. Res. Pap.49, 211–243 (2002).
Foster, G. L. Seawater pH, pCO2 and [CO32-] variations in the Caribbean Sea over the last 130 kyr; a boron isotope and B/Ca study of planktic foraminifera. Earth Planet Sci. Lett.271, 254–266 (2008).
Henehan, M. J. et al. Calibration of the boron isotope proxy in the planktonic foraminifera Globigerinoides ruber for use in palaeo-CO2 reconstruction. Earth Planet Sci. Lett.364, 111–122 (2013).
Curry, W. B. & Oppo, D. Glacial water mass geometry and the distribution of d13C of SCO2 in the western Altantic Ocean. Paleoceanography20, PA1017 (2005).
Yu, J. M., Elderfield, H. & Piotrowski, A. Seawater carbonate ion-d13C systematics and application to glacial-interglacial North Atlantic ocean circulation. Earth Planet Sci. Lett.271, 209–220 (2008).
Pelletier, G., Lewis, E. & Wallace, D. A Calculator for the CO2 System in Seawater for Microsoft Excel/VBA, 1.0 ed. (Washington State Department of Ecology, Olympia; Brookhaven National Laboratory, Upton, 2005).
29.
Sarmiento, J. L. & Gruber, N. Ocean Biogeochemical Dynamics. (Princeton University Press, Princeton, 2006).
30.
Takahashi, T., Broecker, W. & Langer, S. Redfield ratio based on chemical data from isopycnal surfaces. J. Geophys. Res.90, 6907–6924 (1985).
Lynch-Stieglitz, J., van Geen, A. & Fairbanks, R. G. Interocean exchange of Glacial North Atlantic intermediate water: evidence from Subantarctic Cd/Ca and carbon isotope measurements. Paleoceanography11, 191–201 (1996).
Makou, M. C., Oppo, D. W. & Curry, W. B. South Atlantic intermediate water mass geometry for the last glacial maximum from foraminiferal Cd/Ca. Paleoceanography25, PA4101 (2010).
34.
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. Nature428, 834–837 (2004).
Keigwin, L. & Swift, S. A. Carbon isotope evidence for a northern source of deep water in the glacial western North Atlantic. Proc. Natl Acad. Sci. USA114, 2831–2835 (2017).
38.
Johnson, G. C. Quantifying Antarctic bottom water and north atlantic deep water volumes. J. Geophys. Res. Oceans113, C05027 (2008).
39.
Broecker, W. et al. How much deep water is formed in the Southern Ocean? J. Geophys. Reasearch103, 15833–15843 (1998).
Yu, J., Thornalley, D. J. R., Rae, J. & McCave, I. N. Calibration and application of B/Ca, Cd/Ca, and d11B in Neogloboquadrina pachyderma (sinistral) to constrain CO2 uptake in the subpolar North Atlantic during the last deglaciation. Paleoceanography28, 237–252 (2013).
Straub, M. et al. Nutrient conditions in the subpolar North Atlantic during the last glacial period reconstructed from foraminifera-bound nitrogen isotopes. Paleoceanography28, 79–90 (2013).
Marchitto, T. & Broecker, W. Deep water mass geometry in the glacial Atlantic Ocean: a review of constraints from the paleonutrient proxy Cd/Ca. Geochem. Geophys. Geosyst.7, Q12003 (2006).
47.
Yu, J. et al. Sequestration of carbon in the deep Atlantic during the last glaciation. Nat. Geosci.9, 319–324 (2016).
Hoogakker, B. A. A., Elderfield, H., Schmiedl, G., McCave, I. N. & Rickaby, R. E. M. Glacial-interglacial changes in bottom-water oxygen content on the Portuguese margin. Nat. Geosci.8, 40–43 (2015).
Jaccard, S. L., Galbraith, E. D., Martinez-Garcia, A. & Anderson, R. F. Covariation of deep Southern Ocean oxygenation and atmospheric CO2 through the last ice age. Nature530, 207–210 (2016).
Crichton, K. A., Bouttes, N., Roche, D. M., Chappellaz, J. & Krinner, G. Permafrost carbon as a missing link to explain CO2 changes during the last deglaciation. Nat. Geosci.9, 683–686 (2016).
Hain, M. P., Sigman, D. M. & Haug, G. H. Carbon dioxide effects of Antarctic stratification, North Atlantic Intermediate Water formation, and subantarctic nutrient drawdown during the last ice age: diagnosis and synthesis in a geochemical box model. Glob. Biogeochem. Cycle24, GB4023 (2010).
55.
Mehrbach, C., Culberso, C. H., Hawley, J. E. & Pytkowic, R. M. Measurement of apparent dissociation-constants of carbonic-acid in seawater at atmospheric-pressure. Limnol. Oceanogr.18, 897–907 (1973).
Lee, K. et al. The universal ratio of boron to chlorinity for the North Pacific and North Atlantic oceans. Geochim. Cosmochim. Acta74, 1801–1811 (2010).
Menviel, L. et al. Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study. Paleoceanography32, 2–17 (2017).
Schmidt, M. W., Spero, H. J. & Lea, D. W. Links between salinity variation in the Caribbean and North Atlantic thermohaline circulation. Nature428, 160–163 (2004).
Rosenthal, Y., Boyle, E. A. & Slowey, N. Temperature control on the incorporation of magnesium, strontium, fluorine, and cadmium into benthic foraminiferal shells from Little Bahama Bank: prospects for thermocline paleoceanography. Geochim. Cosmochim. Acta61, 3633–3643 (1997).
Barker, S., Greaves, M. & Elderfield, H. A study of cleaning procedures used for foraminiferal Mg/Ca paleothermometry. Geochem. Geophys. Geosyst. 4, 8407 (2003).
62.
Yu, J., Elderfield, H., Greaves, M. & Day, J. Preferential dissolution of benthic foraminiferal calcite during laboratory reductive cleaning. Geochem. Geophys. Geosyst.8, Q06016 (2007).
63.
Manighetti, B., McCave, I. N., Maslin, M. & Shackleton, N. J. Chronology for climate change: developing age models for the Biogeochemical Ocean Flux Study cores. Paleoceanography10, 513–525 (1995).
Barker, S., Kiefer, T. & Elderfield, H. Temporal changes in North Atlantic circulation constrained by planktonic foraminiferal shell weights. Paleoceanography19, PA3008 (2004).
65.
Bertram, C. J., Elderfield, H., Shackleton, N. J. & Macdonald, J. A. Cadmium/calcium and carbon-isotope reconstructions of the glacial northeast Atlantic Ocean. Paleoceanography10, 563–578 (1995).
Yu, J. M. & Elderfield, H. Mg/Ca in the benthic foraminifera Cibicidoides wuellerstorfi and Cibicidoides mundulus: temperature versus carbonate ion saturation. Earth Planet Sci. Lett.276, 129–139 (2008).
Yu, J. M., Day, J., Greaves, M. & Elderfield, H. Determination of multiple element/calcium ratios in foraminiferal calcite by quadrupole ICP-MS. Geochem Geophys Geosyst6, Q08P01 (2005).
Boyle, E. & Keigwin, L. D. Comparison of Atlantic and Pacific paleochemical records for the last 215,000 years: changes in deep ocean circulation and chemical inventories. Earth Planet Sci. Lett.76, 135–150 (1985).
Foster, G. L. & Sexton, P. F. Enhanced carbon dioxide outgassing from the eastern equatorial Atlantic during the last glacial. Geology42, 1003–1006 (2014).
Carlson, A. E. et al. Subtropical Atlantic salinity variability and Atlantic meridional circulation during the last deglaciation. Geology36, 991–994 (2008).
Ziegler, M., Nürnberg, D., Karas, C., Tiedemann, R. & Lourens, L. J. Persistent summer expansion of the Atlantic Warm Pool during glacial abrupt cold events. Nat. Geosci.1, 601 (2008).
Flower, B. P., Hastings, D. W., Hill, H. W. & Quinn, T. M. Phasing of deglacial warming and Laurentide Ice Sheet meltwater in the Gulf of Mexico. Geology32, 597–600 (2004).
Rühlemann, C., Mulitza, S., Müller, P. J., Wefer, G. & Zahn, R. Warming of the tropical Atlantic Ocean and slowdown of thermohaline circulation during the last deglaciation. Nature402, 511 (1999).
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).
Yu, J., Foster, G. L., Elderfield, H., Broecker, W. S. & Clark, E. An evaluation of benthic foraminiferal B/Ca and d11B for deep ocean carbonate ion and pH reconstructions. Earth Planet Sci. Lett.293, 114–120 (2010).
Rae, J. W. B., Foster, G. L., Schmidt, D. N. & Elliott, T. Boron isotopes and B/Ca in benthic foraminifera: proxies for the deep ocean carbonate system. Earth Planet Sci. Lett.302, 403–413 (2011).
Boyle, E. A., Labeyrie, L. & Duplessy, J. C. Calcitic foraminiferal data confirmed by cadmium in aragonitic Hoeglundina—application to the Last Glacial Maximum in the northern Indian Ocean. Paleoceanography10, 881–900 (1995).
Came, R. E., Oppo, D. W., Curry, W. B. & Lynch-Stieglitz, J. Deglacial variability in the surface return flow of the Atlantic meridional overturning circulation. Paleoceanography23, PA1217 (2008).
Jaccard, S. L. & Galbraith, E. D. Large climate-driven changes of oceanic oxygen concentrations during the last deglaciation. Nat. Geosci.5, 151–156 (2012).
Stanford, J. D. et al. Timing of meltwater pulse 1a and climate responses to meltwater injections. Paleoceanography21, PA4103 (2006).
87.
Lambeck, K., Rouby, H., Purcell, A., Sun, Y. Y. & Sambridge, M. Sea level and global ice volumes from the Last Glacial Maximum to the Holocene. P Natl. Acad. Sci. USA111, 15296–15303 (2014).
Rickaby, R. E. M. & Elderfield, H. Evidence from the high-latitude North Atlantic for variations in Antarctic Intermediate water flow during the last deglaciation. Geochem. Geophys. Geosyst.6, Q05001 (2005).
Marchitto, T. M. et al. Improved oxygen isotope temperature calibrations for cosmopolitan benthic foraminifera. Geochim. Cosmochim. Acta130, 1–11 (2014).
Waelbroeck, C. et al. Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records. Quat. Sci. Rev.21, 295–305 (2002).
This contribution is dedicated to W.S. Broecker. We thank R.F. Anderson and Daniel Sigman for the discussions. This work is supported by ARC Future Fellowship (FT140100993) and Discovery Projects (DP140101393 and DP190100894) and NSFC (41676026) to J.Y., DECRA (DE150100107) and Discovery (DP180100048) to L.M., NERC Advanced Fellowship to GLF (NE/D00876X/2), and Australian Laureate Fellowship (FL120100050) to E.J.R.
Author information
Affiliations
Research School of Earth Sciences, The Australian National University, Canberra, ACT, 2601, Australia
J. Yu
, E. J. Rohling
, Y. Dai
& A. P. Roberts
SKLLQG, Institute of Earth Environment, Chinese Academy of Sciences, Xi’an, 710061, China
J. Yu
, Z. D. Jin
& F. Zhang
Climate Change Research Centre, University of New South Wales, Sydney, NSW, 2052, Australia
L. Menviel
Open Studio for Oceanic-Continental Climate and Environment Changes, Qingdao National Laboratory for Marine Science and Technology, Qingdao, 266061, China
Z. D. Jin
Department of Geography, University College London, London, WC1E 6BT, UK
D. J. R. Thornalley
Ocean and Earth Science, University of Southampton, National Oceanography Centre, Southampton, SO14 3ZH, UK
G. L. Foster
& E. J. Rohling
Department of Earth Sciences, University of Cambridge, Cambridge, CB2 3EQ, UK
I. N. McCave
Lamont-Doherty Earth Observatory of Columbia University, 61 Route 9W/PO Box 1000, Palisades, NY, 10964-8000, USA
J. F. McManus
Department of Geosciences, National Taiwan University, Taipei, Taiwan
H. Ren
Center for Climatic Research, Nelson Institute for Environmental Studies, University of Wisconsin-Madison, Madison, WI, 53706, USA
F. He
College of Earth, Ocean, and Atmospheric Sciences, Oregon State University, Corvallis, OR, 97331, USA
F. He
CAS Center for Excellence in Quaternary Science and Global Change, Xi’an, 710061, China
F. Zhang
State Key Laboratory of Marine Geology, Tongji University, Shanghai, 200092, China
P. J. Chen
Contributions
J.Y. conceived the idea and wrote the paper. L.M. assisted with the model data used. Z.J./F.Z. picked the foram shells. E.J.R./Y.D. assisted with the statistics. G.L.F./J.Y. measured the boron isotopes. All authors (L.M., Z.J., D.T., G.F., E.R., N.M., J.M., Y.D., H.R., F.H., F.Z., P.C. and A.R.) contributed to improving the paper.
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/.
