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 Galbraith et al. (2011) and modified as described by Galbraith and de Lavergne (2018). 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; Galbraith et al., 2010). Unlike BLINGv0, BLINGv1.5 allows for variable P : C stoichiometry using the “line of frugality” (Galbraith and Martiny, 2015) and calculates the mass balance of phytoplankton in order to prevent unrealistic bloom magnitudes at high latitudes, reducing the magnitude of disequilibrium O₂, which was very high in BLINGv0 (Duteil et al., 2013; Tagliabue et al., 2016). 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 Galbraith and de Lavergne (2018). A control run was conducted with a radiative forcing equivalent to 270 ppm atmospheric pCO₂ 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 (Laskar et al., 2004), while eccentricity was held constant at 0.03. A range of greenhouse radiative forcings was imposed equivalent to pCO₂ 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⁻², respectively (Myhre et al., 1998).

The biogeochemical component of the model calculates air–sea carbon fluxes using a fixed atmospheric pCO₂ 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 pCO₂ for the carbon cycle means that the DIC_sat is not consistent with the pCO₂ used for the radiative forcing, so that changes in DIC_sat caused by a given pCO₂ 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 pCO₂; 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 Mahowald et al. (2006) as modified by Nickelsen and Oschlies (2015) 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 (Matsumoto et al., 2007); 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 Galbraith and de Lavergne (2018). Runs 43 and 44 and identical to 41 and 42 but the remineralization rate of sinking organic matter is reduced by 25 %.

Download Print Version | Download XLSX

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⁻², 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⁻², 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⁻¹ 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 O₂ 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, DIC_pre, DIC_sat, O₂, and O_2pre, since DIC${}_{\mathrm{dis}}={\mathrm{DIC}}_{\mathrm{pre}}-{\mathrm{DIC}}_{\mathrm{sat}}$, DIC${}_{\mathrm{soft}}=({\mathrm{O}}_{\mathrm{2}\mathrm{pre}}-{\mathrm{O}}_{\mathrm{2}})\cdot r_{\mathrm{C}:{\mathrm{O}}_{\mathrm{2}}}$, DIC${}_{\mathrm{carb}}=\mathrm{DIC}-{\mathrm{DIC}}_{\mathrm{pre}}-{\mathrm{DIC}}_{\mathrm{soft}}$, and O_2dis = O_2pre − O_2sat (O_2sat can be accurately calculated from the conservative tracers temperature and salinity). In the absence of a DIC_sat tracer, it can be estimated from alk_pre, temperature, and salinity. For this large suite of simulations, we only had DIC, alk, O₂, and O_2pre available. Thus, although we could calculate O_2dis and DIC_soft directly, it was necessary to use an indirect method to calculate the other three carbon tracers. Following Bernardello et al. (2014), we first estimate alk_pre using a regression, then calculate DIC_sat using alk_pre, temperature, salinity, and the known atmospheric pCO₂ of 270 ppm. DIC_carb is calculated as $\mathrm{0.5}\cdot (\mathrm{alk}-{\mathrm{alk}}_{\mathrm{pre}})$, and DIC_dis is then calculated as a residual. For more details and an estimate of the error in the method, see the Appendix.

3.1 General climate response to forcings

Differences in the general ocean state among the model simulations are described in detail by Galbraith and de Lavergne (2018). 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 Galbraith and de Lavergne (2018), 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.

Download

3.2 General changes in DIC

Total DIC generally decreases from cold to warm simulations, under the constant pCO₂ of 270 ppm used for air–sea exchange. Changes in DIC_sat drive the largest portion of this trend, decreasing approximately linearly with surface air temperature due to the temperature-dependence of CO₂ solubility, resulting in a difference of 50 µmol kg⁻¹ over this range (see Fig. 3). Because of the constant pCO₂ of 270 ppmv in the biogeochemical module, the simulated range of DIC_sat is larger than would occur if the prescribed biogeochemical pCO₂ were the same as that used to produce the low SAT (Ödalen et al., 2018). DIC_carb is relatively small in magnitude and generally increases with SAT, but has a standard deviation of only 4 µmol kg⁻¹ 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.

Download

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

