Articles | Volume 20, issue 3
Research article
02 Apr 2024
Research article |  | 02 Apr 2024

Assessing transient changes in the ocean carbon cycle during the last deglaciation through carbon isotope modeling

Hidetaka Kobayashi, Akira Oka, Takashi Obase, and Ayako Abe-Ouchi

Atmospheric carbon dioxide concentration (pCO2) has increased by approximately 80 ppm from the Last Glacial Maximum (LGM) to the early Holocene. The change in this atmospheric greenhouse gas is recognized as a climate system response to gradual change in insolation. Previous modeling studies suggested that the deglacial increase in atmospheric pCO2 is primarily attributed to the release of CO2 from the ocean. Additionally, it has been suggested that abrupt change in the Atlantic meridional overturning circulation (AMOC) and associated interhemispheric climate changes are involved in the release of CO2. However, understanding remains limited regarding oceanic circulation changes and the factors responsible for changes in chemical tracers in the ocean during the last deglaciation and their impact on atmospheric pCO2. In this study, we investigate the evolution of the ocean carbon cycle during the last deglaciation (21 to 11 ka BP) using three-dimensional ocean fields from the transient simulation of the MIROC 4m climate model, which exhibits abrupt AMOC changes similar to those observed in reconstructions. We investigate the reliability of simulated changes in the ocean carbon cycle by comparing the simulated carbon isotope ratios with sediment core data, and we examine potential biases and overlooked or underestimated processes in the model. Qualitatively, the modeled changes in atmospheric pCO2 are consistent with ice core records. For example, during Heinrich Stadial 1 (HS1), atmospheric pCO2 increases by 10.2 ppm, followed by a reduction of 7.0 ppm during the Bølling–Allerød (BA) period and then by an increase of 6.8 ppm during the Younger Dryas (YD) period. However, the model underestimates the changes in atmospheric pCO2 during these events compared to values derived from ice core data. Radiocarbon and stable isotope signatures (Δ14C and δ13C) indicate that the model underestimates both the activated deep-ocean ventilation and reduced efficiency of biological carbon export in the Southern Ocean and the active ventilation in the North Pacific Intermediate Water (NPIW) during HS1. The relatively small changes in simulated atmospheric pCO2 during HS1 might be attributable to these underestimations of ocean circulation variation. The changes in Δ14C associated with strengthening and weakening of the AMOC during the BA and YD periods are generally consistent with values derived from sediment core records. However, although the data indicate continuous increase in δ13C in the deep ocean throughout the YD period, the model shows the opposite trend. It suggests that the model either simulates excessive weakening of the AMOC during the YD period or has limited representation of geochemical processes, including marine ecosystem response and terrestrial carbon storage. Decomposing the factors behind the changes in ocean pCO2 reveals that variations in temperature and alkalinity have the greatest impact on change in atmospheric pCO2. Compensation for the effects of temperature and alkalinity suggests that the AMOC changes and the associated bipolar climate changes contribute to the decrease in atmospheric pCO2 during the BA and the increase in atmospheric pCO2 during the YD period.

1 Introduction

Earth's climate has shifted from the colder conditions of the Last Glacial Maximum (LGM) to the warmer conditions of the Holocene. This climatic transition, known as the last deglaciation, occurred approximately 21 to 11 ka BP (thousand years before present). During this period, the atmospheric concentration of carbon dioxide (pCO2) increased by almost 80 ppm (Barnola et al.1987; Petit et al.1999; Siegenthaler et al.2005; Jouzel et al.2007; Lüthi et al.2008). The changes in the carbon cycle that affect the variation in atmospheric pCO2 are closely related to the changes in climate observed during the last deglaciation.

In attempting to elucidate the mechanisms of climate change on the glacial–interglacial scale, previous modeling studies mainly focused on assessing the steady-state difference between the LGM and the preindustrial period. With the development of computational tools, transient climate modeling of the last glacial termination has recently been conducted, using temporal changes in insolation, greenhouse gas concentrations derived from ice core records, and meltwater fluxes from ice sheets (Lunt et al.2006; Timm and Timmermann2007; Liu et al.2009; He et al.2013; Ganopolski and Roche2009; Menviel et al.2011; Ivanovic et al.2016; Obase and Abe-Ouchi2019; Obase et al.2021; Kapsch et al.2022; Bouttes et al.2023). Transient climate modeling has distinct advantages because it avoids unrealistic equilibrium assumptions, and it includes climate responses to internal variability or abrupt change. It also facilitates direct comparison between models and proxies, thereby allowing identification of temporal leads or lags in the process with respect to forcing. In these transient climate modeling studies, changes in atmospheric pCO2, a greenhouse gas, are applied as external forcing. However, understanding the feedback between the climate and the carbon cycle is critical for understanding the long-term changes in climate dynamics. As fundamental research guided by this perspective, earlier modeling studies examined the temporal changes in the ocean carbon cycle during the last deglaciation or late Pleistocene, including glacial–interglacial cycles, using Earth system models of intermediate complexity, e.g., CLIMBER-2 (Bouttes et al.2012a; Brovkin et al.2012; Mariotti et al.2016; Ganopolski and Brovkin2017), Bern3D (Tschumi et al.2011; Menviel et al.2012; Pöppelmeier et al.2023), and LOVECLIM (Menviel et al.2018; Stein et al.2020).

Those earlier studies greatly advanced our understanding of changes in both the climate and the carbon cycle on scales of thousands to tens of thousands of years. However, the mechanisms behind the observed increase in atmospheric pCO2 during the glacial termination are not fully understood. It has been suggested that the release of carbon from the deep Southern Ocean to the atmosphere might have played a role in the rapid increase in atmospheric pCO2 during the last deglaciation (Tschumi et al.2011; Bouttes et al.2012a; Mariotti et al.2016; Menviel et al.2018; Stein et al.2020; Sigman et al.2021; Gray et al.2023; Sikes et al.2023). The release of CO2 can be triggered by disruption of the stratification in the Southern Ocean and changes in the deep-ocean circulation, which are affected by westerly winds in the Southern Hemisphere and brine rejection around Antarctica. In contrast, freshwater input from the Antarctic ice sheet increases stratification and does not lead to enhanced CO2 outgassing. It has also been suggested that variation in North Pacific Intermediate Water (NPIW), associated with changes in the Atlantic meridional overturning circulation (AMOC), might have contributed to the observed increase in atmospheric pCO2 (Okazaki et al.2010; Menviel et al.2014).

To clarify the mechanisms behind long-term changes in ocean dynamics, the concentration or isotopic composition of elements in seawater can provide information on the processes responsible for their distribution. Thus, analysis of geochemical proxies in paleoclimatological archives can help oceanographers gain insight into past ocean variability and its underlying mechanisms.

