Articles | Volume 18, issue 9
Research article
01 Sep 2022
Research article |  | 01 Sep 2022

Glacial state of the global carbon cycle: time-slice simulations for the last glacial maximum with an Earth-system model

Takasumi Kurahashi-Nakamura, André Paul, Ute Merkel, and Michael Schulz

Three time-slice carbon cycle simulations for the last glacial maximum (LGM) constrained by the CO2 concentration in the atmosphere and the increase in the mean concentration of dissolved inorganic carbon in the deep ocean were carried out with a fully coupled comprehensive climate model (the Community Earth System Model version 1.2). The three modelled LGM ocean states yielded different physical features in response to artificial freshwater forcing, and, depending on the physical states, suitable amounts of carbon and alkalinity were added to the ocean to satisfy constraints from paleo-data. In all the simulations, the amount of carbon added was in line with the inferred transfers of carbon among various reservoirs during the evolution from the LGM to the pre-industrial (PI) period, suggesting that the simulated glacial ocean states are compatible with the PI one in terms of the carbon budget. The increase in total alkalinity required to simulate ocean states that were deemed appropriate for the LGM was in broad quantitative accord with the scenario of post-glacial shallow water deposition of calcium carbonate, although a more precise assessment would demand further studies of various processes such as the land chemical weathering and deep-sea burial of calcium carbonates, which have affected the alkalinity budget throughout history since the LGM. On the other hand, comparisons between the simulated distributions of paleoceanographic tracers and corresponding reconstructions clearly highlighted the different water-mass geometries and favoured a shallower Atlantic meridional overturning circulation (AMOC) for the LGM as compared to PI.

1 Introduction

Global climate change during the last 800 kyr was characterized by the periodic variations in ice sheet volume and atmospheric CO2 concentration (pCO2) with a ∼100 kyr period (e.g. Petit et al.1999; EPICA community members2004), i.e. the glacial–interglacial cycles. The synchronous behaviours of those two quantities imply a strong link between them, and it is suggested that the variations in pCO2 at the glacial–interglacial timescale were regulated by the evolution of ice sheets that represented a nonlinear response of the Earth system to the orbital forcing (e.g. Ganopolski and Brovkin2017). This mechanism for the pCO2 variations is also supported by the fact that the periodic changes in ice sheet volume can be numerically simulated even under temporally constant pCO2 if it is sufficiently low (e.g. Ganopolski and Calov2011; Abe-Ouchi et al.2013). Therefore, to comprehend the mechanism of pCO2 variations in the glacial–interglacial cycles, it is fundamental to understand the response of the global carbon cycle to a given change in the states of other climatic components that work as external factors from a carbon cycle viewpoint.

Ideally, to grasp the pCO2 variations over an entire glacial cycle, a time-dependent or transient framework with a fully coupled comprehensive climate model would be beneficial to explicitly deal with a gradually changing background state regulated by slow processes that work at a 104-year timescale such as the changes of the orbital configuration and the growth of ice sheets. From the perspective of the marine carbon cycle, which is considered to have played a major role in the global carbon cycle (e.g. Sigman and Boyle2000), the changes in the inventory of total alkalinity as well as carbon in the entire ocean would become relevant for the dynamics at such a 104-year timescale. The changes in the alkalinity inventory are considered to have occurred as the result of a slowly accumulated imbalance in alkalinity budgets through two different processes as follows.

First, the global sea-level changes accompanied by the growth or retreat of ice sheets affect the area of continental shelves covered by seawater, which is of relevance for the shallow-water ecosystem. For example, the sea level has risen since the LGM and seawater has flooded the continental shelves, so that the modern-day area of shallow water environments for carbonate deposition including coral reef buildup is 3 times as large as that for the LGM (Cartapanis et al.2018). As a result, a substantial amount of alkalinity is considered to have been removed from the LGM ocean by the post-glacial shallow-water deposition of CaCO3, and it is estimated that this process would account for the removal of total alkalinity by as much as 0.44–3.6×1017 equivalent (Milliman1993; Opdyke2000; Ridgwell et al.2003; Vecsei and Berger2004; Husson et al.2018; Köhler and Munhoven2020). Because the removal of alkalinity will cause the increase in the partial pressure of CO2 in the surface ocean, a hypothesis that explains a substantial part of the glacial–interglacial pCO2 variation by the shallow-water deposition of CaCO3 was proposed (the so-called “coral reef hypothesis”; e.g. Berger1982; Opdyke and Walker1992). Although this hypothesis inevitably requires that the pCO2 rise followed the expansion of the area of shallow-water deposition, the actual deglaciation had the opposite sequence; that is to say, a substantial part (∼70 ppm) of the recorded pCO2 increase happened during the first 10 000 years following the LGM (Bereiter et al.2015; Köhler et al.2017), while the expansion of the shallow-water environments is estimated to have mostly occurred afterward (Cartapanis et al.2018). Alternatively, it has been demonstrated with a simplified carbon cycle model that the post-glacial buildup of shallow-water carbonates is consistent with the observed ∼20 ppm rise of pCO2 during the late Holocene (Ridgwell et al.2003).

Second, the total inventory of alkalinity is subject to change also through the deposition or dissolution of CaCO3 in the deep sediments of the open ocean, because the amount of alkalinity in the ocean is affected by a long-term budget of CaCO3 through its input by chemical weathering on land and the removal by sedimentary burial at the ocean floor. In particular, if there is an imbalance between the input and removal, the alkalinity inventory will respond so that the flux balance becomes restored at a timescale of 104 years (e.g. Boudreau et al.2018). For example, a positive imbalance (i.e. more input than removal) would be removed by an increased CaCO3 deposition, because seawater with a higher alkalinity concentration is more favourable to the deposition of CaCO3. This intrinsic negative feedback process by the ocean is known as “carbonate compensation”, which has been considered to be a major mechanism affecting a long-term climate evolution (e.g. Broecker and Peng1987; Archer et al.2000; Brovkin et al.2007; Chikamoto et al.2008; Lord et al.2016; Boudreau et al.2018; Kobayashi et al.2021).

Although transient simulations including those “slow” processes are advantageous to cope with the evolution of the global carbon cycle spanning the glacial–interglacial timescale, the computational cost for such long-term simulations with a comprehensive climate model is still prohibitive. Instead, a more practical way is to conduct time-slice simulations, for example, for the LGM, that are subject to reconstructed boundary conditions such as the greenhouse gas concentrations and orbital configuration. Such time-slice simulations may also serve as an initial state for shorter transient runs that cover a part of the full glacial cycle (e.g. the last deglaciation) instead of a “true” LGM state that results from the evolution of the global carbon cycle since (at least) the last interglacial period.

The main purpose of this study is to carry out fully coupled LGM time-slice simulations including the global carbon cycle. In our model runs, each of which had 2500 model years, we fixed the ice sheet configuration, the orbital configuration, and the inventory of total alkalinity in the ocean as external factors or a background field assuming that the effects of their variations were minor at a millennial timescale. Instead, processes that govern the marine carbon cycle at shorter timescales such as air–sea gas exchange, the biological pump, and the change of tracer distributions by ocean circulation (e.g. Hain et al.2014) were explicitly modelled and simulated. Therefore, the simulations in this study inevitably do not capture the influence of long-term trends in the boundary conditions. For example, the properties of very old water do not explicitly reflect the past state of boundary conditions but only approximate the true properties.

To constrain the carbon cycle state during the LGM, we used two important restrictions in terms of the size of carbon reservoirs: atmospheric pCO2 and an estimated rise of the mean concentration of dissolved inorganic carbon (DIC) in the deep ocean. The LGM pCO2 level is known to be approximately 190 ppm based on ice-core records (Bereiter et al.2015; Köhler et al.2017). The increase in DIC concentrations in the deep ocean has been estimated based on an empirical linear relationship between the 14C ventilation age and modern concentrations of DIC below a depth of 2000 m, suggesting that the average DIC concentrations in the abyss were higher during the LGM by 85–115 µmol kg−1 (corresponding to 730–980 GtC; Sarnthein et al.2013) or 82 µmol kg−1 (687 GtC; Skinner et al.2015) as compared to the modern day. We adjusted the mean concentrations of DIC and total alkalinity in the entire ocean to satisfy the two constraints in addition to adopting the PMIP4 protocol regarding LGM boundary conditions (Kageyama et al.2017).

The increase in deep-water DIC storage during the LGM should be the outcome of two combined factors: the change in the vertical distribution of stored carbon and the change in the mean concentration. The increased vertical contrast contributes to the decrease in pCO2 as the increase in alkalinity does, while the increase in the mean DIC concentration has the opposite effect. The adjustment of DIC concentrations for the entire ocean in this study will indicate whether the net effect of these two factors is in line with the estimated increase in the deep DIC concentrations.

We thereby evaluate the LGM marine carbon cycle from the perspective of the specific change in the pH of seawater that is compatible with the low-pCO2 atmosphere. The carbonate ion concentration [CO32-] in the deep ocean is tightly linked to processes that affect the redistribution of carbon in the ocean and thus is a valuable tracer for investigating the global carbon cycle (e.g. Yu et al.2014). The simulated [CO32-] and the mass accumulation rate (MAR) of CaCO3 at the ocean floor are both compared to reconstructions. We also discuss the compatibility between the modelled LGM and modern states in terms of the carbon and alkalinity budget. That is to say, we examine the likelihood that we can reconcile the total inventory of DIC and alkalinity of the modelled LGM states with the modern counterparts, taking into account the changes in the sizes of several reservoirs during the transition from the LGM to PI.

2 Methods

2.1 Models

