<figure><img src="image_1.png" /></figure>

HAL
open science

# The formation of the ocean's anthropogenic carbon reservoir

Daniele Iudicone, Keith B. Rodgers, Yves Plancherel, Olivier Aumont,
Takamitsu Ito, Robert M. Key, Gurvan Madec, Masao M Ishii

► To cite this version:

Daniele Iudicone, Keith B. Rodgers, Yves Plancherel, Olivier Aumont, Takamitsu Ito, et al. The formation of the ocean's anthropogenic carbon reservoir. Scientific Reports, 2016, 6, pp.35473. 10.1038/srep35473 . hal-01404420

HAL Id: hal-01404420

https://hal.sorbonne-universite.fr/hal-01404420

Submitted on 28 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

# The formation of the ocean's anthropogenic carbon reservoir

Daniele Iudicone¹, Keith B. Rodgers², Yves Plancherel³, Olivier Aumont⁴, Takamitsu Ito⁵, Robert M. Key², Gurvan Madec⁴,⁶ & Masao Ishii⁷

Received: 21 August 2015

Accepted: 28 September 2016

Published: 03 November 2016

The shallow overturning circulation of the oceans transports heat from the tropics to the mid-latitudes. This overturning also influences the uptake and storage of anthropogenic carbon (Cant). We demonstrate this by quantifying the relative importance of ocean thermodynamics, circulation and biogeochemistry in a global biochemistry and circulation model. Almost 2/3 of the Cant ocean uptake enters via gas exchange in waters that are lighter than the base of the ventilated thermocline. However, almost 2/3 of the excess Cant is stored below the thermocline. Our analysis shows that subtropical waters are a dominant component in the formation of subpolar waters and that these water masses essentially form a common Cant reservoir. This new method developed and presented here is intrinsically Lagrangian, as it by construction only considers the velocity or transport of waters across isopycnals. More generally, our approach provides an integral framework for linking ocean thermodynamics with biogeochemistry.

High quality measurements demonstrate that the oceans have absorbed about a third of anthropogenic CO2 (Cant) emissions¹,². However our understanding of the mechanisms governing this uptake are quite elementary. The air-sea exchange of Cant should be expected to play a first-order role in determining the rate of exchange of Cant across the base of the ocean's planetary boundary layer or mixed layer. This mainly occurs in the Equatorial region and mid-to-high latitudes³. Previous efforts such as model-data synthesis using inverse techniques³⁻⁷ have investigated meridional transport of water and Cant in an Eulerian framework. This work tended to emphasize the role of Cant uptake associated with the Southern Ocean divergence. Recent attention has been devoted to the relationship between anthropogenic heat and carbon uptake by the ocean⁸,⁹. None of this research to date has articulated a mechanistic framework linking ocean uptake of carbon to thermodynamic processes associated with the ocean overturning. Here we develop such a framework based on previous efforts to understand oceanic carbon uptake pathways¹⁰⁻¹³, with special focus on the shallow overturning circulation (SOC).

The SOC includes overturning of tropical waters (TW), subtropical mode waters (STMW) and subpolar mode waters (SPMW), primarily Subantarctic Mode Water (SAMW)¹⁴ (see Fig. 1). Observations indicate that the waters with densities characteristic of the SOC ($\sigma$ < 27.0) and intermediate waters (IW; 27.0 $\le$ $\sigma$ < 27.5) contain as much as 63–83% of the global Cant inventory, yet occupy only 27.1% of the global ocean volume (Fig. 2 and Supplementary Fig. 1, Table 1). The range depends on the method used to estimate Cant concentrations at the time of WOCE¹,¹⁵,¹⁶. To understand the higher efficiency of these water masses in retaining Cant we must first understand the formation mechanism(s) of this reservoir. This exercise is particularly important considering the much shorter times scales of re-emergence into the surface layer (thereby inhibiting further Cant uptake) of these water masses as compared to deeper water masses. In addition, they contribute a much lower fraction to the total Cant inventory relative to what would be expected given their the large surface area of their outcrop regions (94.8%; Fig. 2, Table 1). The apparent mismatch between a large area available for air-sea exchange and a relatively small interior inventory for Cant suggests that (a) air-sea Cant fluxes are small in the tropics relative to high latitudes, (b) Cant is exported from the lighter layers to denser water masses, or both.

We present here a novel suite of analysis tools that aim to deconvolve the roles of gas exchange, ocean interior diapycnal transports, and ocean interior diapycnal diffusive exchanges in prescribing how the ocean interior Cant reservoir is formed. We identify a central role for poleward ocean transport of Cant, related to the poleward transport of heat associated with the SOC. This transport sustains a strong exchange of Cant between

¹Stazione Zoologica Anton Dohrn, I-80121 Naples, Italy. ²AOS Program, Princeton University, Princeton, NJ 08544, USA. ³Department of Earth Sciences, University of Oxford, Oxford OX1 3AN, UK. ⁴Sorbonne Universités, UPMC, Univ Paris 06-CNRS-IRD-MNHN, LOCEAN/IPSL, Paris, France. ⁵School of Earth and Atmospheric Sciences, Georgia Institute of Technology, Atlanta, Georgia, USA. ⁶National Oceanography Centre, Southampton SO14 3ZH, UK. ⁷Oceanography and Geochemistry Department, Meteorological Research Institute, JMA, Tsukuba-city, Ibaraki 305-0052, Japan. Correspondence and requests for materials should be addressed to D.I. (email: iudicone@szn.it)
<figure><img src="image_3.png" /></figure>

**Figure 1.** Schematic of the ocean overturning circulation (scheme adapted from ref. 35). The focus of this work is the shallow overturning circulation (SOC), spanning tropical thermocline waters (TW), subtropical mode waters (STMW), and Subantarctic Mode Water (SAMW), with the latter two water masses coupled thermodynamically. Also important in the diagram is the deep overturning, comprising North Atlantic Deep Water (NADW) and Antarctic Bottom Water (AABW). Antarctic Intermediate Water (AAIW) lies between the upper and deeper overturning structures. The Ekman layer (EL) has a significant poleward flow in the tropics which is inverted for higher latitudes. At mid-latitudes the poleward export of water is dominated by the Western Boundary Currents (WBC) and the eddy transport. In the interior, the overturning is delimited by the shallow, ventilated thermocline and the interior pycnocline, mainly set by the Southern Ocean wind-driven flow and by eddies in the Antarctic Circumpolar Current (ACC)³⁵.

<figure><img src="image_4.png" /></figure>

**Figure 2.** The global oceans density layering: outcrop areas, volumes and carbon inventories. Global ocean outcrop surface (grey), volume (black) in density space from monthly density fields from the World Ocean Atlas 2009⁴¹ and global C<sub>ant</sub> inventory in density space (magenta) from GLODAP in red (ΔC*¹). Areas under the curves have been normalised to unity. Important water mass bands are demarcated, distinguishing tropical waters (TW; σ<sub>0</sub> ≤ 24.50), subtropical mode waters (STMW; 24.50 ≤ σ<sub>0</sub> < 26.60), subpolar mode waters (SPMW; 26.60 ≤ σ<sub>0</sub> < 27.05), intermediate waters (IW; 27.05 ≤ σ<sub>0</sub> < 27.55), and deep waters (DW; σ<sub>0</sub> > 27.55).

STMW and SPMW, driven by near-surface diabatic processes. These diabatic processes and subsequent internal down-gradient diffusion of C<sub>ant</sub> play first-order roles in controlling the routing and partitioning of C<sub>ant</sub> in the ocean, and they contribute to the large subduction of C<sub>ant</sub> occurring over the Southern Ocean and North Atlantic¹.

## Basic requirements for global uptake and storage of C<sub>ant</sub>

