Projected reversal of oceanic stable carbon isotope ratio depth gradient with continued anthropogenic carbon emissions

Paleoceanographic records suggest that the present-day vertical gradient in the stable carbon isotopic composition (δ13C) of dissolved inorganic carbon in the ocean was reversed during the Paleocene-Eocene Thermal Maximum, an early period of relatively rapid release of carbon into the climate system. Here we present simulations from an observationally constrained ocean model under various greenhouse gas emissions scenarios. We project a decrease in the globally averaged δ13C of dissolved inorganic carbon in the surface ocean of between −1.8 to −6.3 ‰ by 2100. This reduction is driven by oceanic absorption of anthropogenic carbon dioxide, which is depleted in carbon-13. Our findings suggest an elimination or reversal of the natural vertical gradient in the δ13C of dissolved inorganic carbon by 2100 unless anthropogenic carbon emissions are reduced soon. We conclude that the Paleocene-Eocene Thermal Maximum is a geologic analogue of future global carbon cycle perturbations under continued rapid anthropogenic carbon emissions. Uptake of 13C-depleted anthropogenic carbon dioxide in the surface ocean may eliminate or reverse the natural vertical gradient in the isotopic composition of dissolved inorganic carbon in the ocean, analogous to events at the Paleocene-Eocene Thermal Maximum, suggest ocean model simulations.