We used the Community Earth System Model (CESM; Hurrell et al.2013) version 1.2, which is a coupled atmosphere–ocean–land comprehensive Earth-system model. The ocean component, the Parallel Ocean Program version 2 (POP2), was configured to include the Biogeochemical Elemental Cycling model (BEC; Moore et al.2004, 2013; Lindsay et al.2014). The BEC included a nutrient–phytoplankton–zooplankton–detritus (NPZD) marine ecosystem model and handled sinking processes of biological particles to compute the distribution of various biogeochemical tracers such as phosphate, apparent oxygen utilization (AOU), and carbonate ion concentration. For the sinking of particulate organic matter (POM), the “ballast effect” was taken into account based on Armstrong et al. (2002). The ocean component was extended by the carbon-isotope package developed by Jahn et al. (2015) so that we were able to explicitly simulate the carbon-isotope composition of seawater. The land component (the Community Land Model; CLM4.0) included a prognostic treatment of the terrestrial carbon and nitrogen cycles to calculate the reservoir size of land carbon. A low-resolution configuration of CESM was used in this study to reduce the computational cost (Shields et al.2012); that is to say, the ocean component had a nominal 3 irregular horizontal grid with 60 vertical levels, while the atmosphere component, Community Atmosphere Model version 4 (CAM4), had a T31 spectral dynamical core (horizontal resolution of 3.75) with 26 vertical levels.

In this study, the size of each carbon reservoir including the atmospheric pCO2 is calculated by the carbon cycle model. However, we switched off the feedback from the carbon cycle to the climate system by prescribing a constant pCO2 that determined the radiative forcing. We treated the marine biogeochemical system as being “closed” to exclude the drift of the total inventory from prescribed amounts; that is, there were no fluxes of biogeochemical matter through the top and bottom boundaries of the ocean domain except for the air–sea gas exchange. Particulate matter that reached the ocean floor was assumed to be remineralized in the bottom-most cells. To enhance the numerical stability, we replaced the original pH solver in POP2 with SolveSAPHE (Solver Suite for Alkalinity–pH Equations) developed by Munhoven (2013). We introduced this pH solver originally to avoid numerical instability especially in test runs or in a spin-up phase and continued to use it for consistency. We took into account the alkalinity that arises from the dissociation of the solvent water itself in addition to the carbonate and borate alkalinity, which provided a sufficiently accurate solution (Munhoven2013).

For diagnostic purposes, we adopted the Model of Early Diagenesis in the Upper Sediment of Adjustable complexity (MEDUSA(v.2); Munhoven2021) to explicitly calculate the preservation or dissolution of CaCO3 at the ocean floor that corresponded to each of our CESM simulations. MEDUSA is a one-dimensional advection–diffusion–reaction model that describes the early diagenetic processes involving carbonates and organic matter (OM) in the surface sediment of 10 cm thickness as a function of time. The configuration of MEDUSA was the same as in our previous application (Kurahashi-Nakamura et al.2020) except for diffusive boundary layers (DBLs) to better represent the solute fluxes across the water–sediment boundaries (Munhoven2021) in this study.

2.2 Experiments and analyses

We performed four main time-slice experiments: one control experiment that corresponds to the PI reference period and three experiments with the LGM as a target (Table 1). For the PI run (expPI) we followed the DECK (Diagnostic, Evaluation and Characterization of Klima) PI control protocol (Eyring et al.2016) for the physical forcing to CESM. The control experiment was initialized from the final state of the spin-up run of Kurahashi-Nakamura et al. (2020) and run for 2500 model years. For the model parameters, we followed the default settings of CESM1.2 (, last access: 1 August 2022).

For the LGM runs, we followed the experimental design established by the Paleoclimate Modelling Intercomparison Project Phase 4 (PMIP4; Kageyama et al.2017). We set the radiative and orbital forcing parameters as specified in this protocol and adjusted the ice sheet configuration based on the ICE-6G-C reconstruction (Argus et al.2014; Peltier et al.2015). The dust and other aerosol forcings of the atmosphere and marine biogeochemistry are based on Albani et al. (2014). The mean salinity and nutrient concentrations of the ocean were modified according to the change in ocean volume (Kageyama et al.2017). Namely, we increased the salinity by 1 psu and all other tracer concentrations by 3 %.

Table 1Overview of the experimental designs in this study.

a See also Fig. 1. bEyring et al. (2016). cKageyama et al. (2017).

Download Print Version | Download XLSX

Figure 1Regions allotted for the additional freshwater forcing (shown in red). Additional 0.1 Sv in total was uniformly given to the specified region in the North Atlantic for expLGMws, while, for expLGMss, 0.1 and 0.25 Sv were subtracted from the North Atlantic and Southern Ocean, respectively. Corresponding compensation (i.e. freshwater flux having the opposite sign) was applied to the other regions homogeneously to keep the total volume of seawater constant.

For the baseline LGM experiment (expLGM), we carried out a 2500-year fully coupled carbon cycle run which was branched from a separate physical spin-up without the marine biogeochemistry. We also conducted two sensitivity experiments with the same length for the LGM (expLGMws and expLGMss) with additional freshwater forcing to examine the dependency of the biogeochemical state on the physical ocean state. For the first sensitivity experiment (expLGMws), 0.1 Sv additional freshwater in total was uniformly added to the high-latitude North Atlantic Ocean throughout the simulation. The southern boundary of the freshwater-forcing region was defined at 50 N, and the others by pre-defined CESM1.2 region masks (Fig. 1). For the second sensitivity experiment (expLGMss), we extracted 0.1 Sv from the same high-latitude region in the North Atlantic Ocean and 0.25 Sv in total from the Weddell Sea and Ross Sea regions. The primary motivation for the additional freshwater forcing was to prepare different physical ocean states as in previous studies (e.g. Gu et al.2020; Muglia and Schmittner2021), and the freshwater amount was determined empirically to realize an AMOC with significantly different characteristics compared to the baseline state. Implicitly the additional freshwater forcing takes into account the uncertainty in the glacial freshwater budget in high latitudes due to ice sheet calving and iceberg transport and melting processes (Merino et al.2016) that were not explicitly modelled in our framework. The total seawater volume and hence the global-mean concentrations of tracers were conserved by adding uniform compensating fluxes of the opposite sign over the rest of the global ocean.

The mean concentration of DIC in the LGM simulations was set based on an independent proxy-based estimate of the increase in the concentration in the deep ocean by 100 µmol kg−1 as compared to the modern value by Sarnthein et al. (2013) (Table 1). In fact, the increment of the DIC concentration for our three individual LGM simulations was determined through corresponding preparatory simulations (Fig. 2). To initialize the preparatory runs, we increased the global mean DIC concentration homogeneously by 100 mmol m−3 compared to the PI amount as a first guess (Table 2). Depending on the respective LGM simulations, we also added 75–100 meq m−3 of total alkalinity on top of the ocean-volume effect to tune the model so that a pCO2 level between 180 and 190 ppm was reached. After a 1000-year model integration for each simulation, we calculated the resulting mean DIC anomalies in the ocean deeper than 2000 m in accordance with Sarnthein et al. (2013). To set the initial mean concentration of DIC for the main LGM simulations, we took the difference between each of the posterior DIC anomalies of the preparatory runs and the observation-based estimate and adjusted the first guess (i.e. 100 mmol m−3) by adding the difference. Thereby we obtained second-guess values that were used to initialize the respective main LGM runs (Table 1) so that we were able to achieve the DIC anomalies in the deep ocean in better accordance with the estimate by Sarnthein et al. (2013). For the main LGM runs, we newly adjusted the total alkalinity independently of the preparatory runs. A total of 40 to 80 meq m−3 of alkalinity depending on the experiments was added in addition to the increase due to the change in the seawater volume (Table 1). The increment was chosen similarly such that the model predicted an atmospheric pCO2 level of approximately 190 ppm based on ice-core records (Bereiter et al.2015; Köhler et al.2017).

Figure 2Workflow for the LGM carbon cycle simulations in this study.


Table 2The first-guess adjustment of the global mean concentrations of DIC and total alkalinity for the preparatory runs and the resultant DIC anomaly averaged over the depths of more than 2000 m. In the last column, the differences between each of the LGM runs and expPI that are averaged over the last 50 years are shown.

Download Print Version | Download XLSX

The initial δ13C of DIC was uniformly zero, and we assumed the δ13C value of the atmosphere during the LGM to be −6.46 ‰ (Schmitt et al.2012). Although the treatment of 13C in the current model configuration is not self-consistent because the air–sea gas exchange should have affected the atmospheric δ13C in reality, we took the approach for three practical reasons.

  1. The prescribed atmospheric δ13C set to a reliable observed value was expected to result in a better model representation of δ13CDIC in the ocean that could be compared with observation-based data.

  2. The deviation of the atmospheric pCO2 from the required 190 ppm was very small, so that the simulated state would approximate the consistent LGM state in a reasonable way, including the atmospheric δ13C value.

  3. As far as we recognize, the available carbon isotope package (Jahn et al.2015) does not deal with atmospheric δ13C that evolves interactively and self-consistently with the air–sea gas exchange in the model, and therefore, a further substantial model development will be needed to implement it, which is beyond the scope of this study.

For each experiment, we also diagnosed and analysed upper-sediment properties with regard to CaCO3 by running MEDUSA separately as a stand-alone model, taking the appropriate boundary conditions from the corresponding CESM state. As in Kurahashi-Nakamura et al. (2020), one MEDUSA column was coupled to the deepest grid cell of each POP2 water column, and there was no lateral exchange of information among the MEDUSA columns.

It should be noted that the experiments in this study had a different total amount of carbon stored in the whole (atmosphere–ocean–land) system, reflecting the uncertain physical ocean states. In our experiments, instead of fixing the total amount, we tuned the size of the atmospheric reservoir to a specific amount (i.e. ∼190 ppm), which practically governed the terrestrial reservoir's size, and also tuned the size of the deep-ocean reservoir as well. The size of the shallow-ocean reservoir differed among the different experiments depending on the vertical gradient of DIC concentration in the ocean, hence the ocean circulation.