Formation of the ocean's C<sub>ant</sub> reservoir is due to air-sea fluxes, ocean circulation and small-scale diffusive processes. Here we focus on those processes that maintain overturning and subduction of water masses, thereby contributing to the carbon inventory³. Discrete water masses, here defined as the water contained between isopycnal layers, are a macroscopic expression of oceanic adjustment to air-sea interaction in its various forms. The
<table><caption>**Table 1.** Water mass specific GLODAP ΔC\* and TTD-based C\_ant inventory estimates, volume inventory and outcrop area. Density values, outcrop areas and volumes from monthly WOA09. Winter outcrop areas were computed using the local maximal density over the year. TW = Tropical Waters, STMW = Subtropical Mode Waters, SPMW = Subpolar Mode Waters, IW = Intermediate Waters, DW = Deep Waters. See also Figs S1 and S2. \(1\)Excluding the Arctic and marginal seas. \(2\)[1], \[2]. \(3\)[4]. \(4\)The bottom of the MWML is defined by the deepest depth of the mixed layer found in the monthly climatology of [9]. \(5\)Fraction relative to whole ocean inventory.</caption><thead><tr><th></th><th></th><th>Total</th><th>TW σ<sub>0</sub> &lt; 24.5</th><th>STMW 24.5 &lt; σ<sub>0</sub> &lt; 26.6</th><th>SPWM 26.6 &lt; σ<sub>0</sub> &lt; 27.0</th><th>IW 27.0 &lt; σ<sub>0</sub> &lt; 27.5</th><th>DW 27.5 &lt; σ<sub>0</sub></th></tr></thead><tbody><tr><td colspan="8">(a) WHOLE OCEAN<sup>(1)</sup></td></tr><tr><td colspan="8">Absolute values</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[PgC] - ΔC<sup>*</sup><sup>(2)</sup></td><td>104,7</td><td>7,5</td><td>28,7</td><td>24,2</td><td>26,3</td><td>18,0</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[PgC] - TTD<sup>(3)</sup></td><td>140,6</td><td>8,8</td><td>29,5</td><td>23,3</td><td>27,5</td><td>51,4</td></tr><tr><td>Volume</td><td>[10<sup>18</sup> m<sup>3</sup>]</td><td>1312,3</td><td>15,2</td><td>64,8</td><td>75,1</td><td>200,0</td><td>957,1</td></tr><tr><td>Outcrop area</td><td>[10<sup>12</sup> m<sup>2</sup>]</td><td>348,0</td><td>164,7</td><td>120,9</td><td>22,9</td><td>35,7</td><td>3,7</td></tr><tr><td>Winter Outcrop Area (WOA09)</td><td>[10<sup>12</sup> m<sup>2</sup>]</td><td>348,5</td><td>129,4</td><td>129,1</td><td>34,9</td><td>37,1</td><td>18</td></tr><tr><td colspan="8">Relative fractions</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - ΔC<sup>*</sup></td><td>—</td><td>7,2</td><td>27,4</td><td>23,1</td><td>25,1</td><td>17,2</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - TTD</td><td>—</td><td>6,3</td><td>21,0</td><td>16,6</td><td>19,6</td><td>36,6</td></tr><tr><td>Volume</td><td>[%]</td><td>—</td><td>1,2</td><td>4,9</td><td>5,7</td><td>15,2</td><td>72,9</td></tr><tr><td>Outcrop area</td><td>[%]</td><td>—</td><td>47,3</td><td>34,7</td><td>6,6</td><td>10,3</td><td>1,1</td></tr><tr><td>Winter Outcrop area</td><td>[%]</td><td>—</td><td>37,1</td><td>37,0</td><td>10,0</td><td>10,6</td><td>5,2</td></tr><tr><td colspan="8">Cumulative fractions</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - ΔC<sup>*</sup></td><td>—</td><td>7,2</td><td>34,6</td><td>57,7</td><td>82,8</td><td>100</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - TTD</td><td>—</td><td>6,3</td><td>27,2</td><td>43,8</td><td>63,4</td><td>100</td></tr><tr><td>Volume</td><td>[%]</td><td>—</td><td>1,2</td><td>6,1</td><td>11,8</td><td>27,1</td><td>100</td></tr><tr><td>Outcrop area</td><td>[%]</td><td>—</td><td>47,3</td><td>82,1</td><td>88,6</td><td>98,9</td><td>100</td></tr><tr><td>Winter Outcrop area</td><td>[%]</td><td>—</td><td>37,1</td><td>74,2</td><td>84,2</td><td>94,8</td><td>100</td></tr><tr><td colspan="8">(b) MWML ONLY<sup>(1,4)</sup></td></tr><tr><td colspan="8">Absolute values</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[PgC] - ΔC<sup>*</sup></td><td>18,3</td><td>4,8</td><td>6,9</td><td>2,8</td><td>3,0</td><td>0,9</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[PgC] - TTD</td><td>19,3</td><td>5,7</td><td>7,0</td><td>2,6</td><td>2,9</td><td>1,0</td></tr><tr><td>Volume</td><td>[10<sup>18</sup> m<sup>3</sup>]</td><td>37,3</td><td>9,5</td><td>12,7</td><td>5,4</td><td>6,8</td><td>2,9</td></tr><tr><td colspan="8">Relative fractions</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - ΔC<sup>*</sup></td><td>17.5<sup>(5)</sup></td><td>26,2</td><td>37,7</td><td>15,3</td><td>16,4</td><td>4,9</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - TTD</td><td>13.7<sup>(5)</sup></td><td>29,5</td><td>36,3</td><td>13,5</td><td>15,0</td><td>5,2</td></tr><tr><td>Volume</td><td>[%]</td><td>2.8<sup>(5)</sup></td><td>25,5</td><td>34,0</td><td>14,5</td><td>18,2</td><td>7,8</td></tr><tr><td colspan="8">Cumulative fractions</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - ΔC<sup>*</sup></td><td>—</td><td>26,2</td><td>63,9</td><td>79,2</td><td>95,6</td><td>100</td></tr><tr><td>C<sub>ant</sub> inventory</td><td>[%] - TTD</td><td>—</td><td>29,5</td><td>65,8</td><td>79,3</td><td>94,3</td><td>100</td></tr><tr><td>Volume</td><td>[%]</td><td>—</td><td>25,5</td><td>59,5</td><td>74,0</td><td>92,2</td><td>100</td></tr></tbody></table>

quasi-adiabatic behaviour of the ocean interior makes isopycnal layers the preferred horizons for studying both physics and biogeochemistry.

The complexity of the processes setting the ocean circulation can be simplified by considering the volume budget of waters in density space, i.e. the integral framework developed by Walin¹⁷ and extended through numerous subsequent publications¹⁸⁻²⁰. In contrast to a purely kinematic account through an overturning streamfunction, this framework facilitates an evaluation of diabatic processes acting on different water masses. Access to diabatic processes allows one to link C\_ant overturning across water masses, and finally its routing into the ocean interior to the thermodynamic processes (air-sea buoyancy fluxes, etc) that maintain the circulation. It also warrants mention that the microscale irreversible processes that diffuse C\_ant across density interfaces are also best considered as a component of the thermodynamic controls on C\_ant. This contribution is evaluated separately.

Our goal is thus to develop a framework that links C\_ant uptake rates and uptake pathways to thermodynamic processes and the associated overturning circulation (Fig. 1). Given the rich complexity of the processes involved, we will proceed first using a *conceptual* model for C\_ant inventory formation, i.e., a model that contains only those ingredients essential for representing the C\_ant inventory formation. Our goal is to quantify the relative roles of air-sea exchange, water mass transformations, and diffusive exchanges using a streamlined conceptualization. The model is respectful of the outcrop areas and interior volumes of the observed water masses, and yet able to reproduce the observed inventory. To achieve this, we integrate elements derived from observations with information on diapycnal exchanges derived from a general circulation model.

Consider an ocean layer with volume V\_ρ bounded by two isopycnal surfaces, ρ and ρ + Δρ, that outcrops at the surface and leaves the domain via an open boundary (Fig. 3). A net volume flux Δψ exits the domain while
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematics of the water masses formation in an idealized oceanic basin. Here namely a subdomain of the global ocean with an open boundary. (a) Volume budget. Following standard definitions we denote as $V_ρ$ the volume sandwiched between the isoneutral surfaces $ρ$ and $ρ + Δρ$, the ocean lateral boundary and the ocean surface. $Δψ_ρ$ is the net volume flux of fluid exiting the domain $V_ρ$ across the open boundary; $Ω_ρ$ and $Ω_{ρ+Δρ}$ are the total volume fluxes across the isoneutral surfaces $ρ$ and $ρ + Δρ$, respectively (i.e. the diapycnal volume flux, defined as positive when toward denser waters (or overturning; units are Sverdrup (Sv) where $1 \text{ Sv} = 10^6 \text{ m}^3\text{/s}$)). (b) Tracer budget for $V_ρ$. Red lines indicate the tracer isosurfaces (i.e., surfaces on which the tracer is constant), which are not necessarily aligned with the isopycnals. The rate of change of the tracer inventory C in $V_ρ$ is given by the sum of the diapycnal transports $Φ_ρ$ and $Φ_{ρ+Δρ}$ across the isoneutral surfaces $ρ$ and $ρ + Δρ$, respectively, the air-sea fluxes across the layer outcrop surface $S_ρ$, the diffusive transports across $ρ$ and $ρ + Δρ$ and the tracer transport by the flow $Δψ_ρ$ across the lateral boundary (lateral red arrow). The sign convention is such that all the fluxes are positive if toward denser layers. The dashed line indicates the base of the mixed layer. The fluxes are computed over the entire isopycnal. In the main text the diapycnal fluxes of volume and tracer are discussed separately for the mixed layer (MWML) and for the interior.</figcaption></figure>

the diapycnal fluxes $Ω_ρ$ and $Ω_{ρ+Δρ}$ occur across the isopycnal surfaces $ρ$ and $ρ + Δρ$, respectively. The evolution of the tracer inventory $V_c$ in such a density layer, where $V$ is the constant water mass volume, $c$ its tracer concentration and $Φ_a$ is the air-sea fluxes over the outcrop area, is given by:

$$V_i \frac{dc_i}{dt} = \Phi_a + \Delta_{\text{diffusion}} + \Delta_{\text{overturning}} + c \Delta \Psi \quad (1)$$

The inventory will change if the combination of the air-sea flux of the tracer, the divergence across the bounding isopycnal surfaces due to diffusive fluxes ($Δ_{\text{diffusion}}$) and transport ($Δ_{\text{overturning}}$), and of the exchange $cΔΨ$ with the surrounding basins is non-zero.

We now invoke Occam's Razor in order to construct a conceptual model significantly simplified relative to those used previously to evaluate the $C_{\text{ant}}$ distribution$^{21–25}$. Consider the global ocean to consist of $n$ density layers that are in steady state, internally homogenous, and exposed to an atmospheric tracer whose concentration ($c_{\text{atm}}$) increases exponentially with time starting from a null value (i.e., $c_{\text{atm}} = c_0(e^{t/\tau_{\text{atm}}} - 1)$, similar to $C_{\text{ant}}$) (we do not consider here partial pressures). The simplest form of Eq. 1 for each density layer $i$ can be written as

$$V_i \frac{dc_i}{dt} = \phi_i + [\beta_i(c_{i-1} - c_i) - \beta_{i+1}(c_i - c_{i+1})] + [\omega_i c_{i-1} - \omega_{i+1} c_i] \quad (2)$$
where the subscripts $i-1$ and $i+1$ indicate the layers immediately above and below the layer $i$, $V_i$ is the constant volume of watermass, and $c_i$ its tracer concentration and where the diabatic transport of water (i.e., the water flux across the isopycnal) is $\omega_i \ge 0$ $\forall i$. (See Methods for details on the assumptions). The first term on the right is the air-sea flux integrated over $A_\rho$, the outcrop area. It is defined as $\phi_i = A_\rho k(c_{\text{atm}} - c_i)$, with transfer velocity $k$. Defining $H_i$ as a vertical scale that can be considered equal to the annual mean thickness of the mixed layer, we get a time scale $\tau_i = H_i/k$ (here $\tau_i$ and $H_i$ are considered equal for all the layers). The second term (in square brackets) is the diffusive term $\Delta_{\text{diffusion}}$ of Eq. 1, where $\beta_i$ is the equivalent of a volume transport across the interface between the layer $i$ and the neighboring density layers $i-1$ and $i+1$. The term in square brackets roughly mimics the convergence of a down-gradient flux of tracer $\beta\Delta c$, i.e., the diffusive transport of the tracer from tracer-rich layers to adjacent, tracer-depleted layers. The last terms represent water mass overturning, i.e., the effect of diapycnal volume transports $\omega_i$ and $\omega_{i+1}$ across the bounding layer surfaces. Given that the isopycnals do not necessarily coincide with the tracer isopleths it can happen that net tracer flux across a boundary is not zero even if the net water flux is. On the other hand, even if a tracer has a constant value across adjacent layers ($c_{i-1} = c_i = c_{i+1}$), a water volume convergence would act to change the layer tracer inventory while the microscale turbulence that promotes the tracer diffusion would play no role. These extreme cases illustrate the profound difference between the two diabatic processes. The water volume flux is partially driven by microscale turbulence that acts to change the temperature and the salinity of the water mass, thus changing its density. Since we consider the global oceans there is no transport $c\Delta\Psi$ outside the domain. The equations are solved numerically.