A nthropogenic carbon (C) emissions have led to an accumulation of 13 C-depleted carbon in the atmosphere and the upper ocean, diluting the isotopic composition of carbon dioxide (CO 2 ) and dissolved inorganic carbon (DIC) (referred to as the 13 C Suess effect) 1,2 . The oceanic absorption of anthropogenic CO 2 has disturbed naturally formed vertical gradients of the δ 13 C of DIC (δ 13 C-DIC) where δ 13 C is defined as [( 13 C/ 12 C) sample /( 13 C/ 12 C) standard −1] with the Vienna Pee Dee Belemnite standard. Surface seawater δ 13 C-DIC values are typically more positive relative to deep waters 3,4 , due primarily to the preferential downward transport of 12 C relative to 13 C by marine biota 5 . Vertical gradients of δ 13 C-DIC result from the interplay between the oceanic biological carbon pump, air-sea CO 2 exchange rates, thermodynamic fractionations, and ocean circulation as well as rapid invasion of 13 C depleted atmospheric CO 2 (refs. 3,[6][7][8][9][10][11]. Because all of these controlling processes are sensitive to climate and environmental changes, the δ 13 C-DIC values recorded in marine carbonate fossils have been used to gain valuable insights into past changes of the global carbon cycle and ocean circulation 4,8,10,[12][13][14][15] . However, the preindustrial δ 13 C-DIC distribution and its anthropogenic evolutions are uncertain 1,16 , due to the difficulty of accounting for the imprint of 13 C-depleted emissions from fossil fuel burning and deforestation. These uncertainties in turn hinder projections of future oceanic δ 13 C-DIC, which can provide a long-term perspective of natural and human-induced changes in the global carbon cycle 17 . The best known geologic analogue to the future perturbation of the global carbon cycle is the Paleocene-Eocene Thermal Maximum (PETM;~56 million years ago) 17 , which was characterized by carbon isotope excursions of −(3-5)‰ on the Earth's surface 9,12,13 . These anomalies were accompanied by rising temperatures and widespread ocean acidification, and have been attributed to a rapid carbon release of~10,000 GtC to the climate system 12,13,18 . Foraminifera fossil δ 13 C records also suggest that the pre-PETM depth gradients might have been eliminated or even reversed during the PETM onset 9,13 . However, because the carbon emission rates for the PETM onset were estimated to be an order of magnitude slower than those during industrialization 8,19,20 , the PETM deep ocean might have had sufficient time to fully equilibrate to changing atmospheric CO 2 without necessarily invoking a vertical gradient reversal 8 . Nonetheless, several studies 8,11,21 suggested that relatively rapid carbon emissions could have delayed deep ocean δ 13 C-DIC excursions compared to the surface, causing a temporal reversal of the vertical δ 13 C-DIC gradient. Such sensitivity indicates that the early geochemical evolution of the vertical δ 13 C-DIC gradient can provide an important constraint on the PETM onset duration, which has previously been suggested to be between 3 and 20 kyr 8,19,20 .
Here, we use an observationally constrained global oceanic carbon cycle model 22,23 and four representative concentration pathways (RCPs) 24 to estimate the temporally evolving oceanic δ 13 C-DIC from the preindustrial era to the 21st century (Methods). We find that the 21st century greenhouse gas emissions will induce elimination or even reversal of the natural vertical δ 13 C-DIC gradient by the end of the 21st century and elucidate the geochemical mechanisms underlying such radical changes in surface δ 13 C-DIC. Our 21st century projections combined with multi-millennial simulations for hypothetical futures now allow us to compare the δ 13 C-DIC excursions between the future and the PETM onset, which will benefit from a recent compilation 25 of the PETM foraminifera δ 13 C records.

Results and discussion
Modern ocean 13 C Suess effect. The oceanic δ 13 C-DIC changes from the preindustrial era to 2018 are assessed using a Monte Carlo experiment (Methods), for which uncertainties in various model parameters (Supplementary Table 1) are allowed to propagate in our model into the δ 13 C-DIC simulations. According to the simulations, the globally averaged surface ocean δ 13 C-DIC decline relative to 1780 CE pre-industrial conditions amounts to -(1.1 ± 0.2) ‰ as of 2018 (Fig. 1a), with 84% of the uncertainty arising from uncertainties in air-sea CO 2 exchange rates and the atmospheric 13 C Suess effect (Supplementary Fig. 1a; Supplementary Note 2). A dominance of air-sea CO 2 exchange rates for the uncertainty is in line with a previous study 26 based on an Earth System Model which suggested that air-sea gas exchange rates alone explain 63% of the total variance of simulated oceanic 13 C Suess effects over a time period of 1858-2008.
Our surface estimate is consistent with a previous model-based estimate of Schmittner et al. 3 who suggested a surface-averaged change of −0.67‰ as of the 1990's and an observation-based estimate of −0.76 ± 0.12‰ as of 1992 by Sonnerup et al. 27 . On the other hand, our estimated 13 C Suess effect is more negative than that of Eide et al. 1 who suggested a globally averaged 13 C-depletion of only 0.4‰ at 200 m as of 1994 (Fig. 2c, e); yet would be consistent if the previous estimate is corrected for the uncertainty of 0.15-0.24‰ suggested by the same study 1 and a later independent study 28 . The rate of upper ocean δ 13 C-DIC decline has been accelerating since 1960, keeping up with the atmospheric δ 13 C-CO 2 change (Fig. 1a), which is also consistent with global ocean coral records 29 . Compared to the surface ocean 13 C Suess effect, below a depth of 74 m of the water column (referred to as subsurface hereinafter) the averaged 13 C Suess effect is considerably smaller, attaining an averaged value of only −0.2 ± 0.1‰ due mainly to slow exchange rates between surface and deep waters.
The most pronounced 13 C-depletion in δ 13 C-DIC occurs in the subtropical surface waters, North Atlantic Deep Water, and global mode and intermediate waters (Fig. 2a). These surface regions are characterized by relatively higher sensitivities of surface DIC to increasing atmospheric CO 2 (i.e., the lower Revelle factor; Fig. 3e). The higher sensitivities, together with longer exposure times to the atmosphere, lead to the highest accumulation rate of 13 C-depleted anthropogenic CO 2 (refs. 30,31 ). Relatively weaker depletion occurs in upwelling-dominated areas such as the subpolar surface waters, equatorial surface waters, and the Southern Ocean surface waters where mixing from below dilutes the anthropogenic imprints at the sea surface. As a result, the perturbation ratio (defined as the ratio of deviations from the respective preindustrial values) of surface δ 13 C-DIC to atmospheric δ 13 C-CO 2 exhibits a large latitudinal contrast depending on the oceanic uptake of CO 2 and the air-sea equilibrium states relative to vertical mixing rates 30,31 . For example, as of 2000, the surface to atmospheric δ 13 C perturbation ratio ranges from 0.1 in the highly convective Southern Ocean to 0.7 ± 0.1 in the relatively stable Northern Hemisphere subtropical gyres (Fig. 3b), in agreement with a previous study 1 .
The best estimate for the preindustrial δ 13 C-DIC distributions can be derived by combining our constrained oceanic 13 C Suess effect with the contemporary δ 13 C-DIC observations 4 ( Fig. 2b; Supplementary Fig. 2; Methods). This approach is chosen rather than using simulated preindustrial δ 13 C-DIC distributions due to the large sensitivity of the preindustrial δ 13 C-DIC estimates to a poorly constrained boundary condition of terrestrial carbon inputs 23 (Supplementary Fig. 1a) and the unexplored sensitivity to the biological carbon pump 32 . Compared to a recent study from Eide et al. 33 , our estimate suggests more 13 C-enrichment in the Southern Hemisphere mode and intermediate waters with values attaining up to 2.6 ± 0.3‰, and more 13 C-depletion at 0.5 ± 0.2‰ in the Antarctic Bottom Water (Fig. 2b, d, f). The relatively larger discrepancies in the Southern Hemisphere, compared to the Northern Hemisphere, may arise due to sparse observational coverage, which we argue makes the previous observation-based estimate more prone to uncertainty. Nevertheless, our results support a previous finding that surface water δ 13 C-DIC was more positive during preindustrial times than the present-day observations, revealing higher horizontal and vertical gradients associated with water mass distributions in the upper ocean 33,34 . Our preindustrial estimate of oceanic δ 13 C-DIC can serve as a reference state with which past estimates of oceanic δ 13 C-DIC (e.g., δ 13 C-DIC during the Last Glacial Maximum, 21,000 years ago) can be compared 35 (see data availability statement).
Projected 21st century changes in δ 13 C-DIC. We use four CO 2 emission scenarios 24 and the linear relationship between the atmospheric CO 2 and the δ 13 C of CO 2 estimated based on the last decades of observations 36 (Supplementary Fig. 3) for the 21st century projections (Methods). Compared to the effects of the different CO 2 emission scenarios (Fig. 1b), the effects of potential changes in air-sea CO 2 exchange rates (another primary source of uncertainty identified for the estimated oceanic 13 C Suess effect as of 2018) appear to be an order of magnitude smaller (Supplementary Figs. 4 and 5; Supplementary Note 3). Additional sources of uncertainty, including temporal changes in ocean circulation and the biological carbon pump, also appear to be minor for the 21st century changes in DIC and ocean pH, as shown by the close agreements between our estimates and those projected from multiple climate models 37 (Supplementary Fig. 6b). Likewise, the effects of ocean circulation and biological changes are likely small for the 21st century δ 13 C-DIC change 26 . The modulating effects of climate change on δ 13 C-DIC distributions may emerge on timescales that are longer than the timescales of the geochemical effects from rising atmospheric CO 2 (refs. 11,38 ). Hence, we mainly focus on the effects of different CO 2 emissions scenarios applied to the full model setup with time-invariant ocean circulations and biogeochemical model parameters (Supplementary Table 1) for the 21st century δ 13 C-DIC projections.
All of the RCP scenarios suggest that oceanic δ 13 C-DIC becomes increasingly more 13 C-depleted near the sea surface, as the ocean takes up more 13 C-depleted atmospheric CO 2 (Figs. 1b and 3a; Supplementary Fig. 7a). The globally averaged surface ocean 13 C Suess effect is projected to range from -1.6‰ for RCP2.6 to -2.2‰ for RCP8.5 as of the year 2050. As time progresses, the surface ocean 13 C Suess effect diverges to a range  13 C Suess effect as of 1994. b Our estimate for the preindustrial δ 13 C-DIC obtained by correcting the global database 4 for the oceanic 13 C Suess effect. c The oceanic 13 C Suess effect estimated by Eide et al. 1 . d The preindustrial δ 13 C-DIC estimated by Eide et al. 33 . e Difference between Eide et al. 1 and our estimate. f Difference between Eide et al. 33 and our estimate. g The 95% confidence interval for our estimated 13 C Suess effect. h The 95% confidence interval for our estimated preindustrial δ 13 C-DIC.
of −(1.8-6.3)‰ as of 2100 across the RCP scenarios. The perturbation ratios of surface δ 13 C-DIC to atmospheric δ 13 C-CO 2 also grow in time, reaching a maximum value greater than one at the center of the subtropical gyres under the RCP2.6 scenario that stabilizes atmospheric CO 2 after 2040 ( Fig. 3b). Accordingly, surface δ 13 C-DIC remains nearly constant between 2040 and 2100 ( Fig. 1b). On the other hand, the perturbation ratio increases modestly in the RCP8.5 scenario (Fig. 3b) where the ocean surface does not have a sufficient time to equilibrate with the rapidly changing atmospheric δ 13 C-CO 2 ( Fig. 1b; ref. 11 ).
The most pronounced and robust change towards the end of this century across the RCP scenarios is the appearance of the lowest δ 13 C-DIC values in the thermocline waters (shallower than a depth of 1000 m), which is evident in the global ocean low latitudes and the North Pacific high latitudes (Fig. 3a). The tropical and North Pacific high-latitude thermocline waters coincide with the regions of the highest accumulation rates of the 13 C-depleted DIC remineralized from sinking organic particles, which contributed to the local minimum of the preindustrial δ 13 C-DIC (Fig. 2b). With the invasion of isotopically light anthropogenic carbon, the thermocline waters become further depleted in δ 13 C-DIC, emerging as new global δ 13 C-DIC minima. The thermocline δ 13 C-DIC minima gradually expand to the entire global ocean thermocline as time progresses.
Unlike the well-ventilated thermocline and relatively young North Atlantic Deep Water, the δ 13 C-DIC values in deep waters remain relatively unchanged with a subsurface averaged 13 C Suess effect less than -1‰ as of 2100 (Fig. 1b, c), becoming the waters of the most enriched δ 13 C-DIC. As a result, the naturally formed vertical gradients of δ 13 C-DIC are eliminated or reversed towards the end of the 21st century. For example, the difference between the globally averaged surface and subsurface δ 13 C-DIC becomes 0.2‰ for RCP2.6, -0.7‰ for RCP4.5, -1.6‰ for RCP6.0, and -3.8‰ for RCP8.5 as of 2100, compared to the preindustrial value of 1.6 ± 0.2‰ (Figs. 1c and 3a). Such varying degrees of vertical gradient reversals indicate a large sensitivity of the gradient disruption to the rapidity of CO 2 emissions, given the present-day ocean circulation rates. Furthermore, the regional magnitude of gradient reversal is larger in strongly stratified low latitude ocean than convective high latitudes, suggesting a sensitivity to ocean ventilation state as well (Fig. 3a). These reversed vertical gradients of δ 13 C-DIC are likely to persist at least over the next few centuries due to slow exchange rates between surface and deep waters and also to slow deep ocean circulation.
Projected global warming (e.g., ref. 39 ) is expected to further reduce surface δ 13 C-DIC through increasing air-sea CO 2 exchange rates in high latitudes (mostly due to sea ice melting) and enhanced thermodynamic isotopic fractionations whose effects are most pronounced in low latitudes (Supplementary Fig. 5; Supplementary Note 3). These warming-driven surface δ 13 C-DIC reductions can additionally elevate the ratio of the ocean to atmosphere 13 C Suess effect, and further, enhance the vertical δ 13 C-DIC gradient reversal as of 2100. When globally averaged, the warming effects are small with an additional surface δ 13 C-DIC decline of only -0.1‰ compared to the geochemically driven 13 C Suess effect of -3.7‰ under the RCP6.0 scenario ( Supplementary Fig. 4). Yet, regionally the effects can be as large as 50% in the Weddell Sea ( Supplementary Fig. 5). A more comprehensive assessment of ocean stratification and circulation The ratios of surface δ 13 C-DIC perturbation to atmospheric δ 13 C-CO 2 perturbation where the perturbation is defined as a change since the preindustrial times. c The ratios of surface δ 13 C-DIC perturbation to DIC perturbation where the perturbation is defined as a change since the preindustrial times. d The e-folding timescale for the air-sea δ 13 C equilibrium, estimated as (DIC/ CO 2 aq ) × (k w /z ml ) −1 where CO 2 aq is the aqueous CO 2 concentrations at the sea surface, k w is the piston velocity for air-sea CO 2 exchange, and z ml is the thickness of the mixed layer 40 . e The Revelle factor, defined as (DIC/CO 2 ) × (∂CO 2 /∂DIC) where CO 2 is oceanic pCO 2 . change effects on δ 13 C-DIC projections remains a subject for future studies.
Geochemical mechanisms for the ocean 13 C Suess effect. We elucidate the geochemical mechanisms by which the 21st century CO 2 emissions can drive such radical changes in surface δ 13 C-DIC. The rapid evolution of the oceanic 13 C Suess effect is somewhat surprising given that the ocean's ability to absorb anthropogenic CO 2 in fact diminishes under higher DIC conditions (Fig. 3e). For example, under the RCP8.5 the surfaceaveraged 13 C Suess effect increases eight times from −0.8‰ in 2000 to 6.2‰ in 2100 (Fig. 1b), whereas the surface-averaged anthropogenic DIC increases four times from 40 to 201 μmol kg −1 over the same time period (Supplementary Fig. 8b). The perturbation ratios of the surface δ 13 C-DIC to DIC also increase regionally by up to a factor of two from −15‰ per a 1 mmol kg −1 increase in 2000 to −30‰ per a 1 mmol kg −1 increase in 2100 under the RCP8.5 (Fig. 3c). Similarly, the perturbation ratio of surface δ 13 C-DIC to pH increases as the surface ocean becomes more acidic (Supplementary Fig. 6a, c). The progressively amplified response of surface δ 13 C-DIC, relative to surface DIC or pH, manifests as sharper vertical gradients for the 13 C Suess effect than anthropogenic DIC 30 . As a result, the vertical gradient of oceanic δ 13 C-DIC undergoes more radical changes by the end of the 21st century, while the difference between the surface-and subsurfaceaveraged DIC decreases moderately by 20-70% ( Supplementary  Fig. 8c).
The amplified response of δ 13 C-DIC, compared to 12 C, results from changes in seawater chemistry under higher CO 2 conditions. Elevated surface DIC brings surface δ 13 C-DIC closer to an equilibrium with atmospheric δ 13 C-CO 2 through the dependency of the air-sea exchange rates of 13 C/ 12 C on DIC speciation 40 . The globally averaged e-folding timescale, with which surface δ 13 C-DIC would approach a value in equilibrium with the atmosphere in the absence of other sources and sinks, decreases from 14 years as of the year 2000 to 6-12 years as of the year 2100 (with the smallest reduction in RCP2.6 and the largest in RCP8.5) (Fig. 3d). The enhanced air-sea δ 13 C equilibration rates can, in turn, facilitate the surface ocean Suess effect through either allowing for the surface ocean to catch up with atmospheric δ 13 C-CO 2 changes or alleviating the δ 13 C-DIC disequilibrium effects on a natural component of δ 13 C-DIC. These mechanisms are demonstrated in the following two idealized numerical experiments.
In the idealized fixed atmospheric CO 2 experiment, we fix the atmospheric CO 2 at a preindustrial value of 280 ppm while depleting the atmospheric δ 13 C-CO 2 over time as in the full model setup. In this experiment, the oceanic DIC and air-sea δ 13 C equilibration rates remain fixed at their preindustrial values in all RCP scenarios. This experiment yields the surface-averaged 13 C Suess effects that account for 60-80% of the 13 C Suess effects from the full model setup for the year of 2100 (Fig. 4). The remaining 20-40% (with the lowest value in RCP2.6 and the highest in RCP8.5) can be attributed mainly to the effect of enhanced air-sea equilibrium rates on δ 13 C-DIC, which facilitates the upper ocean depletion of δ 13 C-DIC.
To further explore the effect of enhanced air-sea equilibration rates on the natural component of δ 13 C-DIC, we perform an additional idealized experiment by fixing the atmospheric δ 13 C-CO 2 (or the ratio of 13 C/ 12 C) at a preindustrial value of -6.5‰ while atmospheric CO 2 increases as in the full model setup. The oceanic 13 C Suess effect from this fixed atmospheric δ 13 C-CO 2 setup is nearly negligible when globally averaged, and also small regionally with the estimated 13 C Suess effects that are an order of magnitude smaller than the Suess effects from the full model setup (Supplementary Fig. 7b). Yet, such modulations of the natural component of δ 13 C-DIC are interesting because 13 C-depletion in low latitudes are accompanied by opposing 13 C-enrichments in high latitudes (north of 40°N or south of 40°S) and deep waters, suppressing the subsurface 13 C Suess effect through deep water formation in high latitudes. This geochemically driven enrichment effect, which depends on the magnitude of CO 2 increases ( Supplementary Fig. 7b), can retard deep ocean δ 13 C-DIC responses to atmospheric δ 13 C-CO 2 depletion on millennial timescales (Fig. 5c; See below).
Comparison with PETM depth gradient reversal. The magnitudes of our projected surface δ 13 C-DIC excursion as of 2100 for the intermediate to high emission scenarios are comparable with the PETM carbon isotope excursion of −(3-5)‰ (Fig. 1b). The projected elimination or reversal of the vertical δ 13 C-DIC gradients (Fig. 1c) are also qualitatively similar to the reported changes in the benthic-planktic δ 13 C-DIC contrasts during the PETM onset 9,13 . Beyond the 21st century, future evolutions of oceanic δ 13 C-DIC are highly uncertain due to uncertainty in future atmospheric CO 2 changes and potential feedbacks from the ocean. Yet, our multi-millennial simulations for hypothetical futures, assuming constant atmospheric CO 2 after a year of 2500 and also considering geochemical effects only, reveal that the magnitude and duration of the gradient reversals are sensitive to the atmospheric δ 13 C-CO 2 excursions and the local ventilation states of the ocean. The surface ocean δ 13 C-DIC might drop by up to~2‰ under RCP2.6,~4‰ under RCP4.5, and up to~6‰ under RCP6.0, which outpace deep ocean δ 13 C-DIC decreases of similar magnitudes (Fig. 6a). The duration over which the gradient reversal persists ranges from none under RCP2.6 to a few millennia in poorly ventilated deep North Pacific under RCP6.0 and RCP8.5 ( Fig. 6a; Supplementary Fig. 9).
Analogous to the future projections, Fig. 6b, c shows that reversal or at least elimination of the oceanic vertical δ 13 C-DIC gradients are evident in all pairs of planktic and benthic foraminifera species from two open-ocean sites (ODP Sites 690 and 1209) where high-resolution PETM records are currently available 25 . A precise comparison of the δ 13 C-DIC excursions between the future and the PETM onset is hampered due to the difficulty in reconstructing high temporal resolution PETM δ 13 C records (e.g., refs. 25,41 ) and large uncertainties in age models (e.g., refs. 42,43 ). In particular, benthic foraminiferal extinction and CaCO 3 dissolution that might have led to data gaps during the PETM onset challenge a precise determination of the magnitude and duration of the vertical gradient reversal (e.g., ref. 18 ). Century-scale reversal or elimination events that are likely to occur under some RCP scenarios would not be detected in a typical pelagic deep ocean sediment core. Nevertheless, it is worth exploring whether the observed δ 13 C-DIC gradient reversal is consistent with previously suggested carbon emission rates during the PETM onset. For this, we apply the following two emission estimates to our model: a relatively rapid increase in atmospheric CO 2 over 5 kyr 19 and a slow increase over 20 kyr 8 . In doing so, it is important to consider different ocean ventilation states, because the duration and magnitude of the vertical δ 13 C-DIC gradient reversals are also sensitive to ocean stratification and ventilation rates 11,21 .
Our sensitivity experiments (Methods) suggest that the surface to deep (average over 2-3 km depth) gradient reversal only occurs when an ocean model with poorly ventilated deep waters (an average ventilation age three times longer than the present-day ocean) is forced with rapidly increasing atmospheric CO 2 over 5 kyr (Fig. 5a, b; Supplementary Fig. 9). While the surface and well-ventilated deep waters keep up with the atmospheric δ 13 C-CO 2 change, the δ 13 C-DIC responses of stagnant deep waters (i.e., the North Pacific in the present-day configuration) are substantially delayed, temporally reversing the δ 13 C-DIC depth gradients. The simulated magnitude and duration of the vertical gradient reversal are comparable with those from a PETM record (the lower panel of Fig. 6b), and exhibit strong sensitivity to the atmospheric δ 13 C-CO 2 excursions (Fig. 5b). On the other hand, no discernible gradient reversal is simulated in the presentday ocean circulation state with deep ocean ventilation ages <1400 years (Fig. 5a), except for some localized reversals ( Supplementary Fig. 10).
Under scenarios of slowly increasing atmospheric CO 2 over 20 kyr (ref. 8 ), no reversals of the vertical δ 13 C-DIC gradient are apparent even in poorly ventilated regions (Fig. 5c). Slow increases in atmospheric CO 2 (even slower than the vertical ocean mixing timescales) give the deep ocean sufficient time to fully equilibrate to changing atmospheric CO 2 . Interestingly, the deep ocean δ 13 C-DIC exhibits distinct responses to the same atmospheric δ 13 C-CO 2 excursion of −4‰ between the ×5 scenario where atmospheric CO 2 increases approximately fivefold and the ×2 scenario where atmospheric CO 2 increases approximately twofold (Fig. 5c). The 13 C enrichment effect in highlatitude surface from enhanced air-sea δ 13 C equilibrium tends to offset the depleting effect from the invasion of 13 C-depleted CO 2 . This offset is greater in the ×5 scenario than in the ×2 scenario, such that during the perturbation period the vertical δ 13 C-DIC gradients become eliminated in relatively well-ventilated regions (i.e., the North Atlantic in the present-day configuration) under the ×5 scenario (Fig. 5c).
The evidence of a pronounced vertical gradient reversal at the ODP Site 690 (Fig. 6b), therefore, supports a rapid carbon emission rate towards a 5 kyr perturbation period and also suggests more stagnant deep waters during the PETM onset with slower ventilation rates than the present-day ventilation rates. Our inference is in agreement with previous studies 44,45 suggesting an overall decline in deep ocean ventilation rates during the PETM onset, possibly linked with the reorganization of deep ocean circulation patterns 10,46,47 . There remain important aspects of the observed PETM δ 13 C-DIC excursion that are not captured by our simulation, including much delayed and reduced declines in benthic δ 13 C, slight 13 C-enrichments in benthic δ 13 C prior to declines, and different magnitudes of benthic δ 13 C declines between the two ODP sites (Fig. 6b, c). These disparities imply that concomitant transient changes in ocean circulation, marine biology, land-derived carbon inputs, and marine sedimentary dissolutions of CaCO 3 might have been as important as the geochemical effects for the early PETM δ 13 C excursion (e.g., refs. 10,11,21,44 ). In fact, Kirtland Turner and Ridgwell 11 showed that the CO 2 -climate feedbacks during the PETM onset can delay the time that it takes for the surface δ 13 C-DIC minimum to propagate to the deep ocean δ 13 C-DIC minimum by up to 40%.
Despite the analogy between the future and the PETM onset in terms of large oceanic δ 13 C-DIC excursions and the existence of vertical gradient reversal, the rate at which the 21st century anthropogenic carbon isotope excursion occurs is at least one order of magnitude faster than PETM excursion rates. The 21st century δ 13 C-DIC gradient reversal rates (taking only~3 centuries from the preindustrial era to the maximum surface δ 13 C-DIC excursion) (Fig. 6a) appear to be much faster than those of the PETM (taking at least 3 kyr from the pre-PETM to the maximum surface δ 13 C-DIC excursion) (Fig. 6b, c). Given the fact that the PETM is the best known geological period when the most rapid carbon emissions have occurred over the Cenozoic (e.g., refs. 20,48 ), our comparison suggests that the time rates of 21st century δ 13 C-DIC excursion and associated gradient reversal may be unprecedented over the Cenozoic.

Methods
The ocean carbon cycle model. We use an offline global ocean carbon cycle model 23 where a simple ocean biogeochemistry module is embedded in an observationally constrained ocean circulation inverse model (OCIM) 22 . The model has a horizontal resolution of 2°× 2°and a total of 24 vertical layers where the layer thickness increases from 36 m near the top to 633 m near the bottom. The circulation model uses a simplified version of the primitive equations using hydrostatic, rigid lid, and Boussinesq approximation, and produces ocean transport estimates that are optimally consistent with various observations of temperature, salinity, CFC-11, radiocarbon 22 . Large-scale ocean circulation and ventilation features are shown to be consistent with independent observation-based estimates 49 (See also Supplementary Fig. 10 for the simulated deep ocean ventilation age).
The ocean biogeochemical processes are formulated following the OCMIP2 protocol 50 where ocean productivity is simulated by restoring model surface PO 4 towards the observed PO 4 . Simple parameterizations for the production and remineralization of organic and inorganic carbon are employed. The carbon isotope model uses two prognostic variables of DI 13 C and DI 12 C (the latter approximated as DIC), and the isotopic signature of DIC is estimated as δ 13 C-DIC = [(DI 13 C/DI 12 C) sample /(DI 13 C/DI 12 C) standard −1] with the Vienna Pee Dee Belemnite standard. The model employs a temperature-dependent thermodynamic equilibrium fractionation from Zhang et al. 7 and a CO 2 -dependent photosynthetic fractionation following three different empirical formulations (Supplementary Table 1). The model includes riverine carbon inputs and a simple parameterization for the sedimentary burial of inorganic carbon.
Preindustrial steady-state solutions are obtained using a time-efficient Newton's method 51 whereas industrial changes are simulated by taking time steps with an  24 for a time period of 2019-2100. Our model is based on an annual mean climatology that is invariant over time except for the atmospheric CO 2 and δ 13 C-CO 2 forcing that change over industrial times. The model was previously shown to simulate the oceanic uptake and storage of anthropogenic carbon 22 , as well as the oceanic 13 C Suess effect 23 (Fig. 1a), consistent with previous independent estimates. We refer readers to Kwon et al. 23 and Supplementary Note 1 for details of the carbon isotope model formulations. Although our primary focus is the geochemical effects from changing atmospheric CO 2 , our offline model framework allows us to perform an experiment where some of the biogeochemistry model inputs are taken from results of an Earth System Model (Supplementary Note 3).
A Monte Carlo experiment for historical simulations. In order to explore a wide range of potential sources of uncertainty in simulating the oceanic δ 13 C-DIC until 2018, we perform a Monte Carlo experiment of 1400 ensemble members where the model parameter values and setups are randomly drawn from plausible ranges that are assumed to be uniformly distributed (Supplementary Table 1). Potential sources of uncertainty considered in this study include (A) preindustrial δ 13 C values for atmospheric CO 2 of −(6.3-6.5)‰ 53 , which combined with observed values over 1980-2018 determine the atmospheric 13 C Suess effect over industrial times, (B) temperature-dependent thermodynamic equilibrium fractionation factors for air-sea CO 2 exchange 7 and time-varying sea surface temperatures over industrial times [54][55][56] , (C) the globally uniform δ 13 C values of riverine carbon inputs of -27 ± 2‰ for dissolved organic carbon, −30 ± 2‰ for particulate organic carbon, and −15 ± 2‰ for DIC 57 , (D) the magnitude of non-riverine terrestrial carbon inputs with a δ 13 C value of −26‰ 58 , including uncertainties in groundwater driven fluxes 59 and the carbon export from coastal vegetation 60 , that are assumed to be uniformly distributed along the coastal margins except around the Antarctica, (E) the air-sea CO 2 exchange rates formulated following Wanninkhof 61 and Najjar et al. 50 , which are linearly scaled such that the globally averaged transfer rate ranges between 13 and 17 cm h −1 (refs. 62,63 ), (F) CO 2 dependent fractionation factors for the photosynthetic uptake of carbon, taken from the three different empirical formulations 5,64,65 , (G) present-day ocean circulation states including different mixing parameterizations 22 , and (H) the amount of inorganic carbon buried in marine sediments 66 (See Supplementary Note 1 for the model parameters). Using the Monte Carlo experiment, a median value with a 95% confidence interval (two times one standard deviation) is reported for the estimated 13 C Suess effect and DIC change. We also perform regression analyses using the model input parameter values and the corresponding model estimates to explore the relative contributions of the different sources to uncertainty.
Estimation of preindustrial δ 13 C-DIC. We estimate a preindustrial δ 13 C-DIC distribution by combining the global compilation 4 of observed δ 13 C-DIC over 1972-2016 with our estimated oceanic 13 C Suess effect. We first map the observed δ 13 C-DIC onto our 2°× 2°model grid cells by averaging all data points falling within each grid cell at each year. To obtain a climatological mean distribution, we average the gridded data over time such that each grid cell has an averaged year of data collection and an averaged δ 13 C-DIC value (Supplementary Fig. 2a). Then, the mapped δ 13 C-DIC is corrected for the 13 C Suess effect with our estimate that is taken from the averaged year of data collection for each grid cell ( Supplementary  Fig. 2c). The uncertainty of the resulting preindustrial δ 13 C-DIC estimate (Fig. 2h) includes two standard mapping errors for the present-day δ 13 C-DIC observations, measurement uncertainties of ±0.1‰, and the uncertainty in our estimated 13 C Suess effect (Fig. 2g). This approach gives a preindustrial state of δ 13 C-DIC ( Fig. 2b; Supplementary Fig. 2b) that is similar to our simulated preindustrial δ 13 C-DIC with reduced uncertainty.
Projections for the 21st century and sensitivity experiments. As a primary source of uncertainty in future projection, we consider the geochemical effects of different CO 2 emission scenarios, as represented by the four RCPs (i.e., RCP2.6, Fig. 5 Simulated responses of oceanic δ 13 C-DIC to the different atmospheric CO 2 and δ 13 C-CO 2 forcing previously suggested for the early PETM. Surface δ 13 C-DIC values averaged over the global ocean are shown in solid lines. The deep (2-3 km depths) ocean δ 13 C-DIC values averaged over the North Atlantic (N.Atl., north of 20°N) and North Pacific (N.Pac., north of 20°N) are shown in dotted and dashed lines, respectively. a Results from the rapid PETM carbon emission setup where atmospheric CO 2 linearly increases twofold over 5 kyr according to Penman et al. 69 with maximum atmospheric δ 13 C-CO 2 excursions of −4‰ (red and green) and −5‰ (yellow and gray). We use the present-day ocean circulation state. b Same as a except that we use an ocean circulation state with slower ventilation rates (an average three fold increase in ventilation ages). c Results from the slow PETM carbon emission setup where atmospheric δ 13 C-CO 2 decreases by 4‰ over 20 kyr according to Cui et al. 8 with two different magnitudes of atmospheric CO 2 increases as suggested by Cui et al. 8 . Atmospheric CO 2 is assumed to linearly increase five fold over 20 kyr in red and green lines (×5) and twofold in yellow and gray lines (×2). We use the ocean circulation state of slower ventilation rates, as used in (b).
RCP4.5, RCP6.0, and RCP8.5) 24 . To this end, we use a model configuration, considered as a typical model setup in carbon isotope modeling 3,28,67 , taken from the Monte Carlo experiment (the "full" model setup in Supplementary Table 1). Using the full model setup, we assume that the ocean circulation, air-sea CO 2 exchange rates including sea ice effects, and sea surface temperature and salinity remain unchanged with time throughout the simulation. An uncertainty associated with this assumption is tested in Supplementary Note 3 where we relax the assumption of unchanged air-sea CO 2 exchange rates, and sea surface temperature and salinity using the Community Earth System Model version 2 (CESM2)-based estimates 39,68 . We also assume that the same linear relationship between  Supplementary Fig. 3), holds for the atmosphere from 2019 to 2100.
We also perform two idealized experiments in order to elucidate how changing oceanic carbon chemistry affects the oceanic δ 13 C-DIC changes. In order to suppress the effect of changing oceanic carbonate chemistry, we perform "fixed.CO2" experiment that is identical to the full model setup except that the atmospheric CO 2 is fixed at a preindustrial value of 280 ppm while the δ 13 C of atmospheric CO 2 is changing as in the full model setup. Thus, the difference between the full model setup and the "fixed.CO2" setup indicates the effects of increasing oceanic DIC on simulated oceanic δ 13 C-DIC changes, which are dominated by the effects of enhanced air-sea carbon isotope equilibration rates under higher DIC conditions 40 . In order to explore the effect of changing oceanic carbonate chemistry on the natural component of oceanic δ 13 C-DIC, we also perform "fixed.ratio" experiment where atmospheric CO 2 changes as in the full model setup while the δ 13 C-CO 2 remains fixed at a preindustrial value of −6.5‰. The difference between the full model setup and "fixed.ratio" setup is then interpreted as an effect of changing carbonate chemistry on natural δ 13 C-DIC. In the "fixed.ratio" model, we turn off the kinetic fractionations during air-sea CO 2 exchange and the CO 2 -dependent photosynthetic fractionations in order to focus only on the effect of changing DIC on the redistribution of δ 13 C-DIC. Because the combined effects of kinetic fractionations and the CO 2 -dependent photosynthetic fractionations are negligible with values within ±0.1‰ for all RCP scenarios, whether we turn on or off both effects in the "fixed.ratio" setup do not make discernible differences in our interpretations of the experiment results.
Multi-millennial simulations for hypothetical futures. Although not a primary focus in this study, we also perform multi-millennial simulations to explore how sensitive the magnitude and duration of the vertical δ 13 C-DIC gradient reversal are to the RCP scenarios. In these hypothetical deep future simulations, we extend the RCP-based future projections by assuming that atmospheric CO 2 and its δ 13 C-CO 2 stay constant from the year 2500 to a year 8000 at the values of the year of 2500. Specifically, the RCP scenarios for atmospheric CO 2 until the year of 2500 are taken from Moss et al. 24 and the same linear relationship between atmospheric CO 2 and δ 13 C-CO 2 used for our 21st century projections are assumed for the extended future simulations. Both atmospheric CO 2 and its isotopic signature are kept constant from the year of 2500 onwards until the vertical δ 13 C-DIC gradients are recovered. We use the full model setup without considering changes in ocean circulation, air-sea gas exchange rates, and biogeochemical processes.
Sensitivity experiments with estimated PETM atmospheric CO 2 . Based on the present-day ocean model configurations (See below), we use two estimates of the early PETM atmospheric CO 2 changes following Penman and Zachos 19 and Cui et al. 8 . The two estimates are chosen because they represent both sides of the spectrum for the PETM onset period of 3-20 kyr 8,12,20,69 . In the "rapid" setup, atmospheric CO 2 linearly increases from 750 to 1500 ppm over 5 kyr according to Penman and Zachos 19 . The "rapid" setup is branched into two experiments with two different maximum atmospheric δ 13 C-CO 2 excursions of −4‰ (refs. 8,70 ) and −5‰, the latter corresponding to a typical δ 13 C excursion of land plants during the PETM 70 . Atmospheric δ 13 C-CO 2 is assumed to linearly decrease from −6.5‰ at a model year of zero to either −10.5‰ and −11.5‰ at a model year of 5 kyr, which corresponds to the atmospheric CO 2 perturbation period. In the "slow" setup, atmospheric CO 2 linearly increases from 834 ppm to either 1500 ppm or 4200 ppm over 20 kyr, as denoted as "×2" and "×5" scenarios respectively, following Cui et al. 8 . Both "×2" and "×5" experiments use the same atmospheric δ 13 C-CO 2 excursion of −4‰, by assuming that atmospheric δ 13 C-CO 2 decreases from −6.5‰ at a model year of zero to −10.5‰ at a model year of 20 kyr, which corresponds to the atmospheric CO 2 perturbation period.
In both "rapid" and "slow" setups, we use two ocean circulation states: the present-day ocean circulation as in the full model setup and a circulation where an averaged deep ocean (2-3 km depth) ventilation age is three times larger at 2237 years compared to 758 years for the present-day ocean (See Supplementary Fig. 10 for the distribution of deep ocean ventilation age). The slow ocean ventilation state was previously named as "KL" model in Kwon et al. 71 , and has slower meridional overturning rates of 12 Sv (1 Sv = 10 6 m 3 s −1 ) for the North Atlantic Deep Water and 5 Sv for the Antarctic Bottom Water, compared to the present-day circulation model 49 of 20 Sv and 16 Sv, respectively. Despite the different overturning rates and ventilation states, both circulation models are based on the present-day configurations of ocean bathymetry and atmospheric buoyancy and momentum forcing. Therefore, none of the circulation states might represent a plausible state for the PETM. Nevertheless, the PETM ocean would still have regions of relatively wellventilated and poorly ventilated deep waters 44,45 , although they would not correspond to the North Atlantic and North Pacific, respectively, for the present-day configuration. Hence, the basin-scale contrast in deep ocean δ 13 C-DIC responses should be regarded as a contrast between well vs. poorly ventilated deep waters.
We only explore the geochemical effects from changing atmospheric CO 2 and δ 13 C-CO 2 on δ 13 C-DIC changes, assuming that ocean circulation, ocean temperature, air-sea CO 2 exchange rates, and marine biological pump remain fixed at those from the present-day ocean throughout the multimillennial simulations.
Foraminifera δ 13 C records for the PETM onset. We revisit the early PETM evolutions of surface and deep ocean δ 13 C-DIC, as inferred from planktic and benthic foraminifera δ 13 C records, using the latest compilation of Shaw et al. 25 . Two ODP Sites 690 and 1209 are chosen for the PETM records because the two ODP sites are located in the open ocean and their δ 13 C records have high temporal resolutions, the latter which is crucial for resolving the δ 13 C-DIC evolution during the PETM onset. We assume that the δ 13 C values of mixed layer planktic foraminifera represent the surface δ 13 C-DIC whereas the δ 13 C values of epifaunal benthic foraminifera represent the δ 13 C-DIC at the depths of the ODP sites. Those inferred δ 13 C-DIC values are then compared with our model-based estimates averaged over 0-74 m depths for planktic and over 2-3 km depths for benthic.
The original sources of data presented in Fig. 6b, c are as follows: For the ODP Site 690 (South Atlantic, 65°S, 1°E, 2100 m deep), the δ 13 C records of N.truempyi, Acarinina spp., A.mckannai, A.praepentacamerata are from refs. 13,42,[72][73][74] . For the ODP Site 1209 (North Pacific, 33°N, 159°E, 2387 m deep), the δ 13 C records of N.truempyi are from ref. 41 , the A.soldadoensis data are from ref. 75 , and the M.velascoensis data are from refs. 69,75 . No corrections are applied to the foraminifera δ 13 C records because of the lack of information. Nevertheless, uncertainties are likely small in this single species comparison at the two open ocean sites where a reduced photosymbiosis has been reported for the PETM acarininids and morozovellids (therefore less offset from seawater δ 13 C-DIC) 25 . The age models for the foraminifera δ 13 C records are from Thomas and Shackleton 42 and Bains et al. 43 for ODP Site 690 and Westerhold et al. 41 for ODP Site 1209.

