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.

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 CO2 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 CO2 change such as peatland carbon dynamics and CaCO3
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
CO2, 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.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
CO2 (pCO2) according to the Ocean Carbon-Cycle Model
Intercomparison Project (OCMIP) protocols. The pCO2 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 : ΔO2= 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−1
and is remineralized at a constant rate of 0.02 day−1 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−1, 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 CO2 and O2 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 O2, TALK, and
PO4 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 (CO2,
CH4, N2O). For the experiment at 125 ka (115 ka), the
atmospheric CO2, CH4, and N2O 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−1 (−0.06 PgC yr−1) and
−0.01 Pmol yr−1 (−0.01 Pmol yr−1), 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 pCO2
and TALK budget are small between the two experiments. Such changes would
affect the DIC budget by about 32 PgC. The CO2 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−1
at 125 ka. In addition, the difference in sedimentation rates between the
two experiments (∼6.10-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):

where the DIC at saturation (DICsat) describes the DIC
concentration when the water parcel is in full equilibrium with the
atmospheric CO2 when it is last in contact with the surface.
In our model, changes in DICsat 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 (TALKpre), preformed
phosphate (PO4pre), surface silicate, salinity, and
temperature. In addition, the atmospheric CO2 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 (DICdis) measures
the disequilibrium state of the surface water with respect to the atmosphere.
This parameter is computed from the DICtot (output),
from which the biological and saturated DIC
components (DICbio and DICsat) are subtracted.
Therefore, all components are included in its calculation. A negative
DICdis 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 CO2 consumption at the surface,
which tends to increase the timescale needed for the water parcel to
equilibrate. However, a positive DICdis 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 DICdis and
DICsat 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).

(3)DICbio=DICsoft+DICcarb

The remineralization of soft tissues (hereafter DICsoft)
contributes via phosphate (PO4) remineralization through a carbon
phosphorus stoichiometric ratio rC:P=122. This component is
calculated from the difference between the total and the preformed
PO4 according to

(4)DICsoft=rC:P(PO4tot-PO4pre).

The carbonate pump contributes through the dissolution of CaCO3 hard
shells, calculated as difference between the total and the preformed
alkalinity and PO4 following

(5)DICcarb=0.5[TALKtot-TALKpre+rN:P(PO4tot-PO4pre)],

where rN:P=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
ΔDICtot, ΔDICsat,
ΔDICsoft, ΔDICcarb,
and ΔDICdis 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

(6)PO=O2+rO:PPO4,

where rO:P=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 1Difference 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.

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).

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−3 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 3Zonally 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 (∼60∘ S) compared to 115 ka (∼58∘ 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 4Zonally 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−3).

Table 1Difference in southern-sourced water (ΔSSW) in
the global ocean and for each basin. Row 1 shows the volume
(VSSW) according to our PO ≥0.57 mol O m−3
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
(ΔVSSW, ΔAgeSSW, and
ΔDICSSW) between 125 and 115 ka is the greatest
in the Atlantic (−37 %, −262 years, and −0.92 g C m−3). 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 (PO4) 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
PO4 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 5Difference in (a) export production of carbon at 100 m
depth (ΔEPC) and (b) phosphate at the surface
(ΔPO4) 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
DICtot, DICsat, DICsoft,
DICcarb, and DICdis are averaged over a 500 m
depth interval. The global amount of DICtot is −314.1 PgC in
the ocean under the warmer conditions (Fig. 6a, gray
ΔDICtot). Since the atmospheric CO2 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 ΔDICtot, 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 ΔDICtot 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,
ΔDICtot 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 6DIC differences between 125 and 115 ka (ΔDICx)
in (a) the global ocean and the (b) Atlantic, (c) Indian, and
(d) Pacific basins. The ΔDICx is averaged over
a 500 m depth interval where “x” refers to the different components of
the DIC. The DICtot is represented by the gray bars and is
decomposed into its four components: ΔDICsat (blue),
ΔDICsoft (green),
ΔDICcarb (orange), and
ΔDICdis (purple). The sum throughout the water
column of each components is given by the legend.

The ΔDICtot 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
ΔDICtot are controlled by the changes in
ΔDICsat. 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 ΔDICsat throughout the water column.

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 ΔDICtot (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
ΔDICsoft and
ΔDICcarb, translating to a less effective
soft-tissue and carbonate pump at 125 ka. Positive changes in
ΔDICtot 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 CO2, inducing ΔDICdis to
be more negative. The negative values of ΔDICdis
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 ΔDICdis (Fig. 7d). However, the
upper layers of the Atlantic Ocean are mostly simulated with higher
DICdis (positive ΔDICdis, 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
TALKpre (not shown here), possibly provoked by the melting of the
sea ice.

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 ΔDICtot
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
ΔDICtot 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
ΔDICtot, 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 ΔDICtot 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
ΔDICsoft and Fig. 5a suggest). The strong
positive ΔDICdis 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
ΔDICsat 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 DICsat.

The Pacific basin shows the strongest DICtot 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 ΔDICdis 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.

The ocean plays an important role in storing carbon and, thus, in the long
term regulation of atmospheric CO2 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 δ13C, a water mass
tracer strongly but inversely related to carbon and PO4 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
δ13C decreases in different regions proximal to the Southern
Ocean (Govin et al., 2009; Ninnemann and Charles, 2002). Bottom-water δ13C 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.

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 CO2 is fixed in our model,
preventing a direct assessment of ocean carbon changes in atmospheric
CO2, the decrease in deep carbon storage and shoaling of the DIC pool
during the warm 125 ka interval is generally consistent with higher
atmospheric CO2 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 CO2 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 DICtot 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 CO2 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 CO2 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 CO2 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.

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 DICtot 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 CO2 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.

We thank Victor Brovkin and the two anonymous referees for their positive and
constructive comments, which helped to clarify the manuscript. We also thank
the editor, Laurie Menviel, for the time she dedicated in processing our
manuscript and for her additional comments and suggestions. We thank
Nadine Goris for her valuable feedback on the first draft of the manuscript.
We are grateful to Chuncheng Guo for his technical assistance. This work was
supported by the Research Council of Norway-funded
projects THRESHOLDS (254964) and ORGANIC (239965) and
the Bjerknes Centre for Climate Research project BIGCHANGE. We acknowledge
the Norwegian Metacenter for Computational Science and Storage Infrastructure
(Notur/Norstore) projects nn2345k, ns2345k, nn1002k, and ns1002k for
providing the computing and storing resources essential for this
study.

Menviel, L., Joos, F., and Ritz, S. P.: Simulating atmospheric CO2,
δ13C and the marine carbon cycle during the Last
Glaciale/Interglacial cycle: possible role for a deepening of the mean
remineralization depth and an increase in the oceanic nutrient inventory,
Quaternary Sci. Rev., 56, 46–68,
https://doi.org/10.1016/j.quascirev.2012.09.012, 2012. a

We analyze the changes in oceanic carbon dynamics, using a state-of-the-art Earth system model, by comparing two quasi-equilibrium states: the early, warm Eemian (125 ka) versus the cooler, late Eemian (115 ka). Our results suggest a considerably weaker ocean dissolved inorganic carbon storage at 125 ka, an alteration of the deep-water geometry and ventilation in the South Atlantic, and heterogeneous changes in phosphate availability and carbon export between the Pacific and Atlantic basins.

We analyze the changes in oceanic carbon dynamics, using a state-of-the-art Earth system model,...