The atmospheric bridge communicated the δ¹³C decline during the last deglaciation to the global upper ocean

During the early part of the last glacial termination (17.2–15 ka) and coincident with a ∼35 ppm rise in atmospheric CO₂, a sharp 0.3‰–0.4‰ decline in atmospheric δ¹³CO₂ occurred, potentially constraining the key processes that account for the early deglacial CO₂ rise. A comparable δ¹³C decline has also been documented in numerous marine proxy records from surface and thermocline-dwelling planktic foraminifera. The δ¹³C decline recorded in planktic foraminifera has previously been attributed to the release of respired carbon from the deep ocean that was subsequently transported within the upper ocean to sites where the signal was recorded (and then ultimately transferred to the atmosphere). Benthic δ¹³C records from the global upper ocean, including a new record presented here from the tropical Pacific, also document this distinct early deglacial δ¹³C decline. Here we present modeling evidence to show that rather than respired carbon from the deep ocean propagating directly to the upper ocean prior to reaching the atmosphere, the carbon would have first upwelled to the surface in the Southern Ocean where it would have entered the atmosphere. In this way the transmission of isotopically light carbon to the global upper ocean was analogous to the ongoing ocean invasion of fossil fuel CO₂. The model results suggest that thermocline waters throughout the ocean and 500–2000 m water depths were affected by this atmospheric bridge during the early deglaciation.

1 Introduction

Atmospheric CO₂ increased by 80–100 ppm between the Last Glacial Maximum (LGM) and the Holocene (Marcott et al., 2014; Monnin et al., 2001). During the initial ∼35 ppm rise in CO₂ between 17.2 and 15 ka, ice core records also document a 0.3 ‰ contemporaneous decline in atmospheric δ¹³C (Bauska et al., 2016; Schmitt et al., 2012) (Fig. 1a and b, interval highlighted in grey). Notably, this millennial-scale trend was punctuated by an interval of even more rapid change, with a 12 ppm CO₂ increase (Marcott et al., 2014) and a −0.2 ‰ decrease in δ¹³CO₂ (Bauska et al., 2016) occurring in an interval of just ∼200 years between 16.3 and 16.1 ka (Fig. 1a and b, interval highlighted in red). Hypotheses proposed to explain these observations include increased Southern Ocean ventilation (e.g., Skinner et al., 2010; Burke and Robinson, 2012), a poleward shift and/or enhanced Southern Hemisphere westerlies (Toggweiler et al., 2006; Anderson et al., 2009; Menviel et al., 2018), and reduced iron fertilization (Martínez-García et al., 2014; Lambert et al., 2021). However, the chain of events leading to the atmospheric changes and the location(s) where the isotope signal originated are not yet established.

Figure 1 (a) Ice core records of atmospheric CO₂ (Bereiter et al., 2015; Marcott et al., 2014). (b) δ¹³CO₂ records (Bauska et al., 2016; Schmitt et al., 2012). (c) WOA-18 Pacific zonal mean (120–160^∘ E) salinity; the magenta star marks the GeoB17402-2 site. (d) C. mundulus δ¹³C record for upper-intermediate-depth and mode waters in the western equatorial Pacific. The millennial- and centennial-scale events in these records are highlighted in grey and red, respectively.

Marine proxy records can provide further constraints on the possible mechanisms. For instance, during the early deglaciation, surface and thermocline-dwelling foraminifera around the global ocean also recorded a distinct δ¹³C drop (e.g., Hertzberg et al., 2016; Lund et al., 2019; Spero and Lea, 2002), an observation replicated by shallow benthic records from the tropical–subtropical Atlantic and Indian oceans (Lynch-Stieglitz et al., 2019; Romahn et al., 2014). These observations have been interpreted to reflect a spread of high-nutrient, low-δ¹³C waters originating in the Southern Ocean that were subsequently transported throughout the upper ocean via a so-called intermediate water teleconnection (Martínez-Botí et al., 2015; Pena et al., 2013; Spero and Lea, 2002). According to this hypothesis, formerly isolated carbon from deep waters was upwelled in the Southern Ocean (Anderson et al., 2009) in response to a breakdown of deep-ocean stratification (Basak et al., 2018). This carbon would have then been carried by Antarctic Intermediate Water (AAIW) and Southern Ocean Mode Water (SAMW) to low latitudes where it outgassed to the atmosphere in upwelling regions like the eastern equatorial Pacific (EEP). We term this scenario “bottom-up” transport because ¹³C-depleted carbon passes through the upper ocean globally and is recorded in marine proxy records there before entering the atmosphere (and being recorded in ice cores). The alternative scenario to explain the early deglacial decline in planktic (and shallow benthic) δ¹³C we term “top-down”. This recognizes the importance of air–sea exchange in rapidly (on the order of 1 year) and globally (e.g., Schmittner et al., 2013) conveying an isotopic signal from the atmosphere to the ocean surface, followed by propagation of the δ¹³C signal from surface to upper intermediate depths occurring on a multi-decadal to centennial timescale (Heimann and Maier-Reimer, 1996; Broecker et al., 1985; Eide et al., 2017). Although these timescales allow an atmospheric δ¹³C decline to be propagated throughout the upper ocean, this top-down effect has largely been overlooked in the interpretation of marine planktic and benthic δ¹³C records, at least until recently (Lynch-Stieglitz et al., 2019).

The implications of the top-down scenario differ significantly from those of the bottom-up scenario. First, negative δ¹³C excursions recorded in the upper ocean need not be associated with enhanced influx of nutrients (based on the notion that the extra nutrients came from a previously isolated deep-ocean reservoir along with isotopically depleted respired metabolic carbon). Second, a top-down scenario does not require a specific or even a single initial path of carbon into the atmosphere. Outgassing to the atmosphere could occur anywhere at the ocean surface, with a negative δ¹³C signal that then propagates globally through air–sea gas exchange – akin to the ongoing fossil fuel CO₂ emissions and the propagation of their isotopically depleted signal down through the ocean (Eide et al., 2017).

