Ocean carbon inventory under warmer climate conditions – the case of the Last Interglacial

During the Last Interglacial period (LIG), the transition from 125 to 115 ka provides a case study for assessing the response of the carbon system to different levels of high-latitude warmth. Elucidating the mechanisms responsible for interglacial changes in the ocean carbon inventory provides constraints on natural carbon sources and sinks and their climate sensitivity, which are essential for assessing potential future changes. However, the mechanisms leading to modifications of the ocean's carbon budget during this period remain poorly documented and not well understood. Using a state-of-the-art Earth system model, we analyze the changes in oceanic carbon dynamics by comparing two quasi-equilibrium states: the early, warm Eemian (125 ka) versus the cooler, late Eemian (115 ka). We find considerably reduced ocean dissolved inorganic carbon (DIC; −314.1 PgC) storage in the warm climate state at 125 ka as compared to 115 ka, mainly attributed to changes in the biological pump and ocean DIC disequilibrium components. The biological pump is mainly driven by changes in interior ocean ventilation timescales, but the processes controlling the changes in ocean DIC disequilibrium remain difficult to assess and seem more regionally affected. While the Atlantic bottom-water disequilibrium is affected by the organization of sea-ice-induced southern-sourced water (SSW) and northern-sourced water (NSW), the upper-layer changes remain unexplained. Due to its large size, the Pacific accounts for the largest DIC loss, approximately 57 % of the global decrease. This is largely associated with better ventilation of the interior Pacific water mass. However, the largest simulated DIC differences per unit volume are found in the SSWs of the Atlantic. Our study shows that the deep-water geometry and ventilation in the South Atlantic are altered between the two climate states where warmer climatic conditions cause SSWs to retreat southward and NSWs to extent further south. This process is mainly responsible for the simulated DIC reduction by restricting the extent of DIC-rich SSW, thereby reducing the storage of biological remineralized carbon at depth.

1 Introduction

The Last Interglacial (LIG, or Eemian) is composed of a warm onset around 125 ka characterized by warmer temperature at high latitudes relative to the present and a progressive cooling toward 115 ka when the last glaciation initiated (Masson-Delmotte et al., 2010; Otto-Bliesner et al., 2006). Evidence from land, ice, and ocean records identify the former as the period with the most intense global warming during the last 200 000 years (Capron et al., 2014; Dorthe Dahl-Jensen et al., 2013; Turney and Jones, 2010) mainly due to changes in the orbital configurations. If anthropogenic greenhouse gas emissions continue unabated, a climatically anomalously warm state is expected to occur in the near future with a warming by the end of this century that may be equivalent to the high-latitude reconstructed temperature for the LIG (Otto-Bliesner et al., 2013). The changes in the warm Eemian period may therefore be considered an analog for a future warmer climate.

Few model-based studies examine the carbon cycle dynamics for the LIG period with a particular focus on the ability of models to simulate the transient changes in atmospheric CO₂ concentration, which remains relatively stable around 270–280 ppm without displaying any trends (Lourantou et al., 2010; Schneider et al., 2013). Most of these studies provide a better understanding of the land carbon budget, particularly highlighting the importance of temperature changes on the land vegetation and slow processes of CO₂ change such as peatland carbon dynamics and CaCO₃ shallow-water accumulation (Schurgers et al., 2006; Kleinen et al., 2016; Brovkin et al., 2016). Although there are numerous studies that have analyzed the role of the ocean carbon cycle in regulating the atmospheric CO₂, especially for the interglacial-glacial transition period (Menviel et al., 2012; Ridgwell, 2001; Sigman and Boyle, 2000), to the authors' knowledge there is no study that investigate in details changes in marine carbon and nutrient cycling during the Eemian period of the LIG (125–115 ka).

With respect to changes in large-scale ocean circulation, reconstructions indicate that deep Atlantic circulation patterns and water mass geometries likely change over this interval. While mid-depth ventilation of northern-sourced waters (NSWs) is maintained in the North Atlantic by millennial-scale sustained North Atlantic Deep Water formation throughout the LIG (McManus et al., 2002; Mokeddem et al., 2014), southern-sourced waters (SSWs) expanded at depth toward 115 ka (Govin et al., 2009). In addition to large-scale circulation changes, temperature-induced changes in carbon solubility pumps and biological production are expected to alter the ocean carbon budget, in particular in the interior ocean. Other changes such as sea ice extent and ocean ventilation could also affect ocean carbon sequestration rate during the LIG. Elucidating the mechanisms responsible for changes in the ocean carbon distribution and inventory is of interest as it provides past constraints and context for evaluating the response of natural carbon sources and sinks to future climate change. This study aims to fill this knowledge gap by analyzing and comparing, in terms of ocean carbon dynamic, two opposite states of the LIG: the early and warm Eemian onset (125 ka) versus the cooler and late Eemian (115 ka). Using a state-of-the-art Earth system model, our study addresses the regional differences in the ocean carbon storage and the underlying mechanisms.

The paper is organized as follows: in Sect. 2, we describe the model, the experiment design, and the terms and metrics used to quantify the differences in carbon dynamics during the two periods. Section 3 presents the results of the model simulations, while discussions and comparison with previous studies are presented in Sect. 4. Finally, the study is summarized in Sect. 5.

2 Method

2.1 Model description

The present study uses output of an updated version of the Norwegian Earth System Model (NorESM1-ME), which has recently been further developed to efficiently perform multi-millennial and ensemble simulations (Bentsen et al., 2013; Guo et al., 2018; Luo et al., 2018). This model includes an isopycnal-coordinate ocean general circulation model based on the Miami Isopycnic Coordinate Ocean Model (MICOM; Bleck et al., 1992) and a biogeochemical ocean module adapted from the Hamburg Oceanic Carbon Cycle (HAMOCC5) model (Maier-Reimer, 1993; Maier-Reimer et al., 2005; Tjiputra et al., 2013). The inorganic seawater carbon chemistry in HAMOCC5 includes prognostic partial pressure of CO₂ (pCO₂) according to the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) protocols. The pCO₂ is computed as a function of temperature, salinity, dissolved inorganic carbon (DIC), total alkalinity (TALK), and pressure. This adapted version of HAMOCC5 does not include prognostic weathering fluxes but does employ a 12-layer sediment model following Heinze et al. (1999), which is particularly relevant for long-term transient simulations. The horizontal resolution of the land and atmospheric components is approximately 2^∘, while the ocean and ice components have higher resolutions of approximately 1^∘. In the vertical, the ocean model adopts 53 isopycnal layers.