Simulated changes in DIC_dis are of the same magnitude as the DIC_soft changes, to which much greater attention has been paid. For a global average buffer factor between 8 and 14 (Zeebe and Wolf-Gladrow, 2001), a rough, back-of-the-envelope calculation shows that a 1 µmol kg⁻¹ change in DIC corresponds to a 0.9–1.6 ppm change in atmospheric pCO₂ based on a DIC concentration of 2300 µmol kg⁻¹ and pCO₂ of 270 ppm. Thus, the increase in the global average DIC_dis in these simulations could have contributed more than a 40 ppm change in the atmospheric pCO₂ stored in the ocean during the glacial compared to today. It is important to recognize that the drawdown of CO₂ by disequilibrium storage would have resulted in a decrease of DIC_sat, given the dependence of the saturation concentration on pCO₂, so this estimate should not be interpreted as a straightforward atmospheric pCO₂ 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 DIC_soft

The biogeochemical model used here is relatively complex, with limitation by three nutrients (N, P, and Fe), denitrification, and N₂ 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 DIC_soft, has highly predictable behavior. As shown in Fig. 4, the global DIC_soft 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 DIC_soft is not as well correlated with either of these parameters separately as it is with their product “age × export”.

Figure 4Globally averaged DIC_soft vs. global export, ideal age, and age × export. Globally averaged DIC_soft 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.

Download

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 (Martin et al., 1987). Kwon et al. (2009) 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 (Matsumoto et al., 2007) or ecosystem structure. As a result, the relationship between DIC_soft 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 DIC_soft between two states can be given by a simple linear regression:

or in terms of pCO₂:

Note that m₂ is a function of the buffer factor and the climate state (atmospheric pCO₂ and DIC). Based on the results here, m₁=0.036 and $m_{\mathrm{2}}=\mathrm{0.065},\mathrm{0.042},\mathrm{0.029}$ for modern (405 ppm pCO₂), pre-industrial (270 ppm) and glacial (180 ppm) conditions, respectively. Note that we have not varied pCO₂ 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 CO₂. Note that, as for the disequilibrium estimate above, the soft tissue pump CO₂ 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 DIC_soft 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 DIC_soft concentrations during the glacial then during the Holocene (Galbraith and Jaccard, 2015). In contrast, the model suggests that greater ocean ventilation rates in the glacial state (Fig. 2a) would have led to reduced global DIC_soft. As discussed by Galbraith and de Lavergne (2018), 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⁻¹ DIC_soft were stored in the glacial ocean. This would bring the simulated glacial DIC_soft close to, but still less than, the simulated pre-industrial value.

We propose that the apparent remaining shortfall in simulated glacial DIC_soft 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 DIC_soft 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 DIC_dis

The ocean basins below 1 km depth are largely filled by surface waters subducted to depth in regions of deepwater formation (Gebbie and Huybers, 2011). 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 DIC_dis concentration can be approximated by a simple mass balance:

Here, f_AABW represents the fraction of deepwater originating in the SO, and DIC${}_{{\mathrm{dis}}_{\mathrm{AABW}}}$ and DIC${}_{{\mathrm{dis}}_{\mathrm{NADW}}}$ represent the DIC_dis concentrations at the sites of deepwater formation (see Fig. 5). North Atlantic deep water forms with negative DIC_dis, reflecting surface undersaturation, while the Southern Ocean is supersaturated (DIC_dis>0). These opposing tendencies between NADW and AABW cause a partial cancellation of DIC_dis when globally averaged, which makes the disequilibrium component small in the modern ocean. Theoretically, the simulated DIC_dis could change either due to changes in f_AABW or the end-member compositions. Although the exact values of DIC_dis 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 f_AABW and f_NADW, so that deep DIC_dis is strongly controlled by the global balance of AABW vs. NADW in each simulation (see Fig. 6). Global DIC_dis becomes much larger when f_AABW is larger, similar to the dynamic evoked by Skinner (2009). This is also illustrated by the depth transects of DIC_dis in Fig. 7.

Figure 5(a) Global average fraction of northern- (reddish colors) and southern-sourced (bluish colors) water. (b) Annual average values of DIC_dis 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 60^∘S, are not shown. The size of the symbols corresponds to the SAT.

Download

Figure 6Global average DIC_dis 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.

Download

Figure 7DIC_dis (µmol kg⁻¹) 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).

Download