In this paper we take a two-pronged approach to help elucidate the more likely of these end-member scenarios. First, we present a new benthic δ¹³C record from the western equatorial Pacific (WEP) at a depth of 566 m that fills an important data gap from intermediate water depths in the Pacific basin. The site is located in the pathway of SAMW and AAIW to the upper tropical Pacific (Fig. 1c) and is also shallow enough to be sensitive to δ¹³CO₂ changes in the top-down scenario. Second, the early deglacial section of this record is interpreted with insights gained from analyzing a transient deglacial simulation conducted with the Earth system model LOVECLIM (Menviel et al., 2018). The specific LOVECLIM simulation we utilize starts with a scenario of excess respired carbon accumulated in a more stratified deep ocean with reduced ventilation rates. Although it is not clear if such a glacial carbon scenario is correct (Cliff et al., 2021; Stott et al., 2021), we can still make use of the ability of the model to simulate how the ocean communicates stored carbon and its isotopic composition to the atmosphere during deglaciation (the focus of this paper).

In the transient LOVECLIM simulation, sequestered respired carbon from the deep and intermediate waters is ventilated through the Southern Ocean, leading to a sharp decline in δ¹³CO₂, consistent with ice core records. We evaluate the two different δ¹³C transport scenarios by partitioning the simulated carbon pool and its stable isotope signature into a preformed (DIC_pref, being the carbon that is transported passively by ocean circulation) and a respired (DIC_soft, the accumulated respired carbon since the water parcel was last in contact with the atmosphere) component. Because the LOVECLIM transient experiment does not explicitly simulate either preformed or respired carbon as additional numerical tracers, the respired carbon is instead estimated by apparent oxygen utilization (AOU) – the difference between oxygen saturation and simulated [O₂] (see Sect. 2.4). If the top-down transport scenario was the mechanism responsible for the δ¹³C decline in marine proxy records from the upper 1000 m of depth, the preformed signal should dominate, while a regenerated signal would dominate in the bottom-up scenario. The carbon partitioning framework is not new – previous studies have used this framework to study the mechanisms that lead to lower glacial atmospheric CO₂ (Ito and Follows, 2005; Ödalen et al., 2018; Khatiwala et al., 2019) and processes that control δ¹³CO₂ and marine carbon isotope composition (Menviel et al., 2015; Schmittner et al., 2013). This diagnostic framework has also been applied to study the carbon cycle perturbation in response to a weaker Atlantic Meridional Overturning Circulation (AMOC) (Schmittner and Lund, 2015), albeit in experiments that were performed under constant pre-industrial conditions. However, what is new here is the application of a second Earth system model (cGENIE; Cao et al., 2009) to fully evaluate the AOU-based offline approach against an explicit respired organic matter δ¹³C tracer.

2 Methods

After describing the new foraminiferal δ¹³C record in Sect. 2.1, we summarize the LOVECLIM model and a published deglacial transient simulation in Sect. 2.2. We then summarize the cGENIE Earth system modeling framework and deglacial experiments in Sect. 2.3 before describing the δ¹³C tracer partitioning framework in Sect. 2.4.

2.1 Stable isotope analyses and age model for piston core GeoB17402-2

The WEP piston core GeoB17402-2 (8^∘ N, $\mathrm{126}{}^{\circ }{\mathrm{34}}^{\prime }$ E; 556 m of water depth) (Fig. 1c) was recovered from the expedition SO-228. Planktic foraminiferal samples for ¹⁴C age dating were picked from the greater than 250 µm size fraction of sediment samples and were typically between 2 and 5 mg. All new radiocarbon ages were measured at the University of California Irvine Accelerator laboratory. An age model (Fig. S1) was developed for this core with BChron using the Marine20 calibration curve (Heaton et al., 2020) without any further reservoir age correction.

For benthic foraminiferal δ¹⁸O and δ¹³C measurements approximately four to eight Cibicidoides mundulus (C. mundulus) were picked. These samples were cleaned by first cracking the tests open and then sonicating them in deionized water, after which they were dried at low temperature. The isotope measurements were conducted at the University of Southern California on a GV Instruments Isoprime mass spectrometer equipped with an autocarb device. An in-house calcite standard (Ultissima marble) was run in conjunction with foraminiferal samples to monitor analytical precision. The 1 standard deviation for standards measured during the study was less than 0.1 ‰ for both δ¹⁸O and δ¹³C. The stable isotope data are reported in per mil with respect to Vienna Pee Dee Belemnite (VPDB).

2.2 LOVECLIM deglacial transient simulation

The LOVECLIM model (Goosse et al., 2010) consists of a free-surface primitive equation ocean model ($\mathrm{3}{}^{\circ }\times \mathrm{3}{}^{\circ }$, 20 vertical levels), a dynamic–thermodynamic sea ice model, an atmospheric model based on quasi-geostrophic equations of motion (T21, three vertical levels), a land surface scheme, a dynamic global vegetation model (Brovkin et al., 1997), and a marine carbon cycle model (Menviel et al., 2015). To study the sensitivity of the carbon cycle to different changes in oceanic circulation, a series of transient simulations of the early part of the last deglaciation (19–15 ka) (Menviel et al., 2018) was performed by forcing LOVECLIM with changes in orbital parameters (Berger, 1978), changes in the freshwater surface balance and northern hemispheric ice-sheet geometry and albedo (Abe-Ouchi et al., 2007), and starting from an LGM simulation that best fit oceanic carbon isotopic (¹³C and ¹⁴C) records (Menviel et al., 2017).

