Journal cover Journal topic
Biogeosciences An interactive open-access journal of the European Geosciences Union
Journal topic
Biogeosciences, 15, 3761-3777, 2018
https://doi.org/10.5194/bg-15-3761-2018
Biogeosciences, 15, 3761-3777, 2018
https://doi.org/10.5194/bg-15-3761-2018

Research article 20 Jun 2018

Research article | 20 Jun 2018

# The devil's in the disequilibrium: multi-component analysis of dissolved carbon and oxygen changes under a broad range of forcings in a general circulation model

Disequilibrium DIC
Sarah Eggleston1,a and Eric D. Galbraith1,2,3 Sarah Eggleston and Eric D. Galbraith
• 1Institut de Ciència i Tecnologia Ambientals (ICTA), Universitat Autònoma de Barcelona, 08193 Barcelona, Spain
• 2Institució Catalana de Recerca i Estudis Avançats (ICREA), Pg. Lluís Companys 23, 08010 Barcelona, Spain
• 3Department of Earth and Planetary Science, McGill University, Montréal, Québec H3A 2A7, Canada
• anow at: Laboratory for Air Pollution & Environmental Technology, Empa, Überlandstrasse 129, 8600 Dübendorf, Switzerland
Abstract

The complexity of dissolved gas cycling in the ocean presents a challenge for mechanistic understanding and can hinder model intercomparison. One helpful approach is the conceptualization of dissolved gases as the sum of multiple, strictly defined components. Here we decompose dissolved inorganic carbon (DIC) into four components: saturation (DICsat), disequilibrium (DICdis), carbonate (DICcarb), and soft tissue (DICsoft). The cycling of dissolved oxygen is simpler, but can still be aided by considering O2, ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{sat}}}$, and ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$. We explore changes in these components within a large suite of simulations with a complex coupled climate–biogeochemical model, driven by changes in astronomical parameters, ice sheets, and radiative forcing, in order to explore the potential importance of the different components to ocean carbon storage on long timescales. We find that both DICsoft and DICdis vary over a range of 40 µmol kg−1 in response to the climate forcing, equivalent to changes in atmospheric pCO2 on the order of 50 ppm for each. The most extreme values occur at the coldest and intermediate climate states. We also find significant changes in O2 disequilibrium, with large increases under cold climate states. We find that, despite the broad range of climate states represented, changes in global DICsoft can be quantitatively approximated by the product of deep ocean ideal age and the global export production flux. In contrast, global DICdis is dominantly controlled by the fraction of the ocean filled by Antarctic Bottom Water (AABW). Because the AABW fraction and ideal age are inversely correlated among the simulations, DICdis and DICsoft are also inversely correlated, dampening the overall changes in DIC. This inverse correlation could be decoupled if changes in deep ocean mixing were to alter ideal age independently of AABW fraction, or if independent ecosystem changes were to alter export and remineralization, thereby modifying DICsoft. As an example of the latter, we show that iron fertilization causes both DICsoft and DICdis to increase and that the relationship between these two components depends on the climate state. We propose a simple framework to consider the global contribution of DICsoft+DICdis to ocean carbon storage as a function of the surface preformed nitrate and DICdis of dense water formation regions, the global volume fractions ventilated by these regions, and the global nitrate inventory.

1 Introduction

The controls on ocean carbon storage are not yet fully understood. Although potentially very important, given the large inventory of dissolved inorganic carbon (DIC) the ocean contains (currently 38 000 Pg C vs. 700 Pg C in the pre-industrial atmosphere), the nuances of carbon chemistry, the dependence of air–sea exchange on wind stress and sea ice cover, the intricacies of ocean circulation, and the activity of the marine ecosystem all contribute to making it a very complex problem. The scale of the challenge is such that, despite decades of work, the scientific community has not yet been able to satisfactorily quantify the role of the ocean in the natural variations of pCO2 between 180 and 280 ppm that occurred over ice age cycles. This failure reflects persistent uncertainty that also impacts our ability to accurately forecast future ocean carbon uptake .

In order to help with process understanding, proposed conceptualizing ocean carbon storage as consisting of a baseline surface–ocean average, enhanced by two “pumps” that transfer carbon to depth: the solubility pump, produced by the vertical temperature gradient, and the soft-tissue pump, produced by the sinking and downward transport of organic matter. This conceptualization proved very useful, but it fails to deal explicitly with the role of spatially and temporally variable air–sea exchange, and cannot effectively address changes in ocean circulation. A number of other conceptual systems have been employed (e.g., Broecker et al.1985; Gruber et al.1996), for considering both natural changes in the carbon cycle of the past and the anthropogenic transient input of carbon into the ocean.

Here, we use the decomposition laid out by , with the small change that we consider only DIC, rather than total carbon. This theoretical framework defines four components that, together, add up to the total DIC: saturation (DICsat), disequilibrium (DICdis), carbonate (DICcarb), and soft tissue (DICsoft; ). The first two components are “preformed” quantities (DIC${}_{\mathrm{pre}}={\mathrm{DIC}}_{\mathrm{sat}}+{\mathrm{DIC}}_{\mathrm{dis}}$), i.e., they are defined in the surface layer of the ocean and are carried passively by ocean circulation into the interior. In contrast, the latter two are equal to zero in the surface layer and accumulate in the interior due to biogeochemical activity (see Fig. 1). Note that the four components are diagnostic quantities only, intended to aid in understanding mechanisms and clarifying hypotheses, and do not influence the behavior of the model (although they can be calculated more conveniently by including additional ocean model tracers, as described in Methods).

Figure 1Illustration of the decomposition framework used for DIC in this paper. In the surface ocean, DIC is equal to DIC${}_{\mathrm{pre}}={\mathrm{DIC}}_{\mathrm{sat}}+{\mathrm{DIC}}_{\mathrm{dis}}$. Carbon taken up by biological processes in the surface ocean sinks and remineralizes in the water column to add two additional components at depths: DIC${}_{\mathrm{rem}}={\mathrm{DIC}}_{\mathrm{soft}}+{\mathrm{DIC}}_{\mathrm{carb}}$.

Saturation DIC is simply determined by the atmospheric CO2 concentration and its solubility in seawater, which is a function of ocean temperature, salinity, and alkalinity. For example, cooling the ocean will increase CO2 solubility, thereby leading to an increase in DICsat. Given known changes in temperature, salinity, alkalinity, and atmospheric pCO2, the effective storage of DICsat can be calculated precisely.

At the ocean surface, primary producers take up DIC. The organic carbon that is formed then sinks or is subducted (as dissolved or suspended organic matter) and is transformed into remineralized DIC within the water column (a small fraction is buried at depth). Here we define DICsoft as that which accumulates through the net respiration of organic matter below the top layer of the ocean (in our model, the uppermost 10 m). Thus, DICsoft depends both on the export flux of organic matter, affected by surface ocean conditions including nutrient supply , and the ocean circulation as a whole, including the surface-to-deep export and the flushing rate of the deep ocean, which clears out accumulated DICsoft . The Southern Ocean (SO) is thought to be an important region for such changes on glacial/interglacial timescales, as the ecosystem there is currently iron-limited, and it also plays a major role in deep ocean ventilation . Assuming a constant global oceanic phosphate inventory and constant C : P ratio, DICsoft is stoichiometrically related to the preformed ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ (${\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}$) inventory of the ocean, where ${\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}$ is the concentration of ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ in newly subducted waters, and a passively transported tracer in the interior. The potential to use ${\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}$ as a metric of DICsoft prompted very fruitful efforts to understand how it could change over time , though it has been pointed out that the large variation in C : P of organic matter could weaken the relationship between DICsoft and ${\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}$ . Given that the variability of N : C is significantly smaller than P : C , we use preformed nitrate in the discussion below.

Similar to DICsoft, DICcarb is defined here as the DIC generated by the dissolution of calcium carbonate shells below the ocean surface layer. Note that this does not include the impact that shell production has at the surface; calcification causes alkalinity to decrease in the surface ocean, raising surface pCO2 and shifting carbon to the atmosphere. Rather, within the framework used here, this effect on alkalinity distribution falls under DICsat, since it alters the solubility of DIC. This highlights an important distinction between the four-component framework, which strictly defines subcomponents of DIC, and the “pump” frameworks, which provide looser descriptions of vertical fluxes of carbon and, in some cases, alkalinity. Changes in DICcarb on the timescales of interest are generally thought to be small compared to those of DICsat and DICsoft.

Typically, only these three components are considered as the conceptual drivers behind changes in the air–sea partitioning of pCO2 . However, a fourth component, DICdis, is also potentially significant as discussed by . Defined as the difference between preformed DIC and DICsat, DICdis can be relatively large because of the slow timescale of atmosphere–surface ocean equilibrium of carbon compared to other gases, caused by the reaction of CO2 with seawater . In short, DICdis is a function of all non-equilibrium processes on DIC in the surface ocean, including biological uptake, ocean circulation, and air–sea fluxes of heat, freshwater, and CO2.