We consider three cases of increasing complexity. In each case results are tested against the GLODAP²⁶ C\_ant inventory. In each case we simulate the ocean tracer inventory change over 200 years, with $\tau_{\text{atm}} = 45yr$ and $c_0 = 1$. We choose a value of 0.1 for the density intervals and a value for $k$ that gives $\tau = 360$ days, $H = 60$ m. Density layer volumes and outcropping areas are taken from the observed values for the real ocean (Fig. 2, Table 1). The initial condition is a depleted ocean (i.e., $c_i = 0 \forall i$).

In the first case (AIRSEA) each ocean layer directly absorbs the atmospheric tracer but is isolated from the other layers (i.e., $\beta=0$ and $\omega_i = 0 \forall i$). This model has an inventory (Fig. 4a) that resides mostly in subpolar and intermediate waters, thus capturing one of the main observed features. Furthermore, the highest air-sea fluxes are into intermediate water masses (Fig. 4d).

Insights can be gained through closer inspection of the details of the solution. This case admits an analytical solution (see Methods) that can be approximated, for $e^{\int \tau_{\text{atm}}} \gg 1$, as $c_i = \gamma_i e^{\int \tau_{\text{atm}}}$, where $\gamma_i = \frac{c_0}{\left(1 + \frac{\tau_i}{\tau_{\text{atm}}}\right)}$, with $\tau_i' = \tau_i \frac{V_i}{H_i A_i}$ (i.e., the relaxation time is modulated by the ratio between the total volume and $H_i A_i$ the volume of the boundary layer of the water mass). In the case of the rather shallow subtropical and tropical waters $\frac{V_i}{H_i A_i} \sim 1$ and thus, given that $\frac{\tau_i'}{\tau_{\text{atm}}} \ll 1$, $\gamma_i \sim c_0$. In this case the inventory closely follows the atmosphere. For deep waters, since $\frac{V_i}{H_i A_i} \gg 1$, the opposite occurs, i.e., $\gamma_i \sim 0$ and consequently $c_i \sim 0$.

The air-sea fluxes, $\phi_i = \gamma_i \frac{V_i}{\tau_{\text{atm}}} e^{\int \tau_{\text{atm}}}$, and the inventory, $V_i c_i = V_i \gamma_i e^{\int \tau_{\text{atm}}}$, are both proportional to $V_i \gamma_i$. Since the volume $V_i$ increases with density while $\gamma_i$ decreases with increasing $\frac{V_i}{H_i A_i}$ (i.e., with increasing density) this product has a peak for intermediate waters. The result is a match between the analytical solution for a motionless ocean and GLODAP that is very good up to $\sigma_0 = 27.2$, i.e., to 75% of the observed inventory (Fig. 4d). This peaked distribution is independent of time (its structure, $V_i \gamma_i$, remains constant) while the inventory increases exponentially. In other words, the importance of SPMW and IW relative to global budgets does not depend upon time, thus fixing the internal turnover to their renewal time scales. The present proportion among water masses is thus not a transient despite the strong transient nature of the forcing. Finally, we note that the exponential nature of the solution explains why we found the GCM inventory and its rate of change to be proportional, as presented in the following.

Nevertheless, in the AIRSEA case, deep layers are largely depleted of tracer. By construction, the cumulative (in time) air-sea flux distribution (Fig. 4d) is the same as the inventory when both are viewed cumulatively with respect to density. By adding diffusion ($\beta_i = 10^8 m^3 s^{-1}$ and $\omega_i = 0 \forall i$; AIRSEA + DIFF) the denser water tracer content ($\sigma > 27.4$) has a better match with GLODAP (Fig. 4b), due to the internal transfer of tracer from intermediate waters to the deep layers. Nevertheless, the inventory is largely made of intermediate waters. Further, a difference between cumulative air-sea fluxes and the inventory (Fig. 4e) occurs only for the densest waters (a consequence of the internal transfer described above). Altering the diffusive constant ($\beta_i$) does not significantly alter the outcome, as in contrast to the overturning diffusion acts only between layers with different tracer concentrations. As light waters (i.e., in the tropical and subtropical regions) are nearly saturated and since we ignored the effect of temperature on saturation, concentrations are homogenous for these waters and a gradient occurs only between intermediate and deep waters.

The third case (AIRSEA + DIFF + DIA) incorporates the effect of the diapycnal transport overturning. This case requires evaluation of cross-isopycnal transport by oceanic circulation, here derived from a GCM simulation (see next section). To evaluate the impact of diapycnal transport on inventories and air-sea fluxes, we consider a density layer that has no exchange with the lighter layers yet loses tracer via diapycnal transport toward denser layers (with zero diffusion), i.e., where $\beta_i = \beta_{i+1} = 0$ and $\omega_i = 0$ and $\omega_{i+1} = \text{const}$, the latter thus making the divergence of the diapycnal transport ($\omega_i - \omega_{i+1}$).

In this simple case the tracer increase for each water mass is still exponential, but has a different shape factor:
<figure><img src="image_6.png" /><figcaption>Figure 4. Results from the conceptual model. Here the inventory obtained with conceptual model are compared with GLODAP. First we consider the AIRSEA case (a) with the GLODAP inventory in orange and in blue the inventory predicted by the numerically integrated conceptual model (Eq. 2). Both inventories integrals were normalised to 1 to ease comparison. The analytical solution presented in the text has not been superimposed since it coincides with the numerical one. (b) Same as (a) but for AIRSEA + DIFF case. (c) Same as (a) but for AIRSEA + DIFF + DIA case. The cumulated terms of the Eq. 2 are also presented, with each cumulated variable being normalised to unity. The case AIRSEA is shown in (d) with in blue the rate of change of the inventory (blue) and in red the carbon air-sea fluxes. For each density value these are thus the estimates of the corresponding rates of change occurring for the layers having lower or equal densities. The cumulated GLODAP inventory (orange), WOA outcrop areas (grey), WOA layer volumes (black) are also shown. The color scheme is the same in Fig. 5, where the equivalent GCM estimates are presented. (e) the same as (d) but for the AIRSEA+DIFF case. (f) the same as (d) but for the AIRSEA + DIFF + DIA case.</figcaption></figure>

$$ \gamma_{\omega i} = \frac{1}{\left(1 + \frac{\tau'_{i}}{\tau_{atm}} + \frac{\tau'_{i}}{\tau_{\omega i}}\right)} \qquad (3) $$

where $\tau_{\omega i} = V_i \omega_{i+1}$, is the water mass renewal or turnover time due to the divergence of the diapycnal transport. Evidently, if $\tau_{\omega i} \le \tau_{atm}$, i.e., if the water mass turnover by the overturning is faster or comparable with the atmospheric exponential increase, the inventory formation will be impacted by the former process.

Importantly, in this case

$$ \varphi_{\omega i} = c_0 \frac{\gamma_{\omega i} V_i e^{t/\tau_{atm}}}{\tau_{atm}} \left( 1 + \frac{\tau_{atm}}{\tau_{\omega i}} \right) \qquad (4) $$

Therefore, as a consequence of the term on the right in parenthesis, the relative importance of the flux of a given water mass might be significantly increased in the presence of a fast overturning rate. This result is not surprising. It also anticipates a possible mismatch between the air-sea fluxes and the inventory in correspondence of the water masses having high turnover rates and being impacted by a significant overturning into denser waters, as is the case for the SOC.

It is easy to demonstrate that for $\omega_i > 0$, $\varphi_{\omega i} > \varphi_i$. That is, in the presence of export of water towards denser water classes, the air-sea flux adjusts, causing increased oceanic uptake. Similar arguments can be used to show that the flux will decrease correspondingly in presence of a net overturning convergence.

## Global diapycnal Cant volume fluxes and formation of the Cant inventory

In order to characterize the role of the oceanic overturning, we first introduce a more comprehensive, exact formulation of the carbon inventory time evolution per water mass and evaluate the magnitude of the different terms in this formulation using a General Circulation Model (GCM). We then investigate the limits of validity for the
conceptual model described in the previous section using a GCM-based evaluation of the carbon transport due to overturning and gauge the specific role of SOC in determining the observed Cant inventory.

Consider again the isopycnal layer (Fig. 3) with the same diapycnal fluxes Ωρ and Ωρ+Δρ. (Note that it is also defined as G18 or A19). The net water mass formation rate is defined as the rate of accumulation MΔρ of the water with density ρ < ρ' < ρ + Δρ where M is the water mass formation per unit of density. It is calculated as the time variability of the volume (inflation/deflation) plus the flux entering or exiting the domain Δψ:

$$M\Delta\rho \equiv \frac{\partial V_{\rho}}{\partial t} + \Delta\psi \qquad (5)$$

Conservation of volume requires that MΔρ balance the convergence of the diapycnal fluxes (Fig. 3a), i.e.,

$$M\Delta\rho = \Omega_{\rho} - \Omega_{\rho+\Delta\rho} = -\Delta\rho\frac{\partial\Omega}{\partial\rho}. \qquad (6)$$

The net diapycnal volume transport (Ωρ; the physical overturning or water mass transformation¹⁷; here in the formulation of ref. 20) is given by the sum of the processes acting on density (the diabatic processes related to the ocean thermodynamics) (see Supplementary Material for details on the formulation). The first is the divergence of mixing processes (dρ). This includes microscale diapycnal mixing, - lateral mixing in the mixed layer where isopycnal are nearly vertical. Mixing along the isopycnal, due to mesoscale and submesoscale eddies and filaments in the interior, diabatically impacts the water mass by causing a net transfer of water due to non-linearity of the equation of state (cabbelling, and, when using neutral density, thermobaricity). The second contribution comes from the boundary buoyancy forcing (fρ), i.e., the various surface heat fluxes and the geothermal heating acting on temperature and the freshwater flux due to evaporation, precipitation, ice production and melting and river runoff.

In practice the divergence of mixing processes (dρ) and buoyancy forcings (fρ) are identified by diagnosing and combining the equivalent terms for temperature and salinity tendencies. In fact it can be proven¹⁷,²⁰ (see Supplementary Note) that:

$$\Omega_{\rho} \equiv \frac{\partial}{\partial \rho} \int_{V_{\rho}} d_{\rho} dv + \frac{\partial}{\partial \rho} \int_{V_{\rho}} f_{\rho} dv \qquad (7)$$

where the integrals are evaluated over the volume sandwiched by the isopycnal layer ρ, the ocean bottom and the basin open boundary. Positive transformation values indicate transports toward denser waters¹⁷⁻²¹. Here we use potential density (σ₀) to define the reference frame.