Stable and radioactive carbon isotopes of dissolved inorganic carbon (DIC) in seawater are representative oceanic chemical tracers. In seawater, lighter carbon isotopes are preferentially taken up by phytoplankton during photosynthesis, and they accumulate in the deep ocean as organic matter is degraded (O'Leary1981; Broecker and Maier-Reimer1992; Schmittner et al.2013). Anomalies in the stable carbon isotope signature (δ13C) produced by the biological carbon pump spread globally in association with the deep-ocean circulation. Therefore, the δ13C value differs among water masses within the ocean interior. Radiocarbon (14C) is introduced into seawater through gas exchange at the ocean surface, and it subsequently decreases in concentration through radioactive decay of 14C. Therefore, the radiocarbon isotope signature (Δ14C) serves as an indicator of the deep-water flow rate (Stuiver et al.1983).

Previous modeling studies examined the processes responsible for glacial–interglacial changes in the distribution of δ13C and Δ14C (Kurahashi-Nakamura et al.2017; Menviel et al.2017; Muglia et al.2018; Wilmes et al.2021; Kobayashi et al.2021). It is generally assumed that deep water originating from the Southern Ocean with low δ13C and Δ14C expanded into the deep Atlantic Ocean during the LGM. It is proposed that the carbon isotope distribution during the LGM can be better explained by considering an effective biological pump associated with iron fertilization in the Southern Ocean and a shallower AMOC (Kobayashi et al.2021).

Other proxies used to infer past water mass distribution and deep-ocean circulation include the neodymium isotope ratio (εNd) and the protactinium–thorium ratio (231Pa/230Th). The εNd ratio can be used as a proxy for basin-scale water mass structure because the endmembers differ in each water mass source region (Lippold et al.2012; Howe et al.2016). The sedimentary 231Pa/230Th ratio, used as a proxy for change in flow rate, suggests that the strength of the AMOC might have changed markedly during the last deglaciation (McManus et al.2004; Ng et al.2018). Changes in the AMOC can influence the climatic state by altering not only the meridional interhemispheric heat transport but also the ocean carbon cycle and associated changes in atmospheric pCO2 (Schmittner and Galbraith2008; Menviel et al.2008; Bouttes et al.2012b). Nevertheless, there is ongoing debate regarding the magnitude and the direction of the change in atmospheric pCO2 associated with AMOC variation (Gottschalk et al.2019).

Several earlier modeling studies attempted to estimate AMOC variation by examining changes in the carbon isotope signature (Schmittner and Lund2015; Pöppelmeier et al.2023). By estimating the AMOC that best fits the model and the data for δ13C (Schmittner and Lund2015) or δ13C and multiple chemical tracers in seawater (Pöppelmeier et al.2023), those studies improved our understanding of the relationship between changes in the AMOC and the ocean carbon cycle during the early deglaciation. In recent years, compilation of many sediment core records covering the last deglaciation has deepened our understanding of the spatiotemporal changes in carbon isotope ratios during this period (Zhao et al.2018; Rafter et al.2022; Muglia et al.2023; Skinner et al.2023). By comparing the compiled data with model output, it has been possible to gain valuable insights into deglacial changes in the ocean carbon cycle.

In this study, we conducted transient model experiments of the ocean carbon cycle and compared the model results with recently compiled sediment core records to investigate the mechanisms of deglacial changes in carbon isotope signatures. Additionally, we analyzed the drivers of the changes in atmospheric pCO2 resulting from variations in the ocean carbon cycle. The objectives of this study were to clarify the reproducibility of deglacial carbon cycle changes in the model, understand the mechanisms underlying these changes, and identify the processes that are missing or underestimated in the model.

2 Methods

2.1 Model

In this study, numerical experiments on the ocean carbon cycle were performed using an offline ocean biogeochemical tracer model based on Parekh et al. (2005) within the framework of the CCSR Ocean Component Model version 4.0 (Hasumi2006). This model has an approximate horizontal resolution of 1°, and it includes 44 vertical layers with thicknesses ranging from 5 to 250 m.

The ocean biogeochemical cycle model is forced with monthly averaged output from a climate model simulation designed specifically for the last deglaciation and run with the MIROC 4m atmosphere–ocean general circulation model (AOGCM) (Obase and Abe-Ouchi2019; Obase et al.2021). The boundary conditions include horizontal advection velocity, sea surface height, vertical diffusivity, temperature, salinity, shortwave radiation, wind speed above the sea surface, and sea ice concentration.

Prognostic variables for the ocean biogeochemical cycle model include phosphate, DIC, alkalinity, dissolved organic phosphate, dissolved oxygen, iron, silicate, and carbon isotopes of DIC (13C and 14C). The availability of light, phosphate, and iron is used to determine the rate of phosphate uptake by phytoplankton. Notably, sedimentation processes on the seafloor are not considered, and all particles reaching the seafloor are assumed to dissolve in the deepest layer of the model.

2.2 Experimental design

To evaluate the transient response of the global ocean carbon cycle during the last deglaciation, we conducted offline ocean carbon cycle experiments forced with the outputs of the AOGCM MIROC 4m simulation, covering the period 21 to 11 ka BP. The MIROC 4m simulation focusing on the last deglaciation was performed according to the PMIP protocol (Ivanovic et al.2016) with respect to the changes in orbital parameters and greenhouse gases during this period (Obase and Abe-Ouchi2019). However, the ice sheets were fixed at their 21 ka BP state of the ICE-5G reconstruction. Freshwater inflow from the Northern Hemisphere ice sheets deviates from the PMIP protocol after the latter half of Heinrich Stadial 1 (HS1). This approach seeks to align the simulated AMOC variations with those reconstructed from 231Pa/230Th sediment core records and the associated climate changes that occurred during the Bølling–Allerød (BA) and Younger Dryas (YD) periods.

2.2.1 Steady-state experiment on the ocean carbon cycle for the Last Glacial Maximum

The ocean biogeochemical cycle model was initialized through spin-up under the LGM ocean state (21 ka BP) calculated by the AOGCM. Dust deposition to the ocean surface was taken from a simulation conducted using the SPRINTARS aerosol transport–radiation model computed under LGM climatic conditions (Takemura2005). We assumed that iron deposition at the sea surface accounts for 3.5 wt % of the total dust deposition, with an assumed iron solubility of 1 % (Parekh et al.2005), which was derived from the ratio of wet and dry dust deposition and its solubility. However, it should be noted that some uncertainty is associated with these parameters. The initial distribution of ocean biogeochemical tracers was taken from the climatology of the World Ocean Atlas 2001 (Conkright et al.2002; Locarnini et al.2002) and the Global Ocean Data Analysis Project (Key et al.2004). The initial iron concentration was set to a constant value of 0.6 nmol. The model is initialized with values of atmospheric δ13C and Δ14C of −6.5 ‰ and 0 ‰, respectively.

Notably, the biogeochemical cycle simulation of the LGM ocean performed in this study did not include certain processes such as enhanced Southern Ocean stratification, iron fertilization from glaciogenic dust, and carbonate compensation, as discussed in Kobayashi et al. (2021). These three processes are found to contribute to the glacial reduction in atmospheric pCO2, with values of 294.7 ppm for the preindustrial run (PI_sed in Kobayashi et al.2021) and 217.4 ppm for the LGM run (LGM_all in Kobayashi et al.2021). Therefore, the calculated atmospheric pCO2 during the LGM is expected to be higher than the atmospheric pCO2 reported in Kobayashi et al. (2021). Further experiments with the ocean general circulation model are needed to obtain a physical ocean field that accounts for the enhanced stratification of the Southern Ocean during glacial periods. Therefore, it was not possible to include this process in the model settings adopted for this study.

2.2.2 Transient experiment on the ocean carbon cycle during the last deglaciation

A transient experiment was conducted to investigate the ocean carbon cycle during the last deglaciation, starting from the initial state of the LGM (21 ka BP). Figure 1 illustrates the temporal variations in temperature and AMOC of the AOGCM output imposed as forcing of the ocean biogeochemical cycle model. During the transient experiment, dust deposition was periodically adjusted every 100 years based on the scaling between the LGM and the Holocene, using the reconstructed dust deposition from the Dome Fuji ice core (Dome Fuji Ice Core Project Members2017). Notably, the transient experiment did not account for temporal variations in ocean volume caused by ice sheet changes and associated changes in mean ocean concentrations of biogeochemical tracers (i.e., nutrients, alkalinity, and DIC) (Lhardy et al.2021). This simplification must be revisited in future studies.

Figure 1(a) Deglacial changes in the Atlantic meridional overturning circulation (AMOC; Sverdrup) with 231Pa/230Th reconstructed from Bermuda Rise sediment core data (McManus et al.2004) and compiled sediment core data from the North Atlantic (Ng et al.2018). The strength of the AMOC is defined as the maximum meridional volume transport between 30 and 90° N at depths below 500 m. (b) Deglacial changes in sea surface temperature (SST) in the Southern Ocean (°C) and the difference in SST from the present day (ΔTsource) (Uemura et al.2018). (c) Deglacial changes in global mean ocean temperature (MOT; °C) and the difference in MOT from the present day (“Mix” of Bereiter et al.2018). The model output computed by the AOGCM (OA19) (Obase and Abe-Ouchi2019), shown by blue lines, is compared to reconstructions from geological data, shown by gray lines. The right axes relate to the model output, and the left axes relate to the reconstructions.


Upon completion of the LGM ocean spin-up, we established a restoring term to counteract drifts in δ13C and Δ14C caused by gas exchange between the atmosphere and the ocean. This restoring term includes the exchange of carbon isotopes between the atmosphere and the land, together with the production of 14C in the atmosphere. This restoring term was assumed constant throughout the deglaciation experiment.

3 Results

3.1 Ocean carbon cycle state during the LGM

Figure 2 shows the calculated variations in Δ14C–CO2, δ13C–CO2, and pCO2 in the atmosphere (solid lines), together with their estimated values (dashed lines). Additionally, it illustrates the variations in the average of ΔΔ14C (i.e., the difference in Δ14C between the ocean and the atmosphere), δ13C, and DIC in the middle (500–2000 m) and deep (>2000 m) layers of the Atlantic, Southern (<40° S), and Pacific oceans.

Figure 2(a) Deglacial changes in Δ14C–CO2 (‰) in the atmosphere (solid line) with the reconstruction of IntCal20 (dashed line; Reimer et al.2020). (b–e) Deglacial changes in ΔΔ14C (‰), difference in Δ14C (‰) between the ocean and the atmosphere, averaged in the mid-depth (500–2000 m; solid blue lines) and deep global ocean (2000–5500 m; solid red lines) of the Atlantic Ocean (40° S–90° N), Pacific Ocean (40° S–90° N), and Southern Ocean (90–40° S) with compiled sediment core data (dashed lines) of Rafter et al. (2022). (f) Deglacial changes in δ13C–CO2 (‰) in the atmosphere (solid line) with the reconstruction (dashed line; Schmitt et al.2012). (g–i) Same as (b)(e), respectively, except for δ13C (‰), with compiled sediment core data of Muglia et al. (2023). (k) Deglacial changes in atmospheric pCO2 (ppm; solid line) with ice core data (dashed line; Bauska et al.2021). (l–o) Same as (b)(e), respectively, except for dissolved inorganic carbon (DIC; mmol m−3). LGM: Last Glacial Maximum. HS1: Heinrich Stadial 1. BA: Bølling–Allerød period. YD: Younger Dryas period. The triangles represent the values reported in Kobayashi et al. (2021), with the output of PI_sed plotted at the time of 11 ka BP and the output of LGM_all plotted at the time of 21 ka BP. PI_sed is an ocean carbon cycle model experiment conducted under preindustrial forcing, including carbonate sedimentation processes. LGM_all is an ocean carbon cycle model experiment conducted under LGM forcing, including enhanced salinity stratification in the Southern Ocean, iron fertilization from glaciogenic dust, and carbonate sedimentation processes.


Analysis of the modeled differences between the LGM and the Holocene indicates that the AMOC is weaker (9.0 Sv) during the LGM (21 ka BP) than during the Holocene (11 ka BP), i.e., 17.4 Sv (Figs. 1a and S1 in the Supplement). The basin-wide distributions of ΔΔ14C and δ13C in the Atlantic and Pacific oceans during specific periods are presented in Figs. 3 and 4, respectively. The model results demonstrate that ΔΔ14C, an indicator of ocean ventilation, is lower in the deep Atlantic at 21 ka BP than at 11 ka BP (Figs. 2c and 3a, f), which is in qualitative agreement with reconstructions from sediment core records (Rafter et al.2022). Notably, however, the simulated ΔΔ14C values are less negative than the reconstruction ΔΔ14C values at 21 ka BP in the Atlantic below 3000 m and in the Pacific below 2000 m (Fig. 3a).

Figure 3Oceanic zonal mean distribution of Δ14C (‰), which represents the difference in Δ14C between the ocean and the atmosphere, during key periods of the last deglaciation in the Atlantic and Pacific oceans. The specific periods of interest include (a) the Last Glacial Maximum (21 ka BP), (b) Heinrich Stadial 1 (17 ka BP), (c) just before the Bølling–Allerød (BA) transition (15 ka BP), (d) the BA warm period (13 ka BP), (e) the Younger Dryas period (12 ka BP), and (f) the Holocene (11 ka BP). The contour interval is 40 ‰. The sediment core records used in the figure are compiled in Rafter et al. (2022). Model results are averaged over 200 years, i.e., 100 years before and after each target year. However, for 21 ka BP, the average is taken from results spanning 21.0–20.9 ka BP; for 11 ka BP, the average is taken from results spanning 11.1–11.0 ka BP. The figure also includes a compilation of sediment core records where reconstructed values are plotted for 250 years before and after each target year. The vertical section represents all data within the relevant ocean basins. The abbreviations for the oceans are ATL for Atlantic Ocean, SO for Southern Ocean, and PAC for Pacific Ocean. The top-right notes indicate the number of data points, model–data correlation coefficients, and the root-mean-square error (RMSE) for both the Atlantic and the Pacific basins.

Figure 4Oceanic zonal mean distribution of δ13C (‰) during the last deglaciation in the Atlantic and Pacific oceans. The contour interval is 0.25 ‰. The sediment core records used in the figure are compiled in Muglia et al. (2023). The period over which the model output is averaged is 100 years before and after the year of interest. The reconstructed values are also plotted for 100 years before and after the year of interest.

Regarding δ13C, the model results show a stronger vertical gradient between the surface and the deep ocean in the Atlantic (Fig. 2h), corresponding to a shallower and weaker AMOC at 21 ka BP compared to that at 11 ka BP. This qualitative difference is consistent with the reconstructed δ13C. However, similar to the finding for ΔΔ14C, our model experiment underestimates the reconstruction of δ13C in the deep ocean, particularly in the Southern Ocean (Fig. 2i).

Our previous study, which involved numerical experiments under the climatic conditions of the LGM and accounted for the enhanced stratification in the Southern Ocean and iron fertilization from glaciogenic dust, showed improved quantitative agreement between the model results and the sediment core data for dissolved oxygen, δ13C, and Δ14C (Kobayashi et al.2021). The triangles in Fig. 2 illustrate the changes in the carbon isotope signatures reported in that study. Comparison of the results from the two studies highlights the advances made by Kobayashi et al. (2021) in capturing the dynamics of the Southern Ocean, suggesting that the incorporation of the processes considered in their research could improve model–data agreement. However, their LGM simulation also slightly overestimated the changes in the glacial Pacific, and these discrepancies highlight the difficulty in achieving consistent scenarios that account for all changes in the global ocean within a model.

Atmospheric pCO2 was predicted by running the ocean carbon cycle model. Its value is 278.1 ppm at 21 ka BP and 306.9 ppm at 11 ka BP, i.e., a difference of 28.8 ppm that is relatively small compared to the difference of approximately 80 ppm reconstructed from ice cores (approximately 188–267 ppm in the EPICA Dome C record (Bereiter et al.2015) and approximately 193–273 ppm in the WAIS Divide record (Bauska et al.2021)). This discrepancy could be attributed to several factors, including the relatively small differences in sea surface temperature (SST) observed between the two intervals (Fig. 1b). Using proxy and data assimilation, the global mean SST difference between the LGM and the Holocene has been reported to be 1.7–3.6 °C (MARGO Project Members2009; Tierney et al.2020; Paul et al.2021; Annan et al.2022), whereas the SST difference in our experiment is only 1.6 °C. A small difference in SST leads to a small difference in CO2 solubility between the two periods, which causes underestimation of the magnitude of atmospheric pCO2. Furthermore, it is important to note that this study did not consider specific processes that might have contributed to the reduction in atmospheric pCO2 during the LGM, such as enhanced salinity stratification, iron fertilization from glaciogenic dust, and carbonate compensation, as discussed in Kobayashi et al. (2021).

The differences in the steady-state ocean carbon cycles at 21 and 11 ka BP highlight the difficulty in accurately reproducing the actual changes in atmospheric pCO2 in the transient experiment connecting these periods. Therefore, our analysis focuses on investigation of the impacts of climate change, particularly the notable variations in the AMOC and on the ocean carbon cycle, and reveals the successes and deficiencies of the model through model–data comparison of carbon isotope signatures.

3.2 Carbon isotope changes during the last deglaciation

3.2.1 Radiocarbon isotopes

Here, we present the calculated transient changes in the ocean carbon cycle during the last deglaciation. During the deglaciation, atmospheric Δ14C decreases from approximately 500 ‰ to 0 ‰ (Reimer et al.2020; Fig. 2a). Seawater ΔΔ14C generally increases from the relatively low LGM values but decreases during HS1 and the YD period (Fig. 2b–e). Generally, the state of the AMOC strongly influences the distribution of ΔΔ14C. Both the model and the sediment core records show an increase in ΔΔ14C in the deep ocean during periods when the AMOC is relatively strong, e.g., the BA period and the Holocene (Fig. 3d and f), and present a decrease in ΔΔ14C during periods when the AMOC is relatively weak, e.g., HS1 and the YD period (Fig. 3b and c). The calculated changes in ΔΔ14C are consistent with the pattern observed in the sediment core record during a period characterized by rapid change in the AMOC after the BA transition (14.7 ka BP).

In the Pacific, when the AMOC is strong (i.e., at 13 and 11 ka BP), the ΔΔ14C in the South Pacific is relatively high, with elevated values from the surface to the deep Southern Ocean along the path of the Antarctic Bottom Water (AABW). The sediment core records compiled in Rafter et al. (2022) suggest elevated ΔΔ14C values in the North Pacific intermediate layers during HS1, possibly indicating the influence of the NPIW intrusion. However, the model results do not a provide clear indication of the intrusion of young water masses (Figs. 3 and S4).

The compiled sediment core records globally show a substantial increase in ΔΔ14C during HS1. In contrast, the model experiment does not show such pronounced change (Fig. 2b–e). This discrepancy can be attributed to two main factors: failure of the model experiment to simulate activation of ocean ventilation during HS1 and the greater magnitude of the initial ΔΔ14C values at 21 ka BP relative to the reconstructed values. This is supported by the insufficient carbon sequestration in the ocean calculated during the LGM (Fig. 2k). In other words, considering the glacial–interglacial redistribution of carbon in the atmosphere–ocean system, the relative abundance of 14C to 12C in the atmosphere is higher during the ice age, as manifested in the atmospheric Δ14C–CO2 (Fig. 2a); however, the variation in 12C is not well reproduced in the model experiment (Fig. 2k).

Regarding the latter point of insufficient carbon sequestration during the LGM, the triangles shown in Fig. 2 represent the results of the best LGM simulation (LGM_all) conducted by Kobayashi et al. (2021). That simulation incorporated enhanced salinity stratification and sedimentation processes, which further contribute to accurate reproduction of low ΔΔ14C of deep water during the LGM. As shown in Fig. 2b–e, there is a substantial discrepancy between the model and the reconstruction, particularly in relation to the Southern Ocean during the period of early deglaciation. Incorporation of change in vertical mixing resulting from variation in ocean stratification could potentially improve the simulation of Δ14C in the deglaciation.

3.2.2 Stable carbon isotopes

Next, we focus on the changes in δ13C. For seawater δ13C, the overall trends of change in δ13C and ΔΔ14C are similar and correspond to phases of climatic change; however, δ13C is less sensitive than ΔΔ14C to climate change (Fig. 2). This differential response might be related to biological fractionation of carbon isotopes. A more active AMOC leads to increased biological activity that reduces δ13C in deeper layers, especially in the North Atlantic. This counteracts the influence of lighter carbon transported to the surface by the active AMOC.

During HS1, δ13C decreases gradually in the upper 3000 m of the North Atlantic (Figs. 2h, 4, and S3). The reduction in δ13C can be attributed to several factors (Gu et al.2021) that include increased contribution from southern-sourced deep water with low δ13C endmembers, accumulation of remineralized carbon with low δ13C attributable to a weakened AMOC and reduced ventilation, and potential increase in the δ13C endmember of North Atlantic Deep Water (NADW). The results of this study show no clear change in the NADW endmembers of δ13C (Fig. S5). Therefore, the change in δ13C is attributed to weakened ventilation in the North Atlantic and to expansion of southern-sourced deep water. However, the observed δ13C change is relatively small compared to that derived from sediment core data because the AMOC change is less pronounced than that expected from the 231Pa/230Th reconstruction (McManus et al.2004; Ng et al.2018).

The deep Southern Ocean has its lowest δ13C during the LGM, although the value gradually increases during HS1. However, these observed changes are not reproduced in the model. According to Kobayashi et al. (2021), the low δ13C during the LGM is related to enhanced Southern Ocean stratification and iron fertilization from glaciogenic dust. These processes, which are not considered in this study, contribute to the differences between the model and the observed data.

During the BA period (Fig. 4d) and the Holocene (Fig. 4f), the intensified and deepened AMOC contributes to high δ13C values originating from the North Atlantic penetrating to depths below 2000 m. Basin-averaged δ13C reconstructions also show this increase in δ13C, especially in the Atlantic (Fig. 2g–j). However, in contrast to the calculated change, the sediment core records do not show further reduction in δ13C in the deep ocean during the YD period.

Changes in δ13C in the deep ocean lead to changes in atmospheric δ13C–CO2. There is a sharp drop in atmospheric δ13C–CO2 during HS1, followed by a slight rise during the BA period and then a further decline during the YD period (Schmitt et al.2012; Fig. 2f). However, the model-calculated changes in the ocean carbon cycle do not reproduce this trend in atmospheric δ13C–CO2. The reconstructed δ13C–CO2 increases during the BA period, decreases during the YD period, and then increases again, whereas the model-calculated trend is the opposite. The discrepancy might involve the contribution from changes in vegetation, which is a topic discussed in Sect. 4.3.

Isotope fractionation through temperature-dependent gas exchange and phytoplankton preference for uptake of lighter carbon also plays an important role in δ13C variation. Figure S6 shows the calculated changes in organic carbon export from that in 21 ka BP with qualitative changes in biological flux reconstructed from proxies, specifically opal flux and alkenone flux, in sediment core records (Chase et al.2003; Anderson et al.2009; Bolton et al.2011; Kohfeld and Chase2011; Martínez-García et al.2014; Maier et al.2015; Studer et al.2015; Thiagarajan and McManus2019; Ai et al.2020; Weber2021; Li et al.2022). During HS1, both the model and the proxies show increased biological carbon transport in the polar region of the Southern Hemisphere (Fig. S6a–d). In the polar regions, sea ice is reduced owing to warming, which could result in less light limitation and allow increased biological productivity. However, biological carbon transport is reduced in subpolar regions and in the South Pacific gyres. These changes in southern regions can be attributed to reduced nutrient supply resulting from weakening of the AMOC. Another important factor is the increase in iron limitation associated with the reduced supply of dust-derived iron that affects biological production. During the BA warm period, the enhanced AMOC enhances nutrient transport from the deep ocean to the surface ocean, resulting in increased biological transport in the North Atlantic (Fig. S6e and 6f). These changes in the vertical nutrient transport then propagate to the North Pacific. From the BA period to the YD period, there is an increase in biological export in the Southern Ocean that might be attributable to reduction in sea ice resulting from warming in the Southern Hemisphere. The factors that alter biological production in the model are understood but need to be constrained using additional proxy data with high temporal resolution that can capture millennial-scale variations.

3.3 Deglacial changes in atmospheric pCO2 caused by changes in the ocean carbon cycle

Figure 2k shows the calculated changes in atmospheric pCO2 driven by variations in the climate and the carbon cycle during the last deglaciation. To investigate the factors driving the change in atmospheric pCO2, we decomposed the factors relevant to the partial pressure of CO2 at the sea surface (pCO2os). This parameter controls atmospheric pCO2 through gas exchange between the atmosphere and the ocean. Oceanic pCO2 is affected by temperature, salinity, DIC, and alkalinity, and the influence of those factors on pCO2os can be represented as follows:

(1) p CO 2 os = f ( sDIC , sALK , SST , SSS ) ,

where sDIC is sea surface DIC, sALK is sea surface alkalinity, SST is sea surface temperature, and SSS is sea surface salinity. The function f is determined based on the inorganic chemistry of the carbonate system (Millero1995). We can assess the contribution of each variable to the changes in pCO2os by examining the change in each variable from its original value.

3.3.1 Heinrich Stadial 1

During HS1, the calculated atmospheric pCO2 rises slightly until approximately 17 ka BP, and then it rises sharply to approximately 15 ka BP. From 18 to 15 ka BP, atmospheric pCO2 increases by 10.2 ppm (Fig. 5a), whereas the WAIS Divide ice core record shows an increase of 41.4 ppm during the same period (Bauska et al.2021). The model accounts for approximately one-quarter of the reconstructed changes in atmospheric pCO2. Decomposition analysis of pCO2os reveals that most of the variation in pCO2os is driven by change in SST (Fig. 5a and b). The changes in pCO2os (ΔpCO2os) attributable solely to variations in temperature and salinity (ΔpCO2os(TS)) and in DIC and alkalinity (ΔpCO2os(CA)) are shown in Fig. 6a and b, respectively. It is evident that ΔpCO2os(TS) shows a predominantly positive contribution globally, reflecting the pattern of SST increase (Fig. 6d), because increasing SST reduces CO2 solubility. In other words, the main contributor to the increase in pCO2os during HS1 is warming, especially in the subantarctic region.

Figure 5(a) Temporal changes in the partial pressure of sea surface CO2 (pCO2os, parts per million, gray) during Heinrich Stadial 1 (differences between 18 and 15 ka BP). The contributions of changes in temperature and salinity (purple), temperature (red), salinity (yellow), dissolved inorganic carbon (DIC) and alkalinity (cyan), DIC (green), and alkalinity (blue) to the changes in pCO2os are shown. The thin gray line shows the time series of AMOC strength. (b) Temporal changes in the partial pressure of atmospheric pCO2 and pCO2os (parts per million) during Heinrich Stadial 1. The contributions of changes in temperature and salinity, temperature, salinity, DIC and alkalinity, DIC, and alkalinity to the changes in pCO2os are represented by different-colored bars. (c, d) Similar to panels (a) and (b), respectively, but for the Bølling–Allerød period (differences between 15 and 13 ka BP). (e, f) Similar to panels (a) and (b), respectively, but for the Younger Dryas period (differences between 13 and 12 ka BP).


Figure 6(a) Changes in partial pressure of sea surface CO2 (pCO2os, parts per million) between the early and late Heinrich Stadial 1 (differences between 15 and 18 ka BP). Changes in pCO2os attributable solely to changes in (b) temperature and salinity and (c) dissolved inorganic carbon (DIC) and alkalinity and to changes in (d) sea surface temperature, (e) DIC, and (f) alkalinity between the same periods.

3.3.2 Bølling–Allerød period

At the onset of the BA transition near 14.7 ka BP, atmospheric pCO2 begins to decrease, and this reduction continues until 12.8 ka BP (Fig. 2k). From 15 to 13 ka BP, atmospheric pCO2 decreases by 7.0 ppm (Fig. 5c). During the BA period, the contributions of thermal changes and biogeochemical changes act in opposition to the change in atmospheric pCO2 (Fig. 5c and d). At the onset of the BA period, as the AMOC strengthens (Fig. 1a), both SST and SSS increase in the Northern Hemisphere and decrease in the Southern Hemisphere (Fig. 7d). The net contribution of pCO2os(TS) attributable to changes in CO2 solubility is positive (Fig. 5c and d). However, the enhanced AMOC facilitates the transport of nutrients, carbon, and alkalinity from the deep ocean to the surface ocean, especially in the North Atlantic (Fig. 7e and f). The increase in sDIC leads to an increase in pCO2os, while the increase in sALK leads to a reduction in pCO2os. These opposing effects partially offset each other, resulting in a net reduction in pCO2os (Figs. 5c, d and 7c). During the AMOC overshoot at the BA transition and the subsequent stabilized phase, increased biological production in most of the global ocean contributes to a millennial-scale decrease in sDIC, resulting in a reduction in pCO2os (Fig. 5c and d). In summary, following the recovery of the AMOC, the opposing contributions of ΔpCO2os(TS) and ΔpCO2os(CA) to ΔpCO2os over time control the temporal changes in pCO2os and the subsequent reduction in atmospheric pCO2.

Figure 7Same as Fig. 6 except for the Bølling–Allerød period (differences between 13 and 15 ka BP).

3.3.3 Younger Dryas period

Atmospheric pCO2 rises again at the onset of the YD period (12.8 ka BP), coinciding with the collapse of the AMOC into a weak state (Fig. 1a). From 13 to 12 ka BP, atmospheric pCO2 increases by 6.8 ppm (Fig. 5e). Decomposition analysis of pCO2os reveals that the influences of ΔpCO2os(TS) and ΔpCO2os(CA) on the overall ΔpCO2os are in opposition, and this offset is also observed during the BA period but in the opposite sense. The contribution of ΔpCO2os(TS) increases over time, leading to an increase in pCO2os (Fig. 5e and f). The contribution of ΔpCO2os(CA) is small. As the AMOC weakens, a decrease in sALK contributes to an increase in pCO2os, while a decrease in sDIC contributes to a decrease in pCO2os (Figs. 5e, f and 8c, e, f). However, the net effect of changes in DIC and alkalinity on pCO2os is minimal, resulting in only a slight decrease in ΔpCO2os(CA) during the YD period. From these opposing effects, the overall changes in ΔpCO2os(TS) and ΔpCO2os(CA) indicate an increase in ΔpCO2os during the YD period.

Figure 8Same as Fig. 6 except for the Younger Dryas period (differences between 12 and 13 ka BP).

4 Discussion

The objective of this study was to investigate the transient response of the ocean carbon cycle during the last deglaciation. By comparing the calculated carbon isotope signatures of δ13C and Δ14C with those derived from sediment core records, we can assess the impacts of changes in climate and the AMOC on those signatures. This comparison can also present information to help identify potential biases or missing processes within the model. In addition to changes in atmospheric pCO2, investigating the mechanisms behind changes in carbon isotopes contributes to a more comprehensive understanding of the temporal changes in the global carbon cycle.

4.1 Response of carbon isotope signatures to drastic changes in the deep-ocean circulation

Comparison of ΔΔ14C variations between models and data enables assessment of the accuracy of calculated ocean circulation changes. The reconstructed ΔΔ14C in the deep ocean rises notably during the latter half of HS1 (Rafter et al.2022), in contrast to the less pronounced shift seen in the model experiment (Fig. 2b–e). Two primary factors contribute to this discrepancy. Firstly, the model underestimates the increase in deep-ocean ventilation during the latter half of HS1 period, which is crucial for determining the trend in Δ14C changes in the deep ocean. Secondly, the model calculates higher ΔΔ14C values compared to the reconstructed lower values of ΔΔ14C in the deep ocean during the LGM. These aspects indicate a potential oversight in the model's representation of ocean dynamics, affecting both the simulation of ventilation changes during the latter half of HS1 and ΔΔ14C values during the LGM. The difficulty in reproducing changes in ΔΔ14C might also be related to underestimation of variations in atmospheric pCO2 during this period (Fig. 2k). Processes that might contribute to this problem are discussed in more detail in Sect. 4.2.

After the BA transition, the calculated variations in ΔΔ14C in synchrony with the significant changes in the AMOC are generally consistent with those observed in the reconstruction (Rafter et al.2022). Corresponding to the AMOC change, ΔΔ14C increases in the deep ocean from the Atlantic to the Pacific Ocean (Fig. 2c–e). Subsequently, ΔΔ14C decreases from the Atlantic to the Southern Ocean during the YD period in response to weakening of the AMOC (Fig. 2c–e). This change is consistent with the reconstruction (Rafter et al.2022), but there is an overestimation of the quantitative changes in the deep Atlantic Ocean, as illustrated in Fig. 3e. This overestimation may be related to the challenges in accurately reproducing the deep-ocean circulation fields, which is a topic that is discussed further below.

Assessment of the AMOC changes during the last deglaciation by Pöppelmeier et al. (2023) involved conducting transient model simulations using the Bern3D model. They performed multiple model–data comparisons including carbon isotope ratios, εNd, and 231Pa/230Th. Their research results suggest gradual weakening of the AMOC during HS1, recovery at the BA transition, and subsequent weakening during the YD period, albeit without a complete collapse. The proposed pattern of AMOC change is qualitatively consistent with the ocean modeling of Obase and Abe-Ouchi (2019) (Fig. 1a). The Bern3D study indicates that the proportion of deep water originating from the North Atlantic during the YD period is little different to that during the BA because of the short duration of the YD period. Conversely, this study shows drastic changes in the distribution of Δ14C and δ13C at the basin scale in response to AMOC variations over a period of approximately 1000 years (Figs. 3e and 4e). This extended period of the AMOC stagnation in our model helps to explain the observed discrepancies between the model and the reconstructed data (Rafter et al.2022; Muglia et al.2023) in the deep Atlantic during the YD period.

The carbon isotope ratios calculated by our model suggest that the AMOC during the YD period might be represented as excessively weak or that the duration of the weak AMOC state may be overly extended. However, it is important to acknowledge the inherent differences among models in representing broader deep-ocean circulation patterns, including the AABW and Pacific meridional overturning. These differences result in distinct chemical tracer distributions at the basin scale, underscoring the challenge of conclusively explaining past AMOC variations through a single model or study. Given the systematic biases present in physical and biogeochemical processes within models, a comparative analysis across multiple models is essential for exploring past AMOC variations.

Further information can be obtained by comparing the results of model–data comparisons for ΔΔ14C and δ13C. For δ13C, both the model and the data show similar trends, depicting an increase in deep-water δ13C during the BA period as in ΔΔ14C. However, there is a discrepancy during the subsequent YD period. The model indicates a reduction in deep-water δ13C during the YD period, whereas this feature is absent in the reconstruction of Muglia et al. (2023) (Fig. 2g–j). There are several possible factors that could potentially influence this discrepancy. For example, discrepancies in the directions and magnitudes of changes observed across different sediment cores could reflect inherent variability in environmental signals. Additionally, potential dating inaccuracies within individual sediment core data could result from smoothing effects such as bioturbation and coring artifacts. These complexities in interpreting the sediment core record stem from both natural variability and methodological challenges, highlighting the need for caution when comparing model simulations with sediment core data. Another important factor is the weakening of the simulated AMOC during the YD period. Comparison of the model and sediment data for ΔΔ14C and δ13C suggests that the weakening of the AMOC during the YD period might be overly pronounced in the model (Figs. 3e and 4e). Additionally, the calculated increase in export of biogenic organic matter in the Southern Ocean during the YD period compared to that in the BA period (Fig. S6f and g) contributes to the decrease in δ13C in the deep ocean. This emphasizes the importance of accurately simulating nutrient and iron cycles, especially in iron-limited regions affected by changes in dust-derived iron supply. As the Southern Hemisphere warms and becomes more humid, the supply of iron from dust might decrease (Martin1990; Martínez-García et al.2014). Reproducing changes in δ13C is challenging owing to the intricate interconnections between ocean circulation, biological processes, and atmosphere–ocean gas exchange. Understanding the discrepancies between the model and the data in terms of the δ13C changes will require future sensitivity experiments to clarify their respective contributions and to provide deeper understanding of these factors.

4.2 Insights from carbon isotope ratios: oceanic CO2 release during the deglaciation

Ocean modeling with freshwater-forcing experiments has provided insights into the link between the shutdown and resumption of the AMOC and the changes in atmospheric pCO2. Schmittner and Galbraith (2008) demonstrated that cessation of the AMOC causes increase in atmospheric pCO2 owing to several factors. Firstly, the efficiency of the biological carbon pump in the North Atlantic is relatively high compared to that in the Southern Ocean. Therefore, reduction in the NADW inflow leads to a decrease in biological carbon sequestration in the deep ocean. Secondly, weakening of Southern Ocean stratification associated with shutdown of the AMOC increases the outgassing of CO2 from the ocean to the atmosphere.

The results of this study confirm the gradual increase in atmospheric pCO2 during HS1 and the YD period in parallel with the weakened state of the AMOC; however, such an increase is not directly related to reduction in the regenerated nutrient inventory, as suggested by Schmittner and Galbraith (2008). The reason for this difference is that the contribution of the changes in temperature and alkalinity to the change in pCO2os during the YD period is greater than the contribution of the change in DIC in this study. When the AMOC changes, the non-thermal effects on changes in pCO2os mainly depend on the magnitude of the relative contributions of DIC and alkalinity to pCO2os based on the vertical gradient of DIC and alkalinity between the surface and the deeper ocean (Figs. S7 and S8).

Although our results show the impact of drastic changes in the AMOC on atmospheric pCO2 during the last deglaciation, the model does not fully explain the variations in atmospheric pCO2 during the early deglaciation. Ice core records indicate a rise of approximately 40 ppm in atmospheric pCO2 accompanied by a reduction in atmospheric δ13C–CO2 and Δ14C–CO2 during HS1. However, the calculated variations are insufficient in terms of their amplitude (Fig. 2a, f, and k). A key contributor to this discrepancy is the limited extent of the change in ventilation, evident from the ΔΔ14C changes in the ocean (Fig. 2b–e). However, the current model has difficulty fully reproducing these substantial changes in ventilation.

In addition to the changes in ventilation, biological processes are important in explaining the deglacial carbon cycle changes. Kobayashi et al. (2021) indicated that iron fertilization from glaciogenic dust increases biological production in the subantarctic region, thereby contributing to the reproduction of low δ13C in the deep Southern Ocean during the LGM (triangles in Fig. 2i). However, their study did not reproduce the profoundly old deep water in the deep glacial Southern Ocean, as suggested by radiocarbon data (triangles in Fig. 2d), but it did reproduce the low δ13C values (triangles in Fig. 2i). The change in glaciogenic dust deposition was not considered in this study; therefore, the model might underestimate the glacial–interglacial variation in biological production in the subantarctic region, potentially contributing to the underestimation of the changes in atmospheric and deep-sea δ13C during HS1.

Moreover, while many studies focused primarily on environmental changes in the Atlantic, the contributions from other ocean basins are also important. For example, sediment core records from the North Pacific indicate an increase in ΔΔ14C at depths near 1000 m during HS1 (Okazaki et al.2010; Rae et al.2014; Rafter et al.2022). Analysis of radiocarbon and boron isotopes in sediment cores by Rae et al. (2014) revealed that the extent of the NPIW expanded during HS1. These ventilation changes have the potential to contribute to the rise in atmospheric pCO2. Chikamoto et al. (2012) compared two coupled climate models, i.e., MIROC version 3.1 and LOVECLIM, and showed consistent results for activation of ventilation of NPIW triggered by freshwater inflow into the North Atlantic. The processes of activation of ocean ventilation and subsequent degassing of CO2 in the North Pacific during the early deglaciation could contribute to the currently unexplained increase in atmospheric pCO2.

The high-resolution ice core data obtained from the WAIS Divide record provides valuable insights into the timescale of the changes in the carbon cycle during the last deglaciation. It is suggested that there are two modes of change in relation to atmospheric pCO2: a slow increase on the millennial scale and a rapid increase on the centennial scale (Marcott et al.2014). The rapid increase in atmospheric pCO2 of 10–15 ppm at the end of HS1 (14.8 ka BP) and at the end of the YD period (11.7 ka BP) over a short period of 100–200 years is synchronized with the resumption of the AMOC (McManus et al.2004; Ng et al.2018). However, this study did not reproduce such abrupt changes in atmospheric pCO2. Several factors might contribute to this discrepancy, including insufficient temperature rise in the Southern Hemisphere associated with change in the AMOC, inadequate representation of the vertical concentration gradients of DIC and alkalinity, limitations in capturing atmospheric and oceanic dynamics in the general circulation model, and the influence of small-scale phenomena. Previous modeling studies have suggested that deeper convection in the Southern Ocean and strengthening of westerly winds in the Southern Hemisphere could contribute to the abrupt jump in atmospheric pCO2 during the middle of HS1 (16.3 ka BP) by transporting sequestered carbon from the deep Southern Ocean to the surface (Menviel et al.2018). These processes are related to the challenges in reproducing carbon isotope ratios described above; therefore, this discussion points to the necessity of improving our AOGCM and of refining its experimental setup in future studies, as previewed in Sect. 4.3.

4.3 Improvement of the model and the experimental design: future implications

Model–data comparisons of carbon isotope signatures underscore the importance of refining our climate models to more accurately represent the complex interactions that govern changes in the carbon cycle. Future improvements to the AOGCM used in this study might address the following considerations.

Currently, the AOGCM does not account for temporal changes in ice sheets, and it underestimates the changes in Southern Ocean SST (Obase and Abe-Ouchi2019). Moreover, there is some uncertainty regarding the volume of meltwater flow across the North Atlantic during HS1 (Ivanovic et al.2018; Snoll et al.2023) and the BA period (Kapsch et al.2022; Bouttes et al.2023). It suggests that the problem is integral to the AOGCM because the AMOC response is not consistently and realistically observed, even when realistic freshwater variations are applied. Furthermore, as identified by Obase et al. (2023), the magnitude of ocean warming during the last deglaciation varies depending on the response characteristics of each model, resulting in a range across multiple models. Sherriff-Tadano et al. (2023) demonstrated that the alteration of parameters associated with cloud thermodynamic phase fractions in a climate model reduces the warming bias of SST in the modern Southern Ocean. Using a model with a reduced Southern Ocean warming bias, we expect to obtain different responses in Southern Ocean SST and ocean circulation during the glacial period and in their changes during the deglaciation compared to those derived in this study. Some studies proposed potential alterations in the westerly winds over the Southern Ocean throughout the last deglaciation (Gray et al.2023), but there is a substantial degree of uncertainty concerning the anticipated changes in the atmospheric dynamics. Those uncertainties could have a substantial impact on the results of ocean biogeochemical cycle modeling and, consequently, on atmospheric pCO2. Efforts to reduce bias and to facilitate comprehensive discussion regarding climate model consistency are critical to advancing future climate–carbon cycle modeling. These endeavors are essential to refine our understanding of the complex interactions between climate and the carbon cycle.

In addition to those factors mentioned above, there are several other factors that could contribute to improving the simulation of the ocean carbon cycle during the last deglaciation. One important consideration is the inclusion of critical processes for lowering atmospheric pCO2 during the LGM, which are identified in Kobayashi et al. (2021). Those processes include enhanced stratification of the Southern Ocean, iron fertilization from glaciogenic dust, and carbonate compensation. Understanding the changes in those processes during the deglaciation is critical, and their proper incorporation into future modeling efforts might lead to both improved simulations and better understanding of the dynamics during this period. Studies have also reported that inclusion of the parameterization of vertical mixing, which depends on tidal mixing energy and stratification, could help better reproduce the deep-ocean circulation in the Pacific (Oka and Niwa2013; Kawasaki et al.2022). The introduction of tidal mixing parameterization has also proven effective in reproducing δ13C and Δ14C in the glacial ocean (Wilmes et al.2021). Incorporating the insights from these model developments has the potential to lead to more realistic representation of carbon cycle variations during the glacial period and subsequent deglaciation. Moreover, glacial–interglacial changes in ocean volume due to ice sheet changes also have an impact on the carbon cycle. A recent study analyzing PMIP model outputs highlighted the importance of accurate representation of ocean volume changes and their associated effects on alkalinity adjustments (Lhardy et al.2021). For more accurate simulations, it is critical to perform numerical integration that accounts for temporal changes in ocean volume during the deglaciation.

It is also worth noting that carbon exchange between the atmosphere and the ocean is not the sole driver of deglacial variation in atmospheric pCO2. Changes in terrestrial and soil carbon storage also play important roles in modulating atmospheric δ13C–CO2 and pCO2 during the last deglaciation (Schmitt et al.2012; Jeltsch-Thömmes et al.2019; Bauska et al.2021). Comparison of Δ14C–CO2 and δ13C–CO2 provides a consistent explanation if carbon uptake by vegetation expands during the BA period and declines during the YD period, as suggested by Schmitt et al. (2012). Vegetation growth, prompted by CO2 fertilization, can act as a carbon sink, offsetting the increase in atmospheric pCO2 (Bouttes et al.2012a; Menviel et al.2012). In this study, we applied the same sea surface restoring terms for the carbon isotopes used during the initial spin-up of the LGM throughout the deglaciation experiment. In other words, this approach did not consider changes in vegetation or changes in 14C production in the atmosphere. Future studies using Earth system models that include both terrestrial and oceanic carbon cycle processes would enhance our comprehensive understanding of glacial changes in carbon cycles. A study of the carbon cycle associated with the glacial Dansgaard–Oeschger events, conducted using an Earth system model, revealed that the changes in terrestrial carbon storage at this timescale are as important as those in the oceans (Jochum et al.2022). Moreover, in previous studies using Earth system models of intermediate complexity, temporal variations in atmospheric pCO2 are primarily used to calculate changes in radiative forcing, and several studies have explored the interaction between the carbon cycle and the climate (Bouttes et al.2012a; Ganopolski and Brovkin2017). Although fully coupling a carbon cycle model to a climate model is a more advanced endeavor, we are eager to explore this avenue in future research.

5 Conclusions

To understand the mechanisms of glacial–interglacial variability in the carbon cycle, this study examined the transient response of the ocean carbon cycle to climate change, including the remarkable strengthening and weakening of the AMOC at the BA and YD transitions. This study represents an important step towards comprehensive transient simulations of the carbon cycle using an AOGCM, even though the changes in atmospheric pCO2 are relatively small compared to those derived from ice core reconstructions. The importance of this study lies in its model–data comparisons of carbon isotope ratios that elucidate the impact of AMOC mode changes on the three-dimensional structure of water masses in the Atlantic, Southern, and Pacific oceans.

Our model qualitatively simulates an increase in atmospheric pCO2 during HS1. The calculated increase of approximately 10 ppm in atmospheric pCO2 from 18 to 15 ka BP is mainly caused by an increase in SST. The relatively modest increase in atmospheric pCO2, compared to that of ice core records, might be attributable in part to relatively small increases in SST in the Southern Ocean. Additionally, comparison of carbon isotope signatures between the model and the data highlighted the scope for improvement in the representation of increased ventilation in the deep ocean and the North Pacific. Similarly, there is potential for improvement with respect to changes in surface biological productivity involving the nutrient cycle, including iron. Correction of these elements would substantially improve our understanding of the increase in atmospheric pCO2 during the early deglaciation.

The drastic shifts in the AMOC during the BA and YD periods cause bipolar climate changes. These changes affect not only the temperature and salinity distributions but also the distributions of DIC and alkalinity. Interestingly, the cumulative effects of these changes on atmospheric pCO2 appear to cancel each other out, resulting in only a slight decrease during the BA period and increase during the YD period. It is noticeable that changes in pCO2os due to variations in temperature and alkalinity play a major role in the reduction of atmospheric pCO2 during the BA period.

To simulate transient changes in the carbon cycle, improvements in model accuracy, experimental configurations, and the models themselves are critical for capturing the dynamical and biogeochemical changes in the atmosphere and ocean. Further research is needed to identify the specific processes that influence changes in the ocean carbon cycle over different timescales in individual ocean basins. We emphasize the importance of analyzing carbon isotope variations that can provide valuable insights into past carbon cycle dynamics and contribute to a comprehensive understanding of the glacial–interglacial variations in the ocean carbon cycle.

Code and data availability

The CCSR Ocean Component Model (COCO) is the ocean general circulation model of MIROC, and the code of COCO version 4.0 is included as part of MIROCES2L. The source code of MIROC-ES2L can be obtained from (Ohgaito et al.2021b).


The supplement related to this article is available online at:

Author contributions

HK and AO designed the research. HK conducted the numerical experiments with help from TO. HK performed the analysis and wrote the paper. AO and AAO obtained funding and supervised the study. All authors discussed the results and commented on the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


The authors express sincere gratitude to the two anonymous reviewers for their invaluable and constructive feedback on our paper. We would also like to thank Laurie Menviel for their insightful and helpful comments and for their expert editorial handling. The ocean tracer model simulations in this study were performed at the Information Technology Center of the University of Tokyo. We thank James Buxton, MSc, from Edanz (, last access: 25 March 2024) for editing a draft of this paper. This study was supported by the Cooperative Research Activities of Collaborative Use of Computing Facility of the Atmosphere and Ocean Research Institute, the University of Tokyo.

Financial support

This research was supported by KAKENHI, the Japan Society for the Promotion of Science (grant nos. JP17H06104, JP17H06323, JP19H01963, and JP21K13990), PRESTO, the Japan Science and Technology Agency (grant no. JPMJPR23G4), and the Environment Research and Technology Development Fund (grant no. JPMEERF23S21109) of the Environmental Restoration and Conservation Agency provided by Ministry of the Environment of Japan.

Review statement

This paper was edited by Laurie Menviel and reviewed by two anonymous referees.


Ai, X. E., Studer, A. S., Sigman, D. M., Martínez-García, A., Fripiat, F., Thöle, L. M., Michel, E., Gottschalk, J., Arnold, L., Moretti, S., Schmitt, M., Oleynik, S., Jaccard, S. L., and Haug, G. H.: Southern Ocean upwelling, Earth's obliquity, and glacial-interglacial atmospheric CO2 change, Science, 370, 1348–1352,, 2020. a

Anderson, R. F., Ali, S., Bradtmiller, L. I., Nielsen, S. H. H., Fleisher, M. Q., Anderson, B. E., and Burckle, L. H.: Wind-Driven Upwelling in the Southern Ocean and the Deglacial Rise in Atmospheric CO2, Science, 323, 1443–1448,, 2009. a

Annan, J. D., Hargreaves, J. C., and Mauritsen, T.: A new global surface temperature reconstruction for the Last Glacial Maximum, Clim. Past, 18, 1883–1896,, 2022. a

Barnola, J. M., Raynaud, D., Korotkevich, Y. S., and Lorius, C.: Vostok ice core provides 160,000-year record of atmospheric CO2, Nature, 329, 408–414,, 1987. a

Bauska, T. K., Marcott, S. A., and Brook, E. J.: Abrupt changes in the global carbon cycle during the last glacial period, Nat. Geosci., 14, 91–96,, 2021. a, b, c, d

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549,, 2015. a

Bereiter, B., Shackleton, S., Baggenstos, D., Kawamura, K., and Severinghaus, J.: Mean global ocean temperatures during the last glacial transition, Nature, 553, 39?44,, 2018. a

Bolton, C. T., Lawrence, K. T., Gibbs, S. J., Wilson, P. A., and Herbert, T. D.: Biotic and geochemical evidence for a global latitudinal shift in ocean biogeochemistry and export productivity during the late Pliocene, Earth Planet. Sc. Lett., 308, 200–210,, 2011. a

Bouttes, N., Paillard, D., Roche, D. M., Waelbroeck, C., Kageyama, M., Lourantou, A., Michel, E., and Bopp, L.: Impact of oceanic processes on the carbon cycle during the last termination, Clim. Past, 8, 149–170,, 2012a. a, b, c, d

Bouttes, N., Roche, D. M., and Paillard, D.: Systematic study of the impact of fresh water fluxes on the glacial carbon cycle, Clim. Past, 8, 589–607,, 2012b. a

Bouttes, N., Lhardy, F., Quiquet, A., Paillard, D., Goosse, H., and Roche, D. M.: Deglacial climate changes as forced by different ice sheet reconstructions, Clim. Past, 19, 1027–1042,, 2023. a, b

Broecker, W. S. and Maier-Reimer, E.: The influence of air and sea exchange on the carbon isotope distribution in the sea, Global Biogeochem. Cy., 6, 315–320,, 1992. a

Brovkin, V., Ganopolski, A., Archer, D., and Munhoven, G.: Glacial CO2 cycle as a succession of key physical and biogeochemical processes, Clim. Past, 8, 251–264,, 2012.  a

Chase, Z., Anderson, R. F., Fleisher, M. Q., and Kubik, P. W.: Accumulation of biogenic and lithogenic material in the Pacific sector of the Southern Ocean during the past 40,000 years, Deep-Sea Res. Pt. II, 50, 799–832,, 2003. a

Chikamoto, M. O., Menviel, L., Abe-Ouchi, A., Ohgaito, R., Timmermann, A., Okazaki, Y., Harada, N., Oka, A., and Mouchet, A.: Variability in North Pacific intermediate and deep water ventilation during Heinrich events in two coupled climate models, Deep-Sea Res. Pt. II, 61–64, 114–126,, 2012. a

Conkright, M. E., Garcia, H. E., O'Brien, T. D., Locarnini, R. A., Boyer, T. P., Stephens, C., and Antonov, J. I.: World Ocean Atlas 2001, vol. 4, Nutrients, NOAA Atlas NESDIS, vol. 52, 392 pp., NOAA, Silver Spring, MD, 2002. a

Dome Fuji Ice Core Project Members: State dependence of climatic instability over the past 720,000 years from Antarctic ice cores and climate modeling, Sci. Adv., 3, e1600446,, 2017. a

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716,, 2017. a, b

Ganopolski, A. and Roche, D. M.: On the nature of lead–lag relationships during glacial–interglacial climate transitions, Quaternary Sci. Rev., 28, 3361–3378,, 2009. a

Gottschalk, J., Battaglia, G., Fischer, H., Frölicher, T. L., Jaccard, S. L., Jeltsch-Thömmes, A., Joos, F., Köhler, P., Meissner, K. J., Menviel, L., Nehrbass-Ahles, C., Schmitt, J., Schmittner, A., Skinner, L. C., and Stocker, T. F.: Mechanisms of millennial-scale atmospheric CO2 change in numerical model simulations, Quaternary Sci. Rev., 220, 30–74,, 2019. a

Gray, W. R., de Lavergne, C., Wills, R. C. J., Menviel, L., Spence, P., Holzer, M., Kageyama, M., and Michel, E.: Poleward Shift in the Southern Hemisphere Westerly Winds Synchronous With the Deglacial Rise in CO2, Paleoceanogr. Paleoclimatol., 38, e2023PA004666,, 2023. a, b

Gu, S., Liu, Z., Oppo, D. W., Lynch-Stieglitz, J., Jahn, A., Zhang, J., Lindsay, K., and Wu, L.: Remineralization dominating the δ13C decrease in the mid-depth Atlantic during the last deglaciation, Earth Planet. Sc. Lett., 571, 117106,, 2021. a

Hasumi, H.: CCSR Rep. 25, Chap. CCSR Ocean Component Model (COCO) Version 4.0, Center for Climate System Research, Univ. of Tokyo, Japan, 103 pp., (last access: 25 March 2024), 2006. a

He, F., Shakun, J. D., Clark, P. U., Carlson, A. E., Liu, Z., Otto-Bliesner, B. L., and Kutzbach, J. E.: Northern Hemisphere forcing of Southern Hemisphere climate during the last deglaciation, Nature, 494, 81–85,, 2013. a

Howe, J. N. W., Piotrowski, A. M., Noble, T. L., Mulitza, S., Chiessi, C. M., and Bayon, G.: North Atlantic Deep Water Production during the Last Glacial Maximum, Nat. Commun., 7, 11765,, 2016. a

Ivanovic, R. F., Gregoire, L. J., Kageyama, M., Roche, D. M., Valdes, P. J., Burke, A., Drummond, R., Peltier, W. R., and Tarasov, L.: Transient climate simulations of the deglaciation 21–9 thousand years before present (version 1) – PMIP4 Core experiment design and boundary conditions, Geosci. Model Dev., 9, 2563–2587,, 2016. a, b

Ivanovic, R. F., Gregoire, L. J., Burke, A., Wickert, A. D., Valdes, P. J., Ng, H. C., Robinson, L. F., McManus, J. F., Mitrovica, J. X., Lee, L., and Dentith, J. E.: Acceleration of Northern Ice Sheet Melt Induces AMOC Slowdown and Northern Cooling in Simulations of the Early Last Deglaciation, Paleoceanogr. Paleoclimatol., 33, 807–824,, 2018. a

Jeltsch-Thömmes, A., Battaglia, G., Cartapanis, O., Jaccard, S. L., and Joos, F.: Low terrestrial carbon storage at the Last Glacial Maximum: constraints from multi-proxy data, Clim. Past, 15, 849–879,, 2019. a

Jochum, M., Chase, Z., Nuterman, R., Pedro, J., Rasmussen, S., Vettoretti, G., and Zheng, P.: Carbon Fluxes during Dansgaard–Oeschger Events as Simulated by an Earth System Model, J. Climate, 35, 5745–5758,, 2022. a

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J. M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and Millennial Antarctic Climate Variability over the Past 800,000 Years, Science, 317, 793–796,, 2007. a

Kapsch, M.-L., Mikolajewicz, U., Ziemen, F., and Schannwell, C.: Ocean Response in Transient Simulations of the Last Deglaciation Dominated by Underlying Ice-Sheet Reconstruction and Method of Meltwater Distribution, Geophys. Res. Lett., 49, e2021GL096767,, 2022. a, b

Kawasaki, T., Matsumura, Y., and Hasumi, H.: Deep water pathways in the North Pacific Ocean revealed by Lagrangian particle tracking, Sci. Rep., 12, 6238,, 2022. a

Key, R. M., Kozyr, A., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J. L., Feely, R. A., Millero, F. J., Mordy, C., and Peng, T.-H.: A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP), Global Biogeochem. Cy., 18, GB4031,, 2004. a