Like DICsat, DICdis is a conservative tracer determined in the surface ocean, with no sources or sinks in the ocean interior. Since the majority of the ocean is filled by water originating from small regions of the Southern Ocean and the North Atlantic, the net whole-ocean disequilibrium carbon is approximately determined by the DICdis in these areas weighted by the fraction of the ocean volume filled from each of these sites. Unlike the other three components, DICdis could contribute either additional oceanic carbon storage (DICdis> 0) or reduced oceanic carbon storage (DICdis<0). Although this parameter is implicitly included in most models, studies using preformed nutrients as a metric for biological carbon storage have often ignored the potential importance of DICdis by assuming fast air–sea gas exchange (e.g., Marinov et al.2008a; Ito and Follows2005). In the pre-industrial ocean this is of little importance, given that global DICdis is small because the opposing effects of North Atlantic and Antarctic water masses largely cancel each other. However, showed that DICdis can have a large impact by amplifying changes in DICsoft under constant pre-industrial ocean circulation, and have very recently shown that DICdis can vary significantly in response to changes in ocean circulation states.

The cycling of dissolved oxygen is simpler than DIC. Because O2 does not react with seawater or the dissolved constituents thereof, it has no dependence on alkalinity, and its equilibration with the atmosphere through air–sea exchange occurs approximately 1 order of magnitude faster than DIC (on the order of 1 month, rather than 1 year). Nonetheless, dissolved oxygen can be conceptualized as including a preformed component, which is the sum of saturation and disequilibrium, and an oxygen utilization component, which is given by the difference between the in situ and preformed O2 in the ocean interior. Apparent oxygen utilization (AOU), typically taken as a measure of accumulated respiration, can be misleading if the preformed O2 concentration differed significantly from saturation, i.e., if ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ is significant, as it appears to be in high-latitude regions of dense water formation . If ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ varies with climate state, it might contribute significantly to past or future oxygen concentrations.

Here, we use a complex Earth system model to investigate the potential changes in the constituents of DIC and O2 on long timescales, relevant for past climate states as well as the future. We make use of a large number of equilibrium simulations, conducted over a wide range of radiative, orbital, and ice sheet boundary conditions, as a “library” of contrasting ocean circulations in order to test the response of disequilibrium carbon storage to physically plausible changes in ocean circulation. The basic physical aspects of these simulations were described by . We supplement these with a smaller number of iron fertilization experiments to examine the additional impact of circulation-independent ecosystem changes. In order to simplify the interpretation, we chose to prescribe a constant pCO2 of 270 ppm for the air–sea exchange in all simulations, as also done by . Thus, the changes in DICsat reflect only changes in temperature, salinity, alkalinity, and ocean circulation arising from the climate response, and not changes in pCO2. Nor do they explicitly consider changes in the total carbon or alkalinity inventories driven by changes in outgassing and/or burial ; rather, the alkalinity inventory is fixed, and the carbon inventory varies due to changes in total ocean DIC (since the atmosphere is fixed). As such, the experiments here should be seen as idealized climate-driven changes, and should be further tested with more comprehensive models including interactive CO2.

2 Methods

## 2.1 Model description

The global climate model (GCM) used in this study is CM2Mc, the Geophysical Fluid Dynamics Laboratory's Climate Model version 2 but at lower resolution (3), described in more detail by and modified as described by . This includes the Modular Ocean Model version 5, a sea ice module, static land, and static ice sheets, and a module of Biogeochemistry with Light, Iron, Nutrients and Gases (BLINGv1.5; ). Unlike BLINGv0, BLINGv1.5 allows for variable P : C stoichiometry using the “line of frugality” and calculates the mass balance of phytoplankton in order to prevent unrealistic bloom magnitudes at high latitudes, reducing the magnitude of disequilibrium O2, which was very high in BLINGv0 . In addition, three water mass tracer tags are defined: a southern tracer south of 30 S, a North Atlantic tracer north of 30 N in the Atlantic, and a North Pacific tracer north of 30 N in the Pacific. An “ideal age” tracer is also defined as zero in the global surface layer, and increasing in all other layers at a rate of 1 year per year.

## 2.2 Experimental design

The model runs analyzed here are part of the same suite of simulations discussed by . A control run was conducted with a radiative forcing equivalent to 270 ppm atmospheric pCO2 and the Earth's obliquity and precession set to modern values (23.4 and 102.9, respectively). Experimental simulations were run at values of obliquity (22, 24.5) and precession (90, 270) representing the astronomical extremes encountered over the last 5 Myr , while eccentricity was held constant at 0.03. A range of greenhouse radiative forcings was imposed equivalent to pCO2 levels of 180, 220, 270, 405, 607, or 911 ppm; with reference to a pre-industrial radiative forcing, the radiative forcings are roughly equal to −2.2, −1.1, 0, +2.2, +4.3, and +6.5 W m−2, respectively .

The biogeochemical component of the model calculates air–sea carbon fluxes using a fixed atmospheric pCO2 of 270 ppm throughout all model runs. Note that 270 ppm was chosen to reflect an average interglacial level, rather than specifically focusing on the pre-industrial climate state. This use of constant pCO2 for the carbon cycle means that the DICsat is not consistent with the pCO2 used for the radiative forcing, so that changes in DICsat caused by a given pCO2 change tend to be larger than would be expected in reality. This has a negligible effect on the other carbon components, given that they do not depend directly on pCO2; this has been confirmed by model runs with the University of Victoria Earth System Model using a similar decomposition strategy (Samar Khatiwala, personal communication, 2018).

Eight additional runs were conducted using Last Glacial Maximum (LGM) ice sheets with the lowest two radiative forcings and the same orbital parameters. Iron fertilization simulations calculate the input flux of dissolved iron to the ocean surface assuming a constant solubility and using the glacial atmospheric dust field of as modified by instead of the standard pre-industrial dust field; note that this is not entirely in agreement with more modern reconstructions, which could potentially have an influence on the induced biological blooms, both in magnitude and geographically (e.g., Albani et al.2012, 2016). Four iron fertilization experiments were run with the lowest radiative forcing with LGM ice sheets, as well as one model run similar to the control run. Finally, two simulations were run that were identical to the pre-industrial setup, but the rate of remineralization of sinking organic matter is set to 75 % of the default rate, approximately equivalent to the expected change due to 5 C ocean cooling ; one of these runs also includes iron fertilization. All simulations are summarized in Table 1.

Table 1Simulation overview. A total of 44 simulations were analyzed with varying radiative forcing (RF), obliquity, precession, ice sheets (PI, pre-industrial; LGM, Last Glacial Maximum reconstruction; LGM*, topography of LGM ice sheets but with PI albedo), and with and without iron fertilization. Runs 1–40 are described by . Runs 43 and 44 and identical to 41 and 42 but the remineralization rate of sinking organic matter is reduced by 25 %.

In the following, three particular runs are highlighted for comparison to illustrate cold (CW), moderate (MW), and hot (HW) worlds. These include radiative forcings of −2.2, 0, and +4.3 W m−2, respectively; the former includes LGM ice sheets; and the obliquity and precession is 22 and 90 for glacial-like (GL) and 22 and 270 for PI and WP. These specific runs are distinguished from GL and interglacial-like (IG) scenarios, which refer to averages of four runs each, each with a radiative forcing of −2.2 and 0 W m−2, respectively; the GL runs also have LGM ice sheets.

All simulations were run for 2100–6000 model years beginning with a pre-industrial spinup. While the model years presented here largely reflect runs after having reached steady state, it is important to note that the pre-industrial run (41 in Table 1) still has a drift of 1 µmol kg−1 over the 100 years shown here and thus may not yet be at steady state.

## 2.3 Decompositions

The four-component DIC scheme and three-component O2 scheme described in the introduction can be exactly calculated for any point in an ocean model using five easily implemented ocean model prognostic tracers: DIC, DICpre, DICsat, O2, and O2pre, since DIC${}_{\mathrm{dis}}={\mathrm{DIC}}_{\mathrm{pre}}-{\mathrm{DIC}}_{\mathrm{sat}}$, DIC${}_{\mathrm{soft}}=\left({\mathrm{O}}_{\mathrm{2}\mathrm{pre}}-{\mathrm{O}}_{\mathrm{2}}\right)\cdot {r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}}$, DIC${}_{\mathrm{carb}}=\mathrm{DIC}-{\mathrm{DIC}}_{\mathrm{pre}}-{\mathrm{DIC}}_{\mathrm{soft}}$, and O2dis = O2pre O2sat (O2sat can be accurately calculated from the conservative tracers temperature and salinity). In the absence of a DICsat tracer, it can be estimated from alkpre, temperature, and salinity. For this large suite of simulations, we only had DIC, alk, O2, and O2pre available. Thus, although we could calculate O2dis and DICsoft directly, it was necessary to use an indirect method to calculate the other three carbon tracers. Following , we first estimate alkpre using a regression, then calculate DICsat using alkpre, temperature, salinity, and the known atmospheric pCO2 of 270 ppm. DICcarb is calculated as $\mathrm{0.5}\cdot \left(\mathrm{alk}-{\mathrm{alk}}_{\mathrm{pre}}\right)$, and DICdis is then calculated as a residual. For more details and an estimate of the error in the method, see the Appendix.