The land component in NorESM (CLM4, Community Land Model version 4) is based on version 4 of the CLM family (Lawrence et al., 2012a). The land surface is sub-gridded into three sub-gridded entities: land units, columns, and plant functional types (PFTs). These sub-gridded cells are used to represent large-scale patterns of the landscape, variability in the soil and snow state variables, and the exchanges between land surface and atmosphere, respectively. Each of the sub-grid entities has its own prognostic variables, is independent, and experiences the same atmospheric forcing. Each cell is averaged and weighted with its fractional area.

The marine ecosystem is based on a nutrient–phytoplankton–zooplankton–detritus (NPZD) model that includes dissolved organic carbon (DOC). The inorganic nutrients consist of three macronutrients (phosphate, nitrate, and silicate) and one micronutrient (dissolved iron). A constant Redfield ratio is adopted in the model as P : C : N : ΔO₂ = 1 : 122 : 16 : −172. The phytoplankton growth rate is expressed as a function of temperature, light (Eppley, 1972; Smith, 1936), phosphate, nitrate, and dissolved iron availability, and its loss is regulated by an exudation and mortality rate, in addition to zooplankton grazing. The penetration of light decreases with depth following an exponential function, which responds to a gradual extinction factor formulated as a function of water depth and chlorophyll concentration (Maier-Reimer et al., 2005). The model prescribes a global constant vertical sinking speed of particles produced in the euphotic zone (above 100 m depth). The particulate organic carbon (POC), which comprises dead phytoplankton and zooplankton, sinks through the water column at a speed of 5 m day⁻¹ and is remineralized at a constant rate of 0.02 day⁻¹ when oxygen is available. Other particles such as opal shells and particulate inorganic carbon (PIC) sink at a speed of 60 and 30 m day⁻¹, respectively. Particulates that reach the seafloor without being remineralized interact chemically with the sediment pore water via bioturbation and vertical mixing and advection within the sediment layers. In the model, the air–sea gas exchange of CO₂ and O₂ only occurs between the ocean surface and the atmosphere in ice-free regions and is computed according to the following three components (Wanninkhof, 1992): the gas solubility in seawater, which is computed as a function of surface salinity and temperature according to Weiss (1970, 1974); the gas transfer velocity, which is proportional to the square of the surface wind speed and is computed as a function of the Schmidt number; and finally the air–sea gradient of gas partial pressures. To better elucidate various biogeochemical processes on the carbon cycle, the model is updated to also include preformed O₂, TALK, and PO₄ tracers in the biogeochemical module. At the surface, these preformed tracers are set to their non-preformed value and are advected passively by the ocean circulation in the interior without any other sources and sinks. Finally, in order to provide information of the water mass ages since its last contact with the atmosphere, an idealized age tracer is implemented and simulated in the NorESM model. Hence, the age tracer is set to zero for all water masses at the ocean surface and subsequently transported and mixed passively with circulation in the ocean interior and integrated with the model time step. This tracer is also used to estimate the ventilation rate of different interior water masses.

2.2 Experiment setup