We next extend this density-based volume budget framework by incorporating oceanic tracers. Consider the dynamical budget of a passive tracer within Vρ (bounded by two isopycnal surfaces, ρ and +Δρ) as presented in Eq. 1 but for the full time-varying general case (i.e., a time varying water mass volume; see Supplementary material for the case of active tracers). For a passive tracer c (the same as above), the equation for the rate of inventory change, C, can be written (see Fig. 3b and Supplementary Note²¹):

$$\frac{dC}{dt} \equiv \int_S \phi_a ds + \int_{V_\rho} d_c dv - \Delta\rho \frac{\partial \Phi}{\partial \rho} - \int_A cu \cdot nds \qquad (8)$$

As in Eq. 1 here the first term on the rhs is contribution of the air-sea exchange (the integral over the outcrop of the air-sea fluxes per unit area φa) and the second term is Δdiffusion, the net diffusive flux across the boundaries of Vρ (equivalent to the integral over the volume of the convergence (dρ) of the tracer diffusion). The third term is again Δoverturning, the convergence of the tracer flux (Φρ) due to the diapycnal volume transport the overturning across the bounding isopycnal surface (i.e., by the diapycnal volume transport as defined above). The last term is the tracer transport at the basin boundary, where u is the velocity at the boundary, c is the tracer concentration, and n the unit vector perpendicular to the boundary.

For a given ρ, the diabatic tracer transport across the isopycnal (Φρ) is evaluated²¹ by quantifying the thermodynamic drivers of the overturning (the same as in Eq. 7). Specifically, the diapycnal tracer transport, simplified as the product between the diapycnal transport ω and c in Eq. (2), is due to the combined variation of c and the constituents of the physical diapycnal transport in Eq. (7) along the isopycnals:

$$\Phi_{\rho} \equiv \frac{\partial}{\partial \rho} \int_{V_{\rho}} c d_{\rho} dv + \frac{\partial}{\partial \rho} \int_{V_{\rho}} c f_{\rho} dv. \qquad (9)$$

where the integrals are again evaluated over the isopycnal layer volume. In other words, the tracer transport is expressed as the sum of all the physical processes sustaining the overturning, from solar radiation and evaporation to mixing due to wind or to cabbeling. Eq. (9), and consequently Eq. (8), include a large number of terms derived by using salinity and temperature tendencies to build Eqs 3 and 5. Since the equation is linear, each term can be evaluated separately. As Cant is a passive tracer, these equations and concepts can be used to quantify specific processes that influence formation of the Cant inventory. The same argument holds for any similar biogeochemical tracer. As reported in the Supplementary text, this equation has the same properties as the approach developed by Walin for the volume flux computation. It is intrinsically Lagrangian since it naturally follows isopycnals, i.e., it considers only the velocity of the water relative to the isopycnal). It considers all possible cases, including the Cant exchange across isopycnals that occurs when the sun warms the stratified region below the
<table><caption>**Table 2.** Summary of the main fluxes controlling the model global anthropogenic CO₂ inventory (C<sub>ant</sub>). The values have been estimated using Eq. 4. The diffusive contribution is evaluated as a residual that implies that all kind of diffusive processes are considered, i.e., including any potential numerical diffusion. The coarse-graining of the ocean density structure into five discrete layers is as follows: tropical waters (TW), subtropical mode waters (STMW), subpolar mode waters (SPMW), intermediate water (IW), and deep water (DW). Values are net for each water mass.</caption><thead><tr><th></th><th>Total</th><th>TW</th><th>STMW</th><th>SPWM</th><th>IW</th><th>DW</th></tr></thead><tbody><tr><td>(a) Ocean gain (1995; dC/dt)</td><td></td><td>σ<sub>0</sub> &lt; 24.5</td><td>24.5 &lt; σ<sub>0</sub> &lt; 26.6</td><td>26.6 &lt; σ<sub>0</sub> &lt; 27.05</td><td>27.0 &lt; σ<sub>0</sub> &lt; 27.55</td><td>27.55 &lt; σ<sub>0</sub></td></tr><tr><td>Inventory [PgC/yr]</td><td>2.07</td><td>0.16</td><td>0.54</td><td>0.42</td><td>0.52</td><td>0.43</td></tr><tr><td>Fractions [%]</td><td>—</td><td>7.7</td><td>26.1</td><td>20.3</td><td>25.1</td><td>20.8</td></tr><tr><td>Cumulative fractions [%]</td><td>—</td><td>7.7</td><td>33.8</td><td>54.1</td><td>79.2</td><td>100</td></tr><tr><td>(b) Air-sea</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fluxes [PgC/yr]</td><td>2.04</td><td>0.57</td><td>0.60</td><td>0.36</td><td>0.42</td><td>0.09</td></tr><tr><td>Fractions [%]</td><td>—</td><td>27.9</td><td>29.4</td><td>17.6</td><td>20.6</td><td>4.4</td></tr><tr><td>Cumulative fractions [%]</td><td>—</td><td>27.9</td><td>57.4</td><td>75.0</td><td>95.6</td><td>100</td></tr><tr><td>(c) Diapycnal fluxes divergence (Δ<sub>overturning</sub>)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Global [PgC/yr]</td><td>—</td><td>-0.36</td><td>0.12</td><td>0.12</td><td>-0.07</td><td>0.19</td></tr><tr><td>MWML [PgC/yr]</td><td>—</td><td>-0.47</td><td>0.13</td><td>0.23</td><td>-0.03</td><td>0.16</td></tr><tr><td>Interior [PgC/yr]</td><td>—</td><td>0.11</td><td>0.00</td><td>-0.11</td><td>-0.04</td><td>0.03</td></tr><tr><td>(d) Air-sea + Diapycnal fluxes (φ + Δ<sub>overturning</sub>)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Global [PgC/yr]</td><td>2.04</td><td>0.21</td><td>0.72</td><td>0.48</td><td>0.35</td><td>0.28</td></tr><tr><td>MWML [PgC/yr]</td><td>2.04</td><td>0.10</td><td>0.73</td><td>0.59</td><td>0.39</td><td>0.25</td></tr><tr><td>(e) Residual (a minus d; -Δ<sub>diffusion</sub>)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GLOBAL [PgC/yr]</td><td>—</td><td>-0.05</td><td>-0.18</td><td>-0.06</td><td>0.17</td><td>0.15</td></tr><tr><td>(f) Inventory [PgC]</td><td>94.30</td><td>8.00</td><td>26.00</td><td>19.2</td><td>23.7</td><td>17.4</td></tr></tbody></table>

mixed layer even in absence of any water movement. In this latter case the radiative warming is de facto moving down the isopycnals creating a net diapycnal flux.

By introducing the tracer diapycnal flux due to buoyancy-altering processes Φ<sub>ρ</sub> (Eq. 9), we introduce the concept of tracer transformation associated with water mass transformation. In addition, using Eq. 2, the tracer diapycnal transport associated to the transformation becomes an explicit term in the dynamical balance for the inventory. Note that near-surface buoyancy gradients do not strictly co-vary with C<sub>ant</sub> gradients, therefore the transformation of volume and the associated C<sub>ant</sub> transport may differ in magnitude and by region.

Water mass formation is generally associated with transformations that occur in the mixed layer¹⁹. There, water mass properties are set by air-sea interactions and the largest diapycnal fluxes occur. Here we use a generalization of the concept of water mass formation since here we will evaluate, independently, the transformations occurring in the mixed layer and in the interior. Specifically, to differentiate subsurface from surface processes, we evaluated the change in diapycnal volume and C<sub>ant</sub> fluxes that only occur in the upper maximum winter mixed layer (MWML). Importantly, in this case the net water mass formation coincides with net subduction²⁷⁻²⁹. This is why the MWML is a common reference horizon in physical and biogeochemical studies (e.g., in the concept of preformed nutrients and carbon³⁰).

**Diagnosing connections between the C<sub>ant</sub> inventory, air-sea fluxes and diapycnal volume fluxes.** The physical model used here is the global ORCA2-LIM³¹ ice-ocean model with nominal two-degree global resolution and isoneutral mixing (see Methods). The model was spun up for 1500 years with repeating seasonal forcing fields²⁰,²¹,²⁹. Bi-weekly means of the pertinent physical state variables were saved over the final year for analysis and to drive the biogeochemistry model. A twenty-four compartment biogeochemistry model, PISCES (Pelagic Interaction Scheme for Carbon and Ecosystem Studies)³², was spun up offline using the physical state over 3000 years with pre-anthropogenic boundary conditions. In 1860 the run was split into a control with the same pre-industrial atmospheric concentrations and another with the observed transient atmospheric CO₂. Following standard nomenclature⁶, we refer to the first as “natural” carbon, the second as “contemporary” carbon, and the difference as “anthropogenic carbon”.

By 1995, nearly 80% of the modelled C<sub>ant</sub> had accumulated at densities shallower than σ<sub>0</sub> < 27.55, encompassing the SOC and IW (Table 2, Figs 1 and 2). Up to 34% of the total C<sub>ant</sub> inventory was found in tropical and subtropical waters (σ<sub>0</sub> < 26.6). The model distribution of C<sub>ant</sub> is in good agreement with observational estimates (Tables 1 and 2; Fig. S2) having most accumulation in mode and intermediate waters but these waters are too cool in the model (a shift to denser values is observed) and the deep inventory is generally underestimated. 96% of the total air-sea flux occurs over SOC + IW outcrops and 57% (1.17 PgC/yr, Table 2, red line in Fig. 5) occurs over tropical and subtropical water outcrops. For these shallow water masses the relative importance of air-sea fluxes is greater than the relative fraction of the yearly C<sub>ant</sub> accumulation rate (Table 2, blue line in Fig. 5; 96% vs 80%, or 57% vs 34%). This mismatch between C<sub>ant</sub> uptake via gas exchange and interior accumulation implies a net transfer of C<sub>ant</sub> to higher densities subsequent to initial uptake.
<figure><img src="image_7.png" /><figcaption>**Figure 5. Global carbon budgets.** Model global estimates of the terms in Eq. 8. Integrals in Eq. 8 are performed cumulating from $\sigma_0 = 20.0$. Grey: Integrated outcrop area (i.e., for each density value $\sigma'_0$ the area where $\sigma_0 < \sigma'_0$); Red: Integrated uptake of anthropogenic carbon by gas exchange during 1995; Blue: Integrated rate of interior accumulation during 1995 (Note: in fact, it is very similar to the integrated global model inventory of $C_{ant}$); Black: Integrated volume of the ocean density layers. The difference between the red and the light-blue curve illustrates the importance of processes other than the air-sea fluxes setting the inventory. Additionally, the yellow curve represents the air-sea fluxes minus the diapycnal transports of $C_{ant}$ in the MWML (Fig. 6b; Eq. 8) while the orange curve represents the air-sea fluxes minus the total diapycnal transports (first and last terms in Eq. 8). The difference between the latter and the inventory is due to diffusive $C_{ant}$ fluxes.</figcaption></figure>