3 Results and discussion

## 3.1 General climate response to forcings

Differences in the general ocean state among the model simulations are described in detail by . We provide here a few key points, important for interpreting the dissolved gas simulations. First, the Antarctic Bottom Water (AABW) fraction of the global ocean varies over a wide range among the simulations, with abundant AABW under low and high global average surface air temperature (SAT), and a minimum at intermediate SAT (Fig. 2b). The North Atlantic deep water (NADW) fraction is approximately the inverse of this. The ventilation rate of the global ocean is roughly correlated with the AABW fraction, with rapid ventilation (small average ideal age) when the AABW fraction is high (Fig. 2). The density difference between AABW and NADW is the overall driver of the AABW fraction and ventilation rate in the model simulations. In all simulations colder than the present day, the density difference can be explained by the effect of sea ice cycling in the Southern Ocean on the global distribution of salt: when the Southern Ocean is cold and sea ice formation abundant, there is a large net transport of freshwater out of the circum-Antarctic, causing AABW to become more dense. NADW becomes consequently fresher, because there is less salt left to contribute to NADW. As noted by , the simulated ventilation rates should be viewed with caution, given the poor mechanistic representation of diapycnal mixing in the deep ocean, and the AABW fraction may not be correlated with ventilation rate in the real ocean.

Figure 2Simulated global ocean ventilation and proportion of southern-sourced deepwater. The horizontal axis shows the globally averaged surface air temperature (SAT). (a) The global ocean average ideal age (low age corresponds to rapid ventilation) and (b) the fraction of water in the global ocean originating from the surface south of 30 S are anti-correlated over the range of SAT in these model runs. Orange and blue symbols represent high and low obliquity scenarios, respectively; triangles pointing upward and downward represent greater Northern and Southern hemisphere seasonality (precession 270 and 90); outlines are scenarios with LGM ice sheets; light shading indicates scenarios with LGM ice sheet topography but PI albedo. The size of the symbols corresponds to the SAT.

## 3.2 General changes in DIC

Total DIC generally decreases from cold to warm simulations, under the constant pCO2 of 270 ppm used for air–sea exchange. Changes in DICsat drive the largest portion of this trend, decreasing approximately linearly with surface air temperature due to the temperature-dependence of CO2 solubility, resulting in a difference of 50 µmol kg−1 over this range (see Fig. 3). Because of the constant pCO2 of 270 ppmv in the biogeochemical module, the simulated range of DICsat is larger than would occur if the prescribed biogeochemical pCO2 were the same as that used to produce the low SAT . DICcarb is relatively small in magnitude and generally increases with SAT, but has a standard deviation of only 4 µmol kg−1 over all simulations, so we do not discuss it further.

Figure 3Global average DIC and separate components in simulations 1–36 as a function of globally averaged surface air temperature (SAT). Orange and blue symbols represent high and low obliquity scenarios, respectively; triangles pointing upward and downward represent greater Northern and Southern hemisphere seasonality (precession 270 and 90); outlines are scenarios with LGM ice sheets; light shading indicates scenarios with LGM ice sheet topography but PI albedo.

In contrast to DICsat and DICcarb, DICdis and DICsoft vary nonlinearly with global temperatures, with a clear and shared turning point near the middle of the temperature range. Thus, the extreme values of DICsoft and DICdis occur under the coldest state and at an intermediate state close to the pre-industrial. Both DICdis and DICsoft are strongly correlated with ocean ventilation, quantified here by the global average of the ideal age tracer (r2=0.69 and 0.89, respectively), and thus with each other (r2=0.74). However, whereas found a positive correlation of DICsoft and DICdis under nutrient depletion experiments with constant climate, these experiments indicate that, when driven by the wide range of physical changes explored here, DICsoft and DICdis are negatively correlated.

Simulated changes in DICdis are of the same magnitude as the DICsoft changes, to which much greater attention has been paid. For a global average buffer factor between 8 and 14 , a rough, back-of-the-envelope calculation shows that a 1 µmol kg−1 change in DIC corresponds to a 0.9–1.6 ppm change in atmospheric pCO2 based on a DIC concentration of 2300 µmol kg−1 and pCO2 of 270 ppm. Thus, the increase in the global average DICdis in these simulations could have contributed more than a 40 ppm change in the atmospheric pCO2 stored in the ocean during the glacial compared to today. It is important to recognize that the drawdown of CO2 by disequilibrium storage would have resulted in a decrease of DICsat, given the dependence of the saturation concentration on pCO2, so this estimate should not be interpreted as a straightforward atmospheric pCO2 change. Nonetheless, while this is only a first-order approximation and the model biases are potentially large, it seems very likely that the disequilibrium carbon storage was a significant portion of the net 90 ppm difference.

## 3.3 Climate-driven changes in DICsoft

The biogeochemical model used here is relatively complex, with limitation by three nutrients (N, P, and Fe), denitrification, and N2 fixation, in addition to the temperature- and light-dependence typical of biogeochemical models. The climate model is also complex, including a full atmospheric model, a highly resolved dynamic ocean mixed layer, and many nonlinear sub-grid-scale parameterizations, and uses short (< 3 h) timesteps. The simulations we show span a wide range of behaviors, including major changes in ocean ventilation pathways and patterns of organic matter export.

Thus, it is perhaps surprising that the net global result of the biological pump, as quantified by DICsoft, has highly predictable behavior. As shown in Fig. 4, the global DICsoft varies closely with the product of the global average sinking flux of organic matter at 100 m and the average ideal age of the global ocean. Qualitatively this is not a surprise, given that greater export pumps more organic matter to depth, and a large age provides more time for respired carbon to accumulate within the ocean. However, the quantitative strength of the relationship is striking. As demonstrated in Fig. 4, global DICsoft is not as well correlated with either of these parameters separately as it is with their product “age × export”.

Figure 4Globally averaged DICsoft vs. global export, ideal age, and age × export. Globally averaged DICsoft can be approximated remarkably well by the global export flux of organic carbon at 100 m multiplied by the average age of the ocean. The latter is an ideal age tracer in the model that is set to 0 at the surface and ages by 1 year each model year in the ocean interior. Orange and blue symbols represent high and low obliquity scenarios, respectively; triangles pointing upward and downward represent greater Northern and Southern hemisphere seasonality (precession 270 and 90); outlines are scenarios with LGM ice sheets; light shading indicates scenarios with LGM ice sheet topography but PI albedo. Red boxes indicate Fe fertilization simulations (runs 37–40). The size of the symbols corresponds to the SAT.

It is difficult to assess the likelihood that the real ocean follows this relationship to a similar degree. One reason it might differ is if remineralization rates vary spatially, or with climate state. In the model here, as in most biogeochemical models, organic matter is respired according to a globally uniform power law relationship with depth . showed that ocean carbon storage is sensitive to changes in these remineralization rates, which provides an additional degree of freedom. It is not currently known how much remineralization rates can vary naturally; they may vary as a function of temperature or ecosystem structure. As a result, the relationship between DICsoft and age × export may be stronger in the model than in the real ocean.

Nonetheless, the results suggest that, as a useful first-order approximation, the global change in DICsoft between two states can be given by a simple linear regression:

$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{soft}}\left[\mathrm{µ}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{-\mathrm{1}}\right]={m}_{\mathrm{1}}\cdot \mathrm{\Delta }\left(\mathrm{age}\phantom{\rule{0.125em}{0ex}}\left[\mathrm{yr}\right]×\mathrm{export}\phantom{\rule{0.125em}{0ex}}\left[\mathrm{Pg}\phantom{\rule{0.125em}{0ex}}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{-\mathrm{1}}\right]\right),\end{array}$

or in terms of pCO2:

$\begin{array}{}\text{(2)}& \mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{soft}}\left[\mathrm{ppm}\right]\approx {m}_{\mathrm{2}}\cdot \mathrm{\Delta }\left(\mathrm{age}\phantom{\rule{0.125em}{0ex}}\left[\mathrm{yr}\right]×\mathrm{export}\phantom{\rule{0.125em}{0ex}}\left[\mathrm{Pg}\phantom{\rule{0.125em}{0ex}}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{-\mathrm{1}}\right]\right).\end{array}$