Kobayashi, H., Oka, A., Yamamoto, A., and Abe-Ouchi, A.: Glacial carbon cycle changes by Southern Ocean processes with sedimentary amplification, Sci. Adv., 7, eabg7723,, 2021. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Kohfeld, K. E. and Chase, Z.: Controls on deglacial changes in biogenic fluxes in the North Pacific Ocean, Quaternary Sci. Rev., 30, 3350–3363,, 2011. a

Kurahashi-Nakamura, T., Paul, A., and Losch, M.: Dynamical reconstruction of the global ocean state during the Last Glacial Maximum, Paleoceanography, 32, 326–350,, 2017. a

Lhardy, F., Bouttes, N., Roche, D. M., Abe-Ouchi, A., Chase, Z., Crichton, K. A., Ilyina, T., Ivanovic, R., Jochum, M., Kageyama, M., Kobayashi, H., Liu, B., Menviel, L., Muglia, J., Nuterman, R., Oka, A., Vettoretti, G., and Yamamoto, A.: A First Intercomparison of the Simulated LGM Carbon Results Within PMIP-Carbon: Role of the Ocean Boundary Conditions, Paleoceanogr. Paleoclimatol., 36, e2021PA004302,, 2021. a, b

Li, C., Clementi, V. J., Bova, S. C., Rosenthal, Y., Childress, L. B., Wright, J. D., Jian, Z., and Expedition 379T Scientists: The Sediment Green-Blue Color Ratio as a Proxy for Biogenic Silica Productivity Along the Chilean Margin, Geochem. Geophy. Geosy., 23, e2022GC010350,, 2022. a