The simulation we analyzed for this study is “LH1-SO-SHW” from Menviel et al. (2018). We briefly describe the applied forcing in this simulation: first, a freshwater flux of 0.07 Sv was added to the North Atlantic between 17.6 and 16.2 ka, resulting in an AMOC shutdown. Second, a salt flux was added to the Southern Ocean between 17.2 and 16.0 ka to enhance Antarctic Bottom Water (AABW) formation. Due to its relatively coarse resolution, the model could misrepresent the high-southern-latitude atmospheric or oceanic response to a weaker North Atlantic Deep Water (NADW). Enhanced AABW could have occurred due to a strengthening of the SH westerlies, changes in buoyancy forcing at the surface of the Southern Ocean, an opening of polynyas, or sub-grid processes. Finally, two stages of enhanced Southern Ocean westerlies are prescribed in the simulation at 17.2 and at 16.2 ka; this timing generally corresponds to Southern Ocean warming associated with two phases of NADW weakening during Heinrich Stadial 1 (Hodell et al., 2017). For further details about this experiment, see Menviel et al. (2018).

We chose to focus our analysis on this particular simulation because (1) recent ice core records also suggest enhanced SO westerly winds during Heinrich stadials (Buitzert et al., 2018); (2) LH1-SO-SHW matches some of the important observations (e.g., ice core record of atmospheric CO₂ and δ¹³CO₂) better than the other scenarios presented in Menviel et al. (2018); and (3) the stronger SO wind stress in LH1-SO-SHW leads to an increased transport of AAIW to lower latitudes, which could have impacted the intermediate depths of the global ocean.

2.3 cGENIE simulations

The cGENIE Earth system model is based on a 3-D frictional geostrophic ocean circulation component, plus dynamic and thermodynamic sea ice components, and is configured here at a 36×36 horizontal grid resolution with 16 vertical layers in the ocean. The configuration we employ here lacks a dynamical (GCM) atmosphere, with atmospheric transport fixed and provided via a 2-D energy–moisture balance model (Edwards and Marsh, 2005). The low-resolution ocean component and highly simplified atmospheric component make cGENIE much less computationally expensive to run than LOVECLIM and facilitate multiple sensitivity experiments run to (deep-ocean circulation) steady state to help partition and attribute carbon sources and pathways.

Ocean carbon storage analysis using the cGENIE model has previously utilized a range of preformed tracers, including those of phosphate (P_pref), dissolved inorganic carbon (DIC_pref), dissolved oxygen (O_2pref), and alkalinity (Ödalen et al., 2018). In the model, these are implemented by resetting the current value of the tracer at the ocean surface at each time step to the corresponding “full” tracer; for example, the value of DIC_pref is set to that of surface ocean DIC. Technically, an anomaly is applied to each preformed tracer at the ocean surface at each time step, equal to the difference between the current bulk tracer value and the preformed tracer value (as opposed to simply directly setting the values equal in the code). Because in the numerical scheme, all fluxes, including those induced by ocean circulation and any preformed tracer anomalies, are calculated simultaneously and only summed and applied to update the tracer concentration field at the very end of the model time step, preformed tracer concentrations at the ocean surface and at the end of the time step never exactly equal those of the bulk tracer. Thereafter, these tracers are carried conservatively by ocean circulation, with no loss or gain due to, for instance, organic matter remineralization in the ocean interior.

We expand the diagnostic tracer capabilities of cGENIE here and additionally add DIC_soft, which is the contribution to DIC from respired carbon. This is implemented as a tracer reset to zero at the ocean surface at each time step, but which is incremented by an amount of DIC equal to the remineralization of both particulate and dissolved organic matter and including organic carbon “reflected” (not preserved and buried) from the sediment surface. As for the preformed tracers, ocean circulation also acts on the distribution of DIC_soft in the model. Figure S2 in the Supplement illustrates how DIC is partitioned for the pre-industrial steady state of Cao et al. (2009). Note that we do not explicitly simulate DIC_carb (the contribution to DIC from dissolving CaCO₃, either in the water column or at the sediment surface) as a fourth tracer, but rather simply calculate it as the difference between DIC and DIC_pref+DIC_soft.

We also create a novel addition to the model – preformed and respired ¹³C (δ¹³C_pref and δ¹³C_soft, respectively). These are implemented as DIC_pref and DIC_soft, except for the concentrations of DI¹³C. (In cGENIE, isotopes are carried explicitly as concentrations with delta (δ) values only generated in conjunction with bulk concentrations for output and more convenient input.) Figure S3 illustrates how the δ¹³C signature of DIC is partitioned into explicitly simulated preformed and respired carbon components, with δ¹³C_carb (the contribution to δ¹³C of DIC from dissolved CaCO₃) again calculated by difference.

A full description of the cGENIE tracer scheme, together with δ¹³C tracer decomposition and attribution error analysis for both steady-state carbon cycling and under an idealized perturbation experiment, is available in the Supplement, with the pertinent insights summarized in the Results section.

Finally, we create a transient deglacial-like experiment using cGENIE to approximately mimic some of the key features of a changing climate and carbon cycle simulated by LOVECLIM. Although AOU-based errors in estimating the partitioning of respired vs. preformed δ¹³C are already addressed via the idealized cGENIE steady-state and transient experiments (the Supplement), decoupling in time of atmospheric CO₂ (and δ¹³C), surface climate, biological export, and the large-scale circulation of the ocean (especially the AMOC) across the deglacial transition may induce a more complex evolution of AOU-based error. We address this by then calculating how the AOU-based error changes in a deglacial-like cGENIE experiment. For this, we take a model configuration based on the idealized “glacial” boundary conditions of Rae et al. (2020) (including increased zonal planetary albedo at high Northern Hemisphere latitudes and the orbital configuration at 21 ka). Note that we did not attempt to achieve a glacial-like atmospheric CO₂ value for this spin-up; instead, we prescribed atmospheric CO₂=278 ppm and ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}=-\mathrm{6.5}$ ‰. The spin-up was run for 10 000 years. We then performed a deglacial transient simulation with time-varying salt–freshwater flux into the North Atlantic and the Southern Ocean as well as wind stress forcing over the Southern Ocean (Fig. S8 in the Supplement). We ran this experiment with all the diagnostic tracers described above.