Note that m2 is a function of the buffer factor and the climate state (atmospheric pCO2 and DIC). Based on the results here, m1=0.036 and ${m}_{\mathrm{2}}=\mathrm{0.065},\mathrm{0.042},\mathrm{0.029}$ for modern (405 ppm pCO2), pre-industrial (270 ppm) and glacial (180 ppm) conditions, respectively. Note that we have not varied pCO2 in these simulations, so these equations are only meant to illustrate the mathematical relationship observed in Fig. 4. This simple meta-model may provide a useful substitute for full ocean–ecosystem calculations, and should be further tested against other ocean–ecosystem coupled models with interactive CO2. Note that, as for the disequilibrium estimate above, the soft tissue pump CO2 drawdown is partially compensated for by a decrease in saturation carbon storage, so it would be larger than the net atmospheric effect. In addition, we have not accounted for consequent changes in the surface ocean carbonate chemistry (including changes in the buffer factor).

It is important to point out that the simulated change in DICsoft between interglacial and glacial states appears to be in conflict with reconstructions of the LGM. Proxy records appear to show that LGM dissolved oxygen concentrations were lower throughout the global ocean, with the exception of the North Pacific, implying greater DICsoft concentrations during the glacial then during the Holocene . In contrast, the model suggests that greater ocean ventilation rates in the glacial state (Fig. 2a) would have led to reduced global DICsoft. As discussed by , radiocarbon observations imply that the model ideal age is approximately 200 years too young under glacial conditions compared to the LGM, suggesting a circulation bias that may reflect incorrect diapycnal mixing or non-steady-state conditions. Whatever the cause, if we take this 200-year bias into account, the regression implies that an additional 33 µmol kg−1 DICsoft were stored in the glacial ocean. This would bring the simulated glacial DICsoft close to, but still less than, the simulated pre-industrial value.

We propose that the apparent remaining shortfall in simulated glacial DICsoft could reflect one or more of the following non-exclusive possibilities: (1) the model does not capture changes in remineralization rates caused by ecosystem changes; (2) the model underestimates the glacial increase in the nitrate inventory and/or growth rates, perhaps due to changes in the iron cycle; (3) the ocean was not in steady state during the LGM, and therefore not directly comparable to the GL simulation; (4) the inference of DICsoft from proxy oxygen records is incorrect due to significant changes in preformed oxygen disequilibrium (see below). If either of the first two possibilities is important, it would imply an inaccuracy in the meta-model derived here.

## 3.4 Climate-driven changes in DICdis

The ocean basins below 1 km depth are largely filled by surface waters subducted to depth in regions of deepwater formation . In our simulations, water originating in the surface North Atlantic, termed NADW, and the Southern Ocean, termed AABW, make up 80–96 % of this total deep ocean volume. Thus, to first order, the deep average DICdis concentration can be approximated by a simple mass balance:

$\begin{array}{ll}{\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{deep}}}\approx & \phantom{\rule{0.25em}{0ex}}{f}_{\mathrm{AABW}}\cdot {\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{AABW}}}+\left(\mathrm{1}-{f}_{\mathrm{AABW}}\right)\\ \text{(3)}& & \cdot {\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{NADW}}}.\end{array}$

Here, fAABW represents the fraction of deepwater originating in the SO, and DIC${}_{{\mathrm{dis}}_{\mathrm{AABW}}}$ and DIC${}_{{\mathrm{dis}}_{\mathrm{NADW}}}$ represent the DICdis concentrations at the sites of deepwater formation (see Fig. 5). North Atlantic deep water forms with negative DICdis, reflecting surface undersaturation, while the Southern Ocean is supersaturated (DICdis>0). These opposing tendencies between NADW and AABW cause a partial cancellation of DICdis when globally averaged, which makes the disequilibrium component small in the modern ocean. Theoretically, the simulated DICdis could change either due to changes in fAABW or the end-member compositions. Although the exact values of DICdis in the two polar oceans vary among the simulations in response to climate (the reasons for which are discussed in more detail below), these changes are small relative to the consistent large contrast between fAABW and fNADW, so that deep DICdis is strongly controlled by the global balance of AABW vs. NADW in each simulation (see Fig. 6). Global DICdis becomes much larger when fAABW is larger, similar to the dynamic evoked by . This is also illustrated by the depth transects of DICdis in Fig. 7.

Figure 5(a) Global average fraction of northern- (reddish colors) and southern-sourced (bluish colors) water. (b) Annual average values of DICdis of these water masses determined at 25 m depth in the model during model years and at locations where deep convection occurs. Pink and cyan (red and blue) symbols represent high (low) obliquity scenarios; triangles pointing upward and downward represent greater Northern and Southern hemisphere seasonality (precession 270 and 90); outlines are scenarios with LGM ice sheets; light shading indicates scenarios with LGM ice sheet topography but PI albedo. Five simulations, for which no deep convection events were identified south of 60S, are not shown. The size of the symbols corresponds to the SAT.

Figure 6Global average DICdis as a function of the fraction of the ocean below 1 km derived from the surface Southern Ocean. Orange and blue symbols represent high and low obliquity scenarios, respectively; triangles pointing upward and downward represent greater Northern and Southern hemisphere seasonality (precession 270 and 90); outlines are scenarios with LGM ice sheets; light shading indicates scenarios with LGM ice sheet topography but PI albedo. The size of the symbols corresponds to the SAT.

Figure 7DICdis (µmol kg−1) for simulations (a) cold world; (b) moderate world; (c) hot world (see Sect. 2.2). Depth transects represent the North Atlantic (left), Southern Ocean (center), and North Pacific (right).

We estimated the concentration of DICdis in the regions of AABW and NADW formation, shown in Fig. 5b. The end members vary less significantly than fAABW over the range of simulations, in part due to competing effects of different processes. As discussed in Sect. 3.1, simulations at both the low and high radiative forcing values used show increased AABW production, with a minimum at intermediate values (Fig. 5). The fact that the highest DIC${}_{{\mathrm{dis}}_{\mathrm{AABW}}}$ occurs at low SAT can be attributed to the rapid formation rate of AABW, while the intermediate SAT minimum in AABW volume explains the minimum in global ocean DICdis (Fig. 3). We note that expanded terrestrial ice sheets shift the ratio of AABW to NADW to higher values, due to their impact on NADW temperature and downstream expansion of Southern Ocean sea ice , further increasing DICdis in glacial-like conditions. Sea ice in the Southern Ocean would be expected to exert a further control over DICdis, as this reduces air–sea gas exchange, thus allowing carbon to accumulate beneath the ice. However we did not perform experiments to isolate this effect.

## 3.5 Climate-driven changes in ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$

Similarly to DICdis, ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ is defined as the departure from equilibrium of O2 in the surface ocean with respect to the atmosphere and is advected into the ocean as a conservative tracer. Unlike C, O2 does not react with seawater, and is present only as dissolved diatomic oxygen. Thus, O2 has a much shorter timescale of exchange at the ocean–atmosphere interface, equilibrating 1 order of magnitude faster than CO2. As a result, it is not sensitive to sea ice as long as there remains a fair degree of open water . But as the sea ice concentration approaches complete coverage, O2 equilibration rapidly becomes quite sensitive to sea ice. If there is a significant undersaturation of O2 in upwelling waters, the disequilibrium can become quite large (Fig. 8).