Lippold, J., Luo, Y., Francois, R., Allen, S. E., Gherardi, J., Pichat, S., Hickey, B., and Schulz, H.: Strength and geometry of the glacial Atlantic Meridional Overturning Circulation, Nat. Geosci., 5, 813–816,, 2012. a

Liu, Z., Otto-Bliesner, B. L., He, F., Brady, E. C., Tomas, R., Clark, P. U., Carlson, A. E., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J.: Transient Simulation of Last Deglaciation with a New Mechanism for Bølling-Allerød Warming, Science, 325, 310–314,, 2009. a

Locarnini, R. A., O'Brien, T. D., Garcia, H. E., Antonov, J. I., Boyer, T. P., Conkright, M. E., and Stephens, C.: World Ocean Atlas 2001, vol. 3, Oxygen, NOAA Atlas NESDIS, vol. 51, NOAA, Silver Spring, MD, 286 pp., 2002. a

Lunt, D. J., Williamson, M. S., Valdes, P. J., Lenton, T. M., and Marsh, R.: Comparing transient, accelerated, and equilibrium simulations of the last 30 000 years with the GENIE-1 model, Clim. Past, 2, 221–235,, 2006. a

Lüthi, D., Floch, M. L., Bereiter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U., Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., and Stocker, T. F.: High-resolution carbon dioxide concentration record 650,000–800,000 years before present, Nature, 453, 379–382,, 2008. a