2.4 Separating δ¹³C anomalies into the preformed (Δδ¹³C_pref) and respired (Δδ¹³C_soft) component

The published (Menviel et al., 2018) transient LOVECLIM model experiment that we analyze here does not include the numerical tracers required to explicitly attribute the sources of any given change in δ¹³C in the model ocean. We hence make approximations from AOU calculated in the model experiment but assess the errors inherent in this by means of a set of experiments using a second Earth system model – cGENIE (Cao et al., 2009). This approach is detailed as follows (and expanded upon further in the Supplement).

We assume the following carbon isotopic mass balance:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}\cdot \mathrm{DIC}=& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{pref}}\cdot {\mathrm{DIC}}_{\mathrm{pref}}+{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}\cdot {\mathrm{DIC}}_{\mathrm{soft}}\\ & +{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{carb}}\cdot {\mathrm{DIC}}_{\mathrm{carb}},\end{array}\end{array}$$

where DIC, DIC_pref, DIC_soft, and DIC_carb are the dissolved total inorganic carbon, the preformed and respired organic matter (“Csoft”), and dissolved (calcium) carbonate carbon pools, respectively. δ¹³C_pref, δ¹³C_soft, and δ¹³C_carb are the corresponding isotopic signatures (as ‰) that contribute to the δ¹³C signature of DIC, and it is changes in the δ¹³C of DIC that we assume foraminiferal records reflect.

Any given observed δ¹³C anomaly in the ocean can then be expressed as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=& \mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{pref}}\cdot {\text{DIC}}_{\text{pref}}/\mathrm{DIC})\\ & +\mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}\cdot {\text{DIC}}_{\text{soft}}/\mathrm{DIC})\\ & +\mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{carb}}\cdot {\text{DIC}}_{\text{carb}}/\mathrm{DIC}).\end{array}\end{array}$$

The terms on the right-hand side (RHS) represent the contribution of the preformed, respired, and dissolved (carbonate) components to the overall δ¹³C change, respectively. Since the contribution of CaCO₃ dissolution is small in the upper 1000 m (where GeoB17402-2 is located) in carbon cycle models (see also the Supplement) and since there is no ¹³C fractionation during CaCO₃ formation in the LOVECLIM model, the last term on the RHS can be neglected for the purpose of this study.

We use AOU to estimate respired carbon and its contribution to the δ¹³C changes: $\mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}\cdot {\mathrm{DIC}}_{\mathrm{soft}}/\mathrm{DIC})=\mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}\cdot \mathrm{AOU}\cdot R_{\text{c:-o2}}/\mathrm{DIC})$, where δ¹³C_soft is estimated by the δ¹³C of export particulate organic carbon (POC) in the overlying water column and $R_{\text{c:-o2}}=\mathrm{117}:-\mathrm{170}$.

This leads to

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=& \mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{pref}}\cdot {\text{DIC}}_{\text{pref}}/\mathrm{DIC})\\ & +\mathrm{\Delta }({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}\cdot \mathrm{AOU}\cdot R_{\text{c:-o2}}/\mathrm{DIC}).\end{array}\end{array}$$

The anomaly, defined as the difference between 15 and 17.2 ka, can be expanded as follows.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{15}\mathrm{ka}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{17.2}\mathrm{ka}}\\ & ={\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{pref}}^{\mathrm{15}\mathrm{ka}}\cdot {\mathrm{DIC}}_{\mathrm{pref}}^{\mathrm{15}\mathrm{ka}}/{\mathrm{DIC}}^{\mathrm{15}\mathrm{ka}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{pref}}^{\mathrm{17.2}\mathrm{ka}}\\ & \cdot {\mathrm{DIC}}_{\mathrm{pref}}^{\mathrm{17.2}\mathrm{ka}}/{\mathrm{DIC}}^{\mathrm{17.2}\mathrm{ka}}+{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}^{\mathrm{15}\mathrm{ka}}\cdot {\mathrm{AOU}}^{\mathrm{15}\mathrm{ka}}\\ & \cdot R_{\text{c:-o2}}/{\mathrm{DIC}}^{\mathrm{15}\mathrm{ka}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{soft}}^{\mathrm{17.2}\mathrm{ka}}\\ & \cdot {\mathrm{AOU}}^{\mathrm{17.2}\mathrm{ka}}\cdot R_{\text{c:-o2}}/{\mathrm{DIC}}^{\mathrm{17.2}\mathrm{ka}}\end{array}\end{array}$$

The AOU approach to estimate respired carbon content assumes that the oxygen content of surface waters always reaches equilibrium with the overlying atmosphere. However, studies have shown that this is not always the case, particularly for water masses formed in high latitudes (Bernardello et al., 2014; Ito et al., 2004; Khatiwala et al., 2019; Cliff et al., 2021). As a result, AOU likely overestimates respired carbon content in the deep ocean. Additional errors associated with the AOU approach may result from the nonlinear solubility of O₂ and respiration that does not involve O₂ consumption (i.e., through denitrification or sulfate reduction) (Shiller, 1981; Ito et al., 2004). However, to what extent these biases will affect the relative contribution of the preformed and respired carbon pool to δ¹³C anomalies in a carbon cycle perturbation event has not to our knowledge previously been evaluated. To address this, we performed a deglacial transient simulation with cGENIE (see Sect. 2.3) and then applied Eq. (4) to the output, with the results then compared with the values that are explicitly simulated by cGENIE. We also conducted a simplified (modern-configuration-based) analysis of steady-state and transient error terms (Figs. S2–S7), which we include in full in the Supplement and discuss briefly in the main text.