Figure 8${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ (µmol kg−1) for simulations (a) cold world; (b) moderate world; (c) hot world (see Sect. 2.2). Depth transects represent the North Atlantic (left), Southern Ocean (center), and North Pacific (right).

In the model simulations, the ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ in the Southern Ocean becomes as large as −100µmol kg−1 in the coldest states (note that DICdis and ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ are often anti-correlated). Because the disequilibrium depends on the O2 depletion of waters upwelling at the Southern Ocean surface, this could potentially be even larger if upwelling waters had lower O2. We do not place a large degree of confidence in these values, given the likely sensitivity to poorly resolved details of sea ice dynamics (e.g., ridging, leads) and dense water formation. Nonetheless, the potential for very large disequilibrium oxygen under cold states prompts the hypothesis that very extensive sea ice cover over most of the exposure pathway in the Southern Ocean might have made a significant contribution to the low O2 concentrations reconstructed for the glacial .

## 3.6 Iron fertilization experiments

In addition to the changes driven by ice sheets as well as orbital and radiative forcing, we conducted iron fertilization experiments under glacial and pre-industrial-like conditions, including a simulation with reduced remineralization rates (Fig. 9). As expected, both the global export and DICsoft increase when iron deposition is increased. However, the DICsoft increase is significantly lower in the well-ventilated GL simulations (2.9 µmol kg−1) compared to the PI simulation (7.3 µmol kg−1). This difference is qualitatively in accordance with the age × export relationship (Fig. 4), though with a smaller increase of DICsoft than would be expected from the export increase, compared to the broad spectrum of climate-driven changes. This reduced sensitivity of DICsoft to the global export can be attributed to the fact that the iron-enhanced export occurs in the Southern Ocean, presumably because the remineralized carbon can be quickly returned to the surface by upwelling when ventilation is strong. Thus, the impact of iron fertilization on DICsoft is strongly dependent on Southern Ocean circulation.

Figure 9Iron-fertilized changes in total global average DIC and each of the components (iron fertilization simulation minus associated control run). Simulations were either run under pre-industrial or glacial-like conditions (in the case of the latter, results represent the average of the four GL runs, with error bars showing the standard deviation within the four runs), as well as using 100 and 75 % of the default remineralization rate of organic matter. The close agreement of panels (a) and (b) indicates that the effects of iron fertilization and changes in the remineralization rate are approximately linearly additive in this model.

The iron addition also causes an increase of DICdis of approximately equal magnitude to DICsoft in the PI simulation and of relatively greater proportion in the GL simulations. Because the ocean in the GL simulations is strongly ventilated, the increase in export leads to an increase of DICdis, as remineralized DICsoft is returned to the Southern Ocean surface, where it has a relatively short residence time that limits outgassing to the atmosphere. With rapid Southern Ocean circulation, a good deal of the DIC sequestered by iron fertilization ends up in the form of DICdis, rather than DICsoft as might be assumed. Thus, just as the glacial state has a larger general proportion of DICdis compared to DICsoft, the iron addition under the glacial state produces a larger fraction of DICdis relative to DICsoft. The experiment in which the remineralization rate was reduced by 25 % indicates that the effects of iron fertilization alone on both DICsoft and DICdis are quite insensitive to the remineralization rate (see Fig. 9); for total DIC as well as for each component, the difference between the iron fertilization run and the corresponding control run is similar in panels (a) and (b). While we have not run a simulation under glacial-like conditions with a reduced remineralization rate, this suggests that the effects of iron fertilization and changes in the remineralization rate can be well approximated as being linearly additive in this model.

The tendency to sequester carbon as DICdis vs. DICsoft, in response to iron addition, can be quantified by the global ratio ΔDICdis ∕ ΔDICsoft. Our experiments suggest that this ratio is 0.9 for the pre-industrial state and 3.3 for the glacial-like state. Because of the circulation dependence of this ratio, it is expected that there could be significant variation between models. It is worth noting that found ΔDICdis ∕ ΔDICsoft of 2 in response to iron fertilization, using a modern ocean circulation, as analyzed by . We also note that the quantitative values of DICsoft and DICdis resulting from the altered iron flux should be taken with a grain of salt, given the very large uncertainty in iron cycling models .

These results provide a note of caution for interpreting iron fertilization model experiments, which might be assumed to act primarily on the soft tissue carbon storage. At moderate ventilation rates (the pre-industrial control run), an increase in iron results in an increase both in DICsoft, due to higher biological export, and DICdis, because of the upwelling of C-rich water resulting from higher remineralization, the effect discussed by . However, under a high ventilation state there is only a small increase in DICsoft in response to increased Fe in the surface ocean, as the remineralized carbon is quickly returned to the surface, thus producing a significant increase in DICdis only. Thus, the carbon storage resulting from Southern Ocean iron addition could actually be dominated by DICdis under some climate states, so that the overall impact may be significantly larger than would be predicted from DICsoft and/or O2 utilization.

## 3.7 A unified framework for DICdis and preformed nutrients

The concept of preformed nutrients allowed the production of a very useful body of work, striving for simple predictive principles. This work highlighted the importance of the nutrient concentrations in polar oceans where deep waters form as well as changes in the ventilation fractions of AABW and NADW, given their very different preformed nutrient concentrations . Although the variability of P : C ratios implies significant uncertainty for the utility of ${\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}$ in the ocean, the relative constancy of N : C ratios suggests that ${\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{pre}}}^{-}$ is indeed linked to DICsoft, inasmuch as the global N inventory is fixed .

Figure 10Simple estimation of global DICdis and DICsoft from water mass characteristics. Following Eq. (12) for (a) DICsoft and (b) DICdis separately, this model shows that the sum of the global average DICsoft and DICdis at steady state can be estimated fairly robustly as the result of a simple mass balance of the relevant parameters in the most important ventilating water masses. Here, we take into account upper-ocean water masses (above 1 km) formed in the North Pacific, North Atlantic, and Southern Ocean, and deep water masses formed in the North Atlantic and Southern Ocean. In each plot, the full model output is shown on the x axis and the result of the mass balance approximation on the y axis. Orange and blue symbols represent high and low obliquity scenarios, respectively; triangles pointing upward and downward represent greater Northern and Southern hemisphere seasonality (precession 270 and 90); outlines are scenarios with LGM ice sheets; light shading indicates scenarios with LGM ice sheet topography but PI albedo. Five simulations, for which no deep convection events were identified south of 60S, are not shown. The size of the symbols corresponds to the SAT. The purple square represents the pre-industrial simulation (run 41), and red boxes indicate Fe fertilization simulations (runs 37–40).

However, as shown by the analyses here, DICsoft – reflected by the preformed nutrients – is only half the story. Changes in DICdis can be of equivalent magnitude, and can vary independently of DICsoft as a result of changes in ocean circulation and sea ice. Nonetheless, we find that the same conceptual approach developed for DICsoft can be used to predict DICdis from the end member DICdis and the global volume fractions. The preformed relationships and DICdis can therefore be unified as follows (see Fig. 10):

$\begin{array}{}\text{(4)}& {\mathrm{DIC}}_{\mathrm{soft}}={\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{rem}}}^{-}\cdot {r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}.\end{array}$

Remineralized nitrate can be expressed in terms of the global nitrate inventory and the accumulated nitrate loss due to pelagic and benthic denitrification:

$\begin{array}{ll}\text{(5)}& & {\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{rem}}}^{-}={\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{global}}}^{-}-{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{pre}}}^{-}+{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{den}}}^{-},& {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{upper}}}}^{-}\approx {f}_{\mathrm{SO},\mathrm{upper}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{SO},\mathrm{upper}}}}^{-}+{f}_{\mathrm{NAtl},\mathrm{upper}}\\ \text{(6)}& & \phantom{\rule{1em}{0ex}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NAtl},\mathrm{upper}}}}^{-}+{f}_{\mathrm{NPac},\mathrm{upper}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NPac},\mathrm{upper}}}}^{-},\text{(7)}& & {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{deep}}}}^{-}\approx {f}_{\mathrm{AABW}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{AABW}}}}^{-}+{f}_{\mathrm{NADW}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NADW}}}}^{-}.\end{array}$

Because there is production of intermediate water but no deep convection in the North Pacific, we calculate this mass balance for the upper ocean (above 1 km) and deep ocean separately, dropping the Pacific Ocean term in Eq. (7) for the deep ocean. For brevity, we continue with the derivation for the deep ocean only; the upper ocean follows analogously.

$\begin{array}{ll}& {\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{deep}}}\approx {r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot \left[{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{deep}}}^{-}+{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{den},\mathrm{deep}}}^{-}\\ \text{(8)}& & \phantom{\rule{1em}{0ex}}-\left({f}_{\mathrm{AABW}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{AABW}}}}^{-}+{f}_{\mathrm{NADW}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NADW}}}}^{-}\right)\right]\end{array}$

Combining with Eq. (3),

$\begin{array}{ll}& {\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{deep}}}+{\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{deep}}}\approx {r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot \left({\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{deep}}}^{-}+{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{den},\mathrm{deep}}}^{-}\right)\\ & \phantom{\rule{1em}{0ex}}+{f}_{\mathrm{AABW}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{AABW}}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{AABW}}}}^{-}\right)\\ \text{(9)}& & \phantom{\rule{1em}{0ex}}+{f}_{\mathrm{NADW}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{NADW}}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NADW}}}}^{-}\right).\end{array}$

Finally, the global average is computed by summing the volume-weighted values in the upper and deep ocean:

$\begin{array}{ll}& {\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{global}}}+{\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{global}}}\approx \left[{V}_{\mathrm{upper}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{upper}}}+{\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{upper}}}\right)\\ \text{(10)}& & \phantom{\rule{1em}{0ex}}+{V}_{\mathrm{deep}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{deep}}}+{\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{deep}}}\right)\right]/{V}_{\mathrm{total}}.\end{array}$

Fully expanded, this yields:

$\begin{array}{ll}& {\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{global}}}+{\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{global}}}\approx \frac{{V}_{\mathrm{upper}}}{{V}_{\mathrm{total}}}\left[{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot \left({\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{upper}}}^{-}+{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{den},\mathrm{upper}}}^{-}\right)\\ & \phantom{\rule{1em}{0ex}}+{f}_{\mathrm{SO},\mathrm{upper}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{SO},\mathrm{upper}}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{SO},\mathrm{upper}}}}^{-}\right)\\ & \phantom{\rule{1em}{0ex}}+{f}_{\mathrm{NAtl},\mathrm{upper}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{NAtl},\mathrm{upper}}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NAtl},\mathrm{upper}}}}^{-}\right)\right]\\ & \phantom{\rule{1em}{0ex}}+\frac{{V}_{\mathrm{deep}}}{{V}_{\mathrm{total}}}\left[{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot \left({\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{deep}}}^{-}+{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{den},\mathrm{deep}}}^{-}\right)+{f}_{\mathrm{AABW}}\\ & \phantom{\rule{1em}{0ex}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{AABW}}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{AABW}}}}^{-}\right)\\ \text{(11)}& & \phantom{\rule{1em}{0ex}}+{f}_{\mathrm{NADW}}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{NADW}}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{\mathrm{NADW}}}}^{-}\right)\right],\end{array}$

which can be generalized for any number n of ventilation regions i as

$\begin{array}{ll}& {\mathrm{DIC}}_{{\mathrm{dis}}_{\mathrm{global}}}+{\mathrm{DIC}}_{{\mathrm{soft}}_{\mathrm{global}}}\approx {r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot \left({\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{global}}}^{-}+{\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{den},\mathrm{global}}}^{-}\right)\\ \text{(12)}& & \phantom{\rule{1em}{0ex}}+\sum _{i=\mathrm{1}}^{n}{f}_{i}\cdot \left({\mathrm{DIC}}_{{\mathrm{dis}}_{i}}-{r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{{\mathrm{pre}}_{i}}}^{-}\right).\end{array}$

Thus, total carbon storage as soft and disequilibrium carbon (i.e., everything other than DICsat and DICcarb) varies with the global nitrate inventory, corrected for accumulated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ loss to denitrification, and the difference between DICdis and ${r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}\mathrm{N}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{pre}}}^{-}$ in the polar oceans, modulated by their respective volume fractions. Although this nitrogen-based framework avoids the problem of C : P variability, it is not clear how large the effects of variable C : N might be in the real world. This could be a worthy topic for future exploration.

4 Conclusions

The conceptualization of ocean carbon storage as the sum of the saturation, soft tissue, carbonate, and disequilibrium components can greatly assist in enhancing mechanistic understanding . Our simulations indicate that the disequilibrium component may play a very important role, which has not been broadly appreciated. Changes in the physical climate states, as simulated by our model, tend to drive the soft tissue and the disequilibrium components in opposite directions. However, this is not necessarily true in the real ocean, given that the simulated anti-correlation is not mechanistically required, but instead arises from the fact that fAABW and age × export are anti-correlated in the simulations. On the contrary, the radiocarbon analysis of suggests that the glacial ocean age was significantly greater than in the corresponding simulations, implying that this anti-correlation does not hold in reality. Our iron fertilization experiments explore another aspect of this decoupling, in which age × export increases despite no change in fAABW. There is plenty of scope for these to have varied in additional ways in the real world, not captured by our simulations, including the idealized mechanisms explored by . Although the anti-correlation of DICdis and DICsoft in our simulations results in small overall changes, their magnitudes are sufficient that their total scope for change exceeds that required to explain the glacial/interglacial CO2 change.

Our results also show a surprising capacity for O2 disequilibrium to develop in a cold state. We suggest that this reflects a high sensitivity of O2 to sea ice when sea ice coverage reaches very high fractions. This generally unrecognized potential for sea ice coverage to cause large oxygen undersaturation may have contributed to very low O2 in the Southern Ocean during glacial periods, as suggested by foraminiferal I/Ca measurements .

The results presented here suggest that disequilibrium carbon should be considered as a major component of ocean carbon storage, linked to ocean circulation and biological export in non-linear and interdependent ways. Despite these nonlinearities, the simulations suggest that the resulting global carbon storage can be well approximated by simple relationships. We propose one such relationship, including the global nitrate inventory, and the DICdis and preformed ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in ocean ventilation regions (Eq. 12). As this study represents the result of a single model, prone to bias, it would be very useful to test our results using other GCMs including additional biogeochemical complexity. It would also be useful to consider how disequilibrium carbon can change under future pCO2 levels, including developing observational constraints on its past and present magnitude, and exploring the degree to which inter-model variations in DICdis may contribute to uncertainty in climate projections.

Code availability
Code availability.

All model run scripts, code, and simulation output are freely available from the authors, or can be downloaded from https://earthsystemdynamics.org/cm2mc-simulation-library (Galbraith, 2018).

Appendix A: DIC decomposition

DIC is treated as the sum of four components:

$\begin{array}{}\text{(A1)}& \mathrm{DIC}={\mathrm{DIC}}_{\mathrm{sat}}+{\mathrm{DIC}}_{\mathrm{dis}}+{\mathrm{DIC}}_{\mathrm{soft}}+{\mathrm{DIC}}_{\mathrm{carb}}.\end{array}$

DICsat is the DIC at equilibrium with the atmosphere given the surface ocean temperature, salinity, and alkalinity, and the atmospheric pCO2 calculated following :

$\begin{array}{}\text{(A2)}& {\mathrm{DIC}}_{\mathrm{sat}}=f\left(T,S,{\mathrm{alk}}_{\mathrm{pre}},\mathit{\text{p}}{\mathrm{CO}}_{\mathrm{2}}\right).\end{array}$

In this model, DICsoft is proportional to the utilized O2, which is defined as the difference between preformed and total O2, where the ratio of remineralized C to utilized O2 (${r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}}$) is 106 : 150.

$\begin{array}{}\text{(A3)}& {\mathrm{DIC}}_{\mathrm{soft}}={r}_{\mathrm{C}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}}\cdot \left({\mathrm{O}}_{{\mathrm{2}}_{\mathrm{pre}}}-{\mathrm{O}}_{\mathrm{2}}\right).\end{array}$

DIC derived from CaCO3 dissolution is proportional to the change in alkalinity, correcting for the additional change in alkalinity due to hydrogen ion addition during organic matter remineralization.

$\begin{array}{}\text{(A4)}& {\mathrm{DIC}}_{\mathrm{carb}}=\mathrm{0.5}\cdot \left[\left(\mathrm{alk}-{\mathrm{alk}}_{\mathrm{pre}}\right)+{r}_{\mathrm{N}\phantom{\rule{0.125em}{0ex}}:\phantom{\rule{0.125em}{0ex}}{\mathrm{O}}_{\mathrm{2}}}\cdot \left({\mathrm{O}}_{{\mathrm{2}}_{\mathrm{pre}}}-{\mathrm{O}}_{\mathrm{2}}\right)\right]\end{array}$

Preformed alkalinity, defined as the total alkalinity at the surface and treated as a conservative tracer, is calculated within the model framework but was not written out during the model runs. Therefore, we have reconstructed this parameter a posteriori for each model year through multilinear regressions as a function of century-averaged salinity (S), temperature (T), and preformed O2, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, following the approach of .

$\begin{array}{ll}& {\mathrm{alk}}_{\mathrm{pre}}=\left({a}_{\mathrm{0}}+{a}_{\mathrm{1}}\cdot {S}^{\prime }+{a}_{\mathrm{2}}\cdot {T}^{\prime }+{a}_{\mathrm{3}}\cdot {\mathrm{O}}_{{\mathrm{2}}_{\mathrm{pre}}}\\ & \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}+{a}_{\mathrm{4}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{pre}}}^{-}+{a}_{\mathrm{5}}\cdot {\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}\right)\cdot \mathrm{NAtl}\\ & \phantom{\rule{1em}{0ex}}+\left({b}_{\mathrm{0}}+{b}_{\mathrm{1}}\cdot {S}^{\prime }+{b}_{\mathrm{2}}\cdot {T}^{\prime }+{b}_{\mathrm{3}}\cdot {\mathrm{O}}_{{\mathrm{2}}_{\mathrm{pre}}}\\ & \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}+{b}_{\mathrm{4}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{pre}}}^{-}+{b}_{\mathrm{5}}\cdot {\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}\right)\cdot \mathrm{SO}\\ & \phantom{\rule{1em}{0ex}}+\left({c}_{\mathrm{0}}+{c}_{\mathrm{1}}\cdot {S}^{\prime }+{c}_{\mathrm{2}}\cdot {T}^{\prime }+{c}_{\mathrm{3}}\cdot {\mathrm{O}}_{{\mathrm{2}}_{\mathrm{pre}}}\\ \text{(A5)}& & \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}+{c}_{\mathrm{4}}\cdot {\mathrm{NO}}_{{\mathrm{3}}_{\mathrm{pre}}}^{-}+{c}_{\mathrm{5}}\cdot {\mathrm{PO}}_{{\mathrm{4}}_{\mathrm{pre}}}^{\mathrm{3}-}\right)\cdot \left(\mathrm{1}-\mathrm{SO}-\mathrm{NAtl}\right),\end{array}$