The underlying mechanisms for this transfer are explicitly quantified here using the thermodynamic water mass transformation framework introduced above. Evaluating Eq. 8 for each water mass gives the dynamics of the formation of the corresponding $C_{ant}$ inventory (Table 2). To provide an overview of the role of the different processes we discuss first the cumulative quantities for the different terms. For a given density $\rho$, evaluating Eq. 8 from 0 to $\rho$ we have the rate of change of the inventory for the layers $0 < \rho' < \rho$ given as a function of the air-sea interaction, diapycnal diffusion and transport of $C_{ant}$ across $\rho$ and transport at the basin boundary (here zero since we integrate globally). As shown in Fig. 5, for a given density $\rho$ a large fraction of modelled rate of change of the $C_{ant}$ inventory (left term in Eq. 8; blue line) can be explained by subtracting the diapycnal carbon fluxes $\Phi$ (orange line) $\rho$ from the air-sea fluxes of $C_{ant}$ (first term on right in Eq. 8; red line) using Eq. 8. Note that the diapycnal volume transport generally promotes a transfer of $C_{ant}$ from light to dense waters. The convergence of the diapycnal $C_{ant}$ fluxes in a layer bounded by two isopycnal surfaces (the formation) gives an estimate of the net contribution of the physical overturning to the layer inventory of $C_{ant}$. Quantitatively, the diapycnal fluxes inject 0.12 PgC/yr into STMW, 0.12 PgC/yr into SPMW and 0.19 PgC/yr into deeper waters (Table 2c).

We show here that the mechanisms that promote ocean ventilation are: a) subsurface transport toward extra-tropical latitudes and then subduction (i.e., the net diapycnal convergence of $C_{ant}$ in the MWML); b) the diffusive $C_{ant}$ transport resulting from the diapycnal (mostly vertical) $C_{ant}$ gradients below the MWML. The latter process is the residual between the blue and the orange lines in Fig. 5 (difference of air-sea fluxes and diapycnal formation fluxes Equation 8; Table 2d,e). Diffusion exports $C_{ant}$ away from light, shallow waters into the intermediate and deep waters (Table 2e). The diffusive flux divergence peaks at 0.18 PgC/yr in STMW (Table 2e) and exceeds diapycnal formation fluxes for that layer (0.12 PgC/yr).

Contrasting formation fluxes for the entire ocean and the MWML (Table 2c) emphasises the importance of surface and interior water mass transformations ($\Phi_\rho$) for $C_{ant}$ storage. In fact MWML fluxes are greater than fluxes into the SPMW density classes (by 0.11 PgC/yr; Table 2c). The diapycnal transport of $C_{ant}$ out of TW is also more important than the total fluxes in that density class; -0.47 PgC/yr, or about 26% of the global air-sea flux, vs. -0.36 PgC/yr (Table 2c). Negative values indicate removal of $C_{ant}$ through water mass erosion and a transfer of $C_{ant}$-laden waters into waters of other density classes. Once air-sea fluxes and diapycnal fluxes due to water mass transformations in the MWML account for most of the inventory change for tropical and subtropical water masses (Fig. 2b). In other words, the surface diapycnal fluxes of the SOC and the air-sea fluxes sustain a $C_{ant}$ distribution that is close to the actual inventory change. In the interior, overturning moves $C_{ant}$ into lighter waters while diffusion does the opposite. The two mechanisms nearly cancel for layers that are approximately $C_{ant}$-saturated and thus in near steady-state.

By construction with our cyclo-stationary physical state for the ocean, the water mass volumes do not change, and thus there is no net water mass volume inflation (i.e., net convergence of diapycnal transports). Therefore the interior processes that sustain the overturning of the water masses compensate the diapycnal transport occurring in the MWML (Fig. 1). In the latter layer surface buoyancy fluxes tend to densify the water masses. In fact in the interior the buoyancy diffusion and the penetration of solar radiation below the MWML²⁰ generally move waters from dense layers to lighter ones (driving an upwelling). This diapycnal volume transport (or upwelling) in the interior would have no effect upon $C_{ant}$ if the interior layers were depleted in $C_{ant}$. The presence of $C_{ant}$ in the interior waters is such that they have an effect (by re-entraining $C_{ant}$ into the mixed layer), but their contribution is largely overwhelmed by the diapycnal divergence in the MWML. Finally, the overturning of $C_{ant}$ associated with DW formation is found to be significant (0.19 PgC/yr; mostly NADW formation, see below), and leads to a
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Net MWML diapycnal transports. (a) Volume (Ωρ) and (b) C<sub>ant</sub> (Φ<sub>ρ</sub>) *transformation* fluxes in the MWML due to all processes over the density range 22 ≤ σ<sub>0</sub> ≤ 28. In (a–b), oblique magenta arrows underline the large difference between volume and C<sub>ant</sub> transports associated with the upwelling of CDW. The green arrow points to the density range showing the greatest difference between heat and C<sub>ant</sub> transformation fluxes, the latter having a much smaller (but not negligible) re-emergence at the Equator. Line colors show global (black) and latitude-separated results: North of 30°N (blue), between 30°N and 30°S (orange), south of 30°S (red). In all panels, positive (negative) *transformation* values represent diapycnal transports towards denser (lighter) densities. The slope of the lines represents the *formation*: convergence (negative slope) and divergence (positive slope). Values for σ<sub>0</sub> < 26.4 were smoothed with a 3-points running mean filter.</figcaption></figure>

progressive densification and enrichment of C<sub>ant</sub> in DW source waters while moving from the subtropics to the
high latitude formation sites.

In summary, the MWML component of the overturning dominates the diapycnal fluxes Φρ of Cant and con-
sequently shapes the interior Cant distribution. A closer look shows that in the MWML, transformation-related
diapycnal transport of Cant is associated with the upper poleward branch of the SOC in which light waters are
converted into mode and intermediate waters (Fig. 6). The MWML volume transformation reaches a maximum
across σ₀ = 26.4 (Fig. 6a). This peak densification transfers 0.4 PgC/yr (Fig. 6b) of Cant into denser water masses.
This rate is about 20% of the global Cant uptake for the model (2.07 PgC/yr in 1995, Table 2a). The peak Cant
transformation is not co-located with the peak volume transformation (Fig. 6a), but rather occurs at much lighter
density due to the higher Cant content of those waters.

The diapycnal C<sub>ant</sub> transport towards denser waters within the MWML is quite symmetric for both hemi-
spheres in waters with σ<sub>0</sub> < 26.6 and poleward of 30° (Fig. 6b). Subtropical sources inject C<sub>ant</sub> into denser classes
at a maximum rate of nearly 0.3 PgC/yr. At higher densities in the Northern Hemisphere, a peak for NADW
formation (σ<sub>0</sub> ~ 27.7; 0.2 PgC/yr) is evident (Fig. 6a,b). This process is fed by overturning more than local air-sea
C<sub>ant</sub> fluxes (Fig. 6b). Neither volume nor C<sub>ant</sub> fluxes scale as well in the Southern Ocean. Bottom water formation
in the Antarctic is not well resolved in potential density space<sup>20</sup>, nor does it have a large impact on the model C<sub>ant</sub>
inventory (Fig. 6b). Upwelling and transformation of CDW into SAMW in the Southern Ocean is clearly visible
as negative values in Fig. 6a (~10–15 Sv), but these are not associated with a significant C<sub>ant</sub> transport (magenta
arrows, Fig. 6a,b). The upwelling carbon-depleted CDW increases the air-sea fluxes by increasing the difference
in pCO<sub>2</sub> between the atmosphere and the ocean. Importantly the supply of C<sub>ant</sub> to SAMW by the water mass over-
turning is dominated by the contribution from the subtropics (i.e., by the SOC). In the tropics the diapycnal C<sub>ant</sub>
transport is also very different from the diapycnal volume transport (σ<sub>0</sub> < 24.5, green arrow, Fig. 6b) due to lower
C<sub>ant</sub> content of newly obducted waters.
Overall, the upper branch of the global overturning (that occurring mostly in the MWML) is thus exporting C<sub>ant</sub> to mid and high latitudes. It is important to reconcile this result with the current estimates of the meridional transport of C<sub>ant</sub>. The present consensus is that globally C<sub>ant</sub> converges toward the Equatorward (except for the North Atlantic<sup>2,16</sup>). This contradiction is only apparent. By invoking a water mass framework, we are here evaluat-
ing the C<sub>ant</sub> that is redistributed by the overturning and not on the vertically integrated (meridional) transport. As illustrated with the conceptual model, a large component of C<sub>ant</sub> stays in the isopycnal layers in which it entered via air-sea fluxes. In other words the meridional transport is being mostly sustained by flow *along* isopycnals (i.e., the adiabatic component of the flow). Given that the air-sea fluxes of C<sub>ant</sub> are maximal for subpolar mode and intermediates waters, and as these water masses in the interior exhibit a net flow from the high and mid-latitudes towards the tropics, they provide most of the meridional transport. Importantly, in the North Atlantic the pole-
ward transport of C<sub>ant</sub> in the upper layer by the surface circulation associated with the NADW formation is such that it dominates the meridional transport locally. The role of the overturning considered here is consistent with the view presented in ref. 33, although their kinematic approach did not allow for an explicit connection to ocean thermodynamics.

