Ocean carbon cycle feedbacks in CMIP6 models: contributions from different basins

The ocean response to carbon emissions involves the combined effect of an increase in atmospheric CO₂, acting to enhance the ocean carbon storage, and climate change, acting to decrease the ocean carbon storage. This ocean response can be characterised in terms of a carbon–concentration feedback and a carbon–climate feedback. The contribution from different ocean basins to these feedbacks on centennial timescales is explored using diagnostics of ocean carbonate chemistry, physical ventilation and biological processes in 11 CMIP6 Earth system models. To gain mechanistic insight, the dependence of these feedbacks on the Atlantic Meridional Overturning Circulation (AMOC) is also investigated in an idealised climate model and the CMIP6 models. For the carbon–concentration feedback, the Atlantic, Pacific and Southern oceans provide comparable contributions when estimated in terms of the volume-integrated carbon storage. This large contribution from the Atlantic Ocean relative to its size is due to strong local physical ventilation and an influx of carbon transported from the Southern Ocean. The Southern Ocean has large anthropogenic carbon uptake from the atmosphere, but its contribution to the carbon storage is relatively small due to large carbon transport to the other basins. For the carbon–climate feedback estimated in terms of carbon storage, the Atlantic and Arctic oceans provide the largest contributions relative to their size. In the Atlantic, this large contribution is primarily due to climate change acting to reduce the physical ventilation. In the Arctic, this large contribution is associated with a large warming per unit volume. The Southern Ocean provides a relatively small contribution to the carbon–climate feedback, due to competition between the climate effects of a decrease in solubility and physical ventilation and an increase in accumulation of regenerated carbon. The more poorly ventilated Indo-Pacific Ocean provides a small contribution to the carbon cycle feedbacks relative to its size. In the Atlantic Ocean, the carbon cycle feedbacks strongly depend on the AMOC strength and its weakening with warming. In the Arctic, there is a moderate correlation between the AMOC weakening and the carbon–climate feedback that is related to changes in carbonate chemistry. In the Pacific, Indian and Southern oceans, there is no clear correlation between the AMOC and the carbon cycle feedbacks, suggesting that other processes control the ocean ventilation and carbon storage there.

1 Introduction

Carbon emissions drive an Earth system response via direct changes in the biogeochemical carbon cycle and the physical climate. These changes in the biogeochemical carbon cycle and physical climate further amplify or dampen the Earth system response, with this amplification often referred to as a feedback (Sherwood et al., 2015). The physical climate feedback involves the combined effect from changes in atmospheric water vapour, tropospheric lapse rate, surface albedo and clouds (Ceppi and Gregory, 2017) and from a shift in the regional patterns of ocean heat uptake due to changes in the ocean circulation (Winton et al., 2013). For the carbon cycle, an initial increase in atmospheric CO₂ leads to carbon uptake and storage in land and ocean reservoirs. This response of the carbon cycle to the increase in atmospheric CO₂ is characterised by the carbon–concentration feedback. At the same time, the carbon cycle is further modified by changes in the physical climate, such as warming and an increase in ocean stratification leading to an amplification of the initial increase in atmospheric CO₂. This response of the carbon cycle to changes in the physical climate is characterised by the carbon–climate feedback. These two carbon cycle feedbacks have been extensively used to understand and quantify the response of the global carbon cycle to carbon emissions (Friedlingstein et al., 2003, 2006; Gregory et al., 2009; Boer and Arora, 2009; Arora et al., 2013; Schwinger et al., 2014; Schwinger and Tjiputra, 2018; Williams et al., 2019; Arora et al., 2020). A regional extension of the carbon cycle feedbacks has been also used to explore their geographical distribution and the mechanisms that control the land and ocean carbon uptake and storage in difference regions (Yoshikawa et al., 2008; Boer and Arora, 2010; Tjiputra et al., 2010; Roy et al., 2011).

On a global scale, the carbon–concentration feedback is of comparable strength over the land and ocean, while the carbon–climate feedback is about 3 times stronger over the land than over the ocean on centennial timescales in the CMIP6 Earth system models (Arora et al., 2020). However, there is a substantial geographical variation in the ocean carbon–climate feedback (Tjiputra et al., 2010; Roy et al., 2011) as a result of an interplay between the effect of carbonate chemistry, physical ventilation and biological processes. In the tropics, the carbonate chemistry and the decrease in solubility with warming drives a reduction in the ocean carbon uptake with climate change (Roy et al., 2011; Rodgers et al., 2020). In the North Atlantic, the physical ventilation and its weakening with warming acts to further reduce the ocean carbon uptake with climate change (Yoshikawa et al., 2008; Tjiputra et al., 2010; Roy et al., 2011). In the Southern Ocean, changes in the cycling of biological material with climate change can partly counteract the reduction in the ocean carbon uptake due to the decrease in solubility and physical ventilation with warming (Sarmiento et al., 1998; Bernardello et al., 2014).

The ocean carbon cycle feedbacks can be defined in terms of either the cumulative ocean carbon uptake or the ocean carbon storage (Schwinger et al., 2014; Arora et al., 2020). For the global ocean, these two definition are almost equivalent apart from a small contribution from the land-to-ocean carbon flux from river runoff and the carbon burial in ocean sediments (Arora et al., 2020). However, on a regional scale these two definitions are different, as the ocean carbon storage explicitly includes the convergence of transport of carbon by the ocean circulation. This transport effect leads to different geographical patterns for the ocean carbon storage and the ocean cumulative carbon uptake (Frölicher et al., 2015). This transport effect also leads to a broadly similar geographical distribution for ocean carbon and heat storage, with the “redistribution” of the pre-industrial carbon and heat by changes in the circulation with warming driving a second-order asymmetry between the regional patterns of heat and carbon storage (Winton et al., 2013; Bronselaer and Zanna, 2020; Williams et al., 2021). The combined air–sea transfer and transport effect leads to the Atlantic, Pacific and Southern oceans each storing about 25 %–30 % of the additional heat and carbon in CMIP6 models for a quadrupling of atmospheric CO₂ (Fig. 1a and b), despite their different sizes (see Supplement Fig. S1 for the basins' definition). The Atlantic and Arctic oceans have the largest increase in carbon and heat per unit volume, as given by the dissolved inorganic carbon and temperature (Fig. 1c and d). The Pacific Ocean has the smallest increase in carbon and heat per unit volume (Fig. 1c and d). Our motivation is to explore the mechanisms that lead to these regional variations in carbon storage and carbon cycle feedbacks in the different ocean basins in CMIP6 models.

Figure 1 Carbon and heat storage for the global ocean and different ocean basins in CMIP6 Earth system models: (a) ocean carbon content changes relative to the pre-industrial in petagrams of carbon; (b) ocean heat content changes relative to the pre-industrial in joules; (c) changes in the ocean dissolved inorganic carbon relative to the pre-industrial in moles per cubic metre, expressing the ocean carbon storage changes per volume; and (d) changes in ocean temperature relative to the pre-industrial in degrees Celsius, expressing the ocean heat changes per volume. The solid lines show the model mean and the shading the model range based on the 1 % yr⁻¹ increasing CO₂ experiment over 140 years in 11 CMIP6 models (Table 1). For the definition of the ocean basins see Supplement Fig. S1.

A mechanism that can affect the regional carbon storage is the Atlantic Meridional Overturning Circulation (AMOC). The projected weakening in the AMOC with climate change (Cheng et al., 2013) weakens the ocean physical ventilation and transport of carbon into the ocean interior, which acts to reduce the ocean carbon uptake and storage (Sarmiento and Le Quéré, 1996; Crueger et al., 2008). The weakening in the AMOC with climate change also increases the residence time in the ocean interior and the accumulation of remineralised carbon at depth, which acts to increase the ocean carbon uptake and storage (Sarmiento and Le Quéré, 1996; Joos et al., 1999; Schwinger et al., 2014; Bernardello et al., 2014). Previous studies suggest that the combined effect of these two competing processes leads to a modest reduction in ocean carbon uptake and storage with AMOC weakening and to an ocean carbon–climate feedback that amplifies the increase in atmospheric CO₂ (Sarmiento and Le Quéré, 1996; Joos et al., 1999; Crueger et al., 2008; Schwinger et al., 2014). However, the net effect of AMOC weakening with climate change on the carbon storage is highly uncertain and sensitive to the representation of the vertical carbon gradient and ocean biological processes in Earth system models. This uncertainty motivates us to explore the control of the AMOC on the carbon cycle feedbacks in CMIP6 models, as well as the relative importance of changes in biological processes and physical ventilation for the carbon storage in different ocean basins.

Our aim is to provide insight into the relative contribution of different ocean basins to the ocean carbon cycle feedbacks and the processes that drive this regional partitioning in the CMIP6 models. In Sect. 2, we provide the framework for the carbon cycle feedbacks and explore their geographical distribution in 11 CMIP6 Earth system models (Table 1). In Sect. 3, the ocean carbon cycle feedbacks are separated into contribution from carbonate chemistry, physical ventilation and biological processes, and the controls of the global and regional feedbacks are investigated in diagnostics of CMIP6 models. In Sect. 4, the effect of the AMOC on the global and basin-scale carbon cycle feedbacks is investigated, firstly using an idealised climate model that provides a mechanistic insight and then in diagnostics of CMIP6 models. Section 5 summarises our conclusions and discusses the wider context of our analysis.

Ziehn et al. (2020)Swart et al. (2019g)Swart et al. (2019g)Séférian et al. (2019)Dunne et al. (2020)Boucher et al. (2020)Hajima et al. (2020)Mauritsen et al. (2019)Yukimoto et al. (2019d)Seland et al. (2020)Sellar et al. (2019)

Table 1 List of the 11 CMIP6 Earth system models used in this study along with references for the model description.

2 Ocean carbon cycle feedbacks in CMIP6 models

2.1  Global ocean

In the carbon cycle feedback framework introduced by Friedlingstein et al. (2003, 2006) the ocean carbon gain due to anthropogenic carbon emissions, ΔI_ocean, is expressed as a function, F, of changes in the atmospheric CO₂ and the physical climate:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{I}_{\mathrm{ocean}}=F({\mathrm{CO}}_{\mathrm{2},\mathrm{0}}+\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}},T_{\mathrm{0}}+\mathrm{\Delta }T)-F({\mathrm{CO}}_{\mathrm{2},\mathrm{0}},T_{\mathrm{0}}),\end{array}$$