3 Results

3.1 Physical ocean states

The different forcing factors and boundary conditions used for the PI and LGM time slices yielded four very distinct physical ocean states. As expected, there were noticeable differences in the global ocean circulation and water-mass distributions. The standard LGM experiment (expLGM) had a stronger Atlantic meridional overturning circulation (AMOC) with a similar depth structure compared to that of the pre-industrial run (expPI) (Fig. 3). Both experiments with the additional freshwater forcing resulted in noticeably shallower AMOCs, which shoaled by 500–1000 m in terms of the position of the zero-flux level that separates the upper and lower circulation cells (Fig. 3). On the other hand, the maximum strength of the upper circulation in expLGMss was stronger than in expPI, while that for expLGMws was weaker. In short, expLGMss had a stronger but shallower AMOC compared to PI, while expLGMws had a weaker and shallower AMOC. Both experiments showed stronger bottom circulations than expPI or expLGM.

Figure 3Simulated AMOC for the main experiments in this study (Sv): stream function for (a) expPI, (b) expLGM, (c) expLGMws, and (d) expLGMss. The vertical profiles of the AMOC (e) at 30 N and (f) at 30 S are also shown. The average states of the last 50 years are shown.


The differences in the annually averaged global-mean sea surface temperature (SST) between the respective LGM experiment and expPI were −2.6 K (expLGM), −2.7 K (expLGMss), and −2.0 K (expLGMws). Those values were within the range of various estimates of the mean SST anomaly for the LGM (Annan and Hargreaves2013; Kurahashi-Nakamura et al.2017; Tierney et al.2020; Paul et al.2021), implying that the solubility of CO2 in the model is consistent with the LGM climate. The vertical gradient of modelled LGM salinity was larger than in expPI in general (i.e. comparatively more saline in the very deep ocean) except for the Atlantic Ocean in expLGMws (Fig. 4e–j), which corresponded to more stratified ocean states. However, the degree of stratification differed depending on the magnitude of the additional freshwater forcing. Accordingly, the horizontal distributions of the bottom salinity varied among the LGM experiments (Fig. 4a–g). The best fit to the reconstructed salinity of bottom water at several locations (Adkins et al.2002; Insua et al.2014; Homola et al.2021) was obtained in expLGMss that had the highest bottom salinity, suggesting that the more vigorous penetration of more saline southern-sourced water contributed to the result. In contrast, in expLGMws, the ocean was just slightly more stratified compared to expPI on the whole, and the model–data misfits were the largest.

Figure 4Simulated salinity for the main experiments in this study. The left and middle columns show the global distributions of salinity in the bottom-most grid cells. The absolute values for (a) expPI, (b) expLGM, (c) expLGMws, and (d) expLGMss and the differences between each LGM experiment and expPI are shown (e–g). The coloured dots indicate the reconstructions by Adkins et al. (2002), Insua et al. (2014), and Homola et al. (2021). The meridional sections in the Atlantic (30 W) are also shown for the anomalies (h–j). The average states of the last 50 years are shown.

3.2 Carbon reservoirs

The pCO2 predicted by the carbon cycle module of the model was in a quasi-steady state during at least the last 1000 model years in all the experiments, and the drift of pCO2 in the last 500 years was 0.6 ppm or less (depending on the simulations). In expPI that served as the control case, the predicted pCO2 was 276 ppm (Table 3), which agreed well with the prescribed radiative pCO2 that was used to force the model climate (280 ppm). All the LGM runs reached a pCO2 in the range between 180 and 190 ppm as aimed for, which was approximately 90 ppm lower than in expPI and thus in good agreement with the pCO2 difference obtained from ice cores (Bereiter et al.2015; Köhler et al.2017).

The posterior anomaly of the mean DIC concentration in the deep ocean (>2000 m) as the result of the 2500-year LGM integrations was in the range between 102 and 116 mmol m−3, showing that the second-guess DIC values that were used to initialize the respective LGM runs reasonably captured the target based on the estimate by Sarnthein et al. (2013).

The simulated sizes of the land carbon storage were similar among the three LGM runs and ranged between 1.44×103 and 1.49×103 GtC. They were 340 to 390 GtC smaller than in expPI (Table 3).

Table 3Biogeochemical properties of the model states (the last 50-year average) in this study. For DIC, the anomaly is globally averaged over the depths of more than 2000 m and shows the difference between each of the LGM simulations and expPI.

* Global export production at 100 m depth including particulate organic carbon
(POC) and CaCO3.

Download Print Version | Download XLSX

3.3 Carbon isotopes, phosphate, AOU, and the ideal age

The vertical structure of the distributions of biogeochemical tracers reflected the physical characteristics of each LGM experiment. The two experiments with a shallower AMOC (expLGMws and expLGMss) on the whole showed a similar structure in the meridional sections of the Atlantic Ocean, while those for expLGM revealed markedly different features. For example, δ13CDIC values were higher in the upper ocean by ∼1 ‰ compared to those in expPI and lower in the deeper ocean to a similar degree in expLGMws and expLGMss (Fig. 5b, c), while this contrast was significantly smaller in expLGM (Fig. 5a). As a result, expLGMws and expLGMss demonstrated an excellent model representation of δ13CDIC, especially in terms of the depth that separates the positive and negative anomalies when compared to the data by previous studies (Peterson et al.2014; Yu et al.2020). The clear shallow–deep contrast of δ13CDIC anomalies was also in line with the estimated distribution by Oppo et al. (2018). In the Pacific Ocean, the negative anomalies indicated by data were well reproduced by all the LGM experiments except in the Northern Hemisphere in expLGMws (Fig. 5d–f). It should be noted that δ13CDIC was less equilibrated because of the uniform initial condition and that the drift of δ13CDIC in the last 500 model years in the deep North Pacific (at 30 N, 150 W, 2900 m) was 0.06 ‰ or less depending on the experiments, which was less than the typical magnitude of data uncertainty (e.g. Kurahashi-Nakamura et al.2017).

Figure 5Meridional sections in the Atlantic (a–c; at 30 W) and in the Pacific (d–f; at 150 W) of the simulated δ13C of DIC in the LGM experiments (‰): expLGM (a, d), expLGMws (b, e), and expLGMss (c, f). The differences between each LGM experiment and expPI are shown. The average states of the last 50 years are shown. The dots indicate the reconstructions by Peterson et al. (2014) and Yu et al. (2020).


A notable difference in the vertical structure in the Atlantic Ocean was also found in the distributions of dissolved phosphate (Fig. 6). In expLGMss and expLGMws, the anomaly in phosphate concentration reached more than 1 mmol m−3 in the lower half of the depth range of the Atlantic Ocean, while it was negligible for the expLGM counterpart (Fig. 6a–c). The clear shallow–deep contrast of phosphate anomalies in expLGMss and expLGMws was in agreement with the inversion-based estimate by Oppo et al. (2018). These characteristics of phosphate distribution in the Atlantic Ocean were in accordance with the distribution of the AOU (Fig. 7a–c), which suggested that the strong positive anomalies of phosphate concentration in expLGMss and expLGMws were largely caused by an increase in remineralized phosphate. The distributions of both tracers showed patterns similar to each other also in the Pacific Ocean (Figs. 6d–f and 7d–f), similarly indicating the increase in phosphate concentrations due to amplified remineralization. Moreover, the inferred larger input of biological matter into the deep ocean was in accordance with the negative anomaly of δ13CDIC that covered most parts of the Pacific section. On the contrary, the parts that had smaller positive (or negative) anomalies of AOU corresponded to the positive anomalies of δ13CDIC. Combined with the younger ideal age of seawater in the same regions (Fig. 8d–f), these results suggested that smaller storage of remineralized matter brought the positive anomalies of δ13CDIC, which would have caused the clear mismatch between the simulated δ13CDIC in expLGMws and the data in the North Pacific.

Figure 6Meridional sections in the Atlantic (a–c; at 30 W) and in the Pacific (d–f; at 150 W) of the simulated phosphate concentrations in the LGM experiments (mmol m−3): expLGM (a, d), expLGMws (b, e), and expLGMss (c, f). The differences between each LGM experiment and expPI are shown. The average states of the last 50 years are shown.


Figure 7Meridional sections in the Atlantic (a–c; at 30 W) and in the Pacific (d–f; at 150 W) of the simulated AOU in the LGM experiments (mmol m−3): expLGM (a, d), expLGMws (b, e), and expLGMss (c, f). The differences between each LGM experiment and expPI are shown. The average states of the last 50 years are shown.


Figure 8Meridional sections in the Atlantic (a–c; at 30 W) and in the Pacific (d–f; at 150 W) of the simulated ideal age of seawater in the LGM experiments (years): expLGM (a, d), expLGMws (b, e), and expLGMss (c, f). The differences between each LGM experiment and expPI are shown. The average states of the last 50 years are shown.


The fact that the increased phosphate concentration in the deep water was mainly caused by an increase in remineralized phosphate was also supported by the higher fraction of remineralized phosphate in the total phosphate in the LGM experiments than in expPI (not shown). At the depths of 2000–3500 m in the Atlantic Ocean, expLGMws and expLGMss had an older ideal age than expPI (Fig. 8b, c), suggesting that a more stagnant deep water resulted in the longer-lasting storage of remineralized organic matter, hence the increased phosphate. On the other hand, at the greater depths in the Atlantic and in most parts of the Pacific section below 2000 m, the effect of the increased soft-tissue pump in the mid-latitudes to high latitudes of the Southern Hemisphere (Fig. 10a–c) and the subsequent northward transport of remineralized nutrients by the bottom circulation prevailed over the effect of less-efficient storage of nutrients by the younger water.