An important benefit of a thermodynamic-based water mass framework to C<sub>ant</sub> (Eq. 9) is that it allows for identification and quantification of all mechanisms that control the transfer of C<sub>ant</sub> into the ocean interior. The air-sea exchange of buoyancy sustains a net flux towards higher densities for waters denser than σ<sub>0</sub> = 24.0 (peak of 0.5 PgC/yr at σ<sub>0</sub> = 25.0 in Fig. 7a). Similarly, atmospheric buoyancy fluxes transfer C<sub>ant</sub> (0.6 P gC/yr) toward waters lighter than σ<sub>0</sub> = 22.0. The densification of tropical waters is evenly distributed between the two hemi-spheres (Fig. 7a, σ<sub>0</sub> = 24.5). The densification of light tropical waters during their poleward journey is largely sustained by net evaporation (peak at σ<sub>0</sub> = 23.5 in Fig. 7b). Cooling dominates transformation at higher densities (Fig. 7b), but the both processes are needed to sustain the C<sub>ant</sub> export to STMW and SAMW. The transfer towards light (tropical) waters is due to the combined effects of penetrative solar radiation (up to 0.8 PgC/yr at σ<sub>0</sub> = 23.0) and freshwater input, centered at σ<sub>0</sub> = 21.5 (effectively the Warm Pool). The net transfer of C<sub>ant</sub> from the sub-tropics to the Southern Ocean by buoyancy forcing is balanced by the local net gain of freshwater at the surface (negative peak at σ<sub>0</sub> = 27.3, Fig. 7b). In the Northern Hemisphere, high C<sub>ant</sub> fluxes into NADW are associated with cooling of subtropical waters that have accumulated anthropogenic C<sub>ant</sub> while moving poleward. Categorizing evaporation and cooling as the main drivers of poleward transport of C<sub>ant</sub> holds even when considering mixing, which also influences heat and freshwater budgets (Fig. 7c). However, the surface heat and freshwater fluxes that induce a transport towards the lighter waters across σ<sub>0</sub> = 22.0 (Fig. 7c) are overcompensated by mixing (Fig. 7c). This is more clearly visible in Fig. 7d. The large surface input of heat and freshwater that promotes the forma-
tion of light waters is more than compensated by the (mostly vertical) mixing in the tropical and subtropical part of the SOC. In the SOC, mixing is the thermodynamic counterpart of the Equatorial upwelling and surface divergence. Mixing is, however, subordinate to the release of heat to the atmosphere in further densification of water parcels along the poleward journey at the surface. This holds at least up to the higher water densities, σ<sub>0</sub> = 26.9–27.3, where mixing processes again contribute significantly.

A conceptual model of the global uptake and storage

Finally we return to the conceptual model to see if the addition of diapycnal C<sub>ant</sub> transport from subtropical to
intermediate and dense waters derived from the GCM significantly improves the solution there. If that were to be
the case, then we would have identified the main drivers in forming the reservoir. In other words, the conceptual
model's elements (water mass aspect ratios and the presence of a shallow overturning) are sufficient to capture
the important dynamics.

We assume here that only the shallow component of the overturning is relevant. Specifically, we use the values for ω given from the ocean circulation model (Fig. 6a; Methods) for the export of water from the tropics toward the subtropical and subpolar regions.

The addition of overturning makes a dramatic difference. First, the bimodal distribution of the GLODAP C<sub>ant</sub> inventory is now well reproduced (Fig. 4c; Methods), with a separated peak for the deep reservoir. This is formed by densification of C<sub>ant</sub>-rich water rather than by local air-sea fluxes and it reproduces the importance of NADW formation to the total C<sub>ant</sub> inventory. In addition, the peak for the lighter water masses is now more correctly cen-
tered between SPMW and IW rather than on IW. There is only now a strong separation between the cumulative (in time) air-sea fluxes and inventories, when both are viewed as cumulative quantities with respect to density (Fig. 4f), as the air-sea flux component has now shifted towards lower densities. Stated differently, it is only now with the explicit inclusion of the overturning component that the conceptual model is able to match the results obtained with the ocean circulation model (compare Figs 4f and 5).

The conceptual model here has considerable skill in reproducing the C<sub>ant</sub> inventory and its relation to the air-sea fluxes. It provides mechanistic insight in three ways. First, it illustrates how the aspect ratio of the ocean layers, that is set by the oceanic response to wind and buoyancy fluxes<sup>24</sup>, is the dominant factor, with IW having the optimal geometry. Second, it illustrates the complementary roles played by tracer diffusion across isopycnal surfaces and the overturning circulation. In fact, while diffusion is important, it is the overturning that guarantees a flushing out of the tropical and subtropical regions in favour of mode and IWs, thereby sustaining enhanced uptake via gas exchange at low latitudes. Third, only the dense overturning component (σ<sub>θ</sub> > 27.5) can account for the deep C<sub>ant</sub> reservoir since these waters are weakly influenced by air-sea exchange. This transfer of C<sub>ant</sub> into the densest waters creates an inventory gap between mode/intermediate waters and deep waters.

To summarize, the conceptual model underscores the central role of the shallow overturning circulation (SOC) in determining the interior inventory of C<sub>ant</sub>. The predicted inventory is sensitive to details of the imposed overturning (see also ref. 24). This exercise illustrates the importance of understanding the mechanisms respon-
sible for the overturning circulation and the sensitivity of that circulation to various perturbations, including the
response to climate change.
a

<figure><img src="image_9.png" /></figure>

b

<figure><img src="image_10.png" /></figure>

c

<figure><img src="image_11.png" /></figure>

d

<figure><img src="image_12.png" /></figure>

**Figure 7.** Details of the processes sustaining the carbon overturning. (a) Carbon transformation fluxes ($\Phi_p$) in the MWML due to surface buoyancy fluxes only. Line colors show global (black) and latitude-separated results: North of 30°N (blue), between 30°N and 30°S (orange), south of 30°S (red). (b) Same as (a) (black), showing also the decomposition into the global freshwater (blue) and heat (red) components. (c) Same as (b) (black), showing also the decomposition into the global freshwater (blue) and heat (red) components. (d) Same as (c) (black), showing the decomposition into surface buoyancy fluxes (blue same as (a)) and buoyancy mixing fluxes (solid red line). The dashed red line represents the effect of the lateral mixing processes while the dot-dashed red line indicates the effect of vertical mixing processes. In all panels, positive (negative) transformation values represent diapycnal transports towards denser (lighter) densities. The slope of the lines represents the formation: convergence (negative slope) and divergence (positive slope). Values for $\sigma_0 \le 26.4$ were smoothed with a 3-points running mean filter.

Discussion

The framework presented here provides a means to reconcile the disparate density signatures of $C_{\text{ant}}$ uptake by gas exchange (considered cumulatively in time) and $C_{\text{ant}}$ accumulation in the ocean interior. As such, it is our hope that the methods can be applied more generally to the synthesis of surface carbon products (SOCAT) and ocean interior products (for instance, data comprised by GO-SHIP, or GLODAPv2 ³⁴). In fact the method is quite general, and can be applied to interpret any biogeochemical tracers that are impacted by the overturning circulation of the ocean.

We set out to test the hypothesis that ocean overturning and diffusion strongly modulate the pathways by which $C_{\text{ant}}$ enters the ocean interior, and we developed a schematic view whereby $C_{\text{ant}}$ “rides” or participates in the overturning structures associated with the poleward transport of heat by the ocean³⁵. To test this idea, we applied a thermodynamic water mass framework¹⁷ to the case of $C_{\text{ant}}$. Previous studies have invoked overturning and water mass concepts to interpret $C_{\text{ant}}$ uptake¹,³,¹³,³³ but these have not explicitly considered the ensemble of physical processes involved in regulating the uptake of $C_{\text{ant}}$ by the ocean. Here, the joint use of the simplified
conceptual model and the full three-dimensional OGCM provided a means to deconvolve the mechanistic con-
trols on C<sub>ant</sub> uptake pathways into the ocean.

The ocean C<sub>ant</sub> inventory is strongly influenced by water mass transformations and diffusion, thereby reflect-
ing the large-scale overturning circulation<sup>3</sup>. Beginning with the equatorial divergence, the poleward-flowing
upper branch of the Subtropical Cells absorb C<sub>ant</sub> via gas exchange. In transit these waters become denser in
stages, culminating with a substantial heat loss in western boundary currents and on the poleward flanks of
the Subtropical Cells. This pathway provides a conduit for C<sub>ant</sub> absorbed in the low latitudes to enter the ocean
interior in the subpolar water masses. Both direct overturning and C<sub>ant</sub> diffusion have proven to be important,
consistent with results obtained with kinematic diagnostics of subduction across the mixed layer base<sup>36</sup>.

Although our analysis confirms the importance of the Southern Ocean and the NADW formation region
for C<sub>ant</sub> transport into the ocean interior, our results also highlight the role low latitude accumulation offers in
a critical supporting role. Previous models and kinematic diagnostics<sup>12,37</sup> have led to a local characterization of
subduction of C<sub>ant</sub> within SAMW and AAIW. These studies imply injection rates that are consistent with our
results, namely that 40% of C<sub>ant</sub> accumulates in SAMW and AAIW<sup>1</sup>. The insight provided here is the connection
between mode water formation and large scale water mass pre-conditioning at lower latitudes. The proposed
mechanism, diabatic exchanges between STMW and SAMW, is supported by eddy-permitting models<sup>38</sup> and data
inferences<sup>37,39,40</sup>. Since water mass transformations are directly impacted by eddy parameterizations and model
resolution, the exact numbers may vary but the main processes are represented.

This study has established a framework that links C<sub>ant</sub> uptake to thermodynamics. The climatological poleward
heat transport in the Shallow Overturning Circulation plays a central role. This result motivated our considera-
tion of the more general Shallow Overturning Circulation framework, and the view that more than 50% of SPMW
formation sources are subtropical. More generally, the unified STMW/SPMW grouping is best viewed as a bind-
ing agent between the subtropical and subpolar dynamical regimes. In turn, these regimes may be expected to
experience direct but different impact under climate change. The method developed and applied here is perfectly
general, and it is our hope that this method will also be applied to interpret Earth system models to understand
how the ocean's full carbon reservoir will be perturbed under future climate change.

**Methods**

**Analysis of GLODAP C<sub>ant</sub> from WOCE-era measurements.** The inventory of anthropogenic carbon¹ in density space has been derived using the GLODAP data product¹,²⁶. We consider the partitioning of C<sub>ant</sub> into water masses. We begin with a coarse-grained analysis of the anthropogenic carbon inventory in density space as derived from the ΔC⁻-based GLODAP data product¹,², and from the TTD-based product¹⁵. These estimates have been recently compared⁷. While pattern differences were identified, they remain quite similar in terms of global inventory per water mass in the upper 1000 meters of the water column (see Fig. S2). Since we have chosen to group density layers below intermediate horizons into a single class, we focus only on GLODAP.

We calculated potential density for our reference frame using temperature and salinity climatologies from the World Ocean Atlas 2009 (WOA09)⁴¹. We have chosen five discrete levels: (i) Tropical Waters (TW) with σ₀ ≤ 24.5; (ii) lower thermocline waters (STMW) with 24.5 < σ₀ ≤ 26.6; (iii) sub-arctic and Subantarctic Mode Waters (SAMW) with 26.6 < σ₀ ≤ 27.1; (iv) Intermediate Waters (AAIW, NPIW) with 27.1 < σ₀ ≤ 27.4; and (v) deep waters (DW) with σ₀ > 27.4. The separation between TW and STMW at σ₀ = 24.5 was chosen to match the separation between the equatorially divergent (TW) and convergent (STMW) flows in the Indo-Pacific Subtropical Cells (STCs)⁴².