Maier, E., Méheust, M., Abelmann, A., Gersonde, R., Chapligin, B., Ren, J., Stein, R., Meyer, H., and Tiedemann, R.: Deglacial subarctic Pacific surface water hydrography and nutrient dynamics and links to North Atlantic climate variability and atmospheric CO2, Paleoceanography, 30, 949–968,, 2015. a

Marcott, S. A., Bauska, T. K., Buizert, C., Steig, E. J., Rosen, J. L., Cuffey, K. M., Fudge, T. J., Severinghaus, J. P., Ahn, J., Kalk, M. L., McConnell, J. R., Sowers, T., Taylor, K. C., White, J. W. C., and Brook, E. J.: Centennial-scale changes in the global carbon cycle during the last deglaciation, Nature, 514, 616–619,, 2014. a

MARGO Project Members: Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum, Nat. Geosci., 2, 127–132,, 2009. a

Mariotti, V., Paillard, D., Bopp, L., Roche, D. M., and Bouttes, N.: A coupled model for carbon and radiocarbon evolution during the last deglaciation, Geophys. Res. Lett., 43, 1306–1313,, 2016. a, b

Martin, J. H.: Glacial-interglacial CO2 change: The Iron Hypothesis, Paleoceanography, 5, 1–13,, 1990. a