where the surface air temperature, T, is used as a proxy for the physical climate and subscript 0 denotes the pre-industrial state. By expanding the function F into a Taylor series (Schwinger et al., 2014; Williams et al., 2019), the ocean carbon gain relative to the pre-industrial era, ΔI_ocean, is expressed as

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }{I}_{\mathrm{ocean}}& {\displaystyle }={\left.{\displaystyle \frac{\partial F}{\partial {\mathrm{CO}}_{\mathrm{2}}}}\right|}_{\mathrm{0}}\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}+{\left.{\displaystyle \frac{\partial F}{\partial T}}\right|}_{\mathrm{0}}\mathrm{\Delta }T\\ \text{(2)}& {\displaystyle }& {\displaystyle }+{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}F}{\partial {\mathrm{CO}}_{\mathrm{2}}\partial T}}\right|}_{\mathrm{0}}\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\mathrm{\Delta }T\text{(3)}& {\displaystyle }& {\displaystyle }+{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}F}{\partial {\mathrm{CO}}_{\mathrm{2}}^{\mathrm{2}}}}\right|}_{\mathrm{0}}\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{2}}+{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}F}{\partial T^{\mathrm{2}}}}\right|}_{\mathrm{0}}\mathrm{\Delta }{T}^{\mathrm{2}}+R^{\mathrm{3}},\end{array}$$

where R³ contains the third-order and higher derivatives. By ignoring the second- and higher-order terms but keeping the terms for the non-linear relationship between atmospheric CO₂ and climate change, Eq. (3) is rewritten as

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{I}_{\mathrm{ocean}}=\mathit{\beta }\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}+\mathit{\gamma }\mathrm{\Delta }T+\mathcal{N}\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\mathrm{\Delta }T,\end{array}$$

where the ocean carbon–concentration feedback parameter is defined as $\mathit{\beta }={\left.\frac{\partial F}{\partial {\mathrm{CO}}_{\mathrm{2}}}\right|}_{\mathrm{0}}$, the ocean carbon–climate feedback parameter is defined as $\mathit{\gamma }={\left.\frac{\partial F}{\partial T}\right|}_{\mathrm{0}}$ and the non-linearity of the ocean carbon cycle feedbacks is defined as $\mathcal{N}={\left.\frac{{\partial }^{\mathrm{2}}F}{\partial {\mathrm{CO}}_{\mathrm{2}}\partial T}\right|}_{\mathrm{0}}$.

The carbon cycle feedback parameters, β and γ, are traditionally estimated using Earth system model simulations with the couplings between the carbon cycle and radiative forcing switched either on or off: a fully coupled simulation, a radiatively coupled simulation and a biogeochemically coupled simulation (Friedlingstein et al., 2006; Arora et al., 2013; Jones et al., 2016; Arora et al., 2020). Any combination of these three simulations can be used to estimate the carbon cycle feedback parameters; however, each combination yields somewhat different results due to the non-linearity of the system (Gregory et al., 2009; Zickfeld et al., 2011; Schwinger et al., 2014; Arora et al., 2020). Here, we estimate the carbon cycle feedback parameters using the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment, in which the atmospheric CO₂ concentration increases from its pre-industrial value of around 285 ppm until it quadruples over a 140-year period, following the recommended C⁴MIP protocol of experiments (Jones et al., 2016; Arora et al., 2020). To remove the effect of model drift and reduce model biases, the pre-industrial control simulation (piControl) was used to estimate the pre-industrial state in the Earth system models. For simplicity, we ignore the effect of the air temperature increase in the BGC simulation on the feedbacks, which has a contribution of less than 5 % (Arora et al., 2020), such that

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathit{\beta }={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{ocean}}^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & \mathit{\gamma }={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{ocean}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{ocean}}^{\mathrm{BGC}}}{\mathrm{\Delta }T}},\end{array}\end{array}$$

where ΔCO₂ is the increase in atmospheric CO₂ and ΔT is the increase in air surface temperature in the fully coupled Earth system (i.e. COU simulation) relative to the pre-industrial state. The carbon–climate feedback parameter, γ, in Eq. (5) corresponds to the effect of climate change under rising atmospheric CO₂ and hence includes the effect of the non-linearity, 𝒩ΔCO₂, in Eq. (4) (Schwinger et al., 2014).

The ocean carbon–concentration feedback parameter, β, is positive in all the CMIP6 models (Table 2). The ocean carbon–climate feedback parameter, γ, is negative in all the Earth system models (Table 2), indicating that the ocean takes up less carbon in response to climate change. The variability in β amongst the Earth system models, as described by the coefficient of variation, CV, is relatively small on the global scale (CV = 0.09) when compared with the variability in γ (CV = 0.43) (Table 2). However, for the uncertainty in the ocean carbon gain due to carbon emissions, the carbon–concentration feedback contributes to a spread of 62 PgC, while the carbon–climate feedback contributes only to a spread of 25 PgC amongst the CMIP6 Earth system models on a global scale and to a quadrupling of atmospheric CO₂, where the spread corresponds to 1 standard deviation.

Table 2 Carbon–concentration feedback parameter based on carbon storage, β^* (PgC ppm⁻¹), and carbon–climate feedback parameter based on carbon storage, γ^* (PgC K⁻¹), for the global ocean and different ocean basins in 11 CMIP6 Earth system models, along with the inter-model mean; standard deviation; and coefficient of variation (CV), estimated as the standard deviation divided by the mean. The estimates are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂). Note that for the global ocean, β^* and γ^* are equivalent to β and γ.

2.2 Regional ocean

The carbon cycle feedbacks for the global ocean in Eq. (4) can be further separated into contributions from different ocean regions such that

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{I}_{\mathrm{ocean}}=\sum _{n=\mathrm{1}}^{\mathrm{global}}{\mathit{\beta }}_n\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}+\sum _{n=\mathrm{1}}^{\mathrm{global}}{\mathit{\gamma }}_n\mathrm{\Delta }T+\sum _{n=\mathrm{1}}^{\mathrm{global}}{\mathcal{N}}_n\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}\mathrm{\Delta }T,\end{array}$$

where n denotes the different ocean regions; ΔCO₂ and ΔT are the global changes in atmospheric CO₂ and the surface air temperature, respectively; and β_n, γ_n and 𝒩_n are the carbon cycle feedback parameters and their non-linearity for each ocean region n.