3.4 Carbonate ion concentration

The simulated carbonate ion concentration [CO32-] in all the experiments but expLGM showed a good agreement with the observation-based data in the Atlantic Ocean for the respective time periods, especially at depths deeper than 2500 m, while the distribution suggested a penetration of the northern-sourced water in expLGM that was too deep (Fig. 9a–d). The increased storage of remineralized carbon in the deep water for expLGMss and expLGMws resulted in lower [CO32-] accompanied by a lower pH of seawater as compared to those in expLGM, and then led to a larger vertical gradient of the concentration. However, the differences in the concentration between LGM and PI were larger than the reconstructed differences by Yu et al. (2020) (Fig. 9e–g). The effect of the increased storage of remineralized carbon was superposed onto the general rise of [CO32-] caused by the globally increased alkalinity, so that [CO32-] was higher in all the LGM experiments than in expPI in general in the Atlantic Ocean and the Pacific Ocean as well. The elevated [CO32-] was not well in accord with the estimated change in the deep Pacific that shows only a little change in [CO32-] of no more than 5 mmol m−3 (Yu et al.2013), although the estimated increase in [CO32-] by ∼25 mmol m−3 in the Weddell Sea (Rickaby et al.2010) was consistent with the results of expLGMws and expLGMss.

Figure 9Meridional sections in the Atlantic (a–g; at 30 W) and in the Pacific (h–j; at 150 W) of the simulated carbonate ion concentrations (mmol m−3). The absolute values for (a) expPI, (b) expLGM, (c) expLGMws, and (d) expLGMss and the differences between each LGM experiment and expPI are shown (e–j). The dots indicate the reconstruction by Yu et al. (2020). The average states of the last 50 years are shown.


3.5 Export production

Although the global sum of the export production in the LGM experiments was similar to that in expPI (Table 3), the spatial distributions showed a remarkable difference between them (Fig. 10). On the one hand, in all the LGM experiments, the modelled carbon export exhibited a notable increase compared to the modern counterpart in the so-called high-nutrient, low-chlorophyll (HNLC) regions: the eastern equatorial Pacific and the mid-latitude band (40–60 S) in the Southern Hemisphere, in particular. The increased atmospheric dust input in the LGM experiments boosted the biological production in those regions by iron fertilization (Fig. 10a–c), which was accompanied by more vigorous consumption of macronutrients, resulting in significantly smaller concentrations of nitrate at the surface in those regions (Fig. 10d–f). The increase in carbon export was also seen in the upwelling region along the west coast of the African continent. This noticeable glacial increase in the carbon export in those regions was consistent with the reconstructed relative changes in export production provided by Kohfeld et al. (2005). There was no clear correlation between the distribution of carbon-export anomaly and the sea ice extent (Fig. 10a–c). Reduction of the annual carbon export due to the expanded sea ice as discussed in previous studies (e.g. Kurahashi-Nakamura et al.2007; Sun and Matsumoto2010; Gupta et al.2020) was not significant in the Southern Hemisphere, which suggested that the dust effect of increasing the production overwhelmed the sea-ice effect. On the other hand, in the northern North Atlantic, the distributions of carbon export distinguished expLGMws from the other LGM runs. In expLGMws, the comparatively inactive deep convection in the northern North Atlantic caused a lower supply of macronutrients to the surface water, which resulted in a smaller amount of carbon export (Fig. 10b) and a smaller nitrate concentration (Fig. 10e) in that region. In the other two LGM experiments, on the contrary, the more vigorous vertical mixing stimulated biological production, leading to an increase in carbon export and less depleted nitrate concentrations.

Figure 10Carbon export at a depth of 100 m (a–c; gC m−2 yr−1) and the concentration of nitrate at the surface (d–f; mmol m−3). The differences from those in expPI are shown for each of the LGM runs: expLGM (a, d), expLGMws (b, e), and expLGMss (c, f). In panels (a)(c), the maximum extent of sea ice during the winter season of the respective hemisphere is also shown as solid lines (LGM) and dotted lines (PI). The average states in the last 50 years are shown.

3.6 CaCO3 in the upper sediments

The global sum of the mass accumulation rate (MAR) of CaCO3 for each experiment was 0.094 GtC yr−1 (expPI), 0.14 GtC yr−1 (expLGM), 0.12 GtC yr−1 (expLGMws), and 0.087 GtC yr−1 (expLGMss). Although the simulated modern MAR is ∼25 % smaller than the estimate (∼0.13 GtC yr−1) by Cartapanis et al. (2018), the LGM values except in expLGM approximated the estimate within the 2σ uncertainties for the glacial period (mean: ∼0.11 GtC yr−1, σ: ∼0.024 GtC yr−1) by the same study.

The spatial distributions of the simulated CaCO3 MAR also reasonably agreed with sedimentary data in the Atlantic Ocean, the Indian Ocean, and the Southern Ocean for all experiments (Fig. 11a–d), although the model underestimated the MAR in the eastern Pacific in expPI as in our preceding study (Kurahashi-Nakamura et al.2020). The ratio of the glacial MAR to the modern value showed more distinguishable differences among the LGM experiments (Fig. 11e–g). In the North Atlantic, the reconstruction indicated a significantly lower MAR in the mid-latitude North Atlantic Ocean, which was only well reproduced in expLGMws. The anomalies of the CaCO3 MAR in that region were governed by the supply of CaCO3 to the sediments rather than its preservation and reflected a direct influence of the different magnitude of vertical mixing by local deep convection on the production in line with the results shown in Sect. 3.5. The export of CaCO3 also depended on the sea-ice distribution (see Fig. 10a–c), because, in the BEC model, CaCO3 production is scaled by the difference between local seawater temperature and the freezing point of seawater (Moore et al.2004). In the South Atlantic, the MAR ratio was generally higher than the reconstruction due to the overestimated carbonate ion concentrations that provided more favourable environments for the preservation of CaCO3. Therefore, expLGMss that had the least additional alkalinity showed the smallest ratio of MAR of the three LGM experiments. In the Southern Ocean, the simulated MAR of CaCO3 was very low in all experiments, which resulted from low CaCO3 productivity compared to opal fixation.

Figure 11Mass accumulation rate (MAR) of CaCO3 in the upper sediment simulated with MEDUSA. The left column shows the absolute values (g cm−2 kyr−1) for (a) expPI, (b) expLGM, (c) expLGMws, and (d) expLGMss on a logarithmic scale with overlaid dots indicating the observation-based data by Cartapanis et al. (2018). The ratio of MAR in the respective LGM runs to that in the expPI run is shown in the right panels: (e) expLGM, (f) expLGMws, and (g) expLGMss. For the ratio, the locations with MAR less than 1×10-10 g cm2 kyr−1 were removed from the plot domain to highlight regions of large contrasts.

The spatial distributions of the CaCO3 weight fraction in the upper sediment simulated by MEDUSA did not show distinctive differences among the LGM experiments (not shown). A key proxy-based criterion to assess the model results is the fact that the weight fraction for the LGM is lower in the Atlantic Ocean and higher in the Pacific Ocean than for the PI (Catubig et al.1998), and all three LGM simulations satisfied this requirement despite the quite different global ocean states.

4 Discussion

4.1 Alkalinity and carbon inventories

A prerequisite for an acceptable marine biogeochemical state for the LGM is that the inventories of total alkalinity and DIC are compatible with the PI inventories. In other words, the differences in those ocean inventories between the two time periods need to be consistent with the changes in the inventories of the other reservoirs across the deglaciation. In the three LGM experiments of this study, 40 to 80 meq m−3 of alkalinity was added in addition to the increase due to the change in the seawater volume (Table 1). The excess of alkalinity corresponds to 0.5×1017 to 1.0×1017 equivalent in terms of the inventory in the ocean.

Assuming that those amounts of alkalinity were removed from the ocean by a net deposition of CaCO3, 2.5×1016 to 5.0×1016 mol of CaCO3 needed to be transferred to a different reservoir to account for the alkalinity decrease. Estimated amounts of the post-glacial shallow-water deposition of CaCO3 by coral reef buildup following the sea-level rise range from 2.2×1016 to 18×1016 mol (Milliman1993; Opdyke2000; Ridgwell et al.2003; Vecsei and Berger2004; Husson et al.2018; Köhler and Munhoven2020), which would suggest that the extra alkalinity could be readily removed by that process. This is also supported by another estimation for the modern inorganic carbon burial in shallow-water environments given by Cartapanis et al. (2018). The burial of CaCO3 in deep-sea sediments in the pelagic oceans and the influx from the chemical weathering on land are other processes that affect the alkalinity inventory. Although the total alkalinity flux by these two processes during the past 20 kyr is uncertain, a simple estimate gives a net increase in alkalinity by 3×1017 equivalent, if we assume the input by land weathering is 3.5×1016 eq kyr−1 and the removal by deep burial is 2×1016 eq kyr−1 on average as discussed in Cartapanis et al. (2018). To counteract the net increase in alkalinity, 15×1016 mol of CaCO3 needs to be deposited, and the range of the estimated amount of CaCO3 deposition by coral reef buildup covers the value as shown above (Fig. 12a).

Figure 12Processes discussed in Sect. 4.1 that contributed to the change in the ocean inventories of alkalinity (a) and carbon (b) since the LGM. Positive values indicate addition of alkalinity or carbon to the ocean and negative values removal of them. For the component “coral-reef build up”, the range of estimates by previous studies (Milliman1993; Opdyke2000; Ridgwell et al.2003; Vecsei and Berger2004; Husson et al.2018; Köhler and Munhoven2020) is shown. “Land weathering of carbonate”, “deep burial of carbonate”, and “shallow deposition of organic carbon” are based on the time-average fluxes estimated by Cartapanis et al. (2018). “Change in the atmospheric reservoir”, “change in the terrestrial reservoir”, and “pyrite oxidation” show the amount from ice-core records (Bereiter et al.2015; Köhler et al.2017), the range of estimates by previous studies (Kemppinen et al.2019), and the estimate by Kölling et al. (2019), respectively. The dotted lines correspond to the lowest estimate for the respective components.