3 Results

The new GeoB17402-2 benthic δ¹³C record from the intermediate WEP documents a −0.3 ‰ to −0.4 ‰ decline during the early deglaciation (Fig. 1d). Although the foraminiferal δ¹³C proxy can be complicated by temperature and carbonate ion changes (Bemis et al., 2000; Schmittner et al., 2017) and may thus not solely reflect seawater DIC δ¹³C changes, core-top patterns of benthic foraminiferal δ¹³C are highly correlated with present-day seawater DIC δ¹³C (Schmittner et al., 2017). The apparent lag between the onset of decline in benthic δ¹³C at site GeoB17402-2 (Fig. 1d) and in δ¹³CO₂ appears to be due to the relatively large age model uncertainty below 154 cm in the GeoB17402-2 record (median age ∼16.2 years) of up to 1–2 kyr (2 SD) (Fig. S1). Despite this age uncertainty, the new benthic record from the tropical Pacific captures a similar δ¹³C decline as recorded from similar depth sites in the tropical–subtropical Atlantic and Indian oceans (Lynch-Stieglitz et al., 2014, 2019; Romahn et al., 2014).

Figure 2 Time series from the LOVECLIM transient experiment (Menviel et al., 2018). (a) Freshwater input into the North Atlantic and the Southern Ocean. (b) Southern Hemisphere westerly wind forcing. (c) Simulated NADW, AABW, AAIW, and NPIW maximum stream function in LOVECLIM; 21-year moving averages are shown for the maximum stream function to filter the high-frequency variability. (d) Ice core record of atmospheric CO₂ (Bereiter et al., 2015; Marcott et al., 2014) and LOVECLIM-simulated atmospheric CO₂. (e) The Taylor glacier δ¹³CO₂ record (Bauska et al., 2016) and LOVECLIM-simulated δ¹³CO₂.

To investigate whether the early deglacial δ¹³C decline observed at these sites in the upper ocean is dominated by the preformed or respired component, we carried out an in-depth carbon cycle analysis of the LOVECLIM transient simulation (Menviel et al., 2018). In response to the applied freshwater input to the North Atlantic (Fig. 2a), the AMOC significantly weakens from its glacial state (Fig. 2c). This has only a minor effect on the atmospheric CO₂ and δ¹³CO₂ (Fig. 2d and e). In contrast, enhanced ventilation of AABW and AAIW caused by a combined freshwater- (Fig. 2a) and wind-stress-driven (Fig. 2b) breakdown of stratification leads to an atmospheric CO₂ increase of ∼25 ppm and δ¹³CO₂ decline of −0.35 ‰ between 17.2 and 15 ka (Fig. 2d and e). This is a consequence of stronger upwelling, bringing ¹³C-depleted deep waters to the upper ocean with δ¹³C generally decreasing by 0.2‰–0.3‰ at most locations in the upper 1000 m (Fig. 3a, d, and g). In all sectors of the Southern Ocean below a depth of 400 m, δ¹³C increases by 0.1‰–0.2‰ due to stronger ventilation. Throughout the mid-depth North Atlantic, δ¹³C decreases by more than 0.3‰–0.4‰ due to the weakening AMOC (Fig. 3g). Finally, the LOVECLIM simulates a North Pacific deepwater mass when AMOC slows down (Menviel et al., 2014), and this leads to stronger ventilation and +0.3‰–0.4‰ Δδ¹³C in the North Pacific below depths of 1000 m (Fig. 3a).

Figure 3 Ocean basin zonal mean anomalies (15 ka minus 17.2 ka) as simulated in LOVECLIM. (a–c) Pacific zonal mean anomaly (160^∘ E–140^∘ W). The magenta star marks the GeoB17402-2 site. (d–f) Indian zonal mean anomaly (50–90^∘ E). (g–i) Atlantic zonal mean anomaly (60–10^∘ W). The magenta circles mark the 78GGC and 33GGC sites discussed in Sect. 4.3.

Decomposing the LOVECLIM Δδ¹³C signal into the Δδ¹³C_soft and Δδ¹³C_pref component, we find that the entire water column of the Southern Ocean is characterized by a strong positive Δδ¹³C_soft (indicting a loss of respired carbon) and a strong negative Δδ¹³C_pref (Fig. 3b, c, e, f, h, and g). In the rest of the global upper ocean below depths of 1000 m, Δδ¹³C_soft is negative but of a magnitude smaller than 0.1 ‰, whereas a 0.2‰–0.3‰ decrease in Δδ¹³C_pref accounts for most of the Δδ¹³C signal. In the deep Indo-Pacific, Δδ¹³C_soft and Δδ¹³C_pref show opposite signs, with the positive Δδ¹³C_soft dominating the net Δδ¹³C (Fig. 3a–f). In the deep North Atlantic, Δδ¹³C_soft and Δδ¹³C_pref are both negative (Fig. 3h and i), leading to the largest decrease in Δδ¹³C across the ocean basins (Fig. 3g).

Figure 4 Ocean basin zonal mean anomalies (15 ka minus 17.2 ka), but for the cGENIE deglacial transient simulation. Panels are organized as in Fig. 3.