Inventories and volume-weighted average concentrations were calculated for each density layer. Binning was
performed in two ways. First, density binning was performed over the full ocean volume. Second, the same pro-
cedure was followed, but the maximum winter mixed layer depth was used to separate the ocean interior from the
region directly impacted by seasonal variations. We used a data-derived mixed layer depth climatology⁴³ to define
the needed levels. We used the 2008 product, that includes Argo data and captures, to first-order, the shape of the
maximum wintertime mixed layer in remote regions such as the Southern Ocean.

**The Ocean Carbon Model.** The ice-ocean model we used is a nominal two-degree global version of ORCA2-LIM-PISCES³² with 31 vertical levels. It includes the classical Gent and McWilliams mesoscale parameterization⁴⁴ with lateral mixing evaluated along isoneutral surfaces. The lateral exchange coefficient depends on the baroclinic instability growth rate⁴⁵. The vertical mixing scheme uses turbulent closure⁴⁶, and there is a diffusive bottom boundary layer parameterization⁴⁷. Background vertical diffusivity increases downward to account for decreased stratification and increased small-scale turbulence near the bottom. Values range from 0.12 10⁻⁴ m²s⁻¹ in the top kilometer to 1.2 10⁻⁴ m²s⁻¹ at 5000 m. Vertical mixing due to convection was parameterized by locally enhancing vertical diffusivity. Air-sea fluxes of heat and freshwater (evaporation) are evaluated using bulk formulae⁴⁸. A clear-water penetrative short-wave solar radiation formulation was used and surface salinities were restored to climatology. Climatological European Remote Sensing Satellite-1 and -2 scatterometer monthly mean wind stresses were used for the tropics, and the NCEP-NCAR climatological dataset was used poleward of 50°N and 50°S. At the bottom, geothermal heating was represented as a spatially variable heat source.

The model was spun up for 1500 years with climatological repeating seasonal forcing fields. Natural and anthropogenic climate variability were not considered. These fields are identical to those used in previous studies²⁰,²¹,²⁹. A detailed description of the model and the simulation characteristics are included in those references. From the final spinup year, 15-day means of circulation fields and tendency terms were stored for analysis.

We used the PISCES (Pelagic Inter- action Scheme for Carbon and Ecosystem Studies) model³² biogeochem-
istry with 24 components. Phytoplankton growth can be limited by five different nutrients: nitrate, ammonium,
phosphate, silicate and iron. Four living pools are represented: two phytoplankton size classes/groups (nanophy-
toplankton and diatoms) and two zooplankton size classes (microzooplankton and mesozooplankton). Diatoms
differ from nanophytoplankton by their need for Si, by higher requirements for Fe and by higher half-saturation constants because of their larger mean size. For all living compartments, the ratios between C, N and P are kept constant. The internal Fe contents of phytoplankton groups and the Si of diatoms are a function of the external concentrations in nutrients and the light level. The Chl/C ratio is modelled explicitly. All zooplankton elemental ratios are kept constant. There are three non-living components: semilabile dissolved organic matter (with timescales of several weeks to several years), small, and large sinking particles. As with the living components, constant Redfield ratios are imposed for C/N/P while the iron, silicon and calcite pools of the particles are fully simulated and, their ratios are allowed to vary relative to organic carbon. Nutrients are supplied from three sources: atmospheric dust deposition, rivers, and sediment mobilization. The global PISCES configuration has been used to study past climates, to understand the interannual variability in marine productivity or ocean-atmosphere carbon fluxes, to assess the impact of climate change or ocean acidification on marine ecosystems and air-sea carbon fluxes, to evaluate geo-engineering strategies to mitigate climate change.

The biogeochemical model was run using the same constraints as the physical model including lateral mixing on isoneutral surfaces. The model was spun up offline for 3000 years using state variables from the final year of the spinup of LIM-ORCA2. Pre-anthropogenic boundary conditions were maintained for sea surface atmospheric CO₂ concentrations. This long spinup run was split in 1860. A control run continued through to the end of 1999 with pre-industrial atmospheric CO₂ concentrations. A second run followed the observed atmospheric CO₂ transient. As before, we refer to the first run as “natural” carbon and the second run as “contemporary” carbon. The difference in carbon variables between the runs, our focus here, is “anthropogenic carbon”. The modelled Cant inventory compares rather well with GLODAP and TTD estimates. The model has a deep ocean Cant inventory that is higher than GLODAP, but less than TTD. The meridional transport at 30°S also compares well with published oceanic inversions (see Table 1).

## Theory

**Tracer dynamical budget in a water mass framework.** Detailed theoretical derivation are presented in the Supplementary Note.

**The conceptual model: assumptions and analytical solutions.** We developed a conceptual ocean model having n density layers. This ocean was exposed to an atmospheric tracer boundary condition that grew exponentially with time. The density layers were internally homogeneous.

A steady state ocean circulation implies that the net overturning convergence is zero (volumes neither inflate nor deflate). Here we only considered the subsurface (in the MWML, generally about 100–200 m) expression of the overturning and it is divergent. We discarded the overturning component in the interior since, as illustrated in the GCM analysis, the tracer concentration at depth is much lower than at surface. The drastic assumption of well-mixed layers is at odds with the real ocean, but this simplification contributes to process-attribution. In practice we are assuming that: (1) most of the convergence occurs at the surface, in the mixed layer; (2) the layer specific inventory and air-sea flux properties are controlled not only by the mixed layer but also by the internal volume. This is acceptable especially for layers that are not completely saturated but have a significant vertical exchange (i.e., subduction) and thus for the various mode and intermediate waters which are the focus of this study. We also disregarded the impact of carbonate chemistry, such as solubility effects, and variability in the surface mixed layer thickness.

The resulting equation is solved numerically. For diapycnal transport ω, we used an analytical description of the diapycnal transport that approximately matches the transport in Fig. 4a. For dense water formation (σ₀ > 27.5) we amplified overturning by a factor 7.5 to compensate for underestimation of surface values for water masses with very large volume. This factor was chosen to represent adiapycnal tracer transport (Supplementary Fig. 3) similar to the Cant overturning in Fig. 4b. Our choice has only a minimal impact on mode and intermediate water masses and on the air-sea fluxes. The two peaks in Fig. 4b correspond to SOC and DW formation (essentially NADW).

Finally, valuable insight was achieved by evaluating analytical solutions, which can be easily derived for the AIRSEA case, that is, here we considered a situation where only air-sea exchange was allowed (βᵢ = 0 and ωᵢ = 0 ∀i). We obtained an analytically result by seeking a solution of the form cᵢ = c₀γᵢe^{t/τ-atm}.

This yielded