Besides alkalinity, the compatibility with regard to the carbon inventory is examined as well (Fig. 12b). Out of the added 0.9×1017 to 1.2×1017 mol C, which corresponds to 72 to 94 mmol m−3, 2.5–5.0×1016 mol was to be removed by the shallow-water deposition of CaCO3. The rest (4.0×1016 to 9.5×1016 mol, or 480 to 1140 GtC), therefore, needs to be removed by other processes. The most probable one is the growth of the atmospheric reservoir. The increase in pCO2 by 90 ppm between the LGM and PI corresponds to a 190 GtC expansion in the size of this reservoir. Another potential pool would be the terrestrial carbon reservoir. The modelled terrestrial reservoir in this study indicates a growth of 340 to 390 GtC. The sum of the changes in both the reservoirs is 530–580 GtC, which leaves −100 GtC to 610 GtC (i.e. 480–580 to 1140–530) to be explained by further processes.

A plausible process to fill the gap could be again the post-glacial shallow-water deposition. An estimate of the modern deposition flux of organic carbon on shelves is 50–500 GtC kyr−1 (Cartapanis et al.2018), which might suffice to remove the excess as long as it is a typical amount during a substantial portion of the Holocene, although there are no direct records of the organic carbon deposition on continental margins over time (Cartapanis et al.2018). Presumably, CO2 release through pyrite oxidation on exposed continental shelves during a glacial sea-level lowstand (Kölling et al.2019) would add to the complexity of the change in carbon reservoirs across the deglaciation. It is expected that the sulfuric acid produced during pyrite oxidation is buffered with the dissolution of calcium carbonate minerals to lead to the release of CO2. The associated CO2 release rate is estimated to be between 12 and 36 GtC kyr−1 if the acid is completely buffered by carbonate mineral dissolution. The CO2 release triggered by the erosion and oxidation of pyrite that happened underground during the LGM might be delayed due to the transfer of CO2 to the surface. For the LGM, a model-based estimate of the cumulative CO2 release by buffering acids was 140 GtC (Kölling et al.2019). As an extreme case, even if all of the released CO2 reached the atmosphere during the deglaciation with a delay, that volume might be readily counteracted by the organic carbon deposition referred to above (50–500 GtC kyr−1; Cartapanis et al.2018). Another highly uncertain aspect is the relative change of the terrestrial reservoir size. Various previous studies compiled in Kemppinen et al. (2019) and Jeltsch-Thömmes et al. (2019) suggest that the growth of the terrestrial carbon reservoir from LGM to PI spans −500 to 1500 GtC, although the majority supports a (positive) growth toward PI. The residual carbon excess (i.e. 480–1140 GtC) that needs to be removed apart from the contribution of the CaCO3 deposition is in the range of these additional processes. To sum up, the removal of the whole excess of carbon is consistent with a combined contribution of the mechanisms mentioned above.

4.2 Simulated glacial biogeochemical states

Past studies attempted to constrain the LGM ocean state by utilizing the distribution of the stable carbon isotope ratio δ13CDIC (e.g. Tagliabue et al.2009; Menviel et al.2017; Kurahashi-Nakamura et al.2017; Muglia et al.2018; Gu et al.2020; Kobayashi et al.2021; Muglia and Schmittner2021). It turned out that this tracer is useful in determining the spatial distribution of water masses but much less in terms of the strength of the AMOC. Our study supports this conclusion. The noticeable negative anomaly of δ13CDIC in the deep ocean (deeper than ∼2500 m) found in the observation-based reconstructions (Peterson et al.2014; Oppo et al.2018; Yu et al.2020) was reproduced only in the case of a shallower AMOC, no matter what the strength of the AMOC is (expLGMss and expLGMws), which suggests that an ocean state having a shallower AMOC and northern-source deep water would be much preferable. A similar contrast between these shallower-AMOC LGM states and the deeper-AMOC LGM (expLGM) state is visible in the phosphate distribution, and again the shallower states show a much better correspondence to the reconstruction by Oppo et al. (2018), backing up the argument based on δ13CDIC.

Another difference between the shallower and deeper AMOC states is found in the [CO32-] fields in the North Atlantic. The modelled distribution of [CO32-] in the shallower-AMOC experiments (expLGMss and expLGMws) indicates a vertical gradient between the depths of 1000 and 4000 m that is noticeably larger than in expPI: larger by ∼50 mmol m−3 in expLGMss and by ∼70 mmol m−3 in expLGMws, while only larger by ∼10 mmol m−3 in expLGM. The reconstruction by Yu et al. (2020) shows a ∼50 mmol m−3 larger vertical gradient for the LGM in the North Atlantic, which again supports the modelled states with a shallower AMOC. An increase in the vertical gradient by a similar magnitude is also reported by Chalk et al. (2019). The increased vertical contrast in expLGMss and expLGMws is brought about mainly by the larger storage of carbon in the deep water due to a more sluggish circulation (2000–3500 m) and the larger amount of carbon export in the Southern Hemisphere. Contrary to the vertical gradient, however, the LGM experiments in this study overestimate the values of [CO32-] themselves in the North Atlantic by 20–50 mmol m−3 compared to Yu et al. (2020). This would be, at least partly, due to the uniformly increased alkalinity because such a systematic bias does not appear clearly in expPI. In the Southern Ocean, the same dataset by Yu et al. (2020) shows a slightly higher [CO32-] (by 10–20 mmol m−3) in the bottom water for the LGM than for the modern day, also supporting the results of expLGMss and expLGMws. On the other hand, the negative anomaly of [CO32-] in the South Atlantic caused by the expanded Glacial Pacific Deep Water proposed by Yu et al. (2020) is not reproduced in the experiments in our study.

Although the amounts of additional alkalinity applied in this study are in good agreement with previous estimates from the mass-balance point of view (see Sect. 4.1), the simulated [CO32-] is systematically too high, most likely due to the uniformly increased alkalinity. This issue has two aspects. First, to manage the compatibility of the 190 ppm and more reasonable carbonate ion concentrations, one needs to realize the low pCO2 with a smaller amount of added alkalinity. Supplementary LGM simulations without additional alkalinity resulted in higher pCO2 by 20 to 43 ppm depending on the simulations compared to the respective main simulations. The additional simulations depict ocean states that satisfy the constraint of deep-ocean carbon reservoir but lack the contribution of alkalinity increase. Although the contribution of additional alkalinity to pCO2 drawdown seems to be minor compared to the whole change in pCO2 since the LGM, as suggested by Ridgwell et al. (2003), other mechanisms would be required to achieve 190 ppm without (or with a smaller amount of) additional alkalinity. For example, higher solubility given by lower SST, a larger vertical contrast of DIC concentration by an even more stratified ocean (e.g. Kobayashi et al.2021), and/or larger carbon storage in the deep water by a stronger biological pump (e.g. Morée et al.2021). The more efficient carbon storage given by these processes would relax the problem of too-high carbonate ion concentrations. Second, the compatibility with the post-glacial shallow-water deposition of CaCO3 would need to be satisfied, too. In expLGMss that needed the smallest amount of added alkalinity of the three LGM experiments, the applied alkalinity corresponded to 2.5×1016 mol of CaCO3, which is already close to the lower limit of the independently estimated amounts of the shallow-water deposition (i.e. 2.2×1016 mol). Therefore, to incorporate the likely post-glacial deposition of CaCO3 and accompanying reduction of alkalinity into the evolution of the climate from the LGM to the modern day, another source of alkalinity might be at work. More dissolution of CaCO3 in the deep-ocean sediments or more input of alkalinity from the land weathering would be able to serve as the alkalinity source as discussed in Sect. 4.1. Future studies need to deal with the temporal variations in the global CaCO3 budget since the LGM, which is required for a more accurate discussion of the mass balance during the last 21 kyr.

This study utilized the properties of upper sediments with regard to CaCO3 burial as a further criterion to distinguish among different physical ocean states. The preservation or dissolution of CaCO3 in ocean-floor sediments is strongly affected by the carbonate chemistry of the surrounding seawater. Therefore, the degree of CaCO3 preservation that is observed in the sediments is expected to contain local information on the properties of the seawater above the sediments. Although in this context the solid weight fraction of CaCO3 is a typical measure of the model–data fit, interpretation is not always straightforward because it is also affected by the amount of other solid species in the sediment. For example, even if the amount of CaCO3 does not change, the increased dust input during the LGM would reduce the corresponding CaCO3 weight fraction by “diluting” the contribution of CaCO3. On the contrary, the MAR would not be affected by other solid components. In expLGMws, the CaCO3 MAR in the North Atlantic Ocean (40–65 N) is lower than in expPI (Fig. 11), although the higher [CO32-] in the bottom water is an advantage for the preservation of CaCO3 in the upper sediments (Fig. 9). This fact suggests that the negative anomaly in the CaCO3 MAR in that region is controlled by the supply of CaCO3 to the sediments rather than its preservation, and it reflects a direct influence of the different (i.e. weaker) magnitude of vertical mixing by local deep convection, leading to a lower supply of nutrients to the surface. Although the strength of deep convection is not directly related to the intensity of the meridional transport, a weaker convection may suggest a lower sea-surface density in the North Atlantic Ocean in some way and coincide with a reduced rate of meridional transport, considering the positive correlation between the AMOC strength and the meridional density gradient across the Atlantic Ocean (e.g. Rahmstorf1996).