For comparison, Fig. 4 shows the Δδ¹³C, Δδ¹³C_soft, and Δδ¹³C_pref response in a similar deglacial-like transient simulation conducted with cGENIE (see Sect. 2.3 and Fig. S8) in which the respired and preformed components are explicitly simulated. The Δδ¹³C patterns (Fig. 4a, d, and g) are qualitatively similar to that simulated by LOVECLIM (Fig. 3a, d, and g), although the magnitude of positive Δδ¹³C in the deep Pacific and negative Δδ¹³C in the deep North Atlantic is larger in cGENIE (compare Fig. 3a with Fig. 4a and Fig. 3g with Fig. 4g). cGENIE does not simulate any large positive Δδ¹³C_soft or negative Δδ¹³C_pref in the Southern Ocean above depths of 3000 m (Fig. 4), in contrast to the AOU-based results from LOVECLIM (Fig. 3). In the North Atlantic, the magnitudes of negative Δδ¹³C_soft and Δδ¹³C_pref are both larger in cGENIE compared to LOVECLIM.

Figure 5 cGENIE early deglacial transient AOU error analysis for Δδ¹³C_soft (a, c, e) and Δδ¹³C_pref (b, d, f). The anomalies are defined as 15 ka minus 17.2 ka. The errors are defined as AOU-based anomaly minus explicitly simulated anomaly.

To assess the potential errors associated with the AOU-based approach used to process the LOVECLIM output, we also calculated AOU-derived estimates of Δδ¹³C_soft and Δδ¹³C_pref for the cGENIE deglacial transient simulation (see Sect. 2.4). The results suggest that throughout the mid-depth North Atlantic, the AOU-based Δδ¹³C decomposition may introduce errors up to 0.3‰–0.4‰ under a weakening of the AMOC (Fig. 5). In the Southern Ocean (south of 40^∘ S), the AOU-based approach overestimates the magnitude of the positive Δδ¹³C_soft and negative Δδ¹³C_pref by 0.1‰–0.4‰ (Fig. 5); the largest errors occur in the Pacific sector. Based on these results from cGENIE, we suggest that the apparent Δδ¹³C_soft and Δδ¹³C_pref in the Southern Ocean shown in the LOVECLIM decomposition (Fig. 3) are largely overestimated. Nonetheless, both cGENIE and LOVECLIM (after correcting the errors estimated from the cGENIE deglacial transient experiment; see Fig. 5) show that the preformed component contributes −0.1‰ to −0.2 ‰ to the total Δδ¹³C signal in the upper 1000 m of the Southern Ocean. To the north of 40^∘ S in the upper 1000 m of the global upper ocean (except for the upper North Atlantic), the errors are relatively minor (generally much less than 0.1 ‰ in magnitude) and the AOU-based approach can provide a reasonably good estimate (Fig. 5, also Figs. S5 and S7).

Finally, we further evaluate the errors inherent in the AOU-based approach for the decomposition of the different contributions to the δ¹³C changes by means of a series of idealized steady-state and transient cGENIE experiments, as described in the Supplement. From this we find that errors in estimating δ¹³C_soft arise from both errors in AOU (themselves composed of errors due to assuming air–sea equilibrium and because O₂ solubility increases nonlinearly with decreasing temperature) and the assumption that the isotopic signature of carbon released by the remineralization of organic matter at any location in the ocean reflects that of carbon exported from the directly overlying ocean surface. The latter error turns out to be small in LOVECLIM as a consequence of its relatively small (3 ‰) simulated latitudinal variability in organic matter δ¹³C, leaving the better understood AOU-driven error to dominate the net uncertainty in reconstructing δ¹³C_soft. As a further consequence of this, under idealized transient changes in climate and ocean circulation in cGENIE (see the Supplement), the AOU-induced error in δ¹³C_soft is almost invariant throughout the uppermost ca. 500 m of the ocean, simply because the error in AOU itself is close to zero here. This confirms the conclusions drawn from tracer comparisons made in deglacial cGENIE experiments that at the depth of GeoB17402-2, the AOU-based approach is relatively robust.

4 Discussion

4.1 Atmospheric δ¹³C bridge

In the LOVECLIM model ¹³C-depleted carbon is ultimately sourced from the respired carbon that accumulated in the deep and intermediate waters during the glacial period as a consequence of the imposed weakened deepwater formation (Menviel et al., 2017). We show that in this scenario the isotopic signal is first transmitted to the atmosphere through strong outgassing in the Southern Ocean (Fig. 6). The atmosphere then transmits the δ¹³C signal to the rest of the global surface and subsurface ocean through air–sea gas exchange. An illustrative example is the simulated transient δ¹³C minimum event between 16.2 and 15.8 ka in LOVECLIM (Fig. 2c), which originates from the Southern Hemisphere and specifically from enhanced ventilation of AAIW (Fig. 2a). In the model, if the top-down scenario is true, the upper water masses away from the Southern Hemisphere would show a similar magnitude of δ¹³C_DIC changes as δ¹³CO₂. On the other hand, if the bottom-up scenario is true, a large negative δ¹³C anomaly (of respired nature) should first appear in the South Pacific subtropical gyre (STGSP), as the STGSP lies on the pathway between Southern Ocean water masses and those at lower latitudes. Then the signal would progressively spread to the tropics and finally reach the North Pacific. The negative δ¹³C anomaly may also be gradually diluted along its pathway from the South Pacific to the North Pacific. However, in the LOVECLIM simulation, there is no δ¹³C minimum in the upstream STGSP, while the atmosphere-like negative δ¹³C anomaly appears in the EEP thermocline, the North Pacific subtropical gyre (STGNP), and North Pacific Intermediate Water (NPIW) simultaneously (Fig. 7). In addition, the millennial-scale δ¹³C evolution in these upper-ocean water masses to the north of the Equator exhibits a pattern of change that is similar to that of the atmosphere (Fig. 7). The synchronized δ¹³C changes therefore point to the dominant role of atmospheric communication rather than time-progressive oceanic transport of a low δ¹³C signal in LOVECLIM.