We estimated the concentration of DIC_dis in the regions of AABW and NADW formation, shown in Fig. 5b. The end members vary less significantly than f_AABW 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 DIC_dis (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 (Galbraith and de Lavergne, 2018), further increasing DIC_dis in glacial-like conditions. Sea ice in the Southern Ocean would be expected to exert a further control over DIC_dis, 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 DIC_dis, ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ is defined as the departure from equilibrium of O₂ in the surface ocean with respect to the atmosphere and is advected into the ocean as a conservative tracer. Unlike C, O₂ does not react with seawater, and is present only as dissolved diatomic oxygen. Thus, O₂ has a much shorter timescale of exchange at the ocean–atmosphere interface, equilibrating 1 order of magnitude faster than CO₂. As a result, it is not sensitive to sea ice as long as there remains a fair degree of open water (Stephens and Keeling, 2000). But as the sea ice concentration approaches complete coverage, O₂ equilibration rapidly becomes quite sensitive to sea ice. If there is a significant undersaturation of O₂ in upwelling waters, the disequilibrium can become quite large (Fig. 8).

Figure 8${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ (µmol kg⁻¹) 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).

Download

In the model simulations, the ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ in the Southern Ocean becomes as large as −100 µmol kg⁻¹ in the coldest states (note that DIC_dis and ${\mathrm{O}}_{{\mathrm{2}}_{\mathrm{dis}}}$ are often anti-correlated). Because the disequilibrium depends on the O₂ depletion of waters upwelling at the Southern Ocean surface, this could potentially be even larger if upwelling waters had lower O₂. 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 O₂ concentrations reconstructed for the glacial (Jaccard et al., 2016; Lu et al., 2015).

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 DIC_soft increase when iron deposition is increased. However, the DIC_soft increase is significantly lower in the well-ventilated GL simulations (2.9 µmol kg⁻¹) compared to the PI simulation (7.3 µmol kg⁻¹). This difference is qualitatively in accordance with the age × export relationship (Fig. 4), though with a smaller increase of DIC_soft than would be expected from the export increase, compared to the broad spectrum of climate-driven changes. This reduced sensitivity of DIC_soft 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 DIC_soft 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.

Download

The iron addition also causes an increase of DIC_dis of approximately equal magnitude to DIC_soft 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 DIC_dis, as remineralized DIC_soft 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 DIC_dis, rather than DIC_soft as might be assumed. Thus, just as the glacial state has a larger general proportion of DIC_dis compared to DIC_soft, the iron addition under the glacial state produces a larger fraction of DIC_dis relative to DIC_soft. The experiment in which the remineralization rate was reduced by 25 % indicates that the effects of iron fertilization alone on both DIC_soft and DIC_dis 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 DIC_dis vs. DIC_soft, in response to iron addition, can be quantified by the global ratio ΔDIC_dis ∕ ΔDIC_soft. 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 Parekh et al. (2006) found ΔDIC_dis ∕ ΔDIC_soft of 2 in response to iron fertilization, using a modern ocean circulation, as analyzed by Ito and Follows (2013). We also note that the quantitative values of DIC_soft and DIC_dis resulting from the altered iron flux should be taken with a grain of salt, given the very large uncertainty in iron cycling models (Tagliabue et al., 2016).

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 DIC_soft, due to higher biological export, and DIC_dis, because of the upwelling of C-rich water resulting from higher remineralization, the effect discussed by Ito and Follows (2013). However, under a high ventilation state there is only a small increase in DIC_soft 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 DIC_dis only. Thus, the carbon storage resulting from Southern Ocean iron addition could actually be dominated by DIC_dis under some climate states, so that the overall impact may be significantly larger than would be predicted from DIC_soft and/or O₂ utilization.

3.7 A unified framework for DIC_dis 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 (Sigman and Boyle, 2000; Ito and Follows, 2005; Marinov et al., 2008a) as well as changes in the ventilation fractions of AABW and NADW, given their very different preformed nutrient concentrations (Schmittner and Galbraith, 2008; Marinov et al., 2008b). 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 DIC_soft, inasmuch as the global N inventory is fixed (Galbraith and Martiny, 2015).

Figure 10Simple estimation of global DIC_dis and DIC_soft from water mass characteristics. Following Eq. (12) for (a) DIC_soft and (b) DIC_dis separately, this model shows that the sum of the global average DIC_soft and DIC_dis 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 60^∘S, 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).

Download

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

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

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.

Combining with Eq. (3),

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

Fully expanded, this yields:

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

Thus, total carbon storage as soft and disequilibrium carbon (i.e., everything other than DIC_sat and DIC_carb) varies with the global nitrate inventory, corrected for accumulated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ loss to denitrification, and the difference between DIC_dis and $r_{\mathrm{C}:\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.