Martínez-García, A., Sigman, D. M., Ren, H., Anderson, R. F., Straub, M., Hodell, D. A., Jaccard, S. L., Eglinton, T. I., and Haug, G. H.: Iron Fertilization of the Subantarctic Ocean During the Last Ice Age, Science, 343, 1347–1350,, 2014. a, b

McManus, J. F., Francois, R., Gherardi, J.-M., Keigwin, L. D., and Brown-Leger, S.: Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes, Nature, 428, 834–837,, 2004. a, b, c, d

Menviel, L., Timmermann, A., Mouchet, A., and Timm, O.: Climate and marine carbon cycle response to changes in the strength of the Southern Hemispheric westerlies, Paleoceanography, 23, PA4201,, 2008. a

Menviel, L., Timmermann, A., Timm, O. E., and Mouchet, A.: Deconstructing the Last Glacial termination: the role of millennial and orbital-scale forcings, Quaternary Sci. Rev., 30, 1155–1172,, 2011. a

Menviel, L., Joos, F., and Ritz, S.: Simulating atmospheric CO2, 13C and the marine carbon cycle during the Last Glacial–nterglacial cycle: possible role for a deepening of the mean remineralization depth and an increase in the oceanic nutrient inventory, Quaternary Sci. Rev., 56, 46–68,, 2012. a, b