Figure 6 Changes in the air–sea pCO₂ gradient (15 ka minus 17.2 ka) simulated by LOVECLIM.

Figure 7 LOVECLIM-simulated Δδ¹³C in thermocline the EEP (90–82^∘ W, 5^∘ S–5^∘ N; 77–105 m), South Pacific subtropical gyre (STGSP; 160^∘ E–100^∘ W, 40–22^∘ S; 187–400 m), North Pacific subtropical gyre (STGNP; 110^∘ E–140^∘ W, 22–40^∘ N; 187–400 m), and NPIW (167–170^∘ E, 54–57^∘ N; 660 m). The average of 23.8–20 ka (i.e., LGM) is used as a reference level for the Δδ¹³C calculations. The interval of decreasing δ¹³C is highlighted with a grey bar.

In the LOVECLIM simulation, both millennial- and centennial-scale δ¹³CO₂ declines are the result of enhanced deep-ocean and/or intermediate ocean ventilation originating in the Southern Ocean. Using the UVic Earth system model, Schmittner and Lund (2015) showed that a slowdown of AMOC alone is able to weaken the global biological pump and lead to light carbon accumulation in the upper ocean and the atmosphere, without explicitly prescribing any forcing in the Southern Ocean. Despite the different prescribed forcing, Δδ¹³C_pref also dominates the total Δδ¹³C in the upper 1000 m of the global ocean in the UVic experiment (see Fig. 6 in Schmittner and Lund, 2015). Taken together, simulations by all three models suggest that any process that lowers atmospheric δ¹³CO₂ would influence the global upper-ocean δ¹³C. In fact, the same phenomenon has been recurring since the beginning of the industrial era due to fossil fuel burning – known as the Suess effect (Eide et al., 2017). The top-down scenario is also compatible with the concept of a nutrient teleconnection between the Southern Ocean and low latitudes (Palter et al., 2010; Pasquier and Holzer, 2016; Sarmiento et al., 2004). Figure 8 illustrates that stronger upwelling brings excess nutrients to the surface of the Southern Ocean. Unused nutrients are then transported to low latitudes within the upper-ocean circulation (e.g., through mode waters and thermocline waters). However, a nutrient teleconnection does not, in itself, reflect an enhanced flux of ¹³C-depleted DIC from the deep ocean to low latitudes in a “tunnel-like” fashion (and bottom-up transport).

Figure 8 Ocean basin zonal mean PO₄ anomalies (15 ka minus 17.2 ka) as simulated in LOVECLIM.

In the following sections, we present two cases in which the LOVECLIM transient simulation successfully captures the early deglacial δ¹³C_DIC evolution recorded in marine proxies. The model-based Δδ¹³C partitioning then offers a unique opportunity to investigate the controlling mechanisms of the observed marine δ¹³C variability. We acknowledge that there are also places where models (in both LOVECLIM and cGENIE deglacial transient simulations) fail to simulate the observed δ¹³C trend between 17.2 and 15 ka. For instance, models simulate significant positive Δδ¹³C (above 0.4‰–0.5‰) (Figs. 3a and 4a) in the deep tropical North Pacific, whereas observations record no significant trend (Lund and Mix 1998; Stott et al., 2021). Models also simulate very small Δδ¹³C (∼0.1 ‰) in the deep tropical northern Indian Ocean (Figs. 3d and 4d), whereas proxy records document a distinct +0.3‰–0.4‰ trend (Waelbroeck et al., 2006; Sirocko, 2000). The model–data disagreement in the deep Indo-Pacific warrants future study.

Figure 9 (a) Modern sea surface nitrate concentration from the WOA-18 dataset. The sites of ODP 1238 and GGC17/JPC30 are marked as a purple and blue circle, respectively. (b) Neogloboquadrina dutertrei (N. dutertrei, a shallow thermocline species) δ¹³C data from ODP 1238 (Martínez-Botí et al., 2015) and GGC17/JPC30 (Zhao and Keigwin, 2018), as well as LOVECLIM-simulated δ¹³C of DIC at 100 m (average of 82–90^∘ W, 5^∘ S–5^∘ N). The N. dutertrei data are corrected by −0.5 ‰ to normalize to δ¹³C of DIC (Spero et al., 2003). The grey shaded bars highlight the time period we focus on in this study.

4.2 Revisiting EEP thermocline δ¹³C

Waters at EEP thermocline depths are thought to be connected to the deep ocean through AAIW from the south and NPIW from the north. The EEP is therefore a potential conduit for deep-ocean carbon release to the atmosphere. On the other hand, the EEP thermocline is also shallow enough to record an atmospheric δ¹³C signal, either directly through gas exchange at the surface or indirectly through a preformed signal acquired from other parts of the global surface ocean. We select two EEP thermocline δ¹³C records from different oceanographic settings (Fig. 9a): site GGC17/JPC30 is near the coast, featuring relatively low surface nutrients; site ODP1238 is located in the main upwelling zone, with relatively high surface nutrients. Previous studies suggest that the deglacial histories of deepwater influence at the two sites were also distinctively different. At site ODP1238, strengthened deglacial CO₂ outgassing inferred from boron isotope data has been interpreted to reflect respired carbon transported from the Southern Ocean (Martínez-Botí et al., 2015); at site GGC17/JPC30, wood-constrained constant surface reservoir ages over the last 20 kyr suggest that this site was not influenced by old respired carbon from high latitudes (Zhao and Keigwin, 2018). However, the early deglacial planktic δ¹³C records from the two sites show remarkably similar evolution, which is well captured by the LOVECLIM transient simulation (Fig. 9b). By comparing Fig. 3b and c, it is clear that the simulated δ¹³C anomaly in the EEP thermocline (∼100 m) is dominated by the preformed component. The modeling evidence indicates that even though the EEP is the largest CO₂ outgassing region (in terms of absolute ΔpCO₂; Fig. S9) under an enhanced Southern Ocean upwelling scenario, its thermocline δ¹³C is dominantly controlled by the top-down mechanism rather than the bottom-up mechanism as previously suggested (Martínez-Botí et al., 2015; Spero and Lea, 2002). The apparent conundrum can be explained by the fact that the air–sea balance of carbon isotopes is achieved through gross rather than net CO₂ exchange. Collectively, we make the case that in strong upwelling regions (e.g., the EEP) that are remotely connected to the deep ocean, thermocline δ¹³C is still subjected to strong atmospheric overprint.