Traditionally, the carbon cycle feedbacks are defined based on the cumulative carbon uptake from the atmosphere, Δ𝒮 (Tjiputra et al., 2010; Roy et al., 2011) such that the carbon cycle feedback parameters for each ocean region n are expressed as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\mathit{\beta }}_n={\displaystyle \frac{\mathrm{\Delta }{\mathcal{S}}_n^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & {\mathit{\gamma }}_n={\displaystyle \frac{\mathrm{\Delta }{\mathcal{S}}_n^{\mathrm{COU}}-\mathrm{\Delta }{S}_n^{\mathrm{BGC}}}{\mathrm{\Delta }T}},\end{array}\end{array}$$

where the carbon–climate feedback parameter, γ_n, corresponds to the effect of climate change under rising atmospheric CO₂ and hence includes the effect of the non-linearity, 𝒩_nΔCO₂, in Eq. (6).

Alternatively, the carbon cycle feedbacks can be defined based on the carbon storage, such that the regional carbon cycle feedback parameters for each ocean region n are expressed as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathit{\beta }}_n^{*}={\displaystyle \frac{\mathrm{\Delta }{I}_n^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & {\mathit{\gamma }}_n^{*}={\displaystyle \frac{\mathrm{\Delta }{I}_n^{\mathrm{COU}}-\mathrm{\Delta }{I}_n^{\mathrm{BGC}}}{\mathrm{\Delta }T}},\end{array}\end{array}$$

where ΔI_n is the change in the carbon inventory in region n. The carbon storage includes the combined effect from the local air–sea carbon exchange and the transport of carbon by the ocean circulation, such that

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{I}_n=\mathrm{\Delta }{\mathcal{S}}_n+\mathrm{\Delta }{\mathcal{G}}_n,\end{array}$$

where Δ𝒮_n is the regional cumulative ocean carbon uptake from the atmosphere relative to the pre-industrial era, and Δ𝒢_n is the regional cumulative carbon gain relative to the pre-industrial era due to the ocean transport.

By substituting Eqs. (9) and (7) in Eq. (8) the two definitions for the carbon cycle feedback parameters are related by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{\beta }}_n^{*}={\mathit{\beta }}_n+{\displaystyle \frac{\mathrm{\Delta }{\mathcal{G}}_n^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & {\mathit{\gamma }}_n^{*}={\mathit{\gamma }}_n+\left({\displaystyle \frac{\mathrm{\Delta }{\mathcal{G}}_n^{\mathrm{COU}}-\mathrm{\Delta }{\mathcal{G}}_n^{\mathrm{BGC}}}{\mathrm{\Delta }T}}\right).\end{array}\end{array}$$

Eq. (10) shows that the feedback parameters defined by the regional carbon storage, β^* and γ^*, are proportional to the feedback parameters defined by the regional cumulative carbon uptake, γ and β, but further modified by the ocean carbon transport.

2.2.1 Estimates based on carbon storage versus carbon uptake

On a global scale, the transport effect on the carbon storage integrates to zero, ${\sum }_n^{\mathrm{global}}\mathrm{\Delta }{\mathcal{G}}_n=\mathrm{0}$, such that the feedback parameters estimated from the carbon storage, β^* and γ^*, are equivalent to the feedback parameters estimated from the cumulative carbon uptake, β and γ, when ignoring the small carbon exchange between land and ocean. However, on regional scales the effect of the ocean transport on carbon storage leads to different spatial patterns in these carbon cycle feedback parameters, β and β^* and γ and γ^* (Fig. 2).

Figure 2 Geographical distribution of the carbon cycle feedback parameters normalised by area: (a) β^* (gC ppm⁻¹ m⁻²) and γ^* (gC K⁻¹ m⁻²), estimated based on the regional ocean carbon storage, and (b) β (gC ppm⁻¹ m⁻²) and γ (gC K⁻¹ m⁻²), estimated based on the regional cumulative ocean carbon uptake. Results are shown as the inter-model mean and standard deviation based on 11 CMIP6 Earth system models (Table 1). The estimates are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂).

The carbon–concentration feedback parameter estimated from the cumulative carbon uptake, β, is largest and has more inter-model variability in (i) the Southern Ocean, (ii) the eastern boundary upwelling regions, (iii) the Gulf Stream and its extension into the North Atlantic Current, and (iv) the Kuroshio Extension (Fig. 2b). The inter-model variability in β is also significant along the equatorial Pacific, with this variability related to the inter-model spread in the trade winds and equatorial upwelling. In contrast, in the subtropical gyres, the ocean anthropogenic carbon uptake is limited and β is small (Fig. 2b).

The carbon–concentration feedback parameter estimated from ocean carbon storage, β^*, is again large in the North Atlantic but instead large in the Southern Hemisphere subtropical gyres and small in the Southern Ocean south of 50^∘ S relative to β (Fig. 2a). This difference between β and β^* (Supplement Fig. S2) is due to the northward transport of anthropogenic carbon from the Southern Ocean associated with subduction and transport of mode and intermediate waters. The variability in β^* amongst the models is large in the North Atlantic and extends south along the Atlantic western boundary (Fig. 2a).

The carbon–climate feedback parameter estimated from the cumulative carbon uptake, γ, is large and negative in the North Atlantic and in the Southern Ocean from 50 to 65^∘ S (Fig. 2b). In contrast, γ is large and positive in a narrow band between 40 and 45^∘ S, in the Southern Hemisphere eastern boundary upwelling regions, and in polar regions with sea ice. The regions of large ocean carbon loss or uptake from the atmosphere due to climate change, as shown by the large γ, also experience the largest variability in γ amongst the CMIP6 Earth system models (Fig. 2b).

The effect of carbon transport on γ^* is of opposite sign to the effect of the cumulative carbon uptake in most regions (Supplement Fig. S2). This transport effect leads to a γ^* that is overall less negative in the Southern Ocean and in the high latitudes of the North Atlantic but more negative in the Arctic, the equatorial Pacific and along the Atlantic western boundary relative to γ (Fig. 2a). The spread in γ^* amongst the models is largest in the North Atlantic, in the Arctic, along the Atlantic western boundary and in the Southern Ocean (Fig. 2a).

The carbon cycle feedbacks estimated from the cumulative carbon uptake better describe the atmosphere–ocean interaction. The carbon cycle feedbacks estimated from the ocean carbon storage instead better describe the response of the ocean carbon budget to carbon emissions. Here, we focus on the carbon cycle feedbacks estimated from the regional ocean carbon storage to enable diagnostics in terms of the preformed and regenerated carbon pools and to gain more mechanistic insight.

2.2.2 Basin-scale β^* and γ^*

We define the Southern Ocean as south of 35^∘ S and the Arctic Ocean as north of 65^∘ N, and we exclude semi-enclosed seas from our basins' definition (see Supplement Fig. S1 for a map of the regions). The Pacific, Southern and Atlantic oceans contribute equally to the ocean carbon–concentration feedback parameter, as estimated in terms of carbon storage, with an inter-model mean β^* of 0.24, 0.22 and 0.22 PgC ppm⁻¹, respectively (Table 2). The Indian Ocean contributes less than half than the other three basins to β^*, with an inter-model mean of 0.10 PgC ppm⁻¹ (Table 2). The Arctic Ocean has a β^* of only 0.03 PgC ppm⁻¹. The Pacific, Southern, Atlantic, Indian and Arctic oceans have an inter-model mean carbon–climate feedback parameter, defined in terms of carbon storage, γ^*, of −4.8, −3.2, −4.8, −1.2 and −1.2 PgC K⁻¹, respectively (Table 2). The basin-scale variability in β^* amongst the models, as described by the coefficient of variation, CV, is less than 0.15 (Table 2). The basin-scale variability in γ^* amongst the models varies from CV = 0.56 in the Pacific Ocean to CV = 0.85 in the Southern Ocean (Table 2). The inter-model variability in β^* and γ^* for each basin is larger than that of the global ocean (Table 2), which suggests that variability in different basins compensate for each other. For diagnostics of the separate contribution of the ocean carbon uptake and transport on the basin-scale carbon storage and for feedback parameters see Appendix A.

3 Processes controlling the carbon cycle feedbacks in CMIP6 models

To gain insight into the driving mechanisms of the carbon cycle feedbacks and their uncertainty amongst Earth system models, β^* and γ^* may be separated into contribution from the regenerated, the saturated and the disequilibrium ocean carbon pools following the methodology of Williams et al. (2019) and Arora et al. (2020). The ocean dissolved inorganic carbon, DIC, may be defined in terms of these separate carbon pools (Ito and Follows, 2005; Williams and Follows, 2011; Lauderdale et al., 2013; Bernardello et al., 2014):

$$\begin{array}{}\text{(11)}& \mathrm{DIC}={\mathrm{DIC}}_{\mathrm{pref}}+{\mathrm{DIC}}_{\mathrm{reg}}={\mathrm{DIC}}_{\mathrm{sat}}+{\mathrm{DIC}}_{\mathrm{dis}}+{\mathrm{DIC}}_{\mathrm{reg}},\end{array}$$

where DIC_pref is the part of the DIC transferred from the surface into the ocean interior due to the physical ventilation, involving the circulation, and DIC_reg is the part of the DIC accumulated into the ocean interior due to biological regeneration of organic carbon. Similarly, the DIC_pre can be viewed as the part of the DIC associated with the solubility pump and the DIC_reg as the part of the DIC associated with the biological pump. The DIC_pref can be further split into two idealised carbon pools: (i) the DIC_sat representing the amount of DIC that the ocean would have if the whole ocean reached a full chemical equilibrium with the contemporaneous atmospheric CO₂ concentration and (ii) the DIC_dis representing the extent that the ocean departs from a full chemical equilibrium with the contemporaneous atmospheric CO₂. Assuming the changes in the biological organic carbon inventory are small, the changes in the ocean carbon inventory relative to the pre-industrial era, ΔI_ocean in petagrams of carbon, are related to the volume integral of the changes in each of the DIC pools, ΔDIC in moles per cubic metre, as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathrm{\Delta }{I}_{\mathrm{ocean}}& =m\underset{V}{\int }\left(\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{sat}}+\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{dis}}+\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{reg}}\right)\mathrm{d}V\\ & =\mathrm{\Delta }{I}_{\mathrm{sat}}+\mathrm{\Delta }{I}_{\mathrm{dis}}+\mathrm{\Delta }{I}_{\mathrm{reg}},\end{array}\end{array}$$

where $m=\mathrm{12.01}\times {\mathrm{10}}^{-\mathrm{15}}$ PgC mol⁻¹ is a unit conversion from moles to petagrams of carbon.

By substituting Eq. (12) into Eq. (5), for the global ocean, or into Eq. (8), for the regional ocean, the carbon cycle feedback parameters may be diagnosed in terms of these different ocean carbon pools (Williams et al., 2019; Arora et al., 2020):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\mathit{\beta }}^{*}={\mathit{\beta }}_{\mathrm{sat}}+{\mathit{\beta }}_{\mathrm{dis}}+{\mathit{\beta }}_{\mathrm{reg}}={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{sat}}^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{dis}}^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & {\mathit{\gamma }}^{*}={\mathit{\gamma }}_{\mathrm{sat}}+{\mathit{\gamma }}_{\mathrm{dis}}+{\mathit{\gamma }}_{\mathrm{reg}}={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{sat}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{sat}}^{\mathrm{BGC}}}{\mathrm{\Delta }T}}\\ & +{\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{dis}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{dis}}^{\mathrm{BGC}}}{\mathrm{\Delta }T}}+{\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{BGC}}}{\mathrm{\Delta }T}}.\end{array}\end{array}$$

3.1 Contribution from the saturated carbon pool to β^* and γ^*

The saturated part of β^* and γ^* in Eq. (13) is expressed as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\mathit{\beta }}_{\mathrm{sat}}={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{sat}}^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & {\mathit{\gamma }}_{\mathrm{sat}}={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{sat}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{sat}}^{\mathrm{BGC}}}{\mathrm{\Delta }T}},\end{array}\end{array}$$

The changes in the saturated carbon pool relative to the pre-industrial era in Eq. (14) are diagnosed as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathrm{\Delta }{I}_{\mathrm{sat}}& =m\underset{V}{\int }\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{sat}}\mathrm{d}V\\ & =m\underset{V}{\int }\mathrm{\Delta }f\left({\mathrm{CO}}_{\mathrm{2}},T_{\mathrm{ocean}},S,P,\mathrm{Si},{\mathrm{Alk}}_{\mathrm{pre}}\right)\mathrm{d}V,\end{array}\end{array}$$

where Δ is the change relative to the pre-industrial era, T_ocean is the ocean temperature, S is the ocean salinity, P is the ocean phosphate concentration, Si is the ocean silicate concentration, Alk_pre is the preformed alkalinity and f is a non-linear function representing the solution to the ocean carbonate chemistry which provides DIC_sat for the contemporaneous atmospheric CO₂. Here, f is estimated following the iterative solution for the ocean carbonate chemistry of Follows et al. (2006) and by considering the small contribution of minor species (borate, phosphate, silicate) to the preformed alkalinity. In the limit that the ocean hydrogen ion concentration at a chemical equilibrium with the given atmospheric CO₂, [H⁺]_sat, is known or the preformed alkalinity is assumed equal to the carbonate alkalinity, the f function corresponds to the usual solution of the carbonate system that provides DIC_sat based on two knowns: atmospheric CO₂ and either [H⁺]_sat (see Eq. 18) or carbonate alkalinity. The preformed alkalinity is estimated from a multiple linear regression using salinity and the conservative tracer PO (Gruber et al., 1996), with the coefficients of this regression estimated based on the surface (first 10 m) alkalinity, salinity, oxygen and phosphate in each of the Earth system models.

To understand how the ocean carbonate chemistry operates and the mechanisms that control β_sat, consider an ocean buffer factor, B, where the fractional change in the atmospheric CO₂ and saturated carbon inventory is defined relative to the pre-industrial era (Katavouta et al., 2018):

$$\begin{array}{}\text{(16)}& B={\displaystyle \frac{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}{\mathrm{\Delta }{I}_{\mathrm{sat}}/I_{\mathrm{sat},\mathrm{0}}}},\end{array}$$

where subscript 0 denotes the pre-industrial era.

Substituting Eq. (16) into Eq. (14) the saturated part of β^* can be expressed as

$$\begin{array}{}\text{(17)}& {\mathit{\beta }}_{\mathrm{sat}}={\displaystyle \frac{\mathrm{1}}{B({\mathrm{CO}}_{\mathrm{2}},T_{\mathrm{ocean},\mathrm{0}},S_{\mathrm{0}},{\mathrm{Alk}}_{\mathrm{pre},\mathrm{0}})}}{\displaystyle \frac{I_{\mathrm{sat},\mathrm{0}}}{{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}}},\end{array}$$

where $B({\mathrm{CO}}_{\mathrm{2}},T_{\mathrm{ocean},\mathrm{0}},S_{\mathrm{0}},{\mathrm{Alk}}_{\mathrm{0}})$ is the ocean buffer factor for the increasing atmospheric CO₂ but with no climate change (i.e. for the pre-industrial ocean temperature, salinity and alkalinity), as is the case in the BGC run.