$$ \gamma_i = \frac{e^{t/\tau_{atm}}}{\left(1 + \left(1 + \frac{\tau'_i}{\tau_{atm}}\right) e^{t/\tau_{atm}}\right)} \quad (10) $$

where τ'ᵢ = τᵢVᵢ/HᵢAᵢ is the relaxation time scale rescaled by the ratio of the water mass volume to the upper boundary layer (AᵢHᵢ) volume. For a fast adjustment time scale τᵢ, γᵢ = 1 and the water mass concentration perfectly tracks the atmospheric one (c = c_atm).

The shape factor γᵢ is also a function of time, but after a relatively short transient period (when e^{t/τ-atm} ≫ 1) the solution becomes

$$ \gamma_i = \frac{1}{\left(1 + \frac{\tau'_i}{\tau_{atm}}\right)} \quad (11) $$

Thus, for all the water masses the carbon content increases exponentially with the same doubling time as the atmosphere. Relative repartition of carbon among the water masses is constant in time and equal to Eq. 11.
Similarly, the surface fluxes increase exponentially:

$$ \phi_i = c_0 \frac{\gamma_i V_i}{\tau_{atm}} e^{t/\tau_{atm}} \quad (12) $$

In this instance the repartitioning is not time dependent but rather varies with water masses characteristics (volume, surface area, etc). The surface flux is in fact a function of the product of the water volume, that sharply increases with density, by γᵢ, which is constant for the water masses up to the subtropical densities (τ' < τ_atm) and decrease to zero with increasing density. This decrease occurs because the volume strongly increases with density while the surface area is relatively constant for intermediate water masses and goes nearly to zero for the deep ones; Fig. 2). The combination of an increase of the volume and decrease of γᵢ with increasing density create a maximum of air-sea flux in correspondence of the intermediate water masses that is due only to geometrical factors that in turn depend on oceanic circulation patterns⁴⁹,⁵⁰.

## References

1. Sabine, C. L. et al. The ocean sink for anthropogenic CO₂. *Science* **305**, 367–371 (2004).

2. Khatiwala, S., Primeau, F. & Hall, T. Reconstruction of the history of anthropogenic CO₂ concentrations in the ocean. *Nature* **462**, 346–349, doi: 10.1038/nature08526 (2009).

3. Sarmiento, J. L., Orr, J. C. & Siegenthaler, U. A Perturbation Simulation of CO₂ Uptake in an Ocean General Circulation Model. *J. Geophys. Res.* **97**(C3), 3621–3645 (1992).

4. Gloor, M. et al. A first estimate of present and preindustrial air-sea CO₂ flux patterns based on ocean interior carbon measurements and models. *Geophys. Res. Lett.* **30**(1), 10–1 (2003).

5. Mikaloff Fletcher, S. E. et al. Inverse estimates of anthropogenic CO₂ uptake, transport, and storage by the ocean. *Global Biogeochem. Cycles*, **20**(2), GB2002, doi: 10.1029/2005GB002530 (2006).

6. Gruber, N. et al. Oceanic sources, sinks, and transport of atmospheric CO₂. *Global Biogeochem. Cycles* **23**, GB1005, doi: 10.1029/2008GB003349 (2009).

7. Khatiwala, S. et al. Global ocean storage of anthropogenic carbon, *Biogeosciences* **10**, 2169–2191 (2013).

8. Winton, M., Griffies, S. M., Samuels, B. L., Sarmiento, J. L. & Frölicher, T. L. Connecting Changing Ocean Circulation with Changing Climate. *J. Climate* **26**, 2268–2278, doi: 10.1175/JCLI-D-12-00296.1 (2013).

9. Frölicher, T. L. et al. Dominance of the Southern Ocean in Anthropogenic Carbon and Heat Uptake in CMIP5 Models, *Journal of Climate* **28**(2), 862–886 (2015).

10. Doney, S. C. et al. Evaluating global ocean carbon models: The importance of realistic physics. *Global Biogeochem. Cycles* **18**, GB3017, doi: 10.1029/2003GB002150 (2004).

11. Mignone, B. K., Gnanadesikan, A., Sarmiento, J. L. & Slater, R. D. Central role of Southern Hemisphere winds and eddies in modulating the oceanic uptake of anthropogenic carbon. *Geophys. Res. Lett.* **33**, L01604, doi: 10.1029/2005GL024464 (2006).

12. Ito, T., Woloszyn, M. & Mazloff, M. Anthropogenic carbon dioxide transport in the southern ocean driven by Ekman flow. *Nature* **463**, 80–83, doi: 0.1038/nature08687 (2010).

13. Levine, N. M. et al. The impact of the North Atlantic Oscillation on the uptake and accumulation of anthropogenic CO₂ by North Atlantic Ocean mode waters. *Global Biogeochem. Cycles* **25**(3), GB3022, doi: 10.1029/2010GB003892 (2011).

14. Talley, L. D. Shallow, intermediate, and deep overturning components of the global heat budget. *J. Phys. Oceanogr.* **33**(3), 530–560 (2003).

15. Waugh, D. W., Hall, T. M., McNeil, B. I., Key, R. & Matear, R. J. Anthropogenic CO₂ in the oceans estimated using transit-time distributions. *Tellus*, **58**, B:376-390, doi: 10.1111/j.1600-0889.2006.00222.x (2006).

16. DeVries (2014), The oceanic anthropogenic CO₂ sink: Storage, air-sea fluxes and transports over the industrial era, *Global Biogeochem. Cycles* **28**, doi: 10.1002/2013GB004739 (2014).

17. Walin, G. On the relation between sea-surface heat flow and thermal circulation in the ocean. *Tellus* **34**, 187–195, doi: 10.1111/j.2153-3490.1982.tb01806.x (1982).

18. Nurser, A. J. G., Marsh, R. & Williams, R. G. Diagnosing water mass formation from air-sea fluxes and surface mixing. *J. Phys. Oceanogr.* **29**, 1468–1487 (1999).

19. Marshall, J., Jamous, D. & Nilsson, J. Reconciling thermodynamic and dynamic methods of computation of water-mass transformation rates. *Deep-Sea Res. I* **46**, 545–572 (1999).

20. Iudicone, D., Madec, G. & McDougall, T. J. Water-mass transformations in a neutral density framework and the key role of light penetration. *J. Phys. Oceanogr.* **38**(7), 1357–1376 (2008).

21. Iudicone, D. et al. Water masses as a unifying framework for understanding the Southern Ocean Carbon Cycle. *Biogeosciences* **8**, 1031–1052 (2011).

22. Broecker, W. S. Radioisotopes and large-scale oceanic mixing. In *The Sea* Vol. 2, Hill, M. N. editor, Interscience Pub., NY, 88–108, (1963).

23. Oeschger, H. U., Siegenthaler, U., Schotterer, U. & Gugelman A. A box diffusion model to study the carbon dioxide exchange in nature. *Tellus* **27**, 168–192 (1975).

24. Ito, T., Marshall, J. C. & Follows, M. J. What controls the uptake of transient tracers in the Southern Ocean? *Global Biogeochem. Cycles* **18**, GB2021, doi: 10.1029/2003GB002103 (2004).

25. Sarmiento, J. L. & Gruber, N. *Ocean biogeochemical dynamics*. Princeton University Press, 2006.

26. Key, R. et al. A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP). *Global Biogeochem. Cycles* **18**, GB4031, doi: 10.1029/2004GB002247 (2004).

27. Stommel, H. M. Determination of water mass properties of water pumped down from the Ekman layer to the geostrophic flow below. *Proc. Natl. Acad. Sci. USA* **76**, 3051–3055 (1979).

28. Williams, R. G. The role of the mixed layer in setting the potential vorticity of the main thermocline. *J. Phys. Oceanogr.* **21**, 1802–1814 (1991).

29. Iudicone, D., Madec, G., Blanke, B. & Speich, S. The role of Southern Ocean surface forcings and mixing in the global conveyor. *J. Phys. Oceanogr.* **38**, 1377–1400 (2008).

30. Broecker, W. S. "NO" A conservative water-mass tracer. *Earth Planet. Sci. Lett.* **23**, 100–107 (1974).

31. Madec, G., Delecluse, P., Imbard, M. & Lévy, C. OPA 8.1 Ocean General Circulation Model reference manual. *Note du Pole de modélisation, Institut Pierre-Simon Laplace (IPSL)*, France, 11, p. 91 (1998).

32. Aumont, O. & Bopp, L. Globalizing results from ocean *in situ* iron fertilization studies. *Global Biogeochem. Cycles* **20**, GB2017, doi: 10.1029/2005GB002591 (2006).

33. Nakano, H., Ishii, M., Rodgers, K. B., Tsujino, H. & Yamanaka, G. Oceanic CO₂ uptake, transport, storage, and dynamical controls imposed by the thermocline overturning circulation: a modeling study, *Global Biogeochem. Cycles*, **29**, doi: 10.1002/2015GB005128 (2015).
34. Key, R. M. et al. A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAPv2), CDIAC NDP-093, 2015.

35. Vallis, G. K. *Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation*. Cambridge University Press, (2006).

36. Bopp, L., Lévy, M., Resplandy, L. & Sallée, J.-B. Pathways of anthropogenic carbon subduction in the global ocean. *Geophys. Res. Lett.* **42**, doi: 10.1002/2015GL065073 (2015).

37. Sallée, J. B., Matear, R. J., Rintoul, S. R. & Lenton, A. Localized subduction of anthropogenic carbon dioxide in the Southern Hemisphere oceans. *Nature Geoscience* doi: 10.1038/ngeo1523 (2012).

38. Cerovecki, I., Talley, L. D., Mazloff, M. R. & Maze, G. Subantarctic Mode Water formation, destruction and export in the eddy-permitting Southern Ocean State Estimate. *J. Phys. Oceanogr.* **43**, 1485–1511 (2013).

39. McNeil, B. I., Tilbrook, B. & Matear, R. J. Accumulation and uptake of anthropogenic CO₂ in the Southern Ocean, South of Australia between 1968 and 1996: The SAZ Project. *J. Geophys. Res.* **106**, 31431–31445 (2001).

40. Waugh, D., Primeau, F., DeVries, T. & Holzer, M. Recent changes in the ventilation of the Southern Oceans. *Science* **339**, 568–570, doi: 10.1126/science.1225411 (2013).

41. Locarnini, R. A. et al. *World Ocean Atlas 2009, Volume 1: Temperature*. Levitus, S. Ed. NOAA Atlas NESDIS 68, U.S. Government Printing Office, Washington, DC, 184 pp (2010).

42. Siedler, G., Church, J., Gould, J. eds. *Ocean Circulation and Climate: Observing and Modeling the Global Ocean*. San Diego: Academic. 715 pp (2001).

43. de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A. & Iudicone, D. Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. *J. Geophys. Res.* **109**, C12003, doi: 10.1029/2004JC002378 (2004).

44. Gent, P. R. & McWilliams, J. C. Isopycnal mixing in ocean circulation models. *J. of Phys. Oceanogr.* **20**(1), 150–155 (1990).

45. Treguier, A. M., Held, I. M. & Larichev, V. D. Parameterization of quasigeostrophic eddies in primitive equation ocean models. *J. of Phys. Oceanogr.* **27**(4), 567–580 (1997).

46. Blanke, B. & Delecluse, P. Variability of the tropical Atlantic Ocean simulated by a general circulation model with two different mixed-layer physics. *J. Phys. Oceanogr.* **23**(7), 1363–1388 (1993).

47. Beckmann, A. & Doscher, R. A method for improved representation of dense water spreading over topography in geopotential-coordinate models. *J. Phys. Oceanogr.* **27**(4), 581–591 (1997).

48. Timmermann, R. et al. On the representation of high latitude processes in the ORCA-LIM global coupled sea ice ocean model. *Ocean Modelling* **8**(1), 175–201 (2005).

49. Bolin, B. & Eriksson, E. Changes in the carbon dioxide content of the atmosphere and sea due to fossil fuel combustion, p. 130–142, In *The Atmosphere and the Sea in Motion: Scientific Contributions to the Rossby Memorial Volume* Bert, B. ed., New York, Rockefeller Institute Press, (1958).

50. Khatiwala *et al*. Reconstruction of the history of anthropogenic CO₂ concentrations in the ocean. *Nature*, 346–349 (2009).

## Acknowledgements

D.I. was partially funded by the Flagship Project RITMARE – The Italian Research for the Sea – funded by the Italian Ministry of Education, University, and Research within the National Research Program 2011–2013. This material is based upon work supported by the National Aeronautics and Space Administration under Grant No. NNX14AL85G (K.B.R.). K.B.R. was supported by award NA17RJ2612 and NA08OAR4320752 from the National Oceanic and Atmospheric Administration, US Department of Commerce. The statements, findings, conclusions, and recommendations are those of the author(s) and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration, or the US Department of Commerce. K.B.R. was also supported by the National Science Foundation under Grant No. OCE-1155983. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Y.P. acknowledges support from U.K. NERC IRF NE/M017826/1. R. Key was supported by KeyCrafts grant 2012-001, CICS grants NA08OAR4320752 and NA14OAR4320106, NASA grant NNX12AQ22G, NSF grants OCE-0825163 (with a supplement via WHOI P.O. C119245) and PLR-1425989, and DOE UT-Battelle contract #4000133565 to CDIAC.

## Author Contributions

All authors have contributed extensively to the work. D.I. and K.B.R. conceived the research. D.I. performed the theoretical and diagnostic work. D.I. and K.B.R. wrote the article. O.A and G.M. contributed to the research design and performed the numerical integrations. Y.P. and R.K. contributed to the observational data analysis and the interpretation of the results. T.I. provided numerical results and interpretation. M.I. supervised the data analyses and interpretation. All authors assisted in the writing process and revised the manuscript.

## Additional Information

**Supplementary information** accompanies this paper at http://www.nature.com/srep

**Competing financial interests:** The authors declare no competing financial interests.

**How to cite this article:** Iudicone, D. et al. The formation of the ocean's anthropogenic carbon reservoir. *Sci. Rep.* **6**, 35473; doi: 10.1038/srep35473 (2016).

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_13.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

© The Author(s) 2016