4.3 How deep in the ocean can the negative Δδ¹³C_pref signal from the atmosphere penetrate during the early deglaciation?

We have shown that given the dominant control of the preformed δ¹³C component in the upper ocean, some interpretations of planktic δ¹³C records might need to be re-evaluated. Our simulations also reveal that an atmospheric influence can extend beyond thermocline depths into upper intermediate depths – consistent with what Lynch-Stieglitz et al. (2019) proposed. Below 1000 m, a Δδ¹³C_pref signal from the atmosphere may still exist, but it no longer dominates the total Δδ¹³C as Δδ¹³C_soft becomes increasingly important at depth. (The contribution of δ¹³C_carb also increases at depth and can exceed 10 % of the contribution of δ¹³C_soft; see Fig. S3.)

Figure 10 Observed δ¹³C anomaly of 78GGC and 33GGC from the mid-depth of the Brazil Margin at 27^∘ S (Lund et al., 2015) and the LOVECLIM-simulated δ¹³C anomaly at this location.

It has been suggested that deglacial δ¹³C variability in the waters above a depth of 2000 m in the Atlantic could be driven by air–sea exchange (Lynch-Stieglitz et al., 2019). However, mid-depth (1800–2100 m) benthic δ¹³C records from the Brazil Margin (∼27^∘ S) document a sharp decline of 0.4 ‰ at ∼18 ka (Lund et al., 2019), while atmospheric δ¹³CO₂ did not decrease until ∼17 ka (Bauska et al., 2016; Schmitt et al., 2012). Lund et al. (2019) argued that the lagging atmospheric δ¹³CO₂ decline seemed at odds with the idea that δ¹³C_pref contributed to the early benthic δ¹³C decrease at their site. The observed benthic δ¹³C trend between 20 and 15 ka at these Brazil Margin sites is well simulated by LOVECLIM (Fig. 10), allowing us to explore this question further. Before atmospheric δ¹³CO₂ starts to decline in LOVECLIM at ∼17.2 ka, changes in δ¹³C_DIC at a depth of ∼2000 m at the Brazil Margin are dominantly controlled by excess accumulation of respired carbon (indicated by highly negative Δδ¹³C_soft; Fig. S10b), itself a response to the weakened AMOC, while Δδ¹³C_pref is relatively small (Fig. S10c). This is consistent with what previous studies have suggested (Lacerra et al., 2017; Lund et al., 2019; Schmittner and Lund, 2015). Interestingly, LOVECLIM also reveals a strong negative Δδ¹³C_pref signal between 17.2 and 15 ka when atmospheric δ¹³CO₂ declines (Fig. 3i). However, a positive Δδ¹³C_soft (Fig. 3h) signal originating from a loss of respired carbon due to enhanced ventilation at those depths almost completely compensates for the negative Δδ¹³C_pref, which leads to virtually no net change in δ¹³C_DIC in the simulation (Fig. 3g), consistent with the proxy observations (Fig. 10). These results suggest that, between 17.2 and 15 ka, a negative preformed δ¹³C signal from the atmosphere needs to be considered when interpreting benthic δ¹³C records from the upper 2000 m of the South Atlantic. The complexity associated with interpreting marine δ¹³C records further underscores the urgent need to develop more robust means of estimating respired carbon accumulation or release from water masses.

5 Conclusions

A transient simulation conducted by the LOVECLIM Earth system model is used as a realization of plausible pathways of low-δ¹³C-signal transport under a prevailing deglacial scenario that involves Southern Ocean processes. By applying an AOU-based partitioning of carbon isotopic changes into preformed and respired components – a methodology that we scrutinize via a series of additional cGENIE Earth system model experiments – we show that ocean–atmosphere gas exchange likely dominates the negative δ¹³C anomalies documented in global planktic and intermediate benthic δ¹³C records between 17.2 and 15 ka. Numerical simulations further suggest that enhanced Southern Ocean upwelling can transfer δ¹³C signals from respired carbon in the deep ocean directly to the atmosphere. Consequently, atmospheric δ¹³CO₂ declines and this leaves its imprint on the rest of the global upper ocean through air–sea exchange. The preformed component dominates the upper 1000 m and could account for a 0.3‰–0.4‰ decline in marine δ¹³C records during the early deglaciation, whereas the respired component becomes increasingly important at greater depth. At the same time, the amount of upwelling in the Southern Ocean is a forcing imposed on the model rather than directly constrained. It is therefore possible that there were other sites where excess carbon was ventilated to the atmosphere during the deglaciation, which also would have affected δ¹³CO₂. Our findings imply that planktic and upper intermediate benthic δ¹³C records do not provide strong constraints on the site or the mechanisms through which CO₂ was released from the ocean to the atmosphere. Interpretations of early deglacial upper-intermediate-depth benthic δ¹³C records also need to take into account an atmospheric influence. Whereas in the model simulations the source of the atmospheric signal is a direct response to enhanced Southern Ocean upwelling, our results underscore the need to find a way to fingerprint the actual source(s) of ¹³C-depleted carbon that caused the atmospheric δ¹³CO₂ decline.