Menviel, L., England, M. H., Meissner, K. J., Mouchet, A., and Yu, J.: Atlantic-Pacific seesaw and its role in outgassing CO2 during Heinrich events, Paleoceanography, 29, 58–70,, 2014. a

Menviel, L., Yu, J., Joos, F., Mouchet, A., Meissner, K. J., and England, M. H.: Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study, Paleoceanography, 32, 2–17,, 2017. a

Menviel, L., Spence, P., Yu, J., Chamberlain, M. A., Matear, R. J., Meissner, K. J., and England, M. H.: Southern Hemisphere westerlies as a driver of the early deglacial atmospheric CO2 rise, Nat. Commun., 9, 2503,, 2018. a, b, c

Millero, F. J.: Thermodynamics of the carbon dioxide system in the oceans, Geochim. Cosmochim. Ac., 59, 661–677,, 1995. a

Muglia, J., Skinner, L. C., and Schmittner, A.: Weak overturning circulation and high Southern Ocean nutrient utilization maximized glacial ocean carbon, Earth Planet. Sc. Lett., 496, 47–56,, 2018. a

Muglia, J., Mulitza, S., Repschläger, J., Schmittner, A., Lembke-Jene, L., Lisiecki, L., Mix, A., Saraswat, R., Sikes, E., Waelbroeck, C., Gottschalk, J., Lippold, J., Lund, D., Martinez-Mendez, G., Michel, E., Muschitiello, F., Naik, S., Okazaki, Y., Stott, L., Voelker, A., and Zhao, N.: A global synthesis of high-resolution stable isotope data from benthic foraminifera of the last deglaciation, Sci. Data, 10, 131,, 2023. a, b, c, d, e

Ng, H. C., Robinson, L. F., McManus, J. F., Mohamed, K. J., Jacobel, A. W., Ivanovic, R. F., Gregoire, L. J., and Chen, T.: Coherent deglacial changes in western Atlantic Ocean circulation, Nat. Commun., 9, 2947,, 2018. a, b, c, d

Obase, T. and Abe-Ouchi, A.: Abrupt Bölling-Alleröd Warming Simulated under Gradual Forcing of the Last Deglaciation, Geophys. Res. Lett., 46, 11397–11405,, 2019. a, b, c, d, e, f

Obase, T., Abe-Ouchi, A., and Saito, F.: Abrupt climate changes in the last two deglaciations simulated with different Northern ice sheet discharge and insolation, Sci. Rep., 11, 22359,, 2021. a, b

