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

. During the early part of the last glacial termina-tion (17.2–15 ka) and coincident with a ∼ 35 ppm rise in atmospheric CO 2 , a sharp 0.3‰–0.4‰ decline in atmospheric δ 13 CO 2 occurred, potentially constraining the key processes that account for the early deglacial CO 2 rise. A comparable δ 13 C decline has also been documented in numerous marine proxy records from surface and thermocline-dwelling planktic foraminifera. The δ 13 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 atmo-sphere). Benthic δ 13 C records from the global upper ocean, including a new record presented here from the tropical Paciﬁc, also document this distinct early deglacial δ 13 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 ﬁrst 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 inva-sion of fossil fuel CO 2 . 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.


Introduction
Atmospheric CO 2 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 2 between 17.2 and 15 ka, ice core records also document a 0.3 ‰ contemporaneous decline in atmospheric δ 13 C (Bauska et al., 2016;Schmitt et al., 2012) (Fig. 1a and b, interval highlighted in grey). Notably, this millennialscale trend was punctuated by an interval of even more rapid change, with a 12 ppm CO 2 increase (Marcott et al., 2014) and a −0.2 ‰ decrease in δ 13 CO 2 (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.
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 δ 13 C drop  (Bereiter et al., 2015;Marcott et al., 2014). (b) δ 13 CO 2 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 δ 13 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.
(e.g., Hertzberg et al., 2016;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-δ 13 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 13 Cdepleted 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) δ 13 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 δ 13 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). Al-though these timescales allow an atmospheric δ 13 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 δ 13 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 δ 13 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 δ 13 C signal that then propagates globally through air-sea gas exchange -akin to the ongoing fossil fuel CO 2 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 δ 13 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 δ 13 CO 2 changes in the topdown scenario. Second, the early deglacial section of this record is interpreted with insights gained from analyzing a transient deglacial simulation conducted with the Earth sys-tem 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 δ 13 CO 2 , consistent with ice core records. We evaluate the two different δ 13 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 2 ] (see Sect. 2.4). If the top-down transport scenario was the mechanism responsible for the δ 13 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 bottomup scenario. The carbon partitioning framework is not newprevious studies have used this framework to study the mechanisms that lead to lower glacial atmospheric CO 2 (Ito and Follows, 2005;Ödalen et al., 2018;Khatiwala et al., 2019) and processes that control δ 13 CO 2 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) , 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 δ 13 C tracer.

Methods
After describing the new foraminiferal δ 13 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 δ 13 C tracer partitioning framework in Sect. 2.4.

Stable isotope analyses and age model for piston
core GeoB17402-2 The WEP piston core GeoB17402-2 (8 • N, 126 • 34 E; 556 m of water depth) ( Fig. 1c) was recovered from the expedition SO-228. Planktic foraminiferal samples for 14 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 δ 18 O and δ 13 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 δ 18 O and δ 13 C. The stable isotope data are reported in per mil with respect to Vienna Pee Dee Belemnite (VPDB).