Eq. (17) shows that β_sat is proportional to the ocean capacity to buffer changes in atmospheric CO₂ with no changes in the physical climate, $B({\mathrm{CO}}_{\mathrm{2}},T_{\mathrm{ocean},\mathrm{0}},S_{\mathrm{0}},{\mathrm{Alk}}_{\mathrm{pre},\mathrm{0}}{)}^{-\mathrm{1}}$. The rise in atmospheric CO₂ leads not only to an increase in the saturated ocean carbon inventory, ΔI_sat (Fig. 3a, red shade), but also to a decrease in the ocean capacity to buffer changes in atmospheric CO₂ as the ocean acidifies. Accordingly, the buffer factor, B, increases and β_sat decreases with the rise in atmospheric CO₂ at a global and basin scale in all the Earth system models (Fig. 3b, red shade). The buffer factor, B, and thus β_sat also depend on the pre-industrial ocean state due to the non-linearity of the ocean carbonate system. Hence, there is a spread in β_sat amongst the Earth system models forced by the same increase in atmospheric CO₂ (Fig. 3b, red shade) related to their different pre-industrial ocean temperature, salinity and alkalinity.

Figure 3 Ocean carbon storage and ocean carbon–concentration feedback parameter, β^*, for the global ocean and different ocean basins, along with the contribution from the saturated, disequilibrium and regenerated carbon pools in CMIP6 Earth system models: (a) ocean carbon inventory changes relative to the pre-industrial era in the biogeochemically coupled simulation (BGC) and (b) ocean carbon–concentration feedback parameter based on carbon storage, β^* (PgC ppm⁻¹). The solid lines and the shading show the model mean and the model range, respectively, based on the 1 % yr⁻¹ increasing CO₂ experiment over 140 years in 11 CMIP6 models (Table 1). Note that for the global ocean, β^* is equivalent to β.

To understand the mechanisms that control γ_sat, consider the solution to the saturated part of DIC as dictated by the carbonate chemistry,

$$\begin{array}{ll}{\displaystyle }{\mathrm{DIC}}_{\mathrm{sat}}& {\displaystyle }=[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{sat}}+[{\mathrm{HCO}}_{\mathrm{3}}^{-}{]}_{\mathrm{sat}}+[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}{]}_{\mathrm{sat}}\\ \text{(18)}& {\displaystyle }& {\displaystyle }={\mathrm{CO}}_{\mathrm{2}}\left(K_{\mathrm{o}}+{\displaystyle \frac{K_{\mathrm{o}}{K}_{\mathrm{1}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}}+{\displaystyle \frac{K_{\mathrm{o}}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}^{\mathrm{2}}}}\right),\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{sat}}\mathrm{is}{K}_{\mathrm{o}}{\mathrm{CO}}_{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }[{\mathrm{HCO}}_{\mathrm{3}}^{-}{]}_{\mathrm{sat}}\mathrm{is}{\mathrm{CO}}_{\mathrm{2}}{K}_{\mathrm{o}}{\displaystyle \frac{K_{\mathrm{1}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}},\\ {\displaystyle }& {\displaystyle }[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}{]}_{\mathrm{sat}}\mathrm{is}{\mathrm{CO}}_{\mathrm{2}}{K}_{\mathrm{o}}{\displaystyle \frac{K_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}^{\mathrm{2}}}},\end{array}$$

K_o is the solubility, K₁ and K₂ are the ocean carbon dissociation constants, and [H⁺]_sat is the ocean hydrogen ion concentration at a chemical equilibrium with the contemporaneous atmospheric CO₂. K_o, K₁ and K₂ are a function of the ocean temperature and salinity and so depend on the physical climate change, while [H⁺]_sat depends primarily on the changes in atmospheric CO₂. Combining Eq. (18) with Eqs. (14) and (15) and assuming that the ocean temperature and salinity remain at their pre-industrial value in the BGC run with no climate change, γ_sat can be expressed as

$$\begin{array}{}\text{(19)}& {\mathit{\gamma }}_{\mathrm{sat}}={\displaystyle \frac{m}{\mathrm{\Delta }T}}\underset{V}{\int }{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{\Delta }{K}_{\mathrm{o}}+{\displaystyle \frac{\mathrm{\Delta }\left(K_{\mathrm{o}}{K}_{\mathrm{1}}\right)}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}}+{\displaystyle \frac{\mathrm{\Delta }\left(K_{\mathrm{o}}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\right)}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}^{\mathrm{2}}}}\right)\mathrm{d}V.\end{array}$$

By expanding $[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}=[{\mathrm{H}}^{+}{]}_{\mathrm{sat},\mathrm{0}}+\mathrm{\Delta }[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}$, Eq. (19) can be written as

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\mathit{\gamma }}_{\mathrm{sat}}& ={\displaystyle \frac{m}{\mathrm{\Delta }T}}{\mathrm{CO}}_{\mathrm{2}}\underset{V}{\int }\left(\underset{\text{effect of ocean warming under}{\mathrm{pH}}_{\mathrm{0}}}{\underbrace{\left\{\mathrm{\Delta }{K}_{\mathrm{o}}+{\displaystyle \frac{\mathrm{\Delta }\left(K_{\mathrm{o}}{K}_{\mathrm{1}}\right)}{[{\mathrm{H}}^{+}{]}_{\mathrm{0},\mathrm{sat}}}}+{\displaystyle \frac{\mathrm{\Delta }\left(K_{\mathrm{o}}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\right)}{[{\mathrm{H}}^{+}{]}_{\mathrm{0},\mathrm{sat}}^{\mathrm{2}}}}\right\}}}\right.\\ & \left.-\underset{\text{effect of ocean warming under}\mathrm{\Delta }\mathrm{pH}}{\underbrace{\left\{{\displaystyle \frac{\mathrm{\Delta }\left(K_{\mathrm{o}}{K}_{\mathrm{1}}\right)}{[{\mathrm{H}}^{+}{]}_{\mathrm{0},\mathrm{sat}}}}{\displaystyle \frac{\mathrm{\Delta }[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}}+{\displaystyle \frac{\mathrm{\Delta }\left(K_{\mathrm{o}}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\right)}{[H^{+}{]}_{\mathrm{0},\mathrm{sat}}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{\Delta }\left([{\mathrm{H}}^{+}{]}_{\mathrm{sat}}^{\mathrm{2}}\right)}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}^{\mathrm{2}}}}\right\}}}\right)\mathrm{d}V.\end{array}\end{array}$$

The first term in curly brackets in Eq. (20) contributes to the linear part of γ_sat and is controlled by changes in climate only, specifically by the effect of changes in the ocean temperature to the solubility and ocean carbon dissociation constants. Under warming due to changes in climate, this term is negative. The second term in curly brackets in Eq. (20) contributes to the non-linear part of γ_sat, and it depends on changes in pH, due to the increase in atmospheric CO₂, and changes in climate. Under rising atmospheric CO₂ and warming, this term is negative as [H⁺]_sat increases and K_o decreases, and it is smaller than the linear term as $\frac{\mathrm{\Delta }[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{sat}}}<\mathrm{1}$. The adjustment of γ_sat due to changes in the pH, represented by this non-linear term, depends on the increase in atmospheric CO₂; for example, the non-linear term is about 30 % and 50 % of the linear term for a doubling of atmospheric CO₂ and for a quadrupling of atmospheric CO₂, respectively. Hence, the non-linearity of the carbonate chemistry acts to reduce the magnitude of the negative γ_sat consistent with the carbonate system being less sensitive to change in temperature under higher ocean DIC (Schwinger et al., 2014).

Eq. (20) shows that γ_sat is proportional to the changes in solubility due to climate change and is further modified by changes in the ocean carbon dissociation constants with warming and by the non-linearity of the carbonate chemistry. Ocean warming due to climate change leads to a decrease in the saturated carbon inventory, ΔI_sat, in all basins (Fig. 4a, red shade) primarily driven by a decrease in solubility. This decrease in ΔI_sat with warming drives a nearly constant negative γ_sat (Fig. 4b, red shade) with the deviations from a constant value being associated with the non-linearity of the carbonate system. The spread of γ_sat amongst the Earth system models (Fig. 4b, red shade) is relatively small at a global and basin scale, except in the Arctic, and is associated with different pre-industrial ocean states in the models.

Figure 4 Ocean carbon storage and ocean carbon–climate feedback parameter, γ^*, for the global ocean and different ocean basins, along with the contribution from the saturated, disequilibrium and regenerated carbon pools in CMIP6 Earth system models: (a) ocean carbon inventory changes in the fully coupled simulation (COU) minus the biogeochemically coupled simulation (BGC) and (b) ocean carbon–climate feedback parameter based on carbon storage, γ^* (PgC K⁻¹). The solid lines and the shading show the model mean and the model range, respectively, based on the 1 % yr⁻¹ increasing CO₂ experiment over 140 years in 11 CMIP6 models (Table 1). Note that for the global ocean, γ^* is equivalent to γ.

3.2 Contribution from the regenerated carbon pool to β^* and γ^*

The regenerated part of β^* and γ^* in Eq. (13) is expressed as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& {\mathit{\beta }}_{\mathrm{reg}}={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{BGC}}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}}},\\ & {\mathit{\gamma }}_{\mathrm{reg}}={\displaystyle \frac{\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{BGC}}}{\mathrm{\Delta }T}}.\end{array}\end{array}$$

Assuming that the oxygen concentration is close to saturation at the surface, DIC_reg can be estimated from the apparent oxygen utilisation, AOU, and the contribution of biological calcification to alkalinity, Alk (Ito and Follows, 2005; Williams and Follows, 2011; Lauderdale et al., 2013), such that ΔI_reg in Eq. (21) is diagnosed as

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\mathrm{\Delta }{I}_{\mathrm{reg}}& =m\underset{V}{\int }\mathrm{\Delta }{\mathrm{DIC}}_{\mathrm{reg}}\mathrm{d}V=m\underset{V}{\int }(R_{\mathrm{CO}}\mathrm{\Delta }\mathrm{AOU}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{\Delta }\mathrm{Alk}-\mathrm{\Delta }{\mathrm{Alk}}_{\mathrm{pre}}-R_{\mathrm{NO}}\mathrm{\Delta }\mathrm{AOU}\right))\mathrm{d}V,\end{array}\end{array}$$

where R_CO and R_NO are constant stoichiometric ratios and Alk_pre is the preformed alkalinity such that Alk−Alk_pre gives the contribution to alkalinity from biological calcification.

The regenerated part of β^* is associated with changes in ocean biological processes due to the atmospheric CO₂ increase. $\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{BGC}}$ and β_reg are effectively negligible in the Earth system models (Fig. 3a and b, green shade) as these models do not include an explicit dependence of biological production on an increase in carbon availability or decrease in seawater pH. The regenerated part of γ^* is associated with changes in ocean biological processes due to changes in climate, including the effect of changes in the circulation on the sinking rate of particles, the effect of warming on the solubility of oxygen and the effect of changes in alkalinity on the dissolution of the calcium carbonate shells of calcifying phytoplankton. $\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{COU}}-\mathrm{\Delta }{I}_{\mathrm{reg}}^{\mathrm{BGC}}$ and γ_reg are positive and increase in time on a global and basin scale (Fig. 4a and b, green shade), indicating that γ_reg is dominated by the weakening in the ocean physical ventilation due to climate change. This weakening in the ocean physical ventilation leads to a longer residence time of water masses in the ocean interior and so to an increase in the accumulation of carbon from the regeneration of biologically cycled carbon in the deep ocean (Schwinger et al., 2014; Bernardello et al., 2014).