where ${S}^{\prime }=S-\mathrm{35}$, ${T}^{\prime }=T-\mathrm{20}$C, the ai are determined by a regression in the surface North Atlantic, the bi for the SO, and the ci using the model output elsewhere in the surface. The tracers SO and NAtl are set to 1 in the surface Southern Ocean (south of 30 S) and the North Atlantic (north of 30 N), respectively, and are conservatively mixed into the ocean interior. This parametrization induces an uncertainty on the order of 1 µmol kg−1 in globally averaged DICdis (see Fig. A1). As discussed above, however, this is small compared to the signal seen over all simulations.

Finally, DICdis has been back-calculated from the model output as a residual.

Figure A1Shown is the difference between the exact DICdis surface field in µmol kg−1, where DICsat has been calculated using the surface alkalinity ($\mathrm{alk}\left[z=\mathrm{0}\right]={\mathrm{alk}}_{\mathrm{pre}}\left[z=\mathrm{0}\right]$) and ${\mathrm{DIC}}_{\mathrm{soft}}\left[z=\mathrm{0}\right]={\mathrm{DIC}}_{\mathrm{carb}}\left[z=\mathrm{0}\right]=\mathrm{0}$. Differences are shown for (a) cold world; (b) moderate world; (c) hot world.

Author contributions
Author contributions.

EDG conducted the model simulations, SE performed the analysis, and both contributed to writing the manuscript.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The authors would like to thank Raffaele Bernardello for very helpful discussion. Sarah Eggleston was funded by a fellowship from the Swiss National Science Foundation. Eric D. Galbraith acknowledges computing support from the Canadian Foundation for Innovation and Compute Canada, and financial support from the Spanish Ministry of Economy and Competitiveness, through the María de Maeztu Programme for Centres/Units of Excellence in R&D (MDM-2015-0552).

Edited by: Fortunat Joos
Reviewed by: two anonymous referees

References

Albani, S., Mahowald, N. M., Delmonte, B., Maggi, V., and Winckler, G.: Comparing modele and observed changes in mineral dust transport and deposition to Antarctica between the Last Glacial Maximum and current climates, Clim. Dynam., 38, 1731–1755, https://doi.org/10.1007/s00382-011-1139-5, 2012. a

Albani, S., Mahowald, N. M., Murphy, L. N., Raiswell, R., Moore, J. K., Anderson, R. F., McGee, D., Bradtmiller, L. I., Delmonte, B., Hesse, P. P., and Mayewski, P. A.: Paleodust variability since the Last Glacial Maximum and implications for iron inputs to the ocean, Geophys. Res. Lett., 43, 3944–3954, https://doi.org/10.1002/2016GL067911, 2016. a

Bernardello, R., Marinov, I., Palter, J. B., Sarmiento, J. L., Galbraith, E. D., and Slater, R. D.: Response of the Ocean Natural Carbon Storage to Projected Twenty-First-Century Climate Change, J. Climate, 27, 2033–2053, https://doi.org/10.1175/jcli-d-13-00343.1, 2014. a, b, c

Broecker, W. S. and Peng, T.-H.: Gas exchange rates between air and sea, Tellus, 26, 21–35, 1974. a

Broecker, W. S., Takahashi, T., and Takahashi, T.: Sources and flow patterns of deep-ocean waters as deduced from potential temperature, salinity, and initial phosphate concentration, J. Geophys. Res., 90, 6925–6939, https://doi.org/10.1029/JC090iC04p06925, 1985. a

Duteil, O., Koeve, W., Oschlies, A., Bianchi, D., Galbraith, E., Kriest, I., and Matear, R.: A novel estimate of ocean oxygen utilisation points to a reduced rate of respiration in the ocean interior, Biogeosciences, 10, 7723–7738, https://doi.org/10.5194/bg-10-7723-2013, 2013. a, b

François, R., Altabet, M. A., Yu, E.-F., Sigman, D. M., Bacon, M. P., Frank, M., Bohrmann, G., Bareille, G., and Labeyrie, L. D.: Contribution of Southern Ocean surface-water stratification to low atmospheric CO2 concentrations during the last glacial period, Nature, 389, 929–935, https://doi.org/10.1038/40073, 1997. a

Friedlingstein, P., Meinshausen, M., Arora, V. K., Jones, C. D., Anav, A., Liddicoat, S. K., and Knutti, R.: Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks, J. Climate, 27, 511–526, 2014. a

Galbraith, E. D.: Integrated Earth System Dynamics: CM2Mc Simulation Library, available at: https://earthsystemdynamics.org/cm2mc-simulation-library, last access: 11 June 2018.

Galbraith, E. and de Lavergne, C.: Response of a comprehensive climate model to a broad range of external forcings: relevance for deep ocean ventilation and the development of late Cenozoic ice ages, Clim. Dynam., https://doi.org/10.1007/s00382-018-4157-8, online first, 2018. a, b, c, d, e, f, g, h, i

Galbraith, E. D. and Jaccard, S. L.: Deglacial weakening of the oceanic soft tissue pump: global constraints from sedimentary nitrogen isotopes and oxygenation proxies, Quaternary Sci. Rev., 109, 38–48, https://doi.org/10.1016/j.quascirev.2014.11.012, 2015. a

Galbraith, E. D. and Martiny, A. C.: A simple nutrient-dependent mechanism for predicting the stoichiometry of marine ecosystems, P. Natl. Acad. Sci. USA, 112, 8199–8204, https://doi.org/10.1073/pnas.1423917112, 2015. a, b, c

Galbraith, E. D., Gnanadesikan, A., Dunne, J. P., and Hiscock, M. R.: Regional impacts of iron-light colimitation in a global biogeochemical model, Biogeosciences, 7, 1043–1064, https://doi.org/10.5194/bg-7-1043-2010, 2010. a

Galbraith, E. D., Kwon, E. Y., Gnanadesikan, A., Rodgers, K. B., Griffies, S. M., Bianchi, D., Sarmiento, J. L., Dunne, J. P., Simeon, J., Slater, R. D., Wittenberg, A. T., and Held, I. M.: Climate Variability and Radiocarbon in the CM2Mc Earth System Model, J. Climate, 24, 4230–4254, https://doi.org/10.1175/2011jcli3919.1, 2011. a

Galbraith, E. D., Kwon, E. Y., Bianchi, D., Hain, M. P., and Sarmiento, J. L.: The impact of atmospheric pCO2 on carbon isotope ratios of the atmosphere and ocean, Global Biogeochem. Cy., 29, 307–324, https://doi.org/10.1002/2014GB004929, 2015. a

Gebbie, G. and Huybers, P.: How is the ocean filled?, Geophys. Res. Lett., 38, L06604, https://doi.org/10.1029/2011gl046769, 2011. a

Goodwin, P., Follows, M. J., and Williams, R. G.: Analytical relationships between atmospheric carbon dioxide, carbon emissions, and ocean processes, Global Biogeochem. Cy., 22, GB3030, https://doi.org/10.1029/2008gb003184, 2008. a, b

Gruber, N., Sarmiento, J. L., and Stocker, T. F.: An improved method for detecting anthropogenic CO2 in the oceans, Global Biogeochem. Cy., 10, 809–837, https://doi.org/10.1029/96gb01608, 1996. a

IPCC: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, United Kingdom, 2007. a

Ito, T. and Follows, M. J.: Preformed phosphate, soft tissue pump and atmospheric CO2, J. Mar. Res., 63, 813–839, https://doi.org/10.1357/0022240054663231, 2005. a, b, c

Ito, T. and Follows, M. J.: Air-sea disequilibrium of carbon dioxide enhances the biological carbon sequestration in the Southern Ocean, Global Biogeochem. Cy., 27, 1129–1138, https://doi.org/10.1002/2013gb004682, 2013. a, b, c, d, e, f, g

Ito, T., Follows, M. J., and Boyle, E. A.: Is AOU a good measure of respiration in the oceans?, Geophys. Res. Lett., 31, L17305, https://doi.org/10.1029/2004GL020900, 2004. a

Jaccard, S. L., Galbraith, E. D., Martínez-García, A., and Anderson, R. F.: Covariation of deep Southern Ocean oxygen and atmospheric CO2 through the last ice age, Nature, 530, 207–210, https://doi.org/10.1038/nature16514, 2016. a, b

Kohfeld, K. E. and Ridgwell, A.: Glacial-interglacial variability in atmospheric CO2, in: Surface Ocean-Lower Atmosphere Processes, edited by: Quéré, C. L. and Saltzman, E. S., Geophy. Monog. Series, 251–286, https://doi.org/10.1029/2008gm000845, 2009. a

Köhler, P. and Fischer, H.: Simulating low frequency changes in atmospheric CO2 during the last 740 000 years, Clim. Past, 2, 57–78, https://doi.org/10.5194/cp-2-57-2006, 2006. a