In the LGM simulations of this study, the global carbon export by biological production is broadly similar to that in expPI (Table 3), which is in agreement with the estimates by previous model-based studies (Bopp et al.2003; Tagliabue et al.2009; Oka et al.2011; Schmittner and Somes2016). The regional characteristic features of the modelled LGM carbon export (Sect. 3.3) are well supported by other proxy-based studies. The increased nitrate consumption and carbon export in the Southern Ocean are consistent with the estimate of glacial nutrient consumption based on nitrogen isotope ratios (e.g. Martínez-García et al.2014; Kohfeld and Chase2017; Wang et al.2017). Moreover, the decreased productivity in the northern North Atlantic Ocean obtained in expLGMws is in accordance with a reconstruction based on the analysis of dinocyst assemblages (Radi and de Vernal2008). This reconstruction suggests an annual productivity that is lower by 50–150 gC m−2 during the LGM, which is reproduced reasonably well in the experiment. As discussed in the previous paragraph, the reduced productivity is also consistent with the MAR data. The change in the distribution of biological production naturally influences the distribution of phosphate and AOU. The distinctive increase in carbon export in the Southern Ocean brings larger concentrations of both tracers in the bottom water in spite of its younger ideal age compared to expPI, which contributes to the good model–data agreement with the reconstructed phosphate distribution in the Atlantic Ocean. A similar feature is found in the very deep water in the Pacific section, while the increase in the export production in the Northern Hemisphere together with the older ideal age produces very high anomalies of phosphate and AOU in the depths of 500–1500 m of the North Pacific Ocean.

The younger ideal age in the Southern Ocean of the LGM experiments is caused by the changes in local convective mixing rather than by the changes in Antarctic bottom water (AABW) flow. ExpLGM has a very similar magnitude and geometry of AABW to that in expPI (Fig. 3), but nevertheless the ideal age in the Southern Ocean is significantly younger. On the other hand, in expLGMws and expLGMss, more vigorous AABW does not convey the comparatively old water at the depths of 2000–3500 m in the Atlantic to the south of ∼40 S. In addition, the LGM experiments have a deeper mixed-layer depth in the Southern Ocean, which would contribute to the better ventilated Southern Ocean.

4.3 Relationship between the DIC concentrations and water-mass ages

The estimates of the DIC inventories in the deep ocean during the LGM that we used in this study for model calibration are based on the assumption that there is a linear relationship between the local water age and DIC concentrations in the modern deep ocean below the depth of 2000 m that is also applicable to the glacial ocean (Sarnthein et al.2013). The same assumption is employed in another study for estimating the size of the deep-ocean carbon pool (Skinner et al.2015). We simulated the carbon age by adding reservoir ages estimated for the modern surface water in different basins (Matsumoto2007) and the estimated increase in the surface age in the LGM (Skinner et al.2017) to the explicitly modelled ideal age, and we examined the DIC–age relationship in the model states (Fig. 13). An analogous quasi-linear relationship also appears in our model results from the pre-industrial run (Fig. 13a). It turns out, however, that the LGM model oceans have a different structure of the DIC–age relationship (Fig. 13b–d): in the Pacific and Southern oceans the slope of the linear regression lines is less steep, and the intercept is substantially larger, while in the Atlantic Ocean it clearly has a steeper slope.

Figure 13Scatter plots showing the correlations between DIC concentrations and the age of local water for the Atlantic (red), the Pacific (blue), and the Southern Ocean (green). The age was calculated by adding estimated reservoir ages for the surface water (Matsumoto2007; Skinner et al.2017) to the modelled ideal age. The linear-regression coefficients are also shown at the lower right of each panel, where x denotes the simulated age (years) and y the concentration of DIC (mmol m−3). The upper (lower) four panels show the properties of water that is deeper than 2000 m (3000 m).


The reasons for the different structures are twofold. First, the remarkably different slopes of the regression lines for the Atlantic Ocean and the other two oceans are caused by the DIC enrichment of the southern-sourced deep water due to the reinforcement of the biological soft-tissue pump by the increased dust input to the high-latitude Southern Ocean. The carbon-rich southern-sourced deep water strongly influences the Southern and Pacific oceans, so that these two oceans are generally high in DIC compared to the modern deep water. On the other hand, the properties of the Atlantic water can be explained by the mixing of northern-sourced water and the DIC-enriched southern-sourced water, leading to the apparent steeper gradient for the Atlantic Ocean.

The separation of the Atlantic branch from the others (Fig. 13b–d) becomes invisible by confining the depth domain to deeper ranges in some LGM runs, so that an overall linear relationship with a relatively similar slope as in the modern case is obtained. The mixing effect becomes invisible in similar scatter plots but for water deeper than 3000 m for the LGM simulations with a shallower AMOC (Fig. 13g, h), while the northern-sourced component still emerges in the LGM run with a deeper AMOC (Fig. 13f) because the northern-sourced water reaches the depth of 3000 m. In the pre-industrial setting, a very similar relationship is still valid also for the water deeper than 3000 m (Fig. 13e).

Besides the influence of the enriched southern-sourced water, as a second factor, the prior addition of mean DIC to the entire ocean in the LGM simulations would raise the regression lines uniformly. Considering a possible alteration of the DIC–age correlation induced by these two factors, applying the modern regression line to the glacial condition does not necessarily provide an accurate estimate of the glacial mean DIC concentration in the deep ocean.

In the instances of this study, the DIC concentration according to the modern regression line (Fig. 13a) and those based on the glacial lines (Fig. 13b–d) are different by approximately ∼100 mmol m−3 when the mean age of deep water is 2500 years. More than half of the difference is caused by the modified mean DIC concentration in the entire ocean (see Table 1), and the rest results from the carbon enrichment due to the increased biological pump. Unfortunately, it is not straightforward to estimate a more accurate DIC–age relationship for the LGM in the framework of this study because the magnitude of the two influencing factors, especially the mean DIC offset, is uncertain. We adjusted the mean total alkalinity so that it was compatible with the prescribed mean DIC and the 190 ppm in the atmosphere. Although this means that the mean DIC concentration, in turn, can be tuned to yield the 190 ppm with a prescribed alkalinity inventory, it would also be difficult to have a reliable estimate of the alkalinity inventory during the LGM due to the large uncertainty of the observation-based estimate of the alkalinity budget by the post-glacial CaCO3 deposition (see Sect. 4.1).

5 Conclusions and outlook

Three time-slice simulations for the LGM with a comprehensive Earth-system model including a global carbon cycle module were carried out, so that reasonable biogeochemical states in terms of the fit to paleoclimatological and paleoceanographic records have been achieved by tuning of the total inventories of DIC and alkalinity in the ocean. The simulated ocean states are expected to serve as an initial state for future transient simulations of the last deglaciation, i.e. from the glacial (LGM) to the interglacial (PI) climate state. In terms of biogeochemical tracer distributions in the ocean, our results clearly show that the LGM ocean states with a shallower AMOC are characterized by tracer distributions that can be more easily reconciled with the relevant reconstructions as suggested by various previous studies. Model–data comparisons regarding CaCO3 properties in the upper sediments add further support to a state with a weaker and shallower AMOC rather than a stronger circulation, although they do not directly constrain the magnitude of the volume transport.

The examination of the balance of the bulk ALK and DIC budgets shows that all the LGM simulations in this study are likely to be compatible with the PI state within the uncertainties of the available constraints, at least from a mass balance point of view. The question of whether they would indeed reach the PI state following a realistic trajectory would need to be examined in a transient context, because the trajectory would depend on the timing and magnitude of the deglacial sea-level rise that governs the post-glacial deposition on the shelves. Considering that the sea-level rise is a direct consequence of the ice sheet evolution induced by climate changes, it would be fundamental to analyse and discuss the evolution of the coupled carbon cycle–ice sheet system.

As a connection between these two subsystems, shelf or shallow-water processes including the weathering and deposition of biogeochemical matter and their modelling (e.g. Munhoven and François1996; Kölling et al.2019; Börker et al.2020; Lacroix et al.2020) would take an important role. Those processes would be fundamental to the evolution of the global carbon cycle at the glacial–interglacial timescale not only for the post-glacial evolution but also during the evolution in the time period that preceded the LGM. The excess alkalinity was prescribed and added to the ocean in the LGM experiments of this study to satisfy independent observation-based constraints, but in reality, it would be determined as the cumulative imbalance between incoming and outgoing fluxes during the course of the evolution from the last interglacial period to the LGM, which will be another target of transient simulations in the future.

Code and data availability

The newly developed model source codes to tailor CESM1.2 for experiments with additional freshwater forcing and model output of the main experiments will be available at (Kurahashi-Nakamura et al.2022).

Author contributions

TKN developed the model code for the additional freshwater forcing with input from UM. TKN and AP designed the experiments, and TKN carried them out. TKN interpreted and discussed the results with contributions from all co-authors. AP, UM, and MS conceptualized the overarching research project. TKN prepared the manuscript with contributions from all co-authors.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Climate of the Past. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to thank the reviewers for their insightful comments and suggestions. This research was funded by the project PalMod (, last access: 25 August 2022; FKZ: 01LP1505D) within the framework of Research for Sustainable Development (FONA,, last access: 25 August 2022) by the German Federal Ministry for Education and Research (BMBF).

Financial support

This research has been supported by the Bundesministerium für Bildung und Forschung (FKZ: grant no. 01LP1505D).

The article processing charges for this open-access publication were covered by the University of Bremen.

Review statement

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


Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume, Nature, 500, 190–193,, 2013. a

Adkins, J. F., McIntyre, K., and Schrag, D. P.: The Salinity, Temperature, and δ18O of the Glacial Deep Ocean, Science, 298, 1769–1773,, 2002. a, b