LOVECLIM deglacial transient simulation
The LOVECLIM model (Goosse et al., 2010) consists of a free-surface primitive equation ocean model (3 • ×3 • , 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 ( 13 C and 14 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 2 and δ 13 CO 2 ) 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.

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 3 , 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 13 C (δ 13 C pref and δ 13 C soft , respectively). These are implemented as DIC pref and DIC soft , except for the concentrations of DI 13 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 δ 13 C signature of DIC is partitioned into explicitly simulated preformed and respired carbon components, with δ 13 C carb (the contribution to δ 13 C of DIC from dissolved CaCO 3 ) again calculated by difference.
A full description of the cGENIE tracer scheme, together with δ 13 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 δ 13 C are already addressed via the idealized cGENIE steady-state and transient experiments (the Supplement), decoupling in time of atmospheric CO 2 (and δ 13 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 2 value for this spinup; instead, we prescribed atmospheric CO 2 = 278 ppm and δ 13 CO 2 = −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 δ 13 C anomalies into the preformed (∆δ 13 C pref ) and respired (∆δ 13 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 δ 13 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: 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. δ 13 C pref , δ 13 C soft , and δ 13 C carb are the corresponding isotopic signatures (as ‰) that contribute to the δ 13 C signature of DIC, and it is changes in the δ 13 C of DIC that we assume foraminiferal records reflect. Any given observed δ 13 C anomaly in the ocean can then be expressed as follows: The terms on the right-hand side (RHS) represent the contribution of the preformed, respired, and dissolved (carbonate) components to the overall δ 13 C change, respectively. Since the contribution of CaCO 3 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 13 C fractionation during CaCO 3 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 δ 13 C changes: (δ 13 C soft · DIC soft /DIC) = (δ 13 C soft · AOU · R c:-o2 /DIC), where δ 13 C soft is estimated by the δ 13 C of export particulate organic carbon (POC) in the overlying water column and R c:-o2 = 117 : −170.
This leads to δ 13 C = (δ 13 C pref · DIC pref /DIC) The anomaly, defined as the difference between 15 and 17.2 ka, can be expanded as follows.
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 2 and respiration that does not involve O 2 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 δ 13 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 (modernconfiguration-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.

Results
The new GeoB17402-2 benthic δ 13 C record from the intermediate WEP documents a −0.3 ‰ to −0.4 ‰ decline during the early deglaciation (Fig. 1d). Although the foraminiferal δ 13 C proxy can be complicated by temperature and carbonate ion changes (Bemis et al., 2000; and may thus not solely reflect seawater DIC δ 13 C changes, core-top patterns of benthic foraminiferal δ 13 C are highly correlated with present-day seawater DIC δ 13 C . The apparent lag between the onset of decline in benthic δ 13 C at site GeoB17402-2 (Fig. 1d) and in δ 13 CO 2 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  To investigate whether the early deglacial δ 13 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 2 and δ 13 CO 2 (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 2 increase of ∼ 25 ppm and δ 13 CO 2 decline of −0.35 ‰ between 17.2 and 15 ka ( Fig. 2d and e). This is a consequence of stronger upwelling, bringing 13 C-depleted deep waters to the upper ocean with δ 13 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, δ 13 C increases by 0.1‰-0.2‰ due to stronger ventilation. Throughout the mid-depth North Atlantic, δ 13 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‰ δ 13 C in the North Pacific below depths of 1000 m (Fig. 3a).
Decomposing the LOVECLIM δ 13 C signal into the δ 13 C soft and δ 13 C pref component, we find that the entire water column of the Southern Ocean is characterized by a strong positive δ 13 C soft (indicting a loss of respired carbon) and a strong negative δ 13 C pref (Fig. 3b, c, e, f, h, and g). In the rest of the global upper ocean below depths of 1000 m, δ 13 C soft is negative but of a magnitude smaller than 0.1 ‰, whereas a 0.2‰-0.3‰ decrease in δ 13 C pref accounts for most of the δ 13 C signal. In the deep Indo-Pacific, δ 13 C soft and δ 13 C pref show opposite signs, with the positive δ 13 C soft dominating the net δ 13 C (Fig. 3a-f). In the deep North Atlantic, δ 13 C soft and δ 13 C pref are both negative ( Fig. 3h and i), leading to the largest decrease in δ 13 C across the ocean basins (Fig. 3g).
For comparison, Fig. 4 shows the δ 13 C, δ 13 C soft , and δ 13 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 δ 13 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 δ 13 C in the deep Pacific and negative δ 13 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 δ 13 C soft or negative δ 13 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 δ 13 C soft and δ 13 C pref are both larger in cGENIE compared to LOVECLIM.
To assess the potential errors associated with the AOUbased approach used to process the LOVECLIM output, we also calculated AOU-derived estimates of δ 13 C soft and δ 13 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 δ 13 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 δ 13 C soft and negative δ 13 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 δ 13 C soft and δ 13 C pref in the Southern Ocean shown in the LOVECLIM decomposition (Fig. 3)  Finally, we further evaluate the errors inherent in the AOUbased approach for the decomposition of the different contributions to the δ 13 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 δ 13 C soft arise from both errors in AOU (themselves composed of errors due to assuming air-sea equilibrium and because O 2 solubility increases nonlinearly with decreasing temperature) and the assumption that the isotopic signature of carbon released by the remineralization of organic mat-ter 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 δ 13 C, leaving the better understood AOUdriven error to dominate the net uncertainty in reconstructing δ 13 C soft . As a further consequence of this, under idealized transient changes in climate and ocean circulation in cGE-NIE (see the Supplement), the AOU-induced error in δ 13 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.

Atmospheric δ 13 C bridge
In the LOVECLIM model 13 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 δ 13 C signal to the rest of the global surface and subsurface ocean through air-sea gas exchange. An illustrative example is the simulated transient δ 13 C minimum event between 16.2 and 15.8 ka in LOVE-CLIM (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 δ 13 C DIC changes as δ 13 CO 2 . On the other hand, if the bottom-up scenario is true, a large negative δ 13 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 δ 13 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 δ 13 C minimum in the upstream STGSP, while the atmosphere-like negative δ 13 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 δ 13 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 δ 13 C changes therefore point to the dominant role of atmospheric communication rather than time-progressive oceanic transport of a low δ 13 C signal in LOVECLIM.
In the LOVECLIM simulation, both millennial-and centennial-scale δ 13 CO 2 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, δ 13 C pref also dominates the total δ 13 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 δ 13 CO 2 would influence the global upperocean δ 13 C. In fact, the same phenomenon has been recur-  ring 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 South-ern 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 13 C-depleted DIC from the deep ocean to low latitudes in a "tunnel-like" fashion (and bottom-up transport). In the following sections, we present two cases in which the LOVECLIM transient simulation successfully captures the early deglacial δ 13 C DIC evolution recorded in marine proxies. The model-based δ 13 C partitioning then offers a unique opportunity to investigate the controlling mechanisms of the observed marine δ 13 C variability. We acknowledge that there are also places where models (in both LOVECLIM and cGENIE deglacial transient simulations) fail to simulate the observed δ 13 C trend between 17.2 and 15 ka. For instance, models simulate significant positive δ 13 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 δ 13 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.

Revisiting EEP thermocline δ 13 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 δ 13 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 δ 13 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 2 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 δ 13 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 δ 13 C anomaly in the EEP thermocline (∼ 100 m) is dominated by the preformed compo- dutertrei, a shallow thermocline species) δ 13 C data from ODP 1238 (Martínez-Botí et al., 2015) and GGC17/JPC30 (Zhao and Keigwin, 2018), as well as LOVECLIM-simulated δ 13 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 δ 13 C of DIC (Spero et al., 2003). The grey shaded bars highlight the time period we focus on in this study.
nent. The modeling evidence indicates that even though the EEP is the largest CO 2 outgassing region (in terms of absolute pCO 2 ; Fig. S9) under an enhanced Southern Ocean upwelling scenario, its thermocline δ 13 C is dominantly controlled by the top-down mechanism rather than the bottomup 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 2 exchange. Collectively, we make the case that in strong upwelling regions (e.g., the EEP) that are remotely connected to the deep ocean, thermocline δ 13 C is still subjected to strong atmospheric overprint.
4.3 How deep in the ocean can the negative ∆δ 13 C pref signal from the atmosphere penetrate during the early deglaciation?
We have shown that given the dominant control of the preformed δ 13 C component in the upper ocean, some interpretations of planktic δ 13 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 δ 13 C pref signal from the atmosphere may still exist, but it no longer dominates the total δ 13 C as δ 13 C soft becomes increasingly important at depth. (The contribution of δ 13 C carb also increases at depth and can exceed 10 % of the contribution of δ 13 C soft ; see Fig. S3.) It has been suggested that deglacial δ 13 C variability in the waters above a depth of 2000 m in the Atlantic could However, mid-depth (1800-2100 m) benthic δ 13 C records from the Brazil Margin (∼ 27 • S) document a sharp decline of 0.4 ‰ at ∼ 18 ka , while atmospheric δ 13 CO 2 did not decrease until ∼ 17 ka (Bauska et al., 2016;Schmitt et al., 2012).  argued that the lagging atmospheric δ 13 CO 2 decline seemed at odds with the idea that δ 13 C pref contributed to the early benthic δ 13 C decrease at their site. The observed benthic δ 13 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 δ 13 CO 2 starts to decline in LOVECLIM at ∼ 17.2 ka, changes in δ 13 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 δ 13 C soft ; Fig. S10b), itself a response to the weakened AMOC, while δ 13 C pref is relatively small (Fig. S10c). This is consistent with what previous studies have suggested (Lacerra et al., 2017;Schmittner and Lund, 2015). Interestingly, LOVECLIM also reveals a strong negative δ 13 C pref signal between 17.2 and 15 ka when atmospheric δ 13 CO 2 declines (Fig. 3i). However, a positive δ 13 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 δ 13 C pref , which leads to virtually no net change in δ 13 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 δ 13 C signal from the atmosphere needs to be considered when interpreting benthic δ 13 C records from the upper 2000 m of the South Atlantic. The complexity associated with interpreting marine δ 13 C records further underscores the urgent need to develop more robust means of estimating respired carbon accumulation or release from water masses.

Conclusions
A transient simulation conducted by the LOVECLIM Earth system model is used as a realization of plausible pathways of low-δ 13 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 δ 13 C anomalies documented in global planktic and intermediate benthic δ 13 C records between 17.2 and 15 ka. Numerical simulations further suggest that enhanced Southern Ocean upwelling can transfer δ 13 C signals from respired carbon in the deep ocean directly to the atmosphere. Consequently, atmospheric δ 13 CO 2 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 δ 13 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 δ 13 CO 2 . Our findings imply that planktic and upper intermediate benthic δ 13 C records do not provide strong constraints on the site or the mechanisms through which CO 2 was released from the ocean to the atmosphere. Interpretations of early deglacial upper-intermediate-depth benthic δ 13 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 13 C-depleted carbon that caused the atmospheric δ 13 CO 2 decline.
Code and data availability. The stable isotope and radiocarbon data are archived at the National Climatic Data Center -NOAA: https://doi.org/10.25921/3S96-9C26 (Shao et al., 2021). All modeling data generated or analyzed during this study can be made available upon request to the corresponding author (Jun Shao).
The code for the version of the "muffin" release of the cGENIE Earth system model used in this paper is tagged as v0.9.24 and is assigned a DOI: https://doi.org/10.5281/zenodo.4903423 (Ridgwell et al., 2021b).
Configuration files for the specific experiments presented in the paper can be found in the following directory: genieuserconfigs/MS/shaoetal.2021. Details of the experiments, plus the command line needed to run each one, are given in the readme.txt file in that directory. All other configuration files and boundary conditions are provided as part of the code release.
A manual detailing code installation, basic model configuration, and tutorials covering various aspects of model configuration and experimental design, plus results output and processing, is assigned a DOI: https://doi.org/10.5281/zenodo.4903426 (Ridgwell et al., 2021a).
Author contributions. JS designed the research with input from LS. LM provided the LOVECLIM output. AR implemented the new diagnostic tracers in cGENIE. JS performed the cGENIE simulations with help from AR. JS, LM, and MÖ analyzed the model simulations. MM was the chief scientist of the SO-228 expedition and provided samples from the GeoB17402-2 core. JS wrote the paper with contributions from all co-authors. AR wrote the text for the Supplement.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The funding for the expedition SO-228 is from the BMBF (German Ministry of Education and Research). We thank Fortunat Joos and another anonymous reviewer for their constructive comments and prompting that helped to significantly improve this paper.
Financial support. This research has been supported by the National Science Foundation Division of Ocean Sciences (grant nos. 1558990 to Jun Shao and Lowell D. Stott and 1736771 to Andy Ridgwell), the Australian Research Council (grant nos. FT180100606 and DP180100048), and the German Ministry of Education and Research (grant no. 03G0228A).
Review statement. This paper was edited by Hubertus Fischer and reviewed by Fortunat Joos and one anonymous referee.