3.3 Contribution from the disequilibrium carbon pool to β^* and γ^*

The disequilibrium parts of β^* and γ^* are diagnosed using Eq. (13) as

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& {\mathit{\beta }}_{\mathrm{dis}}={\mathit{\beta }}^{*}-{\mathit{\beta }}_{\mathrm{sat}}-{\mathit{\beta }}_{\mathrm{reg}},\\ & {\mathit{\gamma }}_{\mathrm{dis}}={\mathit{\gamma }}^{*}-{\mathit{\gamma }}_{\mathrm{sat}}-{\mathit{\gamma }}_{\mathrm{reg}}.\end{array}\end{array}$$

The disequilibrium part of the carbon cycle feedback parameters is controlled by the ocean physical ventilation. Specifically, β_dis is a function of the pre-industrial ocean physical ventilation and the rate of transfer of the anthropogenic carbon from the ocean surface into the ocean interior. The rise in atmospheric CO₂ leads to an increase in the magnitude of the negative disequilibrium ocean carbon inventory, ΔI_dis, at a global and basin scale (Fig. 3a, blue shade), as the ocean physical ventilation is relatively slow and the ocean carbon transfer over the ocean interior cannot keep up with the rate of increase in atmospheric CO₂. Hence, β_dis is negative at a global and basin scale in all the Earth system models (Fig. 3b, blue shade). However, the rate of increase in the magnitude of ΔI_dis slows down in time (Fig. 3a, blue shade) and β_dis becomes less negative in time (Fig. 3b, blue shade), as more anthropogenic carbon is transferred into the ocean interior while the buffer capacity of the ocean decreases, which brings the ocean closer to an equilibrium with the contemporaneous atmospheric CO₂.

The disequilibrium part of γ^* depends on the weakening of the ocean physical ventilation with climate change. Here, γ_dis is defined based on the climate change impact under rising atmospheric CO₂ (i.e. COU–BGC runs) and so includes (i) the effect of weakening ventilation on the pre-industrial ocean carbon gradient involving the natural carbon, (ii) the effect of weakening ventilation on the anthropogenic carbon, and (iii) the effect from decreasing sea-ice coverage leading to an increase in the ocean in direct contact with the atmosphere. Overall, the effect of the weakening ventilation on the combined anthropogenic and natural carbon leads to a negative ΔI_dis and a negative γ_dis on a global scale and in the Atlantic, Indian, Pacific and Southern oceans after year 40 (Fig. 4a and b, blue shade). In the Arctic, the effect of the decreasing sea-ice coverage drives a slightly positive γ_dis over the first 80 years, while the effect of the weakening ventilation dominates and drives a negative γ_dis after year 80. The disequilibrium part of γ^* becomes more negative in time on a global and basin scale as the ocean ventilation weakens with warming.

3.4 Combined effect of saturated, disequilibrium and regenerated carbon pools on β^* and γ^*

On a global and basin scale, the ocean carbon–concentration feedback parameter, β^*, is positive in all the Earth system models. This positive β^* is explained by the chemical response involving the rise in ocean saturation, β_sat, opposed by the effect of the relatively slow ocean ventilation, such that the physical uptake of carbon within the ocean is unable to keep pace with the rise in atmospheric CO₂, β_dis (Fig. 3b). There is no significant contribution from biological changes to β^* in the Earth system models. The spread in β_sat amongst the Earth system models is small on a global and basin scale (Fig. 5a), reflecting the use of similar carbonate chemistry schemes and bulk parameterisations of air–sea CO₂ fluxes across marine biogeochemical models in CMIP6 (Séférian et al., 2020). The spread in β_dis is also small on a global and basin scale (Fig. 5a) as all the Earth system models have a broadly similar general circulation and physical ventilation in the pre-industrial era.

Figure 5 Carbon cycle feedback parameters based on carbon storage, along with the contribution from the saturated, disequilibrium and regenerated carbon pools, for the global ocean and the different ocean basins in 11 CMIP6 Earth system models (Table 1): (a) carbon–concentration feedback parameter, β^* (PgC ppm⁻¹), and (b) carbon–climate feedback parameter, γ^* (PgC K⁻¹). The estimates are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂). For the inter-model mean and standard deviation of these estimates, see Supplement Table S1. Note that for the global ocean, β^* and γ^* are equivalent to β and γ.

The decrease in solubility and in the physical ventilation with warming reduces the ocean carbon uptake, leading to negative γ_sat and γ_dis, respectively (Fig. 4b). However, the decrease in the ventilation with warming also acts to increase the residence time in the ocean interior, leading to an increase in the regenerated carbon and a positive γ_reg (Fig. 4b, green shading). The combined γ_sat and γ_dis dominate over the opposing γ_reg, leading to an overall negative γ^* on a global and basin scale. On a global scale, γ_sat, γ_dis and γ_reg are of a similar magnitude (Fig. 5b, black circles). The inter-model spread in the global γ^* is mainly driven by the spread in γ_dis and γ_reg (Fig. 5b, black circles) and is associated with the different response of the ocean ventilation to warming in these models. The inter-model spread in the global γ_reg is larger than the spread in the global γ_dis due to the different parameterisations of ocean biogeochemical processes in the models.

In the Southern Ocean, the contributions from the saturated, disequilibrium and regenerated carbon pools to γ^* are of a similar magnitude (Fig. 5b, blue circles), such that the decrease in solubility, the reduction in the physical ventilation and the increase in the regenerated carbon accumulation in the ocean interior due to climate change are equally important. The inter-model spread in γ^* in the Southern Ocean is dominated by the spread in γ_dis and γ_reg. In the Pacific and Indian oceans, the magnitude of γ^* is primarily controlled by the saturated carbon pool and the decrease in carbon solubility due to warming (Fig. 5b, purple and yellow circles). However, the inter-model spread in γ^* in the Pacific and Indian oceans is dominated by the response of the regenerated carbon pool to climate change, γ_reg (Fig. 5b, purple and yellow circles). In the Atlantic Ocean, γ^* is dominated by the disequilibrium carbon pool (Fig. 5b, red circles) and the reduction in the physical ventilation due to climate change, with the contribution from the saturated carbon pool being relatively small. In the Arctic Ocean, γ^* is primarily controlled by the saturated carbon pool and the decrease in solubility, with the contribution from the regenerated carbon pool being negligible (Fig. 5b, green circles).

On regional scales, the contributions from the saturated, disequilibrium and regenerated carbon pools to β^* and γ^* are further modified by local upwelling, changes in alkalinity and the conversion of regenerated carbon to disequilibrium carbon at the ocean surface, as discussed in Appendix B.

3.5 Processes controlling the contribution from different basins to β^* and γ^*

The Southern and Indian oceans contribute 27 % and 12 %, respectively, to the ocean carbon–concentration feedback parameter, β^*, following their fractional volumes of the global ocean (Fig. 6a). However, the Atlantic and Arctic oceans contribute 26 % and 3 % to β^*, respectively, which is significantly more than their fractional volume of 18 % and 1 % of the global ocean. In contrast, the Pacific Ocean contributes only 30 % to β^* despite its fractional volume of 42 % of the global ocean. By definition, the contribution of each basin to β_sat is approximately proportional to the ocean volume contained in each basin (Fig. 6a). However, β_dis is relatively low in the Atlantic and Arctic oceans and high in the Pacific Ocean compared with their respective volumes (Fig. 6a), which may be understood by the Atlantic and Arctic oceans interior being more ventilated and the Pacific Ocean interior being less ventilated than the rest of the ocean. Specifically, the low β_dis and the relatively large contribution from the Atlantic Ocean to β^* are due to a large transfer of anthropogenic carbon into the ocean interior from strong local physical ventilation and transport of carbon from the Southern Ocean. In the Arctic Ocean, the transfer of anthropogenic carbon from the well-ventilated Atlantic Ocean contributes towards a decrease in β_dis and to a relatively large β^*. The Southern Ocean has large anthropogenic carbon uptake from the atmosphere, but its contribution to β^*, as estimated from carbon storage, is relatively small due to large carbon transport to the other basins (Fig. A1).

Figure 6 The fractional contribution (in %) of different ocean basins to the total volume of the ocean and to the ocean carbon cycle feedback parameters based on carbon storage, along with the contribution from the saturated, disequilibrium and regenerated carbon pools based on the inter-model mean of 11 CMIP6 models (Table 1): (a) carbon–concentration feedback parameter, β^*, and (b) carbon–climate feedback parameter, γ^*. The estimates are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂). The regenerated part of β^* is omitted as its contribution is negligible in all basins (Fig. 5a). The Arctic Ocean has a fractional volume of ∼ 1 % of the global ocean and a fractional contribution of 1.2 %, 0.7 %, 2 % and 1 % to β_sat, β_dis, γ_dis and γ_reg, respectively (not explicitly shown in the figure). The combined fractional contribution of the five ocean basins to β^* and γ^* is less than 100 % due to a small contribution, <2 %, from semi-enclosed seas with a fractional volume of less than 0.5 % of the global ocean (Supplement Fig. S1).

The Pacific and Indian oceans' contributions to γ_sat are slightly smaller than expected from their fractional volumes (Fig. 6b), consistent with a low warming per unit volume in these basins (Fig. 1d). The Pacific and Indian oceans' contributions to γ_dis and γ_reg are significantly smaller than expected from their fractional volumes (Fig. 6b), indicating that there is no significant effect from changes in ventilation in these basins. This absence of any significant effect from changes in the ventilation with warming in the Pacific and the Indian oceans leads to their much smaller contribution to γ^* relative to their volumes (Fig. 6b).

The Atlantic Ocean has a contribution of 31 % to γ^*, which is much larger than expected from its fractional volume of 18 % of the global ocean. This large contribution from the Atlantic Ocean is primarily due to γ_dis (Fig. 6b) and the reduction in the physical ventilation due to climate change. The Atlantic Ocean has a smaller contribution to γ_sat than expected from its fractional volume (Fig. 6b), despite experiencing large warming (Fig. 1d), which suggests that the non-linearity of the carbonate system is important in this basin. Specifically, the Atlantic Ocean has a large increase in DIC (Fig. 1c) which acts to significantly reduce the magnitude of the negative γ_sat driven solely by the effect of warming on solubility (see Eq. 20). The Arctic Ocean has a contribution of 8 % to γ^*, which is much larger than expected from its fractional volume of 1 % of the global ocean. This large contribution from the Arctic Ocean is primarily associated with the saturated part of γ^* (Fig. 6b) and a very large warming per unit volume in this basin (Fig. 1d).