Albani, S., Mahowald, N. M., Perry, A. T., Scanza, R. A., Zender, C. S., Heavens, N. G., Maggi, V., Kok, J. F., and Otto-Bliesner, B. L.: Improved dust representation in the Community Atmosphere Model, J. Adv. Model. Earth Sy., 6, 541–570,, 2014. a

Annan, J. D. and Hargreaves, J. C.: A new global reconstruction of temperature changes at the Last Glacial Maximum, Clim. Past, 9, 367–376,, 2013. a

Archer, D., Winguth, A., Lea, D., and Mahowald, N.: What caused the glacial/interglacial atmospheric pCO2 cycles?, Rev. Geophys., 38, 159–189,, 2000. a

Argus, D. F., Peltier, W. R., Drummond, R., and Moore, A. W.: The Antarctica component of postglacial rebound model ICE-6G_C (VM5a) based on GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories, Geophys. J. Int., 198, 537–563,, 2014. a

Armstrong, R. A., Lee, C., Hedges, J. I., Honjo, S., and Wakeham, S. G.: A new, mechanistic model for organic carbon fluxes in the ocean based on the quantitative association of POC with ballast minerals, Deep-Sea Res. Pt. II, 49, 219–236,, 2002. a

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, b, c, d, e

Berger, W.: Deglacial CO2 buildup: Constraints on the coral-reef model, Palaeogeogr. Palaeocl., 40, 235–253,, 1982. a

Bopp, L., Kohfeld, K. E., Le Quéré, C., and Aumont, O.: Dust impact on marine biota and atmospheric CO2 during glacial periods, Paleoceanography, 18, 1046,, 2003. a

Börker, J., Hartmann, J., Amann, T., Romero-Mujalli, G., Moosdorf, N., and Jenkins, C.: Chemical Weathering of Loess and Its Contribution to Global Alkalinity Fluxes to the Coastal Zone During the Last Glacial Maximum, Mid-Holocene, and Present, Geochem. Geophy. Geosy., 21, e2020GC008922,, 2020. a

Boudreau, B. P., Middelburg, J. J., and Luo, Y.: The role of calcification in carbonate compensation, Nat. Geosci., 11, 894–900,, 2018. a, b

Broecker, W. S. and Peng, T.-H.: The role of CaCO3 compensation in the glacial to interglacial atmospheric CO2 change, Global Biogeochem. Cy., 1, 15–29,, 1987. a

Brovkin, V., Ganopolski, A., Archer, D., and Rahmstorf, S.: Lowering of glacial atmospheric CO2 in response to changes in oceanic circulation and marine biogeochemistry, Paleoceanography, 22, PA4202,, 2007. a

Cartapanis, O., Galbraith, E. D., Bianchi, D., and Jaccard, S. L.: Carbon burial in deep-sea sediment and implications for oceanic inventories of carbon and alkalinity over the last glacial cycle, Clim. Past, 14, 1819–1850,, 2018. a, b, c, d, e, f, g, h, i, j

Catubig, N. R., Archer, D. E., Francois, R., Demenocal, P., Howard, W., and Yu, E.-F.: Global deep-sea burial rate of calcium carbonate during the Last Glacial Maximum, Paleoceanography, 13, 298–310,, 1998. a

Chalk, T. B., Foster, G. L., and Wilson, P. A.: Dynamic storage of glacial CO2 in the Atlantic Ocean revealed by boron [CO32-] and pH records, Earth Planet. Sc. Lett., 510, 1–11,, 2019. a

Chikamoto, M. O., Matsumoto, K., and Ridgwell, A.: Response of deep-sea CaCO3 sedimentation to Atlantic meridional overturning circulation shutdown, J. Geophys. Res., 113, G03017,, 2008. a

EPICA community members: Eight glacial cycles from an Antarctic ice core, Nature, 429, 623–628,, 2004. a

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. a, b

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

Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425,, 2011. a

Gu, S., Liu, Z., Oppo, D. W., Lynch-Stieglitz, J., Jahn, A., Zhang, J., and Wu, L.: Assessing the potential capability of reconstructing glacial Atlantic water masses and AMOC using multiple proxies in CESM, Earth Planet. Sc. Lett., 541, 116294,, 2020. a, b

Gupta, M., Follows, M. J., and Lauderdale, J. M.: The Effect of Antarctic Sea Ice on Southern Ocean Carbon Outgassing: Capping Versus Light Attenuation, Global Biogeochem. Cy., 34, e06489,, 2020. a

Hain, M., Sigman, D., and Haug, G.: 8.18 – The Biological Pump in the Past, in: Treatise on Geochemistry (Second Edition), edited by: Holland, H. D. and Turekian, K. K., Elsevier, 8, 485–517,, 2014. a

Homola, K., Spivack, A. J., Murray, R. W., Pockalny, R., D'Hondt, S., and Robinson, R.: Deep North Atlantic Last Glacial Maximum Salinity Reconstruction, Paleoceanography and Paleoclimatology, 36, e2020PA004088,, 2021. a, b

Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., Lamarque, J.-F., Large, W. G., Lawrence, D., Lindsay, K., Lipscomb, W. H., Long, M. C., Mahowald, N., Marsh, D. R., Neale, R. B., Rasch, P., Vavrus, S., Vertenstein, M., Bader, D., Collins, W. D., Hack, J. J., Kiehl, J., and Marshall, S.: The Community Earth System Model: A Framework for Collaborative Research, B. Am. Meteorol. Soc., 94, 1339–1360,, 2013. a

Husson, L., Pastier, A.-M., Pedoja, K., Elliot, M., Paillard, D., Authemayou, C., Sarr, A.-C., Schmitt, A., and Cahyarini, S. Y.: Reef Carbonate Productivity During Quaternary Sea Level Oscillations, Geochem. Geophy. Geosy., 19, 1148–1164,, 2018. a, b, c

Insua, T. L., Spivack, A. J., Graham, D., D'Hondt, S., and Moran, K.: Reconstruction of Pacific Ocean bottom water salinity during the Last Glacial Maximum, Geophys. Res. Lett., 41, 2914–2920,, 2014. a, b

Jahn, A., Lindsay, K., Giraud, X., Gruber, N., Otto-Bliesner, B. L., Liu, Z., and Brady, E. C.: Carbon isotopes in the ocean model of the Community Earth System Model (CESM1), Geosci. Model Dev., 8, 2419–2434,, 2015. a, b

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

Kageyama, M., Albani, S., Braconnot, P., Harrison, S. P., Hopcroft, P. O., Ivanovic, R. F., Lambert, F., Marti, O., Peltier, W. R., Peterschmitt, J.-Y., Roche, D. M., Tarasov, L., Zhang, X., Brady, E. C., Haywood, A. M., LeGrande, A. N., Lunt, D. J., Mahowald, N. M., Mikolajewicz, U., Nisancioglu, K. H., Otto-Bliesner, B. L., Renssen, H., Tomas, R. A., Zhang, Q., Abe-Ouchi, A., Bartlein, P. J., Cao, J., Li, Q., Lohmann, G., Ohgaito, R., Shi, X., Volodin, E., Yoshida, K., Zhang, X., and Zheng, W.: The PMIP4 contribution to CMIP6 – Part 4: Scientific objectives and experimental design of the PMIP4-CMIP6 Last Glacial Maximum experiments and PMIP4 sensitivity experiments, Geosci. Model Dev., 10, 4035–4055,, 2017. a, b, c, d

Kemppinen, K. M. S., Holden, P. B., Edwards, N. R., Ridgwell, A., and Friend, A. D.: Coupled climate–carbon cycle simulation of the Last Glacial Maximum atmospheric CO2 decrease using a large ensemble of modern plausible parameter sets, Clim. Past, 15, 1039–1062,, 2019. a, b

Kobayashi, H., Oka, A., Yamamoto, A., and Abe-Ouchi, A.: Glacial carbon cycle changes by Southern Ocean processes with sedimentary amplification, Science Advances, 7, eabg7723,, 2021. a, b, c

Kohfeld, K. E. and Chase, Z.: Temporal evolution of mechanisms controlling ocean carbon uptake during the last glacial cycle, Earth Planet. Sc. Lett., 472, 206–215,, 2017. a

Kohfeld, K. E., Le Quéré, C., Harrison, S. P., and Anderson, R. F.: Role of Marine Biology in Glacial-Interglacial CO2 Cycles, Science, 308, 74–78,, 2005. a

Köhler, P. and Munhoven, G.: Late Pleistocene Carbon Cycle Revisited by Considering Solid Earth Processes, Paleoceanography and Paleoclimatology, 35, e2020PA004020,, 2020. a, b, c

Köhler, P., Nehrbass-Ahles, C., Schmitt, J., Stocker, T. F., and Fischer, H.: A 156 kyr smoothed history of the atmospheric greenhouse gases CO2, CH4, and N2O and their radiative forcing, Earth Syst. Sci. Data, 9, 363–387,, 2017. a, b, c, d, e

Kölling, M., Bouimetarhan, I., Bowles, M. W., Felis, T., Goldhammer, T., Hinrichs, K.-U., Schulz, M., and Zabel, M.: Consistent CO2 release by pyrite oxidation on continental shelves prior to glacial terminations, Nat. Geosci., 12, 929–934,, 2019. a, b, c, d

Kurahashi-Nakamura, T., Abe-Ouchi, A., Yamanaka, Y., and Misumi, K.: Compound effects of Antarctic sea ice on atmospheric pCO2 change during glacial-interglacial cycle, Geophys. Res. Lett., 34, L20708,, 2007. 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, b, c

Kurahashi-Nakamura, T., Paul, A., Munhoven, G., Merkel, U., and Schulz, M.: Coupling of a sediment diagenesis model (MEDUSA) and an Earth system model (CESM1.2): a contribution toward enhanced marine biogeochemical modelling and long-term climate simulations, Geosci. Model Dev., 13, 825–840,, 2020. a, b, c, d