Obase, T., Menviel, L., Abe-Ouchi, A., Vadsaria, T., Ivanovic, R., Snoll, B., Sherriff-Tadano, S., Valdes, P., Gregoire, L., Kapsch, M.-L., Mikolajewicz, U., Bouttes, N., Roche, D., Lhardy, F., He, C., Otto-Bliesner, B., Liu, Z., and Chan, W.-L.: Multi-model assessment of the deglacial climatic evolution at high southern latitudes, Clim. Past Discuss. [preprint],, in review, 2023. a

Ohgaito, R., Yamamoto, A., Hajima, T., O'ishi, R., Abe, M., Tatebe, H., Abe-Ouchi, A., and Kawamiya, M.: PMIP4 experiments using MIROC-ES2L Earth system model, Geosci. Model Dev., 14, 1195–1217,, 2021a. 

Ohgaito, R., Yamamoto, A., Hajima, T., O'ishi, R., Abe, M., Tatebe, H., Abe-Ouchi, A., and Kawamiya, M.: Core code of MIROC-ES2L. In Geoscientific Model Development, Zenodo [code],, 2021b. a

Oka, A. and Niwa, Y.: Pacific deep circulation and ventilation controlled by tidal mixing away from the sea bottom, Nat. Commun., 4, 2419,, 2013. a

Okazaki, Y., Timmermann, A., Menviel, L., Harada, N., Abe-Ouchi, A., Chikamoto, M. O., Mouchet, A., and Asahi, H.: Deepwater Formation in the North Pacific During the Last Glacial Termination, Science, 329, 200–204,, 2010. a, b

O'Leary, M. H.: Carbon isotope fractionation in plants, Phytochemistry, 20, 553–567,, 1981. a

Parekh, P., Follows, M. J., and Boyle, E. A.: Decoupling of iron and phosphate in the global ocean, Global Biogeochem. Cy., 19, GB2020,, 2005. a, b

Paul, A., Mulitza, S., Stein, R., and Werner, M.: A global climatology of the ocean surface during the Last Glacial Maximum mapped on a regular grid (GLOMAP), Clim. Past, 17, 805–824,, 2021. a

Petit, J. R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., PÉpin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica, Nature, 399, 429–436,, 1999. a

Pöppelmeier, F., Jeltsch-Thömmes, A., Lippold, J., Joos, F., and Stocker, T. F.: Multi-proxy constraints on Atlantic circulation dynamics since the last ice age, Nat. Geosci., 16, 349–356,, 2023. a, b, c, d

Rae, J. W. B., Sarnthein, M., Foster, G. L., Ridgwell, A., Grootes, P. M., and Elliott, T.: Deep water formation in the North Pacific and deglacial CO2 rise, Paleoceanography, 29, 645–667,, 2014. a, b

Rafter, P. A., Gray, W. R., Hines, S. K., Burke, A., Costa, K. M., Gottschalk, J., Hain, M. P., Rae, J. W., Southon, J. R., Walczak, M. H., Yu, J., Adkins, J. F., and DeVries, T.: Global reorganization of deep-sea circulation and carbon storage after the last ice age, Sci. Adv., 8, eabq5434,, 2022. a, b, c, d, e, f, g, h, i, j

Reimer, P. J., Austin, W. E. N., Bard, E., Bayliss, A., Blackwell, P. G., Ramsey, C. B., Butzin, M., Cheng, H., Edwards, R. L., Friedrich, M., Grootes, P. M., Guilderson, T. P., Hajdas, I., Heaton, T. J., Hogg, A. G., Hughen, K. A., Kromer, B., Manning, S. W., Muscheler, R., Palmer, J. G., Pearson, C., van der Plicht, J., Reimer, R. W., Richards, D. A., Scott, E. M., Southon, J. R., Turney, C. S. M., Wacker, L., Adolphi, F., Büntgen, U., Capano, M., Fahrni, S. M., Fogtmann-Schulz, A., Friedrich, R., Köhler, P., Kudsk, S., Miyake, F., Olsen, J., Reinig, F., Sakamoto, M., Sookdeo, A., and Talamo, S.: The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve (0–55 cal kBP), Radiocarbon, 62, 725–757,, 2020. a, b

Schmitt, J., Schneider, R., Elsig, J., Leuenberger, D., Lourantou, A., Chappellaz, J., Köhler, P., Joos, F., Stocker, T. F., Leuenberger, M., and Fischer, H.: Carbon isotope constraints on the deglacial CO2 rise from ice cores, Science, 336, 711–714,, 2012. a, b, c, d

Schmittner, A. and Galbraith, E. D.: Glacial greenhouse-gas fluctuations controlled by ocean circulation changes, Nature, 456, 373–376,, 2008. a, b, c

Schmittner, A. and Lund, D. C.: Early deglacial Atlantic overturning decline and its role in atmospheric CO2 rise inferred from carbon isotopes (δ13C), Clim. Past, 11, 135–152,, 2015. a, b

Schmittner, A., Gruber, N., Mix, A. C., Key, R. M., Tagliabue, A., and Westberry, T. K.: Biology and air–sea gas exchange controls on the distribution of carbon isotope ratios (δ13C) in the ocean, Biogeosciences, 10, 5793–5816,, 2013. a

Sherriff-Tadano, S., Abe-Ouchi, A., Yoshimori, M., Ohgaito, R., Vadsaria, T., Chan, W.-L., Hotta, H., Kikuchi, M., Kodama, T., Oka, A., and Suzuki, K.: Southern Ocean Surface Temperatures and Cloud Biases in Climate Models Connected to the Representation of Glacial Deep Ocean Circulation, J. Climate, 36, 3849–3866,, 2023. a

Siegenthaler, U., Stocker, T. F., Monnin, E., Lüthi, D., Schwander, J., Stauffer, B., Raynaud, D., Barnola, J.-M., Fischer, H., Masson-Delmotte, V., and Jouzel, J.: Stable Carbon Cycle–Climate Relationship During the Late Pleistocene, Science, 310, 1313–1317,, 2005. a

Sigman, D. M., Fripiat, F., Studer, A. S., Kemeny, P. C., Martínez-García, A., Hain, M. P., Ai, X., Wang, X., Ren, H., and Haug, G. H.: The Southern Ocean during the ice ages: A review of the Antarctic surface isolation hypothesis, with comparison to the North Pacific, Quaternary Sci. Rev., 254, 106 732,, 2021. a

Sikes, E. L., Umling, N. E., Allen, K. A., Ninnemann, U. S., Robinson, R. S., Russell, J. L., and Williams, T. J.: Southern Ocean glacial conditions and their influence on deglacial events, Nat. Rev. Earth Environ., 4, 454–470,, 2023. a

Skinner, L., Primeau, F., Jeltsch-Thömmes, A., Joos, F., Köhler, P., and Bard, E.: Rejuvenating the ocean: mean ocean radiocarbon, CO2 release, and radiocarbon budget closure across the last deglaciation, Clim. Past, 19, 2177–2202,, 2023. a

Snoll, B., Ivanovic, R., Gregoire, L., Sherriff-Tadano, S., Menviel, L., Obase, T., Abe-Ouchi, A., Bouttes, N., He, C., He, F., Kapsch, M., Mikolajewicz, U., Muglia, J., and Valdes, P.: A multi-model assessment of the early last deglaciation (PMIP4 LDv1): The meltwater paradox reigns supreme, EGUsphere [preprint],, 2023. a

Stein, K., Timmermann, A., Kwon, E. Y., and Friedrich, T.: Timing and magnitude of Southern Ocean sea ice/carbon cycle feedbacks, P. Natl. Acad. Sci. USA, 117, 4498–4504,, 2020. a, b

Studer, A. S., Sigman, D. M., Martínez-García, A., Benz, V., Winckler, G., Kuhn, G., Esper, O., Lamy, F., Jaccard, S. L., Wacker, L., Oleynik, S., Gersonde, R., and Haug, G. H.: Antarctic Zone nutrient conditions during the last two glacial cycles, Paleoceanography, 30, 845–862,, 2015. a

Stuiver, M., Quay, P. D., and Ostlund, H. G.: Abyssal Water Carbon-14 Distribution and the Age of the World Oceans, Science, 219, 849–851,, 1983. a

Takemura, T.: Simulation of climate response to aerosol direct and indirect effects with aerosol transport-radiation model, J. Geophys. Res., 110, D02202,, 2005. a

Thiagarajan, N. and McManus, J. F.: Productivity and sediment focusing in the Eastern Equatorial Pacific during the last 30,000 years, Deep-Sea Res. Pt. I, 147, 100–110,, 2019. a

Tierney, J. E., Zhu, J., King, J., Malevich, S. B., Hakim, G. J., and Poulsen, C. J.: Glacial cooling and climate sensitivity revisited, Nature, 584, 569–573,, 2020. a

Timm, O. and Timmermann, A.: Simulation of the Last 21 000 Years Using Accelerated Transient Boundary Conditions, J. Climate, 20, 4377–4401,, 2007. a

Tschumi, T., Joos, F., Gehlen, M., and Heinze, C.: Deep ocean ventilation, carbon isotopes, marine sedimentation and the deglacial CO2 rise, Clim. Past, 7, 771–800,, 2011. a, b

Uemura, R., Motoyama, H., Masson-Delmotte, V., Jouzel, J., Kawamura, K., Goto-Azuma, K., Fujita, S., Kuramoto, T., Hirabayashi, M., Miyake, T., Ohno, H., Fujita, K., Abe-Ouchi, A., Iizuka, Y., Horikawa, S., Igarashi, M., Suzuki, K., Suzuki, T., and Fujii, Y.: Asynchrony between Antarctic temperature and CO2 associated with obliquity over the past 720,000 years, Nat. Commun., 9, 961,, 2018. a

Weber, M. E.: Antiphased dust deposition and productivity in the Antarctic Zone over 1.5 million years, PANGAEA [data set],, 2021. a

Wilmes, S.-B., Green, J. A. M., and Schmittner, A.: Enhanced vertical mixing in the glacial ocean inferred from sedimentary carbon isotopes, Commun. Earth Environ., 2, 166,, 2021. a, b

Zhao, N., Marchal, O., Keigwin, L., Amrhein, D., and Gebbie, G.: A Synthesis of Deglacial Deep-Sea Radiocarbon Records and Their (In)Consistency With Modern Ocean Ventilation, Paleoceanogr. Paleoclimatol., 33, 128–151,, 2018. a

Short summary
This study examines the transient response of the ocean carbon cycle to climate change since the last ice age by using an ocean general circulation model. Our carbon cycle model calculates atmospheric pCO2 changes that are consistent with ice core records but whose magnitude is underestimated. Our analysis of carbon isotopes suggests that improving the expression of activated ocean ventilation and suppressing biological productivity are critical in simulating atmospheric pCO2 changes.