Data availability
The global database for the present-day δ 13 C of DIC is available at https:// www.ncdc.noaa.gov/paleo-search/study/21750. The latest compilations of foraminifera δ 13 C for the PETM are available at https://doi.pangaea.de/10.1594/PANGAEA.922272 and https://doi.pangaea.de/10.1594/PANGAEA.922292. The atmospheric CO 2 and its isotopic composition data used for model forcing are from http://scrippsco2.ucsd.edu/data/ atmospheric_co2/. The sea surface temperature data used in this study are from https:// www.metoffice.gov.uk/hadobs/hadisst/data/download.html, https://www.esrl.noaa.gov/psd/ data/gridded/data.noaa.ersst.v5.html, and https://www.esrl.noaa.gov/psd/data/gridded/ data.cobe2.html. The projected atmospheric CO 2 concentrations for the four RCP scenarios are available at https://tntcat.iiasa.ac.at/RcpDb/. The CESM2 Large Ensemble Simulations are available at https://www.cesm.ucar.edu/projects/community-projects/LENS2/datasets.html. All of the model results presented here are made available at https://climatedata.ibs.re.kr/ data/papers/kwon-et-al-2022-commsenv. Fig. 6 δ 13 C-DIC changes during industrial times and early PETM. a Estimates of surface and deep ocean (2-3 km) δ 13 C-DIC averaged over the global ocean are shown. The Monte Carlo experiment-based estimates from 1800 to 2018 are combined with our projections from 2019 to 5000, assuming that atmospheric δ 13 C-CO 2 remains constant after the year 2500. The time periods over which the globally averaged surface δ 13 C-DIC values fall below the globally averaged deep ocean δ 13 C-DIC values are highlighted as yellow shading. The duration of the vertical gradient reversal is sensitive to the geographic locations, as shown for the North Atlantic (N.Atl., north of 20°N) and the North Pacific (N.Pac., north of 20°N). b Red and blue scatters show planktic and benthic foraminifera δ 13 C records, respectively, from ODP Site 690 (South Atlantic, 65°S, 1°E, 2100 m deep). Red and blue solid lines show linear interpolations of the respective scatters, after averaging the data with a 0.01 mbsf bin interval of sediment depth. The top X-axis labels denote the relative ages since the base of the PETM onset. The relative age outside parentheses is based on Thomas and Shackleton 42 while the relative age based on Bains et al. 43 is shown inside parentheses. The base of the PETM onset is denoted as a vertical dotted line and taken from Nunes and Norris 10 . The periods over which planktic δ 13 C values are lower than benthic δ 13 C values are highlighted as yellow shading. In the bottom panel of (b), we overlay the red and blue dashed lines, showing the model-based estimates of the globally averaged surface δ 13 C-DIC and an averaged δ 13 C-DIC at 2-3 km depths of the North Pacific, taken from Fig. 5b red solid and green dashed lines, respectively. The positions of both dashed lines are shifted such that the positions for initial declining match between the foraminifera records and the simulation. c Same as b except that the data from ODP Site 1209 (North Pacific, 33°N, 159°E, 2387 m deep) is shown. The base of the PETM onset (vertical dotted line) is taken from Petrizzo 76 . The age model is from Westerhold et al. 41 .