Kurahashi-Nakamura, T., Paul, A., Merkel, U., and Schulz, M.: Model results and newly developed model code for time-slice carbon cycle simulations with CESM1.2, PANGAEA [code],, 2022. a

Lacroix, F., Ilyina, T., and Hartmann, J.: Oceanic CO2 outgassing and biological production hotspots induced by pre-industrial river loads of nutrients and carbon in a global modeling approach, Biogeosciences, 17, 55–88,, 2020. a

Lindsay, K., Bonan, G. B., Doney, S. C., Hoffman, F. M., Lawrence, D. M., Long, M. C., Mahowald, N. M., Keith Moore, J., Randerson, J. T., and Thornton, P. E.: Preindustrial-Control and Twentieth-Century Carbon Cycle Experiments with the Earth System Model CESM1(BGC), J. Climate, 27, 8981–9005,, 2014. a

Lord, N. S., Ridgwell, A., Thorne, M. C., and Lunt, D. J.: An impulse response function for the “long tail” of excess atmospheric CO2 in an Earth system model, Global Biogeochem. Cy., 30, 2–17,, 2016. 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

Matsumoto, K.: Radiocarbon-based circulation age of the world oceans, J. Geophys. Res., 112, C09004,, 2007. a, b

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

Merino, N., Le Sommer, J., Durand, G., Jourdain, N. C., Madec, G., Mathiot, P., and Tournadre, J.: Antarctic icebergs melt over the Southern Ocean: Climatology and impact on sea ice, Ocean Model., 104, 99–110,, 2016. a

Milliman, J. D.: Production and accumulation of calcium carbonate in the ocean: Budget of a nonsteady state, Global Biogeochem. Cy., 7, 927–957,, 1993. a, b, c

Moore, J. K., Doney, S. C., and Lindsay, K.: Upper ocean ecosystem dynamics and iron cycling in a global three-dimensional model, Global Biogeochem. Cy., 18, GB4028,, 2004. a, b

Moore, J. K., Lindsay, K., Doney, S. C., Long, M. C., and Misumi, K.: Marine Ecosystem Dynamics and Biogeochemical Cycling in the Community Earth System Model [CESM1(BGC)]: Comparison of the 1990s with the 2090s under the RCP4.5 and RCP8.5 Scenarios, J. Climate, 26, 9291–9312,, 2013. a

Morée, A. L., Schwinger, J., Ninnemann, U. S., Jeltsch-Thömmes, A., Bethke, I., and Heinze, C.: Evaluating the biological pump efficiency of the Last Glacial Maximum ocean using δ13C, Clim. Past, 17, 753–774,, 2021. a

Muglia, J. and Schmittner, A.: Carbon isotope constraints on glacial Atlantic meridional overturning: Strength vs depth, Quaternary Sci. Rev., 257, 106844,, 2021. a, b

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

Munhoven, G.: Mathematics of the total alkalinity–pH equation – pathway to robust and universal solution algorithms: the SolveSAPHE package v1.0.1, Geosci. Model Dev., 6, 1367–1388,, 2013. a, b

Munhoven, G.: Model of Early Diagenesis in the Upper Sediment with Adaptable complexity – MEDUSA (v. 2): a time-dependent biogeochemical sediment module for Earth system models, process analysis and teaching, Geosci. Model Dev., 14, 3603–3631,, 2021. a, b

Munhoven, G. and François, L. M.: Glacial-interglacial variability of atmospheric CO2 due to changing continental silicate rock weathering: A model study, J. Geophys. Res.-Atmos., 101, 21423–21437,, 1996. a

Oka, A., Abe-Ouchi, A., Chikamoto, M. O., and Ide, T.: Mechanisms controlling export production at the LGM: Effects of changes in oceanic physical fields and atmospheric dust deposition, Global Biogeochem. Cy., 25, GB2009,, 2011. a

Opdyke, B. N.: Shallow Water Carbonate Deposition and Its Effect on the Carbon Cycle, Cambridge University Press, 161–168,, 2000. a, b, c

Opdyke, B. N. and Walker, J. C.: Return of the coral reef hypothesis: Basin to shelf partitioning of CaCO3 and its effect on atmospheric CO2, Geology, 20, 733–736,<0733:ROTCRH>2.3.CO;2, 1992. a

Oppo, D. W., Gebbie, G., Huang, K.-F., Curry, W. B., Marchitto, T. M., and Pietro, K. R.: Data Constraints on Glacial Atlantic Water Mass Geometry and Properties, Paleoceanography and Paleoclimatology, 33, 1013–1034,, 2018. a, b, c, d

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

Peltier, W. R., Argus, D. F., and Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model, J. Geophys. Res. Sol. Ea., 120, 450–487,, 2015. a

Peterson, C. D., Lisiecki, L. E., and Stern, J. V.: Deglacial whole-ocean δ13C change estimated from 480 benthic foraminiferal records, Paleoceanography, 29, 549–563,, 2014. a, b, c

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

Radi, T. and de Vernal, A.: Last glacial maximum (LGM) primary productivity in the northern North Atlantic Ocean, Can. J. Earth Sci., 45, 1299–1316,, 2008. a

Rahmstorf, S.: On the freshwater forcing and transport of the Atlantic thermohaline circulation, Clim. Dynam., 12, 799–811,, 1996. a

Rickaby, R. E. M., Elderfield, H., Roberts, N., Hillenbrand, C.-D., and Mackensen, A.: Evidence for elevated alkalinity in the glacial Southern Ocean, Paleoceanography, 25, PA1209,, 2010. a

Ridgwell, A. J., Watson, A. J., Maslin, M. A., and Kaplan, J. O.: Implications of coral reef buildup for the controls on atmospheric CO2 since the Last Glacial Maximum, Paleoceanography, 18, 1083,, 2003. a, b, c, d, e

Sarnthein, M., Schneider, B., and Grootes, P. M.: Peak glacial 14C ventilation ages suggest major draw-down of carbon into the abyssal ocean, Clim. Past, 9, 2595–2614,, 2013. a, b, c, d, e, f

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

Schmittner, A. and Somes, C. J.: Complementary constraints from carbon (13C) and nitrogen (15N) isotopes on the glacial ocean's soft-tissue biological pump, Paleoceanography, 31, 669–693,, 2016. a

Shields, C. A., Bailey, D. A., Danabasoglu, G., Jochum, M., Kiehl, J. T., Levis, S., and Park, S.: The Low-Resolution CCSM4, J. Climate, 25, 3993–4014,, 2012. a

Sigman, D. M. and Boyle, E. A.: Glacial/interglacial variations in atmospheric carbon dioxide, Nature, 407, 859–869, 2000. a

Skinner, L., McCave, I., Carter, L., Fallon, S., Scrivner, A., and Primeau, F.: Reduced ventilation and enhanced magnitude of the deep Pacific carbon pool during the last glacial period, Earth Planet. Sc. Lett., 411, 45–52,, 2015. a, b

Skinner, L., Primeau, F., Freeman, E., de La Fuente, M., Goodwin, P. A., Gottschalk, J., Huang, E., McCave, I. N., Noble, T. L., and Scrivner, A. E.: Radiocarbon constraints on the glacial ocean circulation and its impact on atmospheric CO2, Nat. Commun., 8, 16010,, 2017. a, b

Sun, X. and Matsumoto, K.: Effects of sea ice on atmospheric pCO2: A revised view and implications for glacial and future climates, J. Geophys. Res., 115, G02015,, 2010. a

Tagliabue, A., Bopp, L., Roche, D. M., Bouttes, N., Dutay, J.-C., Alkama, R., Kageyama, M., Michel, E., and Paillard, D.: Quantifying the roles of ocean circulation and biogeochemistry in governing ocean carbon-13 and atmospheric carbon dioxide at the last glacial maximum, Clim. Past, 5, 695–706,, 2009.  a, b

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

Vecsei, A. and Berger, W. H.: Increase of atmospheric CO2 during deglaciation: Constraints on the coral reef hypothesis from patterns of deposition, Global Biogeochem. Cy., 18, GB1035,, 2004. a, b, c

Wang, X. T., Sigman, D. M., Prokopenko, M. G., Adkins, J. F., Robinson, L. F., Hines, S. K., Chai, J., Studer, A. S., Martínez-García, A., Chen, T., and Haug, G. H.: Deep-sea coral evidence for lower Southern Ocean surface nitrate concentrations during the last ice age, P. Natl. Acad. Sci. USA, 114, 3352–3357,, 2017. a

Yu, J., Anderson, R. F., Jin, Z., Rae, J. W. B., Opdyke, B. N., and Eggins, S. M.: Responses of the deep ocean carbonate system to carbon reorganization during the Last Glacial-interglacial cycle, Quaternary Sci. Rev., 76, 39–52,, 2013. a

Yu, J., Anderson, R., and Rohling, E.: Deep Ocean Carbonate Chemistry and Glacial-Interglacial Atmospheric CO2 Changes, Oceanography, 27, 16–25,, 2014. a

Yu, J., Menviel, L., Jin, Z. D., Anderson, R. F., Jian, Z., Piotrowski, A. M., Ma, X., Rohling, E. J., Zhang, F., Marino, G., and McManus, J. F.: Last glacial atmospheric CO2 decline due to widespread Pacific deep-water expansion, Nat. Geosci., 13, 628–633,, 2020. a, b, c, d, e, f, g, h, i

Short summary
With a comprehensive Earth-system model including the global carbon cycle, we simulated the climate state during the last glacial maximum. We demonstrated that the CO2 concentration in the atmosphere both in the modern (pre-industrial) age (~280 ppm) and in the glacial age (~190 ppm) can be reproduced by the model with a common configuration by giving reasonable model forcing and total ocean inventories of carbon and other biogeochemical matter for the respective ages.