Kwon, E. Y., Primeau, F., and Sarmiento, J. L.: The impact of remineralization depth on the air-sea carbon balance, Nat. Geosci., 2, 630–635, https://doi.org/10.1038/NGEO612, 2009. a

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, https://doi.org/10.1051/0004-6361:20041335, 2004. a

Le Quéré, C., Rödenbeck, C., Buitenhuis, E. T., Conway, T. J., Langenfelds, R., Gomez, A., Labuschagne, C., Ramonet, M., Nakazawa, T., Metzl, N., Gillett, N., and Heimann, M.: Saturation of the Southern Ocean CO2 Sink Due to Recent Climate Change, Science, 316, 1735–1738, https://doi.org/10.1126/science.1136188, 2007. a

Lu, Z., Hoogakker, B. A. A., Hillenbrand, C.-D., Zhou, X., Thomas, E., Gutchess, K. M., Lu, W., Jones, L., and Rickaby, R. E. M.: Oxygen depletion recorded in upper waters of the glacial Southern Ocean, Nat. Commun., 7, 11146, https://doi.org/10.1038/ncomms11146, 2015. a, b

Mahowald, N. M., Muhs, D. R., Levis, S., Rasch, P. J., Yoshioka, M., Zender, C. S., and Luo, C.: Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates, J. Geophys. Res.-Atmos., 111, d10202, https://doi.org/10.1029/2005JD006653, 2006. a

Marinov, I., Follows, M., Gnanadesikan, A., Sarmiento, J. L., and Slater, R. D.: How does ocean biology affect atmospheric pCO2? Theory and models, J. Geophys. Res., 113, C07032, https://doi.org/10.1029/2007jc004598, 2008a. a, b, c, d

Marinov, I., Gnanadesikan, A., Sarmiento, J. L., Toggweiler, J. R., Follows, M., and Mignone, B. K.: Impact of oceanic circulation on biological carbon storage in the ocean and atmospheric pCO2, Global Biogeochem. Cy., 22, GB3007, https://doi.org/10.1029/2007gb002958, 2008b. a

Martin, J. H.: Glacial-interglacial CO2 change: The Iron Hypothesis, Paleoceanography, 5, 1–13, https://doi.org/10.1029/pa005i001p00001, 1990. a, b

Martin, J. H., Knauer, G. A., Karl, D. M., and Broenkow, W. W.: VERTEX: carbon cycling in the northeast Pacific, Deep-Sea Res. Pt. A, 34, 267–285, https://doi.org/10.1016/0198-0149(87)90086-0, 1987. a

Martiny, A. C., Vrugt, J. A., Primeau, F. W., and Lomas, M. W.: Regional variation in the particulate organic carbon to nitrogen ratio in the surface ocean, Global Biogeochem. Cy., 27, 723–731, https://doi.org/10.1002/gbc.20061, 2013. a

Matsumoto, K., Hashioka, T., and Yamanaka, Y.: Effect of temperature-dependent organic carbon decay on atmospheric pCO2, J. Geophys. Res., 112, G02007, https://doi.org/10.1029/2006JG000187, 2007. a, b

Moore, C. M., Mills, M. M., Arrigo, K. R., Berman-Frank, I., Bopp, L., Boyd, P. W., Galbraith, E. D., Geider, R. J., Guieu, C., Jaccard, S. L., Jickells, T. D., La Roche, J., Lenton, T. M., Mahowald, N. M., Marañón, E., Marinov, I., Moore, J. K., Nakatsuka, T., Oschilles, A., Saito, M. A., Thingstad, T. F., Tsuda, A., and Ulloa, O.: Processes and patterns of oceanic nutrient limitation, Nat. Geosci., 6, 701–710, https://doi.org/10.1038/ngeo1765, 2013. a

Myhre, G., Highwood, E. J., Shine, K. P., and Stordal, F.: New estimates of radiative forcing due to well mixed greenhouse gases, Geophys. Res. Lett., 25, 2715–2718, https://doi.org/10.1029/98GL01908, 1998. a

Nickelsen, L. and Oschlies, A.: Enhanced sensitivity of oceanic CO2 uptake to dust deposition by iron-light colimitation, Geophys. Res. Lett., 42, 492–499, https://doi.org/10.1002/2014gl062969, 2015. a

Ödalen, M., Nycander, J., Oliver, K. I. C., Brodeau, L., and Ridgwell, A.: The influence of the ocean circulation state on ocean carbon storage and CO2 drawdown potential in an Earth system model, Biogeosciences, 15, 1367–1393, https://doi.org/10.5194/bg-15-1367-2018, 2018. a, b, c, d, e, f

Parekh, P., Dutkiewicz, S., Follows, M. J., and Ito, T.: Atmospheric carbon dioxide in a less dusty world, Geophys. Res. Lett., 33, L03610, https://doi.org/10.1029/2005gl025098, 2006. a

Roth, R., Ritz, S. P., and Joos, F.: Burial-nutrient feedbacks amplify the sensitivity of atmospheric carbon dioxide to changes in organic matter remineralisation, Earth Syst. Dynam., 5, 321–343, https://doi.org/10.5194/esd-5-321-2014, 2014. a

Russell, J. L. and Dickson, A. G.: Variability in oxygen and nutrients in South Pacific Antarctic Intermediate Water, Global Biogeochem. Cy., 17, 1033, https://doi.org/10.1029/2000GB001317, 2003. a

Schmittner, A. and Galbraith, E. D.: Glacial greenhouse-gas fluctuations controlled by ocean circulation changes, Nature, 456, 373–376, https://doi.org/10.1038/nature07531, 2008. a

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

Sigman, D. M., Hain, M. P., and Haug, G. H.: The polar ocean and glacial cycles in atmospheric CO2 concentration, Nature, 466, 47–55, https://doi.org/10.1038/nature09149, 2010. a

Skinner, L. C.: Glacial-interglacial atmospheric CO2 change: a possible “standing volume” effect on deep-ocean carbon sequestration, Clim. Past, 5, 537–550, https://doi.org/10.5194/cp-5-537-2009, 2009. a

Stephens, B. B. and Keeling, R. F.: The influence of Antarctic sea ice on glacial-interglacial CO2 variations, Nature, 404, 171–174, https://doi.org/10.1038/35004556, 2000. a

Tagliabue, A., Aumont, O., DeAth, R., Dunne, J. P., Dutkiewicz, S., Galbraith, E., Misumi, K., Moore, J. K., Ridgwell, A., Sherman, E., Stock, C., Vichi, M., Völker, C., and Yool, A.: How well do global ocean biogeochemistry models simulate dissolved iron distributions?, Global Biogeochem. Cy., 30, 149–174, https://doi.org/10.1002/2015GB005289, 2016. a, b

Takahashi, T. Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Watson, A., Bakker, D. C. E., Schuster, U., Metzl, N., Yoshikawa-Inoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Sabine, C., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Körtzinger, A., Steinhoff, T., Hoppema, M., de Baar, H. J. W., Wong, C. S., Delille, B., and Bates, N. R.: Climatological mean and decadal changes in surface ocean pCO2, and net sea-air CO2 flux over the global oceans, Deep-Sea Res. Pt. II, 56, 554–577, https://doi.org/10.1016/j.dsr2.2008.12.009, 2009. a

Toggweiler, J. R., Murnane, R., Carson, S., Gnanadesikan, A., and Sarmiento, J. L.: Representation of the carbon cycle in box models and GCMs: 2. Organic pump, Global Biogeochem. Cy., 17, 1027, https://doi.org/10.1029/2001gb001841, 2003. a, b

Tschumi, T., Joos, F., Gehlen, M., and Heinze, C.: Deep ocean ventilation, carbon isotopes, marine sedimentation and the deglacial CO2 rise, Clim. Past, 7, 771–800, https://doi.org/10.5194/cp-7-771-2011, 2011. a

Volk, T. and Hoffert, M. I.: The Carbon Cycle and Atmospheric CO2: Natural Variations Archean to Present, chap. Ocean Carbon Pumps: Analysis of Relative Strengths and Efficiencies in Ocean-Driven Atmospheric CO2 Changes, American Geophysical Union, 99–110, https://doi.org/10.1029/GM032p0099, 1985.  a

Watson, A. J., Vallis, G. K., and Nikurashin, M.: Southern Ocean buoyancy forcing of ocean ventilation and glacial atmospheric CO2, Nat. Geosci., 8, 861–865, https://doi.org/10.1038/ngeo2538, 2015. a

Williams, R. G. and Follows, M. J.: Ocean Dynamics and the Carbon Cycle: Principles and Mechanisms, Cambridge University Press, 2011. a, b

Zeebe, R. E. and Wolf-Gladrow, D.: CO2 in Seawater: Equilibrium, Kinetics, Isotopes, vol. 65, Elsevier Oceanography Series, Amsterdam, 2001. a, b, c