Two equilibrium experiments are performed over the Eemian (model data available at https://doi.org/10.11582/2018.00038; Kessler, 2018), one near the onset (warmer than today, 125 ka) and one at the end (colder, 115 ka) of the Last Interglacial. Both experimental configurations follow the standard protocols of the third phase of the Paleoclimate Modelling Intercomparison Project (PMIP3; https://pmip3.lsce.ipsl.fr/, last access: 5 December 2018), with a fixed vegetation coverage from the preindustrial boundary conditions. The only differences from the preindustrial configurations are the orbital parameters and the greenhouse gases concentrations (CO₂, CH₄, N₂O). For the experiment at 125 ka (115 ka), the atmospheric CO₂, CH₄, and N₂O levels are prescribed to be 276 ppmv (273 ppmv), 640 ppb (472 ppb), and 263 ppb (251 ppb), respectively. The two experiments are branched off from 1000 years of spin-up with a preindustrial setup and forced with their respective interglacial boundary conditions for 4000 simulation years.

In the last 50 years of Eemian forcing simulations the ocean is close to equilibrium. Only small drifts remain, mainly in the Pacific basin where the equilibrium is still not fully established. Therefore, the global ocean DIC and TALK slightly decrease in the experiment at 125 ka (115 ka) by approximately −0.15 PgC yr⁻¹ (−0.06 PgC yr⁻¹) and −0.01 Pmol yr⁻¹ (−0.01 Pmol yr⁻¹), respectively. However, these drifts are small compared to the absolute ocean budget in DIC (37391 and 37 705 PgC) and TALK (3291 and 3303 Pmol) for the 125 and 115 ka experiments, respectively. The differences in the pCO₂ and TALK budget are small between the two experiments. Such changes would affect the DIC budget by about 32 PgC. The CO₂ flux is relatively constant and depicts the ocean as a weak source to the atmosphere with an outgassing of 0.12±0.06 at 115 ka and 0.15±0.06 PgC yr⁻¹ at 125 ka. In addition, the difference in sedimentation rates between the two experiments ($\sim {\mathrm{6.10}}^{-\mathrm{4}}$ PgC) appears to be negligible compared to the difference in DIC budget.

2.3 DIC decomposition

In order to analyze the oceanic carbon cycle, dissolved inorganic carbon (DIC) is decomposed into its preformed and biological components (Eq. 1), following Bernadello et al. (2014):

$$\begin{array}{}\text{(1)}& {\mathrm{DIC}}^{\mathrm{tot}}={\mathrm{DIC}}^{\mathrm{pre}}+{\mathrm{DIC}}^{\mathrm{bio}}.\end{array}$$

The preformed component of DIC (DIC^pre) comprises saturated and disequilibrium parts (Eq. 2).

$$\begin{array}{}\text{(2)}& {\mathrm{DIC}}^{\mathrm{pre}}={\mathrm{DIC}}^{\mathrm{sat}}+{\mathrm{DIC}}^{\mathrm{dis}},\end{array}$$

where the DIC at saturation (DIC^sat) describes the DIC concentration when the water parcel is in full equilibrium with the atmospheric CO₂ when it is last in contact with the surface. In our model, changes in DIC^sat are mainly affected by changes in preformed TALK. We compute this variable offline with the inorganic carbon chemistry program CO2SYS developed in Matlab (van Heuven et al., 2011), including other parameters such as the model output of preformed alkalinity (TALK^pre), preformed phosphate (PO₄^pre), surface silicate, salinity, and temperature. In addition, the atmospheric CO₂ concentration from each experiment is used. To complete the CO2SYS input, we applied the dissociation constants K1 and K2 introduced by Mehrbach et al. (1973) and refitted by Dickson and Millero (1987). The disequilibrium part of DIC (DIC^dis) measures the disequilibrium state of the surface water with respect to the atmosphere. This parameter is computed from the DIC^tot (output), from which the biological and saturated DIC components (DIC^bio and DIC^sat) are subtracted. Therefore, all components are included in its calculation. A negative DIC^dis occurs when the water parcel sinks into the ocean interior before a full equilibration with the atmosphere is obtained, which leads to an undersaturation of the water parcel. This undersaturation can also be reinforced by biological CO₂ consumption at the surface, which tends to increase the timescale needed for the water parcel to equilibrate. However, a positive DIC^dis translates into a supersaturation. The latter occurs when deep waters, which contain high concentration of DIC because of remineralization processes, upwell or mix vertically with the surface waters (Follows and Williams, 2004). Both DIC^dis and DIC^sat are transported by ocean circulation into the interior ocean.

The biological component of DIC comprises (1) the interior remineralization of organic matter (expressed in carbon), which is produced in the euphotic layer via photosynthesis (also referred to as soft-tissue pump), and (2) the remineralization of planktonic calcium carbonate shells (expressed in carbon; calcium carbonate pump). These two remineralization components are added to form the biological component of DIC, as shown in Eq. (3).

$$\begin{array}{}\text{(3)}& {\mathrm{DIC}}^{\mathrm{bio}}={\mathrm{DIC}}^{\mathrm{soft}}+{\mathrm{DIC}}^{\mathrm{carb}}\end{array}$$

The remineralization of soft tissues (hereafter DIC^soft) contributes via phosphate (PO₄) remineralization through a carbon phosphorus stoichiometric ratio $r_{\mathrm{C}:\mathrm{P}}=\mathrm{122}$. This component is calculated from the difference between the total and the preformed PO₄ according to

$$\begin{array}{}\text{(4)}& {\mathrm{DIC}}^{\mathrm{soft}}=r_{\mathrm{C}:\mathrm{P}}({{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{tot}}-{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{pre}}).\end{array}$$

The carbonate pump contributes through the dissolution of CaCO₃ hard shells, calculated as difference between the total and the preformed alkalinity and PO₄ following

$$\begin{array}{}\text{(5)}& {\mathrm{DIC}}^{\mathrm{carb}}=\mathrm{0.5}[{\mathrm{TALK}}^{\mathrm{tot}}-{\mathrm{TALK}}^{\mathrm{pre}}+r_{\mathrm{N}:\mathrm{P}}({{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{tot}}-{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{pre}}\left)\right],\end{array}$$

where $r_{\mathrm{N}:\mathrm{P}}=\mathrm{16}$ is the Redfield ratio adopted by the model and the phosphate term accounts for the alkalinity changes owing to the soft-tissue pump.

In our analysis, we show differences between the warmer 125 ka and the colder 115 ka experiments. We therefore use the delta notations ΔDIC^tot, ΔDIC^sat, ΔDIC^soft, ΔDIC^carb, and ΔDIC^dis to refer to changes between the warmer and the colder periods.

2.4 Water mass analysis

In order to identify water mass sources, we apply the “PO” tracer as defined by Broecker (1974). It is computed using phosphate and oxygen fields following

$$\begin{array}{}\text{(6)}& \mathrm{PO}={\mathrm{O}}_{\mathrm{2}}+r_{\mathrm{O}:\mathrm{P}}{\mathrm{PO}}_{\mathrm{4}},\end{array}$$

where $r_{\mathrm{O}:\mathrm{P}}=\mathrm{172}$ is the phosphorus-to-oxygen stoichiometric ratio used in the model. This tracer is presumed to be nearly constant for a specific water mass. It is based on the principle that phosphate is released while oxygen is used during remineralization, and vice versa during biological production. The distinction of water masses using PO is useful for contrasting water masses with very different surface PO values. Here, we mainly use PO to identify NSW and SSW masses in the deep ocean below 1000 m depth characterized by low and high PO values, respectively.

Figure 1 Difference in sea surface temperature (ΔSST) between 125 and 115 ka. Only significant differences (i.e., with absolute value greater than the interannual standard deviation over the last 50 years in both 125 and 115 ka) are shown. The green and purple lines correspond to 50 % sea ice extent at 115 and 125 ka, respectively. The black (blue) thick lines depict regions affected by a deepening (shallowing) of the mixed-layer depth (with difference greater than 100 m depth) at 125 ka compared to 115 ka.

3 Results

Section 3.1 describes the sea surface temperature and sea ice changes, while changes in water mass properties are discussed in Sect. 3.2. The near-surface changes that particularly influence the biological pump are addressed in Sect. 3.3, and global to regional oceanic DIC storage changes are presented in Sect. 3.4. The analysis is performed over the average of the last 50 years of the simulations. In addition, the global ocean is divided into three main basins: Atlantic, Indian, and Pacific.

3.1 Sea surface temperature and sea ice

Our model simulates a global and annual increase of sea surface temperature (SST; +0.27 ^∘C) at 125 ka relative to the 115 ka experiment. This warming is mainly simulated at high latitudes (Fig. 1a–d) where higher SSTs are simulated throughout the year in the Southern Ocean (SO), south of Greenland, the Norwegian Sea, and the northern part of the Pacific Ocean. This persistent warming at 125 ka induces the sea ice to melt (Fig. 1, green and purple lines). Consequently, in the areas where sea ice extent is reduced, there is more air–sea gas exchange and an increase in mixed-layer depth (MLD; Fig. 1, black lines in the Southern Ocean). In the Labrador Sea the mixed layer depth at 125 ka is deeper by more than 100 m than at 115 ka. This is due to higher salinity simulated in this region at 125 ka. At lower latitudes, the SSTs vary more spatially and seasonally. For example, in some parts of the Atlantic Ocean (North Atlantic Drift, the equatorial region, and some sections of the subantarctic near the 45^∘ S band) cooler SSTs are simulated over several seasons (ΔSST<0, Fig. 1a–b). While the North Atlantic drift depicts cooler SSTs during boreal winter and spring (Fig. 1a–b), colder SSTs last until the boreal summer in the equatorial region of the Atlantic Ocean (Fig. 1a–c). The 45^∘ S latitude band remains cooler over the four seasons. Here, the MLD seems to be more controlled by changes in salinity instead of temperature distribution.

We note that relative to the preindustrial, based on proxy (Hoffman et al., 2017), our model simulates consistent spatial feature of annual SST anomalies at 125 ka. It simulates, among other things, (i) strongest warming at high latitudes, specifically in parts of the Southern Ocean; (ii) weak cooling at low latitudes; (iii) cooler SST in most of the Indian Ocean; and (iv) a warmer eastern North Pacific. Nevertheless, the amplitude of SST warming and cooling at specific sites tends to be weaker in our model. This feature appears to be common in other global models and could be attributed to their low spatial resolutions (Hoffman et al., 2017).

Figure 2 Southern Ocean overturning stream function in Sverdrup (Sv) for the experiment at 115 ka (a) and 125 ka (b). Red colors represent water masses moving clockwise, while blue colors represent water masses moving counterclockwise.

3.2 Water mass properties

In order to analyze the water mass properties, it is useful to examine the changes in the overturning circulation. Figure 2 shows the global overturning stream function in the Southern Ocean for both experiments. The Antarctic circumpolar current is simulated stronger and deeper at 125 ka compared to 115 ka (Fig. 2). This strengthening mainly occurs in the Pacific section of the Southern Ocean (near 50^∘ S), suggesting an increase of the ventilation rate of the intermediate waters formed in this region. The Atlantic section of the Southern Ocean remains weakly modified. Indeed, using the same model simulations as the present study, Luo et al. (2018, Supplement Fig. S8) show that the surface wind speed in the eastern and western South Atlantic are relatively similar at 125 and 115 ka. In addition, they also show that the simulated Atlantic Meridional Overturning Circulation (AMOC) at 125 ka is as vigorous as at 115 ka but deepen by about 300 m depth. This suggests that the mid-depth and bottom water in the North Atlantic Ocean should be better ventilated at 125 ka than at 115 ka.

These changes in the overturning circulation affect the water mass age tracer. Analyzing this parameter allows us to examine the interior ocean ventilation rate. A reduction in water mass age translates to a faster ventilation rate and vice versa for an increase in age. The differences in water mass age between 125 and 115 ka are presented in Fig. 3, depicting the zonally averaged sections for each ocean basin. The water mass ages in the Atlantic and the Indian oceans show similar patterns, with mean older water masses in the upper layers at 125 ka (roughly +100 years) and younger water masses below 1000 m depth (as much as 500 years younger). This is consistent with a deeper and slightly stronger AMOC at 125 ka as shown in Luo et al. (2018). The Southern Ocean (south of 50^∘ S) contains younger water masses throughout the entire water column in both basins at 125 ka, translating to a stronger ventilation rate. This is likely due to changes in the distribution of SSW and NSW. Figure 4 show that there is a clear distinction between interior water mass structure in the Atlantic basin between 115 ka (Fig. 4a) and 125 ka (Fig. 4c). It shows that below 2000 m depth the SSW retreats further southward at 125 ka relative to 115 ka in the Atlantic basin. This confinement in SSW is induced by the change in the Antarctic sea ice cover (Fig. 1), as explained by Ferrari et al. (2014). However, the model also simulates a modification in SSW density (−0.2 kg m⁻³ at 125 ka compared to 115 ka; Fig. 4a, c). This reduction in water density is mainly driven by the input of low-salinity freshwater from the melting of the Antarctic sea ice and may have an additional impact on determining the Atlantic distributions of NSW and SSW. As a result of sea ice retreat, the winter mixed layer depth in this area increases, resulting in younger water mass age in the Southern Ocean. Near the surface (∼ 800 m depth) the SSW seems to enter further north in the Atlantic basin at 125 ka (Fig. 4c) relative to 115 ka (Fig. 4a).

Figure 3 Zonally averaged section of the difference in water mass age (ΔAge) between 125 and 115 ka in (a) the Atlantic Ocean, (b) the Indian Ocean, and (c) the Pacific Ocean. The dashed lines display ΔAge=0.

While such large redistributions of northern- and southern-origin deep waters only occur in the Atlantic, these changes also influence water properties in the Indian Ocean due to the “downstream” advection of younger deep water into the interior during the warmest period (125 ka – Fig. 3b). In addition to simple advection of younger water northward in the Indian basin, the residence time of the Indian Ocean's deep water must also decrease since the ventilation ages decrease northward at depth. By contrast, in the Pacific Ocean, the zonally averaged bottom-water mass ages are simulated to be older at 125 ka (Fig. 3c). However, this basin can be divided between the western and eastern side. While the western side is also influenced by the younger water masses created in the Atlantic Ocean, the eastern side waters of the basin are simulated to be older at 125 ka by as much as 300 years. These older water masses are created in the Pacific SO and may be affected by a flattening of the isopycnals south of 60^∘ S at 125 ka (Fig. 4d, gray lines). This flattening of the isopycnals is influenced by both sea ice melting and higher SSTs, and suggests stronger stratification and a weaker subduction rate. In addition, the SSW boundary (Fig. 4b, d; light purple dashed lines) seems to be slightly poleward at 125 ka ($\sim \mathrm{60}{}^{\circ }$ S) compared to 115 ka ($\sim {\mathrm{58}}^{\circ }$ S). This may suggest that more surface water originating from the subtropical gyre (more depleted in phosphate) enters the Pacific section of the Southern Ocean. Finally, the Pacific intermediate waters are simulated as younger at 125 ka when compared to 115 ka (Fig. 3c), which is consistent with the strengthening of the upper cell of the overturning circulation previously mentioned for this region.

Figure 4 Zonally averaged section of PO as defined by Broecker (1974) in (a, c) the Atlantic Ocean at 115 and 125 ka, respectively, and (b, d) the Pacific Ocean at 115 and 125 ka, respectively. The light purple dashed lines delimit the water influenced by the SSW from water influenced by the NSW. The dark gray solid lines depict the neutral density (units: kg m⁻³).

Table 1 Difference in southern-sourced water (Δ_SSW) in the global ocean and for each basin. Row 1 shows the volume (V_SSW) according to our PO ≥0.57 mol O m⁻³ criteria. Rows 2 and 3 summarizes the DIC and water mass age mean value for the two period of study. The changes relative to 115 ka are given as a percentage in parenthesis.

The southern-sourced waters are particularly affected in terms of geometry distribution. We therefore divide the changes from these SSWs into the three basins. Table 1 summarizes the changes occurring below 1000 m depth in the SSWs in terms of volume, DIC, and water mass age for each basin and reveals the Atlantic as the most affected area under warmer climate conditions. The relative difference in all those three characteristics (ΔV_SSW, ΔAge_SSW, and ΔDIC_SSW) between 125 and 115 ka is the greatest in the Atlantic (−37 %, −262 years, and −0.92 g C m⁻³). This demonstrates that the ventilation mechanism in the Atlantic sector of the SO is likely to be more sensitive (than in other basins) to climate change.

In response to the different forcings between 125 and 115 ka, significant changes are simulated in deep-water ventilation rates and water mass distribution in the three basins. While the responses in the Atlantic and Indian oceans have some similarity (better deep-water ventilation), the Atlantic basin seems to be the most sensitive and also simulates water mass geometry changes. However, the ventilation rate in the Pacific Ocean depicts an opposite sign of change to that of the Atlantic and Indian oceans. In the next section we discuss how these modification impact the near-surface productivity in our model.

3.3 Near-surface productivity

Figure 5 shows the differences in carbon export production (EPC) and phosphate (PO₄) at the surface of the ocean between 125 and 115 ka. When comparing the 125 and 115 ka experiments, our model simulates two major features characterized by (i) an increase in EPC in the equatorial region of the Atlantic Ocean (Fig. 5a, turquoise rectangle) and (ii) a reduction in EPC of similar magnitude in the equatorial region of the Pacific Ocean (Fig. 5a, purple rectangle). Although there are changes in surface wind speeds and vertical mass fluxes (not shown here) in the equatorial regions, they do not consistently explain the simulated changes in export production (i.e., stronger wind and upwelling in the eastern equatorial Pacific at 125 ka compared to 115 ka, where a decrease in export production is simulated; in parts of the equatorial Atlantic some increase in export production may be related to the simulated stronger upwelling at 125 ka). However, nutrient changes in the Atlantic and Pacific sections of the Southern Ocean are consistent with the changes in export production, and therefore we suggest the following mechanism. In the Atlantic Ocean, the subantarctic water (45^∘ S latitude band) corresponds to the southern-sourced intermediate water formation in the model. This water mass sinks and reemerges along the Equator (Fig. 5, turquoise rectangles). This pattern expected from models and modern observations shows that the intermediate and mode waters formed at high southern latitudes feed the subtropical thermocline and act as a predominant source of nutrients important for sustaining low-latitude biological productivity (Gu and Philander, 1997; Sarmiento et al., 2004b). We acknowledge that there is uncertainty in the complex pathways of the subantarctic water masses toward the Equator simulated in the model. Nevertheless, due to higher preformed and remineralized phosphate in the subantarctic water sinking region, more PO₄ is advected through this “ocean tunnel” to the Equator and therefore leads to an increase in EPC (Fig. 5a, turquoise rectangle). A similar ocean tunnel (Fig. 5, purple rectangle) connects the high and low latitudes in the Pacific Ocean but results in the opposite sign of change. Here, as pointed out in Sect. 3.2, more subtropical waters seem to enter the Pacific Southern Ocean at 125 ka than at 115 ka. These waters are depleted in phosphate compared to southern-sourced waters. As a result, less phosphate is subducted toward the equatorial Pacific (Fig. 5b, purple rectangle), leading to a reduction in the EPC near the equatorial upwelling (Fig. 5a, purple rectangle).

Figure 5 Difference in (a) export production of carbon at 100 m depth (ΔEPC) and (b) phosphate at the surface (ΔPO₄) between 125 and 115 ka. When the absolute value of the difference is below the standard deviation over the last 50 years at 125 and 115 ka, the value returns a NaN. Purple and turquoise rectangles highlight the two “ocean tunnels” linking the sinking/upwelling regions of southern-sourced intermediate waters in our model.

There are no significant changes in the biological activity in the Indian Ocean or the remaining Southern Ocean. Only a weaker carbon export at 125 ka relative to 115 ka is simulated around 40^∘ S and the Arabian Sea.

Thus, the simulations reveal that changes in Southern Hemisphere thermocline ventilation regions modulate basin scale productivity and export production even within an interglacial period with modest changes in external forcing. This result is broadly consistent with previous studies suggesting that the upper limb of the biogeochemical divide is critical for setting biological export production and is sensitive to climate changes (Marinov et al., 2006; Moore et al., 2018; Sarmiento et al., 2004a). Despite zonally homogeneous forcing we find a basinally heterogeneous response both in subantarctic ventilation and in low-latitude productivity which is similar to, albeit more extreme than, the basin specific response simulated for future warming and stratification (Moore et al., 2018).

These changes in surface physical and biological activities could also have implications for the exchanges of carbon between near-surface and interior water masses, and therefore the interior carbon budget. In the next section, ventilation changes at 125 ka are compared to 115 ka by analyzing the simulated water mass properties.

3.4 Global and regional carbon budgets

Figure 6 shows the difference in the carbon inventory vertical profiles between 125 and 115 ka as simulated by our model. The changes in DIC^tot, DIC^sat, DIC^soft, DIC^carb, and DIC^dis are averaged over a 500 m depth interval. The global amount of DIC^tot is −314.1 PgC in the ocean under the warmer conditions (Fig. 6a, gray ΔDIC^tot). Since the atmospheric CO₂ is kept constant in each experiment, this carbon loss at 125 ka compared to 115 ka is implicitly assumed to be balanced out by the changes in the land carbon reservoir. Here, the Atlantic Ocean accounts for 15 % (−49.1 PgC) of that global decrease, while the Indian and the Pacific basins contribute 28 % (−87.2 PgC) and 57 % (−179.0 PgC), respectively. Only in the near-surface layers does the model simulate a positive ΔDIC^tot, which translates to higher surface DIC concentration at 125 ka relative to 115 ka. Most of the ocean interior has lower DIC concentration at 125 ka, with the strongest difference in ΔDIC^tot simulated between 2000 and 3000 m depths for each basin, except for the Atlantic (Fig. 6b–d). The soft-tissue pump and the disequilibrium effect are the main contributors for the reduced carbon inventory depicted at 125 ka at the global scale (Fig. 6a, green and purple bars) – each accounting for a third of the total −314.1 PgC decrease. Similarly, ΔDIC^tot in the Atlantic basin is also predominantly controlled by the biological pump, i.e., the soft-tissue plus carbonate pump throughout the entire water column with a decrease at depth and an increase in the near-surface (Fig. 6b). Except for the upper ocean, the contribution from the saturation component related to temperature and salinity change is generally negligible.

Figure 6 DIC differences between 125 and 115 ka (ΔDIC^x) in (a) the global ocean and the (b) Atlantic, (c) Indian, and (d) Pacific basins. The ΔDIC^x is averaged over a 500 m depth interval where “x” refers to the different components of the DIC. The DIC^tot is represented by the gray bars and is decomposed into its four components: ΔDIC^sat (blue), ΔDIC^soft (green), ΔDIC^carb (orange), and ΔDIC^dis (purple). The sum throughout the water column of each components is given by the legend.

The ΔDIC^tot of the Indian basin also depicts a strong reduction at 125 ka relative to 115 ka. Here the soft-tissue and saturation components simulate the strongest decrease (Fig. 6c, green and blue bars). In addition, near-surface changes in ΔDIC^tot are controlled by the changes in ΔDIC^sat. Similarly, the near-surface layer of the Pacific Ocean is also controlled by the changes in the saturation component (Fig. 6d), simulating a positive difference of about +18 PgC. Changes in the deeper layers are mainly attributed to the disequilibrium effect and the soft-tissue pump, accounting for a decrease of −83.6 and −44.0 PgC at 125 ka relative to 115 ka, respectively. However, the saturation component also has a considerable influence on the carbon storage with persistent negative ΔDIC^sat throughout the water column.

Figure 7 Atlantic zonally averaged section of the difference in (a) ΔDIC^tot, (b) ΔDIC^soft, (c) ΔDIC^carb, (d) ΔDIC^dis, and (e) ΔDIC^sat between 125 and 115 ka. The black dashed lines represent the zero values.

Figure 8 Same as Fig. 7 but for the Indian Ocean.

In order to address the regional changes, we analyze the differences in each DIC components further by calculating the zonally averaged values in each basin. Figures 7, 8, and 9 depict these differences for the Atlantic, Indian, and Pacific basins, respectively. As shown in Fig. 6b, the carbon inventory of the Atlantic is reduced mainly below 1500 m depth. Here the Southern Hemisphere is the most affected region, which depicts the strongest differences in ΔDIC^tot (Fig. 7a, blue shading). This pattern corresponds well to the changes in soft-tissue pump (Fig. 7b). Near the surface, higher carbon export (Fig. 5a) increases the remineralization of organic matter, leading to higher DIC concentration at 125 ka. At depth, the slightly deeper AMOC at 125 ka (compared to 115 ka), leads to better-ventilated mid-depth to bottom waters in the Northern Hemisphere, leading to less remineralized organic matter by reducing the water mass age. This is reflected by the negative ΔDIC^soft and ΔDIC^carb, translating to a less effective soft-tissue and carbonate pump at 125 ka. Positive changes in ΔDIC^tot in near-surface waters and bottom water at 20^∘ N also arise from the soft-tissue and carbonate signal due to the increase of the alkalinity (not shown here) and slightly older water masses along the African coast. The bottom waters in the Southern Hemisphere are mainly controlled by a stronger disequilibrium effect, i.e., negative change in comparison to 115 ka. This change in disequilibrium is due to the sea-ice-induced retreat of the SSW and the inflow of more NSW between 50^∘ S and the Equator at 125 ka. The NSW mass, formed in the North Atlantic, is generally more subject to biological production during its near-surface northward transport before sinking into the interior than the SSW (Duteil et al., 2012). The biological production consumes DIC during photosynthesis and pushes the water mass further out of the equilibrium with the atmospheric CO₂, inducing ΔDIC^dis to be more negative. The negative values of ΔDIC^dis are conserved when the water parcel flows southward into the deep ocean. For this reason, the regions that are no longer influenced by SSW at 125 ka depict a negative ΔDIC^dis (Fig. 7d). However, the upper layers of the Atlantic Ocean are mostly simulated with higher DIC^dis (positive ΔDIC^dis, or less disequilibrium at 125 ka). This could be induced by more SSW (less disequilibrated than NSW) entering the Atlantic basin further north near the surface as Fig. 4c suggests. Finally, the loss of carbon in the Southern Ocean is shown to be mainly attributable to a decrease of the saturation component at 125 ka. This decrease is mainly attributed to lower TALK^pre (not shown here), possibly provoked by the melting of the sea ice.

Figure 9 Same as Fig. 7 but for the Pacific Ocean.

The DIC storage in the Indian Ocean generally shows a decrease at 125 ka, with the strongest changes occurring at depth north of 30^∘ S (Fig. 8a, dark blue shading). Positive ΔDIC^tot are nevertheless simulated in the region that may correspond to the Antarctic Intermediate Water (AAIW) (Fig. 8a, red shading). Similar to the Atlantic basin, the pattern of the soft-tissue pump changes corresponds to the ΔDIC^tot pattern throughout most of the Indian Ocean. This decrease in biological remineralization is in agreement with the water mass age changes seen in Sect. 3.2 (Fig. 3b): younger water masses account for less biologically induced DIC content. However, the bottom and the surface waters show opposite signs in the ΔDIC^tot, which suggests that other processes are acting in these regions. The differences in the carbonate pump remain small and roughly follow the pattern of the soft-tissue pump (Fig. 8c). Changes in the bottom-water ΔDIC^tot can be attributed primarily to the difference in the disequilibrium effect, which is probably affected by larger influence of NSW at 125 ka relative to 115 ka as described for the Atlantic basin, and to a slight reduction in the saturation component. In addition, the disequilibrium component might also be influenced by stronger carbon export in the Southern Ocean (as the positive ΔDIC^soft and Fig. 5a suggest). The strong positive ΔDIC^dis simulated near the surface along the Indian coast is in good agreement with lower SSTs (Fig. 1a–c) allowing more DIC to be absorbed at 125 ka and a strong reduction in carbon export production (Fig. 5a). Finally the negative ΔDIC^sat depicted in the top layers in the north of the Indian Ocean is mainly attributed to a change in water mass origin from 115 to 125 ka. During 115 ka the SSW upwells from the deep ocean into the Arabian Sea. By contrast, at 125 ka, the waters coming from the Indonesian region mix with SSW. These Indonesian throughflow water masses initially coming from the Pacific Ocean are affected by strong precipitation in the Indonesian basin, which reduces the TALK at the surface and therefore DIC^sat.

The Pacific basin shows the strongest DIC^tot decrease in the Northern Hemisphere, mainly due to the reduction in soft-tissue pump (Fig. 9b). This decrease in biogenic carbon is induced by a better-ventilated water mass around 30^∘ N (Fig. 3c), which may come from the increase of the overturning circulation in the upper cell of the Pacific Ocean. In contrast, the DIC inventory of the high-latitude Southern Hemisphere bottom and near-surface waters is larger at 125 ka relative to 115 ka. This is attributed to the changes in soft-tissue pump, which is more effective at 125 ka due to longer residence time of the water masses (Fig. 3c). The carbonate pump has a relatively low impact on the total ΔDIC inventory in the Pacific basin but follows the same pattern as the changes in the organic carbon remineralization. The disequilibrium effect accounts for the strongest decrease of DIC throughout the basin as seen by the negative ΔDIC^dis in almost all regions (Fig. 9d). However, the bottom waters seem to be the most affected and become more undersaturated at 125 ka relative to 115 ka. This may be influenced by the slowing down of the subduction process in the Southern Ocean, induced by the flattening of the isopycnals. The possible higher carbon export production in the SO south of 60^∘ S (Fig. 5a) may also push the water mass further out of equilibrium. Finally, the saturation component is controlling the changes occurring in the near-surface waters (Fig. 9e). Lower saturations are attributed to higher export of calcium carbonate south of 60^∘ S, which lower the alkalinity at the surface and thereby the buffering capacity. On the other hand, the calcium carbonate formation decreases north of 60^∘ S, resulting in a higher alkalinity and buffering capacity.

4 Discussion

The ocean plays an important role in storing carbon and, thus, in the long term regulation of atmospheric CO₂ levels. The processes involved in regulating the ocean carbon inventory are likely to change under warmer future conditions. In this study, we simulate two equilibrium states of the penultimate interglacial period using a state-of-the-art Earth system model and make a first attempt at quantifying the biogeochemical and physical processes responsible for carbon storage changes caused by different (interglacial) orbital configurations and background climates. Significant decreases in the ocean carbon storage capacity are found in a warmer climate. More than 48 % of this decrease is induced by the reduction of the biological pump. This decrease is found to be mainly driven by the shorter residence time of interior deep water masses arising from changes in Southern Ocean sea ice extent that influence the NSW and SSW mass geometry, in addition to changes in overturning circulation in the Atlantic (deeper but almost equally vigorous) and Pacific (stronger upper cell) basins.

Using the available proxy reconstructions during the LIG period allows us to assess the validity of important features in our model results. We assess the validity of the simulated 115–125 ka water mass geometry change using LIG proxy reconstructions of bottom-water δ¹³C, a water mass tracer strongly but inversely related to carbon and PO₄ contents (Eide et al., 2017). Similar to our results, expanded SSW in the late compared to early LIG is inferred from such reconstructions to explain bottom-water δ¹³C decreases in different regions proximal to the Southern Ocean (Govin et al., 2009; Ninnemann and Charles, 2002). Bottom-water δ¹³C reconstructions indicate less influence of NSW in the deep South Atlantic at 115 ka than at 125 ka (Govin et al., 2009; Ninnemann et al., 1999). In addition, persistent (millennial-scale) mid-depth NSW ventilation extending from the LIG and into the subsequent glacial inception is also suggested based on proxy reconstructions (McManus et al., 2002; Mokeddem et al., 2014) and is consistent with model simulations (Born et al., 2011; Wang and Mysak, 2002). In our study, even though the AMOC is simulated slightly stronger and deeper at 125 ka, vigorous AMOC persists until 115 ka, ventilating the North Atlantic mid-depth. We also note that our model may not properly represent North Atlantic overflows due to its sparse resolution. This can further add uncertainties to North Atlantic water ventilation.

Also consistent with our results, ice core proxies indicate that Southern Ocean sea ice extent is greater at 115 ka than at 125 ka (Röthlisberger et al., 2008; Wolff et al., 2006), while our model reproduces the volumetric SSW expansion in response to this increase in Southern Ocean sea ice extent as suggested for glacial climates (e.g., Ferrari et al., 2014). Our model results suggest that similar sea ice (Fig. 1) and SSW expansions (Fig. 4a), albeit muted compared to glacial changes, occurred in response to LIG orbital configuration changes and without continental ice sheet growth (not included in the model), indicating a relatively tight coupling between Antarctic climate, sea ice, and the deep Atlantic water mass geometry changes influencing ocean carbon storage.

The changes in ocean carbon storage simulated by our model are significant and demonstrate that warm (interglacial) ocean carbon content changes with climate forcing. While atmospheric CO₂ is fixed in our model, preventing a direct assessment of ocean carbon changes in atmospheric CO₂, the decrease in deep carbon storage and shoaling of the DIC pool during the warm 125 ka interval is generally consistent with higher atmospheric CO₂ levels at this time. Our estimated changes in ocean carbon budget are in the range of previous modeling studies that also suggest weaker ocean carbon storage during the beginning of the LIG (125 ka) relative to the glacial inception (115 ka). Schurgers et al. (2006) obtain a difference in atmospheric and terrestrial carbon storage of about 40 and 350 PgC, respectively, between the onset and end of the LIG. This potentially translates to a weaker ocean carbon storage of 310 PgC at the onset compared to the late LIG, which corresponds well in absolute magnitude with our findings of 314.1 PgC decrease at 125 ka. However, their simulated atmospheric CO₂ steadily increases over this period, which potentially points toward more carbon needing to be stored in land or ocean toward the end of the LIG.

In contrast, Brovkin et al. (2016) simulate more ocean DIC^tot at 125 ka than 115 ka using simpler EMIC models. This increase is mainly attributed to the shallow-water carbonate precipitation implemented in their model. This process is not included in our model and could explain the differences in the results. However, their simulated change in atmospheric CO₂ after 121 ka is in the opposite direction (an increase of roughly 20 ppm) of the atmospheric trend (a decrease of roughly 3 ppm) observed in the ice core data. This suggests that either (or both) the terrestrial or the oceanic carbon reservoir does not take up enough carbon toward the end of the LIG in order to simulate a decrease in atmospheric CO₂ content.

Concerning the modification in the upper-ocean productivity under warmer climatic conditions, our model study shows a heterogeneous response in phosphate availability and carbon export production especially between the Atlantic and Pacific basins. Moore et al. (2018) also highlight such a biogeochemical divide response for future projections under warmer climatic conditions. This implies that future anthropogenic CO₂ forcings may have a similar impact on the biogeochemical divide to that of past forcings. Therefore, reconstructing and understanding the large-scale productivity responses to past climate forcing are critical for assessing both global and regional sensitivity of the ocean carbon dynamic to climate change.

There are limitations to our study. Factors that could influence ocean carbon storage including sea level, riverine input of nutrients, and atmospheric dust loading are all set to preindustrial levels in our simulations but may be different in the LIG. Global sea level, for example, may be as much as 6–9 m above present levels (Kopp et al., 2009). In addition, our model does not include weathering fluxes, which might influence the carbon budget on such a long timescale. Further, we compare two quasi-equilibrated states, which is unrealistic and ignores transient forcings and shorter-term variability. This may explain differences between our model results and some proxy reconstructions. For example, proxy reconstructions suggest that both NSW and SSW ventilation may have varied considerable near 125 ka (Galaasen et al., 2014; Hayes et al., 2014). The changes suggested by these studies include reductions of NSW and expansions of SSW similar to our modeled 115–125 ka equilibrium difference, albeit occurring as short-lived (centennial-scale) transient events associated with freshwater input episodes during the final phase of northern deglaciation (Galaasen et al., 2014). Our quasi-equilibrated model simulations for 115 and 125 ka also lack ice sheet and the corresponding freshwater input variability, nor do they address such shorter-term changes that could affect the ocean carbon inventory (Stocker and Schmittner, 1997). However, short-lived changes would likely have less impact on the ocean carbon inventory than the longer-term (millennial-scale) changes we address, the latter allowing all carbon system components and ocean dynamics to adjust. Thus, we still expect our model simulations to provide insight into baseline changes and redistribution of ocean DIC forced by the different LIG orbital configurations, supported by the important role of deep Atlantic water mass geometry changes coupled with its similarity to the long-term (millennial-scale) evolution inferred from proxy reconstructions.

5 Conclusions

Ongoing anthropogenic warming raises questions about the oceanic carbon sink and its efficiency under warmer climate conditions. In this study, we use the fully coupled NorESM model to simulate two quasi-equilibrium states of the Last Interglacial: one period is globally colder (115 ka) and one is globally warmer (125 ka) than today. We focus on the differences in ocean carbon cycle that occur at 125 ka in comparison to 115 ka, specifically the differences at global and basin scales. We provide a detailed description of the biogeochemical and physical processes that are responsible for the ocean carbon inventory changes under warmer climate conditions during the LIG at the temporal and spatial scales discussed here.

We find that the global ocean carbon budget decreases during the warm (125 ka) period by 314.1 PgC and is related mainly to better ventilation in the interior ocean. The Pacific Ocean has the largest reduction and accounts for 57 % of the global DIC loss. The response of the Pacific ventilation in a warmer climate shown in this study is consistent with previous studies (Menviel et al., 2014, 2015). The Indian and Atlantic basins account for 28 % and 15 %, respectively. These quantities mostly reflect basin volumes. The SSWs are revealed to play an instrumental role for the DIC changes in the interior below 1000 m depth. In these waters, the Atlantic is highlighted to be the region where the strongest DIC loss occurs per unit volume and is characterized by a stronger ventilation and a DIC^tot decrease of about 37 % compared to its respective value at 115 ka.

The reduced DIC budget at 125 ka occurs mostly in the interior ocean, while there is a weak increase in the top 1000 to 1500 m depths. Two factors that contribute most to the drop in the DIC budget in the interior ocean are (1) a weaker biological component from both the soft-tissue and the carbonate pumps that dominates at the depth between 1000 and 3000 m, and (2) a stronger disequilibrium effect (i.e., more negative) of DIC in the bottom waters. However, the processes affecting the disequilibrium component can arise from different factors such as changes in the physical pump, overturning circulation, or biological pump. No general process could be attributed to its variation, which seems to be regionally affected. While the SSWs seem to become more undersaturated at 125 ka, the NSWs seem to be more saturated. Further experiments with for instance fixed biological productivity or overturning circulation could help to identify the sensitivity of this component to such factors, but they remain too expensive to perform with our model.

The weakening of the biological component at depth is driven by younger water masses simulated in the interior ocean. This decrease in residence time of the water masses is caused by the strong SST modifications that affect the ventilation at 125 ka as compared to 115 ka. Higher SST, especially at high latitudes, induces strong summer sea ice retreat in the Atlantic sector of the Southern Ocean and stratification in the Pacific Ocean. In the Atlantic basin, this results in a more southerly confined SSW and southward expansion of NSW in the deep ocean. These water masses are advected by the Antarctic circumpolar current into the Indian and the western Pacific basins. The eastern Pacific Ocean is influenced by water masses coming from the Pacific sector of the Southern Ocean, with a warmer SST that hinders the ventilation and increases the residence time of the interior water masses on the eastern side of the basin.

Concerning the modification in the upper-ocean productivity under warmer climatic conditions, our model study reveals clear yet heterogeneous changes in phosphate availability and carbon export production especially between the Atlantic and Pacific basins. Such inter-basinal response in the biogeochemical divide is also highlighted by Moore et al. (2018) for future projections under warmer climatic conditions. This implies that changes in the biogeochemical divide could be somewhat similarly impacted from past and future anthropogenic CO₂ forcings, although the basin-specific responses suggest that it may not be a priori simple to predict the pattern or sign of the response of large-scale productivity to a given common forcing. Given the economic importance of basin scale productivity and the sensitivity found in past and future simulations, reconstructing and understanding the pattern and validating the sign and (model) response of large scale productivity to climate forcing are therefore critical for assessing not only the sign but also the sensitivity of global productivity to climate change.

The remaining uncertainties include, among others, the use of preindustrial states for some boundary conditions and the absence of freshwater input, which could modify the spatial response, particularly during the early interglacial period, which might include the final episodes of continental deglaciation. This is due to the lack of knowledge on such forcing during past climate conditions. Additional model-based studies using different Earth system models would be useful to confirm the robustness of our finding and further improve our understanding of the carbon dynamics and the feedback in the ocean under warmer climate conditions. Finally, our model-based study suggests that past warm periods experienced considerable carbon cycle and ocean DIC changes, linked to the response of the interior-ocean ventilation and biological productivity to high-latitude warming and interglacial background climate differences. It also suggests that the Atlantic part of the Southern Ocean is most sensitive to past climate change and hence could be a potential indicator of similar large-scale circulation changes in the future. Close monitoring of this region could be critical to better understand climate feedback on the carbon cycle under future warmer climate conditions.