The Southern Ocean has a contribution of 38 % to γ_dis and of 53 % to γ_reg, which is much larger than expected from its fractional volume (Fig. 6b). These large contributions indicate that changes in the physical ventilation and the accumulation of regenerated carbon due to climate change have a large effect on the Southern Ocean carbon storage. The Southern Ocean contribution to γ_sat is also larger than expected from its fractional volume (Fig. 6b), consistent with a large warming per unit volume in this basin (Fig. 1d). Hence, the apparent small contribution of the Southern Ocean to γ^*, of 21 %, is due to compensation between (i) the large decrease in carbon storage associated with the combined decrease in solubility and physical ventilation and (ii) the large increase in carbon storage associated with a longer residence time and accumulation of regenerated carbon in the Southern Ocean interior.

4 Dependence of the carbon cycle feedbacks on the Atlantic Meridional Overturning Circulation

The ocean carbon cycle feedbacks are controlled by ocean ventilation and the transfer of carbon from the mixed layer to the thermocline and deep ocean. The ocean ventilation involves the seasonal cycle of the mixed layer; the subduction process; and the effects of the eddy, gyre and overturning circulations. Despite this complexity, the strength of the Atlantic Meridional Overturning Circulation (AMOC) and its weakening due to climate change is often used as a proxy for the large-scale ocean ventilation. Here we investigate the dependence of the carbon cycle feedbacks on the AMOC strength and its weakening with climate change.

4.1 Insight from an idealised climate model with a meridional overturning

The idealised climate model of Katavouta et al. (2019) is used to investigate the control of the carbon cycle feedbacks by the AMOC. This idealised model consists of a slab atmosphere, two upper ocean boxes for the southern and northern high latitudes, two boxes for the mixed layer and the thermocline in the low and middle latitudes, and one box for the deep ocean (Fig. 7a). The model solves for the thermocline thickness from a volumetric balance between the surface cooling conversion of light to dense waters in the North Atlantic, the diapycnal transfer of dense to light waters in low and middle latitudes, and the conversion of dense to light waters in the Southern Ocean involving Ekman transport partially compensated for by poleward mesoscale eddy transport (Gnanadesikan, 1999; Johnson et al., 2007; Marshall and Zanna, 2014). The model also accounts for the rate of subduction occurring in the Southern Ocean versus the tropics and subtropics through an isolation fraction for water remaining below the mixed layer and spreading northwards in the thermocline. The model solves for the ocean carbon cycle including physical and chemical transfers but ignores biological transfers and sediment and weathering interactions involving changes in the cycling of organic carbon or calcium carbonate. The ocean carbonate system is solved using the iterative algorithm of Follows et al. (2006). For the model closures and an explicit description of the model budgets and equations, see Katavouta et al. (2019).

Figure 7 (a) A simplified climate model with overturning circulation, including a slab atmosphere, an upper layer of light water consisting of a thermocline layer and a surface mixed layer in the low and middle latitudes, two upper layers at southern and northern high latitudes, and a lower layer of dense water. (b) Meridional overturning in the three experiments forced by a 1 % yr⁻¹ increase in atmospheric CO₂: (i) control experiment (solid line), (ii) weaker pre-industrial overturning (dotted line) and (iii) larger reduction in the overturning with climate change (dashed line). (c) The ocean carbon–concentration feedback parameter, β, and (d) the ocean carbon–climate feedback parameter, γ, along with their saturated and disequilibrium parts, in the experiments with different meridional overturning. The estimates of β and γ are based on the fully coupled simulations (COU) and the biogeochemically coupled simulations (BGC). Note that for the “global” ocean of the simplified model, β and γ are equivalent to β^* and γ^*.

The model is first integrated to a pre-industrial steady state, with the distribution of temperature and DIC depending on the pre-industrial strength of the overturning. The model is then forced by a 1 % yr⁻¹ increase in atmospheric CO₂ concentration from a pre-industrial value of 280 ppm until atmospheric CO₂ quadruples over a 140-year period. This increase in atmospheric CO₂ drives a radiative forcing: $R=a\ln \left({\mathrm{CO}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2},\mathrm{0}}\right)$, where a is 5.35 W m⁻² (Myhre et al., 1998) and subscript ₀ denotes the pre-industrial state. This radiative forcing then drives a radiative response, λΔT_air, and net planetary heat uptake given by the downward heat flux entering the system at the top of the atmosphere, N_TOA (Gregory et al., 2004) such that $R=\mathit{\lambda }\mathrm{\Delta }{T}_{\mathrm{air}}+N_{\mathrm{TOA}}$, where T_air is the temperature of the slab atmosphere and λ is the climate feedback parameter that is assumed constant and equal to 1 W m⁻² K⁻¹ for simplicity. The ocean heat uptake, N, is estimated as the planetary heat uptake minus the atmospheric heat uptake, $N=N_{\mathrm{TOA}}-c(\mathrm{\Delta }{T}_{\mathrm{surf}}-\mathrm{\Delta }{T}_{\mathrm{air}})$, where c=50 W m⁻² K⁻¹ is an air–sea heat transfer parameter and T_surf is the ocean temperature at the surface. The ocean heat uptake is distributed equally over the ocean surface. For this model closure, the ocean heat uptake is more than 95 % of the net planetary heat uptake.

This additional ocean heat uptake reduces the conversion of light to dense waters in the North Atlantic, q_NA, and leads to an overturning weakening, following

$$\begin{array}{}\text{(24)}& \mathrm{\Delta }{q}_{\mathrm{NA}}=-{\displaystyle \frac{AN}{\mathit{\rho }{C}_p{T}_{\mathrm{contrast}}}},\end{array}$$

where A is the model area covered by the low and middle latitudes; ρ is a referenced ocean density; C_p is the specific heat capacity for the ocean; and T_contrast is the temperature contrast between light waters in the low and middle latitudes, involving the mixed layer and thermocline, and the dense waters in the deep ocean.

Three experiments were conducted using this idealised model forced by a 1 % yr⁻¹ increase in atmospheric CO₂: (i) a control experiment with a strong pre-industrial AMOC that experiences a weakening with warming as described by Eq. (24) (solid line in Fig. 7b); (ii) an experiment with a weak pre-industrial AMOC that experiences the same weakening with warming as in the control experiment (dotted line in Fig. 7b); and (iii) an experiment with a strong pre-industrial AMOC that experiences a doubled weakening with warming compared to that in the control experiment, such that the right hand of Eq. (24) is multiplied by a factor of 2 (dashed line in Fig. 7b). In all the experiments, a fully coupled simulation (COU) and a biogeochemically coupled simulation (BGC) are used to estimate the carbon cycle feedback parameters β and γ.

In the idealised model, a weaker pre-industrial meridional overturning leads to a smaller carbon–concentration feedback parameter, β, during the centennial transient response to the increase in atmospheric CO₂ (Fig. 7c, black lines). The saturated part of β does not directly depend on the AMOC strength (Fig. 7c, red lines). In contrast, a weaker pre-industrial AMOC leads to a more negative disequilibrium part of β (Fig. 7c, blue lines) due to a weaker and slower transfer of anthropogenic carbon below the ocean surface. Hence, the pre-industrial AMOC controls β through its effect on the physical ventilation via the disequilibrium carbon pool.

The strength of the pre-industrial meridional overturning has only a small impact on the carbon–climate feedback parameter, γ, that is associated with the saturated carbon pool (Fig. 7d, solid and dotted lines). In the idealised model, a different pre-industrial overturning is associated with a different pre-industrial ocean state, in terms of temperature and DIC, which has a small effect on the saturated part of γ due to the non-linearity of the carbonate chemistry. The carbon–climate feedback parameter, γ, is primarily controlled by the changes in the AMOC with warming and specifically by the dependence of the negative γ_dis on the AMOC weakening with warming (Fig. 7d, solid and dashed lines). Hence, a more pronounced weakening in the AMOC with warming drives a more negative γ, as there is a more pronounced decrease in the physical transfer of carbon into the ocean interior.

Figure 8 The Atlantic Meridional Overturning Circulation (AMOC; in Sv), in the pre-industrial era: (a–k) 11 CMIP6 Earth system models and (l) their inter-model mean. The estimates are based on pre-industrial control simulation (piControl) for years 121 to 140.

Figure 9 Weakening of the Atlantic Meridional Overturning Circulation (AMOC; Sv), relative to the pre-industrial era: (a–k) 11 CMIP6 Earth system models and (l) their inter-model mean. The estimates are based on the fully coupled simulation (COU) under the 1 % yr⁻¹ increasing CO₂ experiment minus the pre-industrial control simulation (piControl). Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂).

4.2 Results from CMIP6 Earth system models

The strength and weakening of the AMOC are diagnosed as the maximum pre-industrial AMOC between 30 and 50^∘ N (Fig. 8) and the maximum AMOC change between 30 and 50^∘ N (Fig. 9), respectively, in 11 CMIP6 Earth system models (Table 3). There is a large spread in the pre-industrial strength of the AMOC amongst the Earth system models (Fig. 8), with a range of 11.3 sverdrups (Sv) in IPSL-CM6A-LR to 23.2 Sv in NorESM2-LM (Table 3). There is also a large spread in the response of the AMOC to climate change amongst the models (Fig. 9), with the magnitude of the AMOC weakening ranging from −5.6 Sv in IPSL-CM6A-LR to −18.7 Sv in MRI-ESM2 (Table 3). Generally, the models with a stronger pre-industrial AMOC simulate a larger weakening in the AMOC with warming. Linear correlations are used to relate the AMOC variability amongst the models with the variability in the carbon cycle feedbacks for the global ocean and the different ocean basins (Table 4).

Table 3 Carbon–concentration feedback parameter based on carbon storage, β^* (PgC ppm⁻¹), and carbon–climate feedback parameter based on carbon storage, γ^* (PgC K⁻¹), in the Atlantic Ocean, and the contribution from the saturated, disequilibrium and regenerated carbon pools, along with the strength of the Atlantic Meridional Overturning Circulation (AMOC) in sverdrups, in the pre-industrial era and its weakening relative to the pre-industrial era due to climate change. The strength and weakening of the AMOC are estimated as the maximum pre-industrial AMOC and maximum AMOC change between 30 and 50^∘ N, respectively (Figs. 8 and 9). Results are shown for 11 CMIP6 Earth system models, along with the inter-model mean, standard deviation and coefficient of variation (CV), estimated as the standard deviation divided by the mean. The estimates are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂).

Table 4 Correlation between the ocean carbon–concentration feedback parameter based on carbon storage, β^*, and the pre-industrial strength of the Atlantic Meridional Overturning Circulation (AMOC) and correlation between the ocean carbon–climate feedback parameter based on carbon storage, γ^*, and the AMOC weakening due to climate change, along with the contribution from the saturated, disequilibrium and regenerated carbon pools based on 11 CMIP6 models (Table 1). For the regenerated part of γ^*, correlations estimated excluding the CNRM-ESM2-1 model, which has substantially larger changes in the regenerated carbon with warming than the rest of the Earth system models, are presented as ${\widehat{\mathit{\gamma }}}_{\mathrm{reg}}$. The correlation is expressed as a correlation coefficient, r, and a p value. The level of significance is assumed to be 0.05, and the statistically significant correlations with p<0.05 are highlighted by bold text. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂). Note that for the global ocean, β^* and γ^* are equivalent to β and γ.

For the global ocean, there is a moderate positive correlation between the pre-industrial AMOC and the ocean carbon–concentration feedback parameter, β (r=0.72; Table 4), such that models with a stronger pre-industrial AMOC have a more positive β. However, there is no significant correlation between the pre-industrial AMOC and the saturated or disequilibrium part of β (Table 4) on the global scale. In the Atlantic Ocean, β^* is strongly correlated with the pre-industrial AMOC in the Earth system models (r= 0.87), and this correlation is due to the contribution from β_dis (Table 4 and Fig. 10a), which suggests that a stronger pre-industrial AMOC leads to a larger β^* via a stronger physical transfer of anthropogenic carbon below the ocean surface there, which is similar to the behaviour of the idealised model (Fig. 7c).

In the Southern Ocean, there is a statistically significant correlation between the pre-industrial AMOC and β^* (r=0.82) that is not associated with the saturated, the disequilibrium or the regenerated carbon pools (Table 4). A stronger pre-industrial AMOC is associated with a stronger return flow of interior water into the surface in the Southern Ocean upwelling branch of the overturning circulation (Marshall and Speer, 2012) and hence with stronger local anthropogenic carbon uptake from the atmosphere in the BGC run there. However, in the Southern Ocean both the anthropogenic carbon uptake from the atmosphere and its transport to the other basins are also controlled by subduction (Sabine et al., 2004; Sallée et al., 2012, 2013) and the residual circulation involving the wind-driven Ekman transport and eddy fluxes (Lauderdale et al., 2013). Hence, there is no significant correlation between the pre-industrial AMOC and β_dis in the Southern Ocean as other processes involving subduction of mode and intermediate waters, as well as the residual circulation are important in setting the pre-industrial physical ventilation and transport of carbon to other basins. There is also no significant correlation between the pre-industrial AMOC and β^* in the Pacific, Indian and Arctic oceans, which suggests that other processes dominate the ocean ventilation there.

Figure 10 (a) Dependence of the carbon–concentration feedback parameter based on carbon storage, β^* (PgC ppm⁻¹), on the pre-industrial strength of the Atlantic Meridional Overturning Circulation (AMOC; Sv), in the Atlantic Ocean, and (b) dependence of the carbon–climate feedback parameter based on carbon storage, γ^* (PgC K⁻¹), on the AMOC weakening with climate change in the Atlantic Ocean, in 11 CMIP6 models (Table 1), along with the contribution from the saturated, disequilibrium and regenerated carbon pools. The black lines correspond to the regression line based on an ordinary least-squares regression, with the corresponding correlation coefficient, r, and p value shown in each panel. The dashed red line and r and p value in red for γ_reg correspond to the regression line and the correlation coefficient estimated by excluding CNRM-ESM2-1, the model with the largest deviation from the rest of the Earth system models in terms of the regenerated carbon pool. The strength and weakening of the AMOC are estimated as the maximum pre-industrial AMOC and maximum AMOC change between 30 and 50^∘ N, respectively (Table 3). Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂).

For the global ocean, the carbon–climate feedback parameter, γ, and its disequilibrium part, γ_dis, are positively correlated with the AMOC weakening due to climate change in the Earth system models (r=0.62 and r=0.71, respectively, Table 4). This correlation for the global ocean is to first order driven by the Atlantic Ocean, as there is no significant correlation between AMOC weakening and γ^* in the Pacific, Indian and Southern oceans or between AMOC weakening and γ_dis in the Arctic Ocean (Table 4). In the Atlantic Ocean, a more pronounced AMOC weakening due to climate change is associated with a larger reduction in the physical transfer of carbon below the ocean surface and so with a more negative γ^* and γ_dis (Fig. 10b), consistent with the inferences for the idealised model (Fig. 7d). In the Arctic Ocean, there is a moderate correlation between γ^* and AMOC weakening (r=0.67; Table 4). This correlation between γ^* and AMOC weakening in the Arctic Ocean is associated with the carbonate chemistry and the saturated carbon pool, γ_sat (Table 4) and it appears to be related to changes in alkalinity (Supplement Fig. S3).

In the Atlantic Ocean, the correlation between γ_reg and AMOC weakening is negative ($r=-\mathrm{0.64}$; Table 4), when excluding the CNRM-ESM2-1 model which has a much larger γ_reg than the rest of the Earth system models (Table 3). A more pronounced AMOC weakening leads to a longer residence time in the Atlantic Ocean interior and so to a larger accumulation of regenerated carbon there and a more positive γ_reg (Fig. 10b). The substantially larger change in the regenerated carbon pool with warming in CNRM-ESM2-1 compared to the rest of the models is probably due to a revised parameterisation for organic matter remineralisation (in PISCESv2-gas), its parameterisation of sedimentation and its interactive riverine input; see Séférian et al. (2020) for a discussion of differences in ocean biogeochemistry amongst the Earth system models.

5 Discussion and summary

Our study reveals the contribution of the Atlantic, Pacific, Southern, Indian and Arctic oceans to the ocean carbon cycle feedbacks in a set of CMIP6 Earth system models. This basin-scale contribution is explained in terms of the effect of the carbonate chemistry, physical ventilation and biological processes. Experiments using an idealised climate model and diagnostics of the CMIP6 models suggest a dependence of the ocean carbon cycle feedbacks on the strength of the pre-industrial Atlantic Meridional Overturning Circulation (AMOC) and on the magnitude of the AMOC weakening with warming.

The global ocean carbon–concentration feedback parameter, β, in CMIP6 models is controlled by competition between opposing contributions from the saturated and the disequilibrium carbon pools, such that (i) the carbonate chemistry drives a positive β that decreases with an increase in atmospheric CO₂ as the ocean acidifies and its capacity to buffer atmospheric CO₂ decreases and (ii) the physical ventilation drives a negative β that becomes less negative in time as the anthropogenic carbon is transferred from the ocean surface into the ocean interior. The global ocean carbon–climate feedback parameter, γ, in CMIP6 models is controlled by competition between the combined decrease in the saturated and disequilibrium carbon pools and the increase in the regenerated carbon pool, such that (i) the decrease in solubility and weakening in physical ventilation with warming drive a negative γ and (ii) the increase in the accumulation of regenerated carbon, associated with a longer residence time in the ocean interior with climate change, drives a positive γ.

5.1  Regional ocean carbon cycle feedbacks

The regional ocean carbon storage is controlled by the local air–sea carbon exchange and the transport of carbon by the ocean circulation. Here, the regional carbon cycle feedbacks are estimated from the regional ocean carbon storage so that the effect of the ocean transport of carbon is included in the feedbacks. This transport effect acts to decrease the carbon–concentration feedback parameter, β^*, in the Southern Ocean and increase β^* in the Southern Hemisphere subtropical gyres and the western boundary of the Atlantic Ocean in CMIP6 models. The transport effect on β^* is consistent with the view that while almost half of the anthropogenic carbon enters through the Southern Ocean, much of this carbon is transferred and stored in the Atlantic, Pacific and Indian oceans (Sabine et al., 2004; Khatiwala et al., 2009; Frölicher et al., 2015). The transport effect also leads to a more spatially uniform and negative carbon–climate feedback parameter estimated in terms of carbon storage, γ^*, compared with the carbon–climate feedback parameter estimated in terms of cumulative carbon uptake.

The Atlantic, Pacific and Southern oceans provide comparable contributions, of between 26 % and 30 %, to the global ocean carbon–concentration feedback, as estimated from the carbon storage, despite their different size. This large contribution from the Atlantic Ocean relative to its volume is associated with an enhanced physical transfer of anthropogenic carbon into the ocean interior due to strong local ventilation and transport of carbon from the Southern Ocean. The Arctic Ocean also provides a larger contribution to β^* relative to its volume, consistent with observation-based estimates of anthropogenic carbon storage (Tanhua et al., 2009). The inter-model variability in the carbon–concentration feedback parameter, β^*, as estimated by the coefficient of variation, is relatively small in all the ocean basins, reflecting the use of a similar carbonate chemistry scheme in these models (Séférian et al., 2020) and a similar general large-scale circulation.

The Atlantic Ocean has a large contribution to the ocean carbon–climate feedback relative to its size, as estimated from the carbon storage, in CMIP6 models, which is mainly due to the disequilibrium carbon pool. This result indicates that the Atlantic Ocean experiences a strong weakening in physical ventilation with climate change. Specifically, the negative γ^* per unit area is largest in the high latitudes of the North Atlantic in CMIP6 models, consistent with γ diagnosed from the air–sea carbon fluxes in previous-generation Earth system models (Roy et al., 2011; Ciais et al., 2013). This response is in accord with previous studies suggesting the effect of climate-driven changes in the circulation on the carbon storage is significant in the North Atlantic (Sarmiento et al., 1998; Winton et al., 2013; Bernardello et al., 2014). The Arctic Ocean also has a large contribution to γ^* relative to its size, primarily attributed to the saturated carbon pool and the large warming per unit volume in this basin.

The relative small contribution from the Southern Ocean to γ^* in CMIP6 is due to partial compensation between the combined decrease in solubility and physical ventilation with warming, driving a negative γ^*, and the large increase in accumulation of regenerated carbon in the ocean interior with climate change, driving a positive γ^*. This compensation between the decrease in solubility and physical ventilation with warming and the increase in accumulation of regenerated carbon with climate change in the Southern Ocean is also supported by sensitivity model experiments (Sarmiento et al., 1998; Bernardello et al., 2014). A large increase in the regenerated carbon pool of the Southern Ocean with climate change was also found in CMIP5 models (Schwinger et al., 2014; Ito et al., 2015). In the Indo-Pacific Ocean, the carbon–climate feedback in CMIP6 is primarily attributed to the saturated carbon pool and the decrease in solubility with warming. This response in the Indo-Pacific Ocean is consistent with a feedback triggered by warming and the reduction in solubility dominating in low and middle latitudes, as revealed in a CMIP5 model (Rodgers et al., 2020) and diagnostics in previous-generation Earth system models (Roy et al., 2011; Ciais et al., 2013). The inter-model variability in γ^* is mainly driven by the spread in the response of physical ventilation and circulation to climate change and in the parameterisation of ocean biogeochemical processes amongst the models.

5.2 Effect of the Atlantic Meridional Overturning Circulation

Our sensitivity experiments with an idealised climate model reveal that a weaker pre-industrial Atlantic Meridional Overturning Circulation (AMOC) leads to a smaller ocean carbon–concentration feedback parameter, β. This dependence of β on the pre-industrial AMOC is controlled by the disequilibrium carbon pool and the physical transfer of anthropogenic carbon from the surface to the ocean interior. The ocean carbon–climate feedback parameter, γ, is instead primarily controlled by the weakening in the AMOC with warming, which reduces the physical ventilation and leads to the disequilibrium part of γ becoming more negative in the idealised model.

Turning to the CMIP6 models, in the Atlantic Ocean, the carbon–concentration feedback parameter, β^*, is strongly correlated to the strength of the pre-industrial AMOC. This correlation of β^* to the AMOC is associated with the disequilibrium carbon pool and the effect of physical ventilation, such that Earth system models with a stronger AMOC provide more efficient transport of anthropogenic carbon into the ocean interior in the Atlantic Ocean, in accord with the behaviour in the idealised climate model. However, there is only a moderate correlation between AMOC strength and β on a global scale, which suggests that other processes related to ocean ventilation, involving the wind-driven gyre circulation, mode water formation and subduction, and the Southern Ocean residual circulation, are also important for controlling the ocean carbon–concentration feedback. In contrast with our results based on the idealised model and the CMIP6 models, Roy et al. (2011) found no direct dependence of β on the AMOC strength in older-generation Earth system models, which suggests that either (i) the sample of four models used in Roy et al. (2011) is too small to reveal the link between the AMOC and β or (ii) the dependence of β on AMOC strength is more pronounced when considering estimates at the quadrupling of atmospheric CO₂ rather than estimates from years 2010 to 2100 in a high-emission scenario as in Roy et al. (2011).

In the Atlantic Ocean, the carbon–climate feedback parameter, γ^*, is strongly correlated to the AMOC weakening with climate change in the CMIP6 models, similarly to the behaviour of the idealised climate model. This correlation of γ^* to AMOC weakening is (i) primarily due to the disequilibrium carbon pool, with a more pronounced AMOC weakening driving a larger reduction in the physical transfer of carbon into the ocean interior and a more negative γ^*, and (ii) to a lesser extent due to the regenerated carbon pool, with a more pronounced AMOC weakening driving a larger accumulation of regenerated carbon in the ocean interior and a more positive γ^*. The dependence of the carbon–climate feedback on the AMOC weakening in the Atlantic Ocean is consistent with experiments using an older-generation Earth system model (Crueger et al., 2008; Yoshikawa et al., 2008). In the Arctic Ocean, there is moderate correlation between γ^* and AMOC weakening in the CMIP6 models, with this correlation being associated with the saturated carbon pool and changes in carbonate chemistry. There is no significant correlation between γ^* and AMOC weakening in the Pacific, Indian and Southern oceans. Hence, on a global and centennial scale, the weakening in the AMOC due to climate change has only a modest impact on ocean carbon uptake, consistent with previous studies (Sarmiento et al., 1998; Joos et al., 1999; Zickfeld et al., 2008).

A caveat in our analysis is that it does not reveal the control of the carbon cycle feedbacks by other ventilation processes beyond the AMOC, such as the formation and subduction of mode and intermediate waters involving the seasonal mixed-layer cycle, as well as the horizontal and vertical circulation. An analysis for the ocean carbon uptake and storage in dynamically based density space (Sallée et al., 2013; Meijers, 2014; Iudicone et al., 2016) is necessary to reveal the contribution of different water masses to β and γ (Tilla Roy, personal communication, 2020; Roy et al., 2021), and to better constrain the effect of local versus remote processes on the regional ocean carbon uptake and storage. Although our diagnostics reveal the relative contribution from the air–sea fluxes and the ocean transport of carbon to the regional carbon storage, further analysis is necessary to identify the effect of the ocean transport of carbon and other traces (e.g. temperature, salinity and nutrients) to the local air–sea fluxes.

In summary, in the 100-year time frame of our analysis, the Atlantic and Arctic oceans provide a large contribution to the carbon cycle feedbacks, as estimated in terms of carbon storage, relative to their size in CMIP6. This large contribution from the Atlantic Ocean is associated with the physical ventilation. The Southern Ocean has a relatively small contribution to the carbon–concentration feedback as estimated in terms of carbon storage, despite its large anthropogenic carbon uptake, due to large carbon transport to the other basins. The Southern Ocean also has a relatively small contribution to the carbon–climate feedback due to competing processes affecting the carbon storage in this basin. The more poorly ventilated Indo-Pacific Ocean has a small contribution to the carbon cycle feedbacks relative to its size. On a global and centennial scale, the carbon cycle feedbacks depend on the strength of the pre-industrial AMOC and its weakening with warming, such that (i) a stronger pre-industrial AMOC leads to a more positive carbon–concentration feedback and (ii) a more pronounced AMOC weakening due to climate change leads to a more negative carbon–climate feedback. Finally, for a quadrupling of atmospheric CO₂ and on a centennial timescale, as considered in this study, the carbon–concentration feedback is substantially larger than the carbon–climate feedback; however, this will not necessarily be the case after the emissions cease and the system adjusts towards equilibrium.

Appendix A: Contribution from ocean carbon uptake and transport to the basin-scale β^* and γ^*

The Southern Ocean accounts on average for about half of the ocean anthropogenic carbon uptake (Fig. A1a red line) in CMIP6 Earth system models despite covering only 27 % of the global ocean, which is consistent with results from CMIP5 model during the historical period (Frölicher et al., 2015), as well as with observation-based estimates (Mikaloff-Fletcher et al., 2006). However, a large portion of the anthropogenic carbon entering through the Southern Ocean is transported into the Atlantic, Indian and Pacific oceans, with the Atlantic Ocean receiving most of this carbon (Fig. A1a) in CMIP6 models, which is consistent with observation-based estimates (Sabine et al., 2004; Khatiwala et al., 2009). Consequently, the carbon–concentration feedback parameter estimated based on the regional cumulative carbon uptake, β (Fig. A1b, red lines) is substantially larger in the Southern Ocean but smaller in the Atlantic, Pacific and Indian oceans than when estimated based on the regional carbon storage, β^* (Fig. A1b, black lines). The feedback parameter β is also smaller than β^* in the Arctic Ocean indicating that ocean transport supplies anthropogenic carbon to the Arctic in CMIP6 models.

Figure A1 Ocean carbon storage and cumulative air–sea carbon flux in the Pacific, Southern, Atlantic, Indian and Arctic oceans in CMIP6 Earth system models, along with the ocean carbon cycle feedback parameters: (a) changes in the ocean carbon storage, ΔI, and cumulative ocean carbon uptake, Δ𝒮, in petagrams of carbon, in the fully coupled simulations (COU) relative to the pre-industrial era; (b) carbon–concentration feedback parameter based on carbon storage, β^* (PgC ppm⁻¹), and based on the cumulative carbon uptake, β (PgC ppm⁻¹); and (c) carbon–climate feedback parameter based on carbon storage, γ^* (PgC K⁻¹), and based on the cumulative carbon uptake, γ (PgC K⁻¹). The solid lines and the shading show the model mean and the model range, respectively, based on the 1 % yr⁻¹ increasing CO₂ experiment over 140 years in 11 CMIP6 models (Table 1). The estimates for the feedback parameters are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC).

The net effect of climate change is a carbon loss from the ocean to the atmosphere and a negative γ in the Atlantic, Indian, Pacific and Southern oceans (Fig. A1c, red lines). The ocean transport acts to reduce the negative γ in the Atlantic, Indian and Southern oceans but to increase the negative γ in the Pacific Ocean (Fig. A1c, compare the red with the black lines). In contrast, in the Arctic Ocean, the decrease in sea-ice coverage due to climate change leads to an increase in the ocean carbon uptake from the atmosphere and a positive γ. The ocean transport opposes this ocean gain in carbon from the atmosphere in the Arctic and leads to a reduction in the ocean carbon inventory and a negative γ^* (Fig. A1c, black lines).

Appendix B: Contribution from saturated, disequilibrium and regenerated carbon pools to the regional β^* and γ^*

The saturated part of β^* is positive everywhere in the ocean (Fig. B1a). By definition, β_sat depends on the ocean volume, so it is larger in deeper regions. The disequilibrium part of β^* is small on the shelves and in the well-ventilated North Atlantic and large in the more poorly ventilated Pacific Ocean and in deep-ocean regions (Fig. B1a). This vertically integrated regional view for the carbon storage reflects that in the BGC simulation, ΔDIC_dis is small in the upper ocean but large and negative in deeper water outside the North Atlantic (see Figs. 10 and 11 in Arora et al., 2020). The regenerated part of β^* is negligible in the CMIP6 models.

The saturated part of γ^* is negative and large in the Arctic Ocean (Fig. B1b), reflecting a large decrease in solubility due to large warming per unit volume (Fig. 1d). Along the western boundary of the North Atlantic Ocean, γ_sat is slightly positive (Fig. B1b), which is related to a large regional gain in carbon (Fig. 2a), acting to reduce the carbonate system sensitivity to changes in temperature (see Eq. 20), combined with a regional increase in alkalinity due to climate change. The disequilibrium part of γ^* is negative and large in the North Atlantic, along the western boundary of the Atlantic Ocean and in the Southern Ocean, indicating that the effect of weakening in the ventilation is significant in these regions. In the Arctic Ocean, the decrease in sea-ice coverage drives an increase in the ocean carbon uptake and leads to areas with a large positive γ_dis. The regenerated part of γ^* has a similar spatial pattern but is of opposite sign to γ_dis (Fig. B1b), which suggests that γ_reg is dominated by changes in the ocean ventilation. Along the Equator and in the eastern boundary upwelling regions, γ_dis and γ_reg are slightly positive and negative, respectively. This negative γ_reg and positive γ_dis in regions of upwelling is consistent with a coupling between the regenerated and the disequilibrium carbon pools, such that a positive ΔDIC_reg brought into the surface waters from depth leads to an excess of carbon at the surface that is converted into a positive ΔDIC_dis (i.e. oversaturation) (Ito and Follows, 2013).

Figure B1 Geographical distribution of the CMIP6 inter-model mean contribution from the saturated, disequilibrium and regenerated carbon pools to the carbon cycle feedback parameters estimated based on carbon storage, β^* and γ^*, normalised by area: (a) β_sat, β_dis and β_reg (gC ppm⁻¹ m⁻²) and (b) γ_sat, γ_dis and γ_reg (gC K⁻¹ m⁻²). The estimates are based on the fully coupled simulation (COU) and the biogeochemically coupled simulation (BGC) under the 1 % yr⁻¹ increasing CO₂ experiment. Diagnostics are from years 121 to 140 (the 20 years up to quadrupling of atmospheric CO₂). The inter-model mean is based on the models in Table 1.