Drivers of decadal trends in the ocean carbon sink in the past, present, and future in Earth system models

The ocean and the land biosphere are the two major sinks of anthropogenic carbon at present. When anthropogenic carbon emissions become zero and temperatures stabilize, the ocean is projected to become the dominant and only global natural sink of carbon. Despite the ocean's importance for the carbon cycle and hence the climate, uncertainties about the decadal variability in this carbon sink and the underlying drivers of this decadal variability remain large because observing the ocean carbon sink and detecting anthropogenic changes over time remain challenging. The main tools that are used to provide annually resolved estimates of the ocean carbon sink over the last decades are global observation-based pCO₂ products that extrapolate sparse pCO₂ observations in space and time and global ocean biogeochemical models forced with atmospheric reanalysis data. However, these tools (i) are limited in time over the last 3 to 7 decades, which hinders statistical analyses of the drivers of decadal trends; (ii) are all based on the same internal climate state, which makes it impossible to separate externally and internally forced contributions to decadal trends; and (iii) cannot assess the robustness of the drivers in the future, especially when carbon emissions decline or cease entirely. Here, I use an ensemble of 12 Earth system models (ESMs) from phase 6 of the Coupled Model Intercomparison Project (CMIP6) to understand drivers of decadal trends in the past, present, and future ocean carbon sink. The simulations by these ESMs span the period from 1850 to 2100 and include four different future Shared Socioeconomic Pathways (SSPs), from low emissions and high mitigation to high emissions and low mitigation. Using this ensemble, I show that 80 % of decadal trends in the ocean carbon sink can be explained by changes in decadal trends in atmospheric CO₂ as long as the ocean carbon sink remains smaller than 4.5 Pg C yr⁻¹. The remaining 20 % are due to internal climate variability and ocean heat uptake, which result in a loss of carbon from the ocean. When the carbon sink exceeds 4.5 Pg C yr⁻¹, which only occurs in the high-emission SSP3-7.0 and SSP5-8.5, atmospheric CO₂ rises faster, climate change accelerates, and the ocean overturning and the chemical capacity to take up carbon from the atmosphere reduce, so that decadal trends in the ocean carbon sink become substantially smaller than estimated based on changes in atmospheric CO₂ trends. The breakdown of this relationship in both high-emission pathways also implies that the decadal increase in the ocean carbon sink is effectively limited to ∼1 Pg C yr⁻¹ dec⁻¹ in these pathways, even if the trend in atmospheric CO₂ continues to increase. Previously proposed drivers, such as the atmospheric CO₂ or the growth rate of atmospheric CO₂, can explain trends in the ocean carbon sink for specific time periods, for example, during exponential atmospheric CO₂ growth, but fail when emissions start to decrease again. The robust relationship over an ensemble of 12 different ESMs also suggests that very large positive and negative decadal trends in the ocean carbon sink by some pCO₂ products are highly unlikely and that the change in the decadal trends in the ocean carbon sink around 2000 is likely substantially smaller than estimated by these pCO₂ products.

1 Introduction

The ocean has taken up around one-quarter of all anthropogenic CO₂ emissions from land-use change and fossil fuels since the beginning of the Industrial Revolution (Friedlingstein et al., 2023; Gruber et al., 2023; Terhaar et al., 2022b). As such, in addition to the land biosphere, it is one of the two major natural sinks of carbon in the Earth system. If land-use change emissions are regarded as part of the land carbon sink, the land becomes almost neutral and the ocean carbon sink becomes the only major natural carbon sink (Friedlingstein et al., 2023). Once temperatures stabilize, the ocean will become the dominant global natural sink of carbon (Silvy et al., 2024) and will store more than half of the anthropogenically emitted CO₂ in around 1000 years (Joos et al., 2013). By taking up carbon from the atmosphere, the ocean effectively slows down global warming (IPCC, 2021) and will contribute to stabilizing global temperatures over the next centuries if emissions reach near-zero (Terhaar et al., 2023; MacDougall et al., 2020). Here, I define the ocean carbon sink as in the Global Carbon Budget as the change in air–sea CO₂ flux due to anthropogenic carbon emissions and anthropogenic climate change in comparison to a relatively stable pre-industrial state (Friedlingstein et al., 2023). Consequently, “anthropogenic” refers to direct effects from anthropogenic emissions and indirect effects from anthropogenically caused climate change.

The overall magnitude of the ocean carbon sink is mainly determined by the ocean overturning circulation, i.e., the rate at which surface waters with increased anthropogenic carbon content can be transported to the deep ocean and be replaced by waters with low anthropogenic carbon content (Sarmiento et al., 1992; Caldeira and Duffy, 2000; Orr et al., 2001). Furthermore, the magnitude of the ocean carbon sink is influenced by the surface ocean capacity to take up more anthropogenic carbon, which itself is determined by the surface ocean carbonate chemistry and especially the alkalinity (Broecker et al., 1979; Terhaar et al., 2022b). Over the historical period, the change in atmospheric CO₂ has been the main driver of changes in the ocean carbon sink and is assumed to be approximately proportional to the strength of the ocean carbon sink (Mikaloff Fletcher et al., 2006; Gruber et al., 2009; Terhaar et al., 2021b). However, this linear relationship between the strength of the ocean carbon sink and atmospheric CO₂ is only assumed to work under exponential atmospheric CO₂ growth (Raupach, 2013; Raupach et al., 2014; Ridge and McKinley, 2021).

Over the last decades, the relatively steady growth of the ocean carbon sink has been weakened by outgassing of natural carbon due to warming and climate change (Joos et al., 1999; McNeil and Matear, 2013; Frölicher et al., 2015) and is superimposed by decadal variability and trends in the ocean carbon sink, i.e., a reduction in the 1990s and an increase since 2000 (Lovenduski et al., 2008, 2007; Le Quéré et al., 2007; Landschützer et al., 2015, 2016). A consensus of the drivers of these trends is still not reached, and possible explanations for these different trends are changes in wind and pressure systems (Le Quéré et al., 2007; Keppler and Landschützer, 2019), variability in the ocean circulation and ventilation (DeVries et al., 2017), and changes in the atmospheric CO₂ growth rate and ocean surface temperature due to the eruption of Mount Pinatubo (McKinley et al., 2020; Frölicher et al., 2011). Moreover, recent studies suggest that the observation-based decadal variability over the last decades might be overestimated (Gloege et al., 2021; Hauck et al., 2023).

Despite the importance of the ocean carbon sink for the global climate and carbon cycle, observing or simulating the ocean carbon sink is still challenging. The two main tools to estimate the annually resolved ocean carbon sink over the past 3 to 7 decades, to provide an annual update every year within the Global Carbon Budget and to understand the drivers of the magnitude and trends in the ocean carbon sink, are observation-based pCO₂ products and global ocean biogeochemical models (GOBMs) (Friedlingstein et al., 2023; Hauck et al., 2023). Observation-based pCO₂ products extrapolate sparse observations of surface ocean pCO₂ using statistical methods and/or machine learning to create global monthly maps of surface ocean pCO₂ (Fay et al., 2021; Gregor and Gruber, 2021; Chau et al., 2022; Rödenbeck et al., 2015; Watson et al., 2020; Landschützer et al., 2015; Gloege et al., 2022; Bennington et al., 2022a, b). These monthly pCO₂ maps are then used to estimate the global ocean carbon uptake. In addition, an estimate of pre-industrial natural outgassing of CO₂ due to the difference in riverine carbon input and carbon sequestration in sediments (Lacroix et al., 2020; Regnier et al., 2022) has to be added to estimate the change in the air–sea CO₂ flux with respect to pre-industrial conditions. Global ocean biogeochemical models (GOBMs) (Orr et al., 2001; Hauck et al., 2020; Terhaar et al., 2024) simulate the ocean carbon sink while being forced with past observed atmospheric CO₂ and observation-based reanalysis data, such as wind, humidity, precipitation, and temperature (Hersbach et al., 2020; Tsujino et al., 2018).

The estimates of both product classes vary in magnitude and decadal trends, with pCO₂ products estimating a larger magnitude of the ocean carbon sink and also generally larger decadal trends over the last 2 decades (DeVries et al., 2023; Friedlingstein et al., 2023). One reason for the low carbon sink in GOBMs might be that the starting year of these simulations is often later than the beginning of the increase in atmospheric CO₂; thus the atmospheric CO₂ in the pre-industrial reference period and the pCO₂ in the ocean are too high (Terhaar et al., 2019, 2024; Bronselaer et al., 2017). Another reason for the low carbon sink is the existing biases in GOBMs in their simulated ocean circulation, especially Southern Ocean and Atlantic Ocean overturning, and the biases in the surface ocean carbonate chemistry (Terhaar et al., 2024). Similar biases were previously identified in Earth system model (ESM) ensembles (Terhaar et al., 2022b, 2021b; Goris et al., 2018). As opposed to the magnitude, the differences in the decadal trends between pCO₂ products and GOBMs might be due to uneven sampling of observations in space and in time, e.g., few observations in the 1980s and 1990s and few observations in the Southern Hemisphere, as demonstrated with a subset of pCO₂ products evaluated with output from a GOBM (Hauck et al., 2023) and ESMs (Gloege et al., 2021). In addition to differences in trends between pCO₂ products and GOBMs, no consensus has yet been made with respect to the underlying drivers of the decadal trends in the ocean carbon sink (Friedlingstein et al., 2023; DeVries et al., 2023; Gruber et al., 2023). The detection of these drivers with pCO₂ products, GOBMs, and other methods, such as data assimilation models (DeVries et al., 2017), is difficult due to the relatively short time period over which enough pCO₂ observations and atmospheric reanalysis data exist, due to the relative homogeneity of drivers over this period, e.g., constantly increasing atmospheric CO₂, and the absence of an alternative climate state against which these ocean carbon sink estimates can be compared.

Here, I use an ensemble of 12 ESMs from phase 6 of the Coupled Model Intercomparison Project (CMIP6) (Table 1) to provide a new perspective on potential drivers of the decadal trends in the ocean carbon sink, i.e., the underlying causes of its multi-decadal variability. Fully coupled ESMs are another tool to quantify and understand the ocean carbon sinks (e.g., Joos et al., 1999; McNeil and Matear, 2013; Frölicher et al., 2015; Goris et al., 2018; Terhaar et al., 2022b, 2021b). As ESMs are fully coupled and not forced with atmospheric reanalysis data, they do not simulate the same inter-annual internal climate variability as pCO₂ products and GOBMs do, and their biases of the surface ocean physics and biogeochemistry are thus larger than surface ocean biases of GOBMs (Terhaar et al., 2022b, 2024). However, ESMs have distinctive advantages compared to pCO₂ products and GOBMs for the analysis of decadal drivers of the ocean carbon sink because (1) they cover a period of 251 years from 1850 to 2100, (2) they cover at least four different future scenarios, and (3) they all have a different internal climate state. The long time period with different climate states in each model gives ample material to perform statistical analyses, and the different future scenarios allow testing to see how robustly potential drivers predict the decadal variability in the ocean carbon sink under continuously rising and strongly decreasing carbon emission trajectories. Furthermore, using different ESMs (e.g., Goris et al., 2018; Terhaar et al., 2022b, 2021b; McKinley et al., 2023) in comparison to large ensembles of one ESM (Fay et al., 2023; McKinley et al., 2016) avoids the risk of having a common bias in that one ESM, which might wrongly influence the analysis. Using the ESM ensemble from CMIP6, I will present how potential drivers of the ocean carbon sink (i.e., the atmospheric CO₂ and its growth rate, ocean heat uptake, and climate variability) drive trends in the ocean carbon sink from 1850 to 2100 in these models.

2 Methods and datasets

2.1 Earth system model ensemble

In this study, I use an ensemble of 12 ESMs from CMIP6 (Table 1). All ESMs from CMIP6 that provide the necessary model output for the following analysis were chosen. For each ESM, only the first ensemble member is used, as averaging over multiple ensemble members would have removed variability and using different numbers of ensemble members per ESM would have biased results towards the ESMs with more ensemble members.

Table 1 CMIP6 models used in this study and the corresponding model groups.

2.2 Calculating the ocean carbon sink

The annually averaged ocean carbon sink was calculated from concentration-driven historical simulations from CMIP (1850–2014) and four different concentration-driven Shared Socioeconomic Pathways (SSPs) (2015–2100): the low-emission high-mitigation SSP1-2.6, the high-emission low-mitigation SSP5-8.5, and the two intermediate pathways SSP2-4.5 and SSP3-7.0 (Riahi et al., 2017). To account for drifts in the historical and SSP simulations, a linear fit was calculated over the annual carbon sink over the years of the pre-industrial control run that correspond to the years 1850 to 2100 in the historical and SSP simulations. The linear change in the carbon sink in the pre-industrial simulations since 1850 was then subtracted from the historical and SSP simulations.

Furthermore, ESMs have biases in the magnitude of the ocean carbon sink due to biases in their respective circulation and surface ocean carbonate chemistry that also affect the size of the decadal trends; i.e., a negative bias in the magnitude of the carbon sink also introduces a negative bias in the decadal trends. To statistically compare the decadal trends in the carbon sink over the ESM ensemble used here, the global estimate of the ocean carbon sink was adjusted for each ESM with respect to biases in its circulation and surface ocean carbonate chemistry following Terhaar et al. (2022). Firstly, the Revelle factor, the inter-frontal Southern Ocean sea surface salinity, and the Atlantic Meridional Overturning Circulation (AMOC) strength were calculated for each model. Afterwards, a multi-linear fit was performed with the three observation-based quantities as predictors and the average ocean carbon sink from 1850 to 2100 as target variables (the period from 2015 to 2100 was used four times for each of the four SSPs). Finally, the biases in each predictor with respect to observation-based estimates of these predictors were calculated and used to adjust the simulated ocean carbon sink based on the determined constants from the multi-linear fit. Overall, this result in an adjustment of 10±7 % (i.e., increased ocean carbon uptake) for the model ensemble used here. The adjustment corrects for known biases in the models' circulations and surface ocean carbonate chemistry and hence reduced differences in the overall magnitude of the simulated carbon sink between ESMs (Terhaar et al., 2022b). This reduction in the difference in the magnitude of the carbon sink also reduces differences between the magnitude of trends and slightly improves the relationships found here (r² in Fig. 3 would have been 0.83 without adjustment instead of 0.91 with adjustment). Nevertheless, the results are quantitatively and qualitatively almost identical with and without that adjustment.

The ocean carbon sink was also calculated for each of the five major ocean basins (Atlantic Ocean, Pacific Ocean, Indian Ocean, Arctic Ocean, and Southern Ocean) using the RECCAP2 biome mask (DeVries et al., 2023; Terhaar et al., 2024), which is a slightly adapted version of a previously developed biome mask (Fay and McKinley, 2014). Regionally, no bias adjustments were performed, as it still remains largely unclear, apart from the Southern Ocean (Terhaar et al., 2021b), how biases in circulation and carbonate chemistry affect the regional ocean carbon sink estimates.

2.3 Atmospheric CO₂ and growth rate

The annually averaged atmospheric CO₂ over the historical period and for each SSP was taken from the CMIP6 forcing files (Meinshausen et al., 2020, 2017). The atmospheric CO₂ growth rate in each year was calculated as the difference in atmospheric CO₂ between that year and the year before.

2.4 Estimating the effect of climate change and ocean heat uptake on the ocean carbon sink

The effect of climate change and ocean heat uptake on the ocean carbon sink in ESMs was calculated based on additional idealized simulations provided by 5 of the 12 ESMs in the ensemble (ACCESS-ESM1-5, CanESM5, MRI-ESM2-0, NorESM2-LM, and UKESM1-0-LL) within the CMIP6 framework. These five ESMs made historical simulations, called “bgc”, where the change in atmospheric CO₂ has no effect on climate change but the carbon cycle still “sees” the increase in atmospheric CO₂. However, other non-CO₂ radiative agents (aerosols, CH₄, N₂O, etc.) still affect the climate in these simulations. These “bgc” simulations were only made for SSP5-8.5 (“ssp585-bgc”) and not for the other SSPs. The difference between the normal historical and SSP5-8.5 simulations (including effects from CO₂ and non-CO₂ radiative agents) and the respective additional “bgc” simulations quantifies the heat and carbon fluxes that are solely due to CO₂-induced climate change and warming.

2.5 Climate modes

To assess the climate variability across the ensemble of the ESM, annual averages of three climate modes were calculated for each ESM over the 251 years of the pre-industrial control simulation: (1) the Atlantic Multi-decadal Oscillation (AMO), (2) the Niño 3.4 index, and (3) the Marshall Southern Annular Mode (SAM) index. The AMO was calculated based on SST anomalies in the North Atlantic between 0 and 80° N. The Niño 3.4 index was calculated based on SST anomalies in the tropical Pacific region from 5° S to 5° N and from 170 to 120° W. The Marshall SAM index was calculated as anomalies of the zonal pressure difference between the latitudes of 40 and 65° S. Anomalies for each index in ESMs were calculated by removing a linear fit over the 251 years of the pre-industrial control simulation.

In addition, observation-based estimates of each climate mode were used. The observation-based AMO index (https://climatedataguide.ucar.edu/sites/default/files/2022-03/amo_monthly.txt, last access: 10 June 2024) and the Niño 3.4 index (https://psl.noaa.gov/gcos_wgsp/Timeseries/Data/nino34.long.data, last access: 10 June 2024) are based on HadISST1 (Rayner et al., 2003). The Marshall SAM index was calculated based on 12 stations: 6 stations at ∼40° S and 6 stations at ∼65° S (Marshall, 2003). To compare each observation-based index to the simulated index in the pre-industrial control simulations, the time series of the observation-based indexes have been detrended by subtracting a linear trend over the respective observation-based index estimate.

2.6 Decadal trends

Decadal trends in different variables are defined here as the slope of linear fits over 10 years.

2.7 Coefficient of determination, p-values, and Bayes factor

To determine the strength of correlations, the coefficient of correlation was calculated throughout this study (r²). In addition, the p-value was calculated to test the hypothesis that the trend change in atmospheric CO₂ from 1 decade to another decade is a significant driver of trends in the global and regional ocean carbon sink. A p-value larger than 0.1 indicates little or no evidence for that hypothesis, a p-value from 0.1 to 0.05 indicates weak evidence or a suggestion of evidence, a p-value from 0.05 to 0.01 indicates evidence or modest evidence, and a p-value from 0.01 to 0.001 indicates strong evidence (Held and Ott, 2018). In addition, an upper bound for the Bayes factor can be calculated following Halsey (2019). Throughout the paper, the p-values are never larger than $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{89}}$, resulting in Bayes factors that are at least 1×10⁸⁶. Based on the Bayes factor, the hypothesis that the trend change in atmospheric CO₂ from 1 decade to another decade is a significant driver of trends in the global and regional ocean carbon sink is 1×10⁸⁶ more likely than the hypothesis that the trend change in atmospheric CO₂ from 1 decade to another decade is not a significant driver of trends in the global and regional ocean carbon sink. As p-values are that small and Bayes factors are that high, I simply refrain from reporting that the p-values are smaller than 0.001.

3 The influence of atmospheric CO₂ on the ocean carbon sink

3.1 Atmospheric CO₂

Over the historical period of CMIP6 simulations from 1850 to 2014, the annually averaged global ocean carbon sink has increased approximately proportionally to the rise in atmospheric CO₂ (Fig. 1a, b, c). Due to the exponential rise in atmospheric CO₂, the cumulative ocean carbon sink is also approximately proportional to the rise in atmospheric CO₂. However, these quasi-linear relationships did not hold from 1920 (atmospheric CO₂ of 304 ppm) to 1960 (317 ppm) and from 1990 (354 ppm) to 1995 (360 ppm), when the ocean carbon sink did not increase while the atmospheric CO₂ continued to increase. These periods manifest themselves as “jumps” in the linear relationship between the atmospheric CO₂ and the cumulative ocean carbon sink (indicated by light-gray shading in Fig. 1d).

After 2014, when the historical period in CMIP6 ends and SSPs start, the link between atmospheric CO₂ and the ocean carbon sink depends strongly on the future scenario of atmospheric CO₂. The linear relationship between the annually averaged carbon sink and atmospheric CO₂ breaks down under all scenarios. Under SSP5-8.5, a pathway with a continuous increase in emissions (Riahi et al., 2017) and exponentially growing atmospheric CO₂ (Fig. 1a), the increase in the ocean carbon sink per increase in atmospheric CO₂ reduces until the ocean carbon sink reaches a maximum just above 6 Pg C yr⁻¹ (Fig. 1b), which is not exceeded even if atmospheric CO₂ rises (Fig. 1e). Under SSP3-7.0, a pathway with slightly smaller emissions and atmospheric CO₂ than SSP5-8.5, the ocean carbon sink also converges to a maximum but at around 5 Pg C yr⁻¹. Under SSP2-4.5, CO₂ emissions start to decline around 2050 (Riahi et al., 2017) and atmospheric CO₂ stabilizes around 600 ppm by 2100 (Fig. 1a). Although atmospheric CO₂ stabilizes, the ocean carbon sink reduces strongly (Fig. 1b). Under SSP1-2.6, atmospheric CO₂ not only stabilizes but starts to reduce by 2080, leading to a strong reduction in the ocean carbon sink (Fig. 1a). In comparison, the relationship between the cumulative ocean carbon sink and atmospheric CO₂ remains almost linear in the two high-emission pathways (SSP3-7.0 and SSP5-8.5), although the slope reduces with warming (Fig. 1f). For the two low-emission pathways (SSP1-2.6 and SSP2-4.5), the relationship breaks down as the ocean continuously takes up carbon, even when atmospheric CO₂ stabilizes and decreases (Fig. 1e, f).

Figure 1 The relationship between atmospheric CO₂ and the global ocean carbon sink. (a) The annually averaged atmospheric CO₂ that was used to force the ESMs from CMIP6 based on observation-based estimates from 1850 to 2014 (black) and based on four different SSPs (SSP1-2.6 in blue, SSP2-4.5 in orange, SSP3-7.0 in red, and SSP5-8.5 in brown) from 2015 to 2100. (b) The resulting ocean carbon sink as simulated by 12 ESMs (Table 1) after being adjusted for biases in circulation and surface ocean carbonate chemistry following Terhaar et al. (2022). The thick lines indicate multi-model means, and the shading indicates the 1σ standard deviation across the model ensemble. Relationships between atmospheric CO₂ and the annually averaged ocean carbon sink (c) for the historical period until 2014 and (e) for the 21st century from 2015 onwards, as well as between atmospheric CO₂ and the cumulative ocean carbon sink (d) for the historical period until 2014 and (f) for the 21st century from 2015 onwards. The light-gray shading in (a)–(d) indicates the time periods from 1920 to 1960 and from 1990 to 1995.

3.2 Atmospheric CO₂ growth rate

As an alternative to the atmospheric CO₂, the growth rate of atmospheric CO₂ was proposed as a key driver for the strength of the ocean carbon sink (McKinley et al., 2017, 2020). Over the historical period, the atmospheric CO₂ growth rate appears to be weakly linearly related to the strength of the ocean carbon sink (Fig. 2c). This relationship weakens after ∼1960, when the prescribed atmospheric CO₂ growth rate is based on direct atmospheric CO₂ observations and not, as before, on relatively smooth observation-based estimates from proxies (light-gray shading in Fig. 2a and gray dots in Fig. 2c, d). The direct observations capture the strong inter-annual variability in the atmospheric CO₂ growth rate that cannot be reconstructed by observation-based estimates from proxies. However, even this relationship between the atmospheric CO₂ growth rate and the strength of the ocean carbon sink breaks down in the 1920s and 1940s (pale-green dots in Fig. 2c, d), when the growth rate is around zero over around 1 decade each time (pale-green shading in Fig. 2a), but the ocean carbon sink does not go back close to zero and instead remains almost stable (pale-green shading in Fig. 2b).

Over the 21st century, the relationship between the ocean carbon sink and the atmospheric CO₂ growth rate breaks down (Fig. 2c). As long as CO₂ emissions and the atmospheric CO₂ growth rate rise, as they do under SSP3-7.0 and SSP5-8.5 (Fig. 2a), the relationship flattens (Fig. 2a). However, the strength of the relationship varies between SSP3-7.0 and SSP5-8.5. Under SSP5-8.5, the relationship also breaks down in the last 2 decades, when the atmospheric CO₂ growth stabilizes but the ocean carbon sink weakens. Under SSP1-2.6 and SSP2-4.5, with declining emissions or declining and even negative atmospheric CO₂ growth rates, the ocean carbon sink reduces – but not along the same path on which it increased over the historical period (Fig. 2c).

Figure 2 The relationship between the atmospheric CO₂ growth rate and the global ocean carbon sink. (a) The annually averaged atmospheric CO₂ growth rate based on atmospheric CO₂ forcing files from CMIP6, which are based on observation-based estimates from 1850 to 2014 (black) and based on four different SSPs (SSP1-2.6 in blue, SSP2-4.5 in orange, SSP3-7.0 in red, and SSP5-8.5 in brown) from 2015 to 2100. (b) The ocean carbon sink as simulated by 12 ESMs (Table 1) after being adjusted for biases in circulation and surface ocean carbonate chemistry following Terhaar et al. (2022). The thick lines indicate multi-model means, and the shading indicates the 1σ standard deviation across the model ensemble. Relationships between atmospheric CO₂ growth rate and the annually averaged ocean carbon sink (c) for the entire period from 1850 to 2100 and (d) only for historical period until 2014. The light-gray shading in panel (a) indicates the period where direction atmospheric CO₂ observations are available, and the pale-green shading in panels (a) and (b) and the pale-green dots in panels (c) and (d) indicate the 1920s and 1940s. The zero growth rate and ocean carbon sink in panels (a) and (b) are shown as dashed black lines.

3.3 Changes in atmospheric CO₂ growth rate determine changes in decadal trends in the ocean carbon sink

3.3.1 Global relationship

Although neither the atmospheric CO₂ nor its growth rate can quantify the strength of the ocean carbon sink over various time periods and different trajectories of atmospheric CO₂, the atmospheric CO₂ growth rate can nevertheless be used to understand changes in the ocean carbon sink on decadal timescales, i.e., decadal trends in the ocean carbon sink. For the period from 1980 to 2018, it has been shown that a slowing of the growth rate in comparison to a linear trend has led to a stagnation of the increase in the ocean carbon sink and that an accelerated increase in the growth rate has led to a strongly increasing carbon sink (McKinley et al., 2020).

Over longer time periods and different future SSPs, ESMs provide more such examples, where changes in the growth rate of atmospheric CO₂ led to changes in decadal trends in the simulated ocean carbon sink (Fig. 2). Around 1915, the atmospheric CO₂ growth rate changes from an increase to a decrease, and, at the same time, the ocean carbon sink stops increasing and starts to decrease. Then, in 1930, the atmospheric CO₂ growth rate increases and the ocean carbon sink also starts to increase simultaneously. Then, in 1940, the atmospheric CO₂ growth rate decreases again and the ocean carbon sink also decreases at the same time. Similarly, the atmospheric CO₂ growth rate changes from a positive trend to a negative trend in 1990, exactly when the ocean carbon sink also starts to slow down. When the atmospheric CO₂ growth rate increases again, the ocean carbon sink also increases. Over the 21st century, the same relationship continues. Under SSP2-4.5, the atmospheric CO₂ growth rate slows down until 2050 and the positive trend in the ocean carbon sink weakens. Once the atmospheric CO₂ growth declines, the trend in the ocean carbon sink becomes negative.

As a slowing or acceleration in the growth rate in comparison to a theoretical linear trend, as in McKinley et al. (2020), is not possible anymore over longer time periods of exponential growth or when atmospheric CO₂ peaks, I generalize the idea of McKinley et al. (2020) that a slowing or acceleration in the atmospheric CO₂ growth rate drives the trends in the ocean carbon sink by defining such slowing or acceleration as the difference in the growth rate in a given decade with respect to the preceding decade. When defining slowing or acceleration in the atmospheric CO₂ growth rate that way, a clear relationship (r²=0.91) emerges over the entire historical period and all four future scenarios over the 21st century (excluding years where the ocean carbon sink exceeds 4.5 Pg C yr⁻¹) between changes in the atmospheric CO₂ growth rate and the decadal trend in the multi-model average of the ocean carbon sink (Fig. 3). Changes in the atmospheric CO₂ growth rate are defined here as the change in trends (linear fit over 1 decade) in atmospheric CO₂ from 1 decade to the next. Trends in the ocean carbon sink are a linear fit over annual values of the global ocean carbon sink in the second decade. It thus appears that it is the change in the growth rate in comparison to the previous decade that drives the decadal trends in the ocean carbon sink and not the difference to an expected linear trend. If, for example, the growth rate strongly reduced from 1 decade to another, the ocean carbon sink would show a negative trend. If the growth rate then stays at that lower level, the carbon sink would not decline further but would stabilize at its new level. This relationship even holds when CO₂ emissions decline strongly, as under SSP1-2.6.

Figure 3 The relationship between changes in the atmospheric CO₂ growth rate and decadal trends in the global ocean carbon sink for the multi-model mean. (a) Decadal trends in the multi-model mean ocean carbon sink compared to changes in decadal trends in atmospheric CO₂, which represent the decadal averaged growth rate of atmospheric CO₂. The dark-blue to yellow circles without a surrounding black line show multi-model averages for all years of the historical period from 1850 to 2014 and for all years from 2015 to 2100 for all four SSPs. All decades from 1850 to 2100 are shown, i.e., 2000–2009, 2001–2010, 2002–2011, etc. The brown line shows a linear fit for all years when the global ocean carbon sink is smaller than 4.5 Pg C yr⁻¹, and the brown shading is the 1σ projection uncertainty. The dots with black lines around them show values from the respective ensemble means of the pCO₂ products (pink) and GOBMs (orange) from the Global Carbon Budget 2023 (Friedlingstein et al., 2023) for the 3 decades between 1990 and 2020. Small deviations from the relationship in SSP1-2.6 are marked by “SSP1-2.6”. (b) The simulated ocean carbon sink in comparison to the expected ocean carbon sink based on the relationship in panel (a) and the prescribed trend change in atmospheric CO₂ in the simulations.

However, this relationship between changes in the atmospheric CO₂ growth rate and the decadal trend in the multi-model average of the ocean carbon sink breaks down if the ocean carbon sink is larger than 4.5 Pg C yr⁻¹ (Fig. 3; r² starts to reduce if years with an ocean carbon sink larger than 4.5 Pg C yr⁻¹ are included). The breakdown likely occurs because climate change and the associated ocean heat uptake and circulation changes become so large that effects on the natural carbon sink reduce the trend in the ocean carbon sink substantially enough. Thus, it is not the carbon uptake of 4.5 Pg C yr⁻¹ itself that causes the breakdown in the relationship but the combined impact of an increasing Revelle factor (Revelle and Suess, 1957) and climate change (Joos et al., 1999; McNeil and Matear, 2013; Frölicher et al., 2015). In SSPs from CMIP6, the combined impact is large enough to affect the relationship identified here when the ocean carbon sink is around 4.5 Pg C yr⁻¹. The breakdown in the relationship also implies that the decadal trends in the ocean carbon sink cannot exceed 0.78±0.10 Pg C yr⁻¹ dec⁻¹ (the uncertainty is the 1σ standard deviation across the ESM ensemble in the decade when the multi-model mean decadal trend is largest). Thus, if the ocean carbon sink is below 4.5 Pg C yr⁻¹, and if its magnitude 10 years ago and the change in the decadal trends in atmospheric CO₂ between the last 2 decades (20 to 10 years ago and 10 years ago to now) is known, the absolute ocean carbon sink this year can be determined (Fig. 3b).

3.3.2 Regional relationships

The relationship between changes in the atmospheric CO₂ growth rate and the global ocean carbon sink holds in all five large ocean basins (Fig. 4a–j) as it has also done from 1980 to 2018 (McKinley et al., 2020). The coefficient of determination (r²) is larger than 0.84 in the Atlantic, Pacific, Indian, and Southern oceans. Only the Arctic Ocean has a smaller correlation coefficient of 0.66.

In the Arctic Ocean, the carbon sink has already been shown to be substantially more affected by climate change than in any other ocean basin (Yasunaka et al., 2023). In the future, when sea ice disappears and the Arctic continues to warm faster than any other region, the importance of climate change for the Arctic Ocean carbon sink will likely remain relatively large, for example, through freshening (Terhaar et al., 2021a) and a change in the seasonal cycle of pCO₂ (Orr et al., 2022), and hence reduce the importance of changes in the atmospheric CO₂ for trends in the ocean carbon sink.

In the Southern Ocean, the simulated trends in the ocean carbon also slightly differ from the expected trends based on changes in trends in atmospheric CO₂ in three brief periods (Fig. 4i). From 1995 to 2005 over the historical period and from 2030 to 2050 under SSP1-2.6, the decadal trend in the ocean carbon sink is larger than expected, whereas it is smaller than expected from 2080 to 2100 under SSP1-2.6. The differences under SSP-1.2.6 are even visible for the global carbon sink (Fig. 3a). As the difference occurs in the multi-model mean, it appears to be a forced response and not a response that is linked to the particular state of the climate in one of the models. The time periods where the differences are visible globally (2030–2050 and 2080 to 2100 under SSP1-2.6) are the times when the growth in atmospheric CO₂ stops and when it starts to decrease in that scenario (Fig. 1c). As the atmospheric CO₂ growth rate changes quickly in these periods (Fig. 2a), first by changing into a decreasing phase and then by transitioning into a stabilizing phase, it appears that a fast transition of the trend change in atmospheric CO₂ temporarily leads to differences in the expected relationship. If the trend change in atmospheric CO₂ decreases quickly, the trend in the ocean carbon sink remains larger than expected, and, if the trend change in atmospheric CO₂ increases quickly, the trend in the ocean carbon sink remains smaller than expected. However, the drivers behind the divergence from the expected decadal trend in the multi-model mean from 1995 to 2005 in the Southern Ocean remain unclear and should be analyzed in future research.

Despite these small differences, the overall relationship between changes in decadal trends in the atmospheric CO₂ and decadal trends in the local and global ocean carbon sink is very strong (r²>0.84, apart from in the Arctic Ocean) and demonstrates how atmospheric CO₂ is the main driver of the externally forced decadal trends in the ocean carbon sink.

Figure 4 The relationship between changes in the atmospheric CO₂ growth rate and decadal trends in the ocean carbon sink in the five major ocean basins. Decadal trends in the ocean carbon sink in the (a) Atlantic Ocean, (c) Pacific Ocean, (e) Indian Ocean, (g) Arctic Ocean, and (i) Southern Ocean compared to changes in decadal trends in atmospheric CO₂, which represent the decadal averaged growth rate of atmospheric CO₂. The dark-blue to yellow circles without a surrounding black line show multi-model averages for all years of the historical period from 1850 to 2014 and for all years from 2015 to 2100 for all four SSPs. The brown line shows a linear fit for all years when the global ocean carbon sink is smaller than 4.5 Pg C yr⁻¹, and the brown shading is the 1σ projection uncertainty. The dots with black lines around them show values from pCO₂ products (pink) and GOBMs (orange) from the Global Carbon Budget 2023 (Friedlingstein et al., 2023) for the 3 decades between 1990 and 2020. The simulated ocean carbon sink in the (b) Atlantic Ocean, (d) Pacific Ocean, (f) Indian Ocean, (h) Arctic Ocean, and (j) Southern Ocean in comparison to the expected ocean carbon sink based on the respective relationships in panels (a), (c), (e), (g), and (i) and the prescribed trend change in atmospheric CO₂ in the simulations. The p-value for each regional relationship is smaller than 0.001.

4 The importance of climate variability on decadal trends on the ocean carbon sink

Internal climate and ocean variability in ESMs can strongly affect the decadal trends in ESMs (Li and Ilyina, 2018) and hence reduces the strength of the relationship between changes in decadal trends in the atmospheric CO₂ and decadal trends in the ocean carbon sink. To quantify the importance of climate variability, I calculated the relationship between changes in decadal trends in the atmospheric CO₂ and decadal trends in the ocean carbon sink not for the multi-model mean but for the individual ESMs. When calculating the linear fit over the results from the individual ESMs, the correlation factor only slightly reduces from r²=0.91 to r²=0.80 (Fig. 5). Although Li and Ilyina (2018) showed that 50–70 ensemble members are needed to remove internal variability for 1 given decade, the relationship remains strong here because it is based on at least 34 decades for each ESM (16 historical plus 18 future decades (SSP1-2.6 and SSP2-4.5)). Thus, two ESMs alone already provide the 50–70 decades that were shown to be necessary (Li and Ilyina, 2018). The 1σ prediction interval around the linear fit is 0.16 Pg C yr⁻¹ dec⁻¹, meaning that 68 % of all trends will be within ±0.16 Pg C yr⁻¹ dec⁻¹ of the predicted trend based on decadal trends in the atmospheric CO₂, 95 % will be within ±0.31 Pg C yr⁻¹ dec⁻¹ of the predicted trend, and virtually all trends (99.7 %) will be within ±0.47 Pg C yr⁻¹ dec⁻¹ of the predicted trend. The largest simulated trend in the ocean carbon sink in one of the ESMs is 0.96 Pg C yr⁻¹ dec⁻¹. This is within the 2σ range of the largest trend as diagnosed by the multi-model mean 0.78±0.10 Pg C yr⁻¹ dec⁻¹.

Figure 5 The relationship between changes in the atmospheric CO₂ growth rate and decadal trends in the global ocean carbon sink for individual ESMs. (a) Decadal trends in the ocean carbon sink for all ESMs individually compared to changes in decadal trends in atmospheric CO₂, which represent the decadal average growth rate of atmospheric CO₂. The dark-blue to yellow circles without a surrounding black line show multi-model averages for all years of the historical period from 1850 to 2014 and for all years from 2015 to 2100 for all four SSPs. The brown line shows a linear fit for all years when the global ocean carbon sink is smaller than 4.5 Pg C yr⁻¹, and the brown shading is the 1σ projection uncertainty. The dots with black lines around them show values from individual pCO₂ products (pink) and GOBMs (orange) from the Global Carbon Budget 2023 (Friedlingstein et al., 2023) for the 3 decades between 1990 and 2020. (b) The simulated ocean carbon sink in comparison to the expected ocean carbon sink based on the relationship in panel (a) and the prescribed trend change in atmospheric CO₂ in the simulations.

The range of simulated trends in ocean carbon sink with different internal climate variability encompasses the ocean carbon sink trend estimates of GOBMs from the Global Carbon Budget 2023 (Friedlingstein et al., 2023), but the trend estimates of the pCO₂ products exceed the range that is simulated by ESMs. For the decade from 1990 to 1999, 7 out of 10 GOBMs fall within the ±1σ range and the remaining 3 GOBMs fall within the ±2σ range. In comparison, only 3 out of 8 pCO₂ products fall within the ±1σ range, 2 pCO₂ products fall within the ±3σ range, 1 pCO₂ product falls within the ±4σ range, and 2 pCO₂ products fall within the ±5σ range. For the decade from 2000 to 2009, 5 GOBMs fall within the ±1σ range, 4 GOBMs fall within the ±2σ range, and the final GOBM falls within the ±3σ range. In comparison, only 2 pCO₂ products fall within the ±1σ range, 1 pCO₂ product falls within the ±2σ range, 2 pCO₂ products fall within the ±3σ range, 1 pCO₂ product falls within the ±5σ range, 1 pCO₂ product falls within the ±6σ range, and 1 pCO₂ product falls within the ±7σ range. For the decade from 2010 to 2019, 9 GOBMs fall within the ±1σ range and the remaining one falls within the ±2σ range. In comparison, only 1 pCO₂ product falls within the ±1σ range, 2 pCO₂ products fall within the ±2σ range, 3 pCO₂ products fall within the ±3σ range, 1 pCO₂ product falls within the ±4σ range, and 1 pCO₂ product falls within the ±5σ range. In general, the pCO₂ product estimates of the decadal trends are not randomly distributed across the possible range that the ESMs suggest. Instead, pCO₂ products systematically overestimate the magnitude of the respective trends that is suggested by ESMs, i.e., a negative trend in the 1990s that is too low and a positive trend in the 2000s and 2010s that is too high.

5 Imprint of climate change and ocean heat uptake on the ocean carbon sink

In addition to atmospheric CO₂ and internal climate variability, climate change and ocean heat uptake also affect the ocean carbon sink and potentially its decadal trends. The ocean heat uptake, for example, causes changes in the ocean circulation, such as stratification and outgassing of natural carbon from the ocean, due to increasing temperatures and reduced solubility (Fig. 6). Across the five ESMs that performed the simulations to quantify the effect of ocean heat uptake on the natural carbon in the ocean (see Methods), the loss of natural carbon from the ocean to the atmosphere is related to the ocean heat uptake via a second-degree polynomic function under strong radiative forcing (SSP5-8.5) (Fig. 6a). Although annual variability hides part of this relationship, the relationship emerges strongly for decadal averages (Fig. 6b).

Figure 6 The relationship between natural carbon loss and ocean heat uptake and its link to the Atlantic Meridional Overturning Circulation (AMOC). (a) Annual and (b) decadal average air–sea CO₂ flux solely caused by ocean heat uptake (for details see Methods) vs. the annual average ocean heat uptake in five different ESMs (ACCESS-ESM1-5 in blue, CanESM5 in orange, MRI-ESM2-0 in green, NorESM2-LM in purple, and UKESM1-0-LL in brown). A second-degree polynomic function (colored lines) was fitted over the decadal average values of each ESM. (c) The CO₂ flux per ocean heat uptake, calculated for each model from the fitted second-degree polynomic function at an ocean heat uptake of 35 ZJ yr⁻¹ against the AMOC at 26.5° N calculated in each ESM from 2004 to 2018 (historical plus SSP5-8.5 simulations). The vertical dashed black line and the gray shading indicate the average observed AMOC at the RAPID array from 2004 to 2018 (McCarthy et al., 2020).

While each of the five ESMs suggests that the loss of natural carbon from the ocean to the atmosphere is related to the ocean heat uptake via a second-degree polynomic function, the amount of carbon loss due to ocean heat uptake varies across ESMs (Fig. 5b). The main reason for these differences is likely the different changes in ocean circulation and stratification due to ocean heat uptake in each ESM. One aspect of the ocean overturning circulation that is expected to change a great deal with climate change and ocean warming in ESMs is the Atlantic Meridional Overturning Circulation (AMOC; Weijer et al., 2020). Across a large ensemble of ESMs from CMIP6, it has been shown that ESMs with an already stronger AMOC also show a larger reduction in the AMOC (Weijer et al., 2020). The larger overturning reduction thus causes the models with a higher AMOC to lose more natural carbon for the same heat uptake (Fig. 5c). Based on this linear relationship, it would be possible to constrain the loss of carbon per heat uptake with observations of the present-day AMOC. However, using only five ESMs to quantify a linear relationship likely does not yield a robust relationship, so I abstain from constraining the loss of carbon per heat uptake. Nevertheless, the observed AMOC at 26.5° N is close to the average of the simulated AMOC at 26.5° N in ESMs, suggesting that the multi-model average sensitivity of air–sea CO₂ fluxes to heat uptake is a good approximation of the real-world sensitivity.

Unfortunately, CMIP6 only provides simulations that allow the quantification of the ocean's natural carbon sink response to ocean warming for SSP5-8.5 and not for other scenarios where the ocean warming slows down or even stabilizes. Thus, it remains impossible for now to quantify the effect of ocean heat uptake for other scenarios and to test if the relationship identified here is robust. However, as differences in the decline in the AMOC are similar across all these scenarios, although the ocean heat uptake is much smaller in the low-emission scenarios (Weijer et al., 2020), the sensitivity of carbon loss to heat uptake might be greater in low-emission scenarios. As such, changes in the ocean heat uptake and its trend might well cause changes in the anthropogenic ocean carbon sink via outgassing of natural carbon. Although these changes are likely small, as decadal average ocean heat uptake does not change quickly. However, these changes might still be partly responsible for differences between the decadal trend in the ocean carbon sink that were expected based on changes in the trend in atmospheric CO₂ and in the decadal trend in the simulated ocean carbon sink, especially those in SSP1-2.6 globally (Fig. 3) and in the Southern Ocean (Fig. 4i). To verify this hypothesis, CMIP simulations allowing the quantification of the ocean's natural carbon sink response to ocean warming would have to be run for other scenarios than SSP5-8.5.

6 Potential caveats and limitations

The strong dependence of decadal trends in ocean carbon sinks on the change in the atmospheric CO₂ growth rate from 1 decade to the other was identified here across an ensemble of state-of-the-art ESMs from CMIP6. The robustness of this relationship depends on the model's ability to represent the internal climate variability and might also be biased if the entire model ensemble is biased, for example, due to relatively coarse resolution or a common unrealistic representation of the physics or biogeochemistry.

If, for example, the internal climate variability on decadal timescales was underestimated by the ESMs used here, the predictability of the decadal trends in the ocean carbon sink by changes in the growth rate of CO₂ would be overestimated. A prerequisite for ESMs to simulate the contribution of the natural variability to decadal trends in the ocean carbon sink is that they also simulate the size of the decadal trends in internal climate modes that are known to most affect the variability of the ocean carbon sink. The major climate modes that are known to influence the decadal variability in the ocean carbon sink are the Niño 3.4 index (Feely et al., 1999; Ishii et al., 2014; McKinley et al., 2004), the Atlantic Multi-decadal Oscillation (Breeden and McKinley, 2016; Keppler et al., 2023), and the Marshall Southern Annular Mode (Le Quéré et al., 2007; Gruber et al., 2019b; Landschützer et al., 2015; Lovenduski et al., 2008; Thompson and Solomon, 2002; Lenton and Matear, 2007; Hauck et al., 2013). The decadal trends in the Niño 3.4 index in ESMs are 17 % (−19 % to 53 %) larger than the decadal trends in the observation-based estimates of the Niño 3.4 index (the numbers in parentheses indicate the standard deviation across ESMs) (Fig. 7b), the decadal trends in the Atlantic Multi-decadal Oscillation in ESMs are 13 % (−7 % to 33 %) larger than the decadal trends in the observation-based estimates (Fig. 7c), and the decadal trends in the Marshall Southern Annular Mode in ESMs are 52 % (30 % to 75 %) larger than the decadal trends in the observation-based estimates (Marshall, 2003) (Fig. 7d). The relatively large decadal trends in climate modes in ESMs suggest that the ESMs are indeed capable of simulating the internal climate variability on decadal timescales. As the decadal trends in climate mode are larger than or equal to the observed ones, there is no indication that the decadal variability in the ocean carbon sink in ESMs (Fig. 7a) might be too small because of an internal climate variability in ESMs that is too small, as previously hypothesized by Gruber et al. (2023) based on small carbon–climate feedbacks in idealized scenarios without variability in the atmospheric CO₂ growth (Arora et al., 2020).

Figure 7 Time series and decadal trends in the ocean carbon sink and climate modes in Earth system models compared to observations. (a) The globally integrated annual air–sea CO₂ flux in the pre-industrial control simulations for 12 ESMs (thin blue lines) and the multi-model average (thick blue line). The same is shown for (b) the Niño 3.4 index, (c) the Atlantic Multi-decadal Oscillation, and (d) the Marshall Southern Annular Mode. For the three climate modes, observation-based estimates are shown based on HadISST1 for panels (b) and (c) and based on Marshall (2003) for panel (d). The decadal trends in these observation-based estimates (orange numbers) are compared to the decadal trends in ESM estimates (blue numbers indicating average and standard deviation across the ESM ensemble).

However, in addition to physical climate-driven variability, which is the dominant driver of variability in the ocean carbon sink (Doney et al., 2009), there is also biology- and biogeochemical-driven climate-related variability in the air–sea CO₂ fluxes (Ostle et al., 2022; Doney et al., 2009; Keller et al., 2012), for example, due to changes in net primary production or remineralization caused by changes in nutrient supply, temperature, and oxygen. Over the North Atlantic, it has been shown that biogeochemical variability is also strongly influenced by climate modes, such as the Atlantic Multi-decadal Oscillation (Ostle et al., 2022). Nevertheless, GOBMs suggest that the influence of physical variability exceeds the influence of biogeochemical variability (Doney et al., 2009; DeVries et al., 2023). Despite different representations of the biogeochemistry and biology across the models from RECCAP2 (DeVries et al., 2023; Rodgers et al., 2023), they all simulate a similar inter-annual and decadal variability and trends in the ocean carbon sink (DeVries et al., 2023; Terhaar et al., 2024), as they are forced with historical atmospheric reanalysis products and share the same internal climate modes. Although the similarity of all GOBMs when forced with historical reanalysis strongly suggest that the physical impact on decadal variability exceeds the biogeochemical impact, detailed regional analyses of the biogeochemical-driven climate-related variability in the air–sea CO₂ fluxes (Ostle et al., 2022; Keller et al., 2015), which exceed the scope of this paper, are necessary. Overall, the dominance of physical variability over biogeochemical variability and the larger decadal trends in climate modes in ESMs than in the real world suggest that the ESMs do not underestimate the natural variability in the ocean carbon sink, although it always remains possible that the pCO₂ products capture real signals that are not yet simulated.

Although the ESMs from CMIP6 used here simulate even larger decadal trends in important climate modes, they might still underestimate decadal trends in the ocean carbon sink driven by climate variability because of their resolution that has increased over the past decades but is still too coarse to explicitly resolve mesoscale ocean eddies. Higher-resolved ESMs are still computationally too expensive to be run within the CMIP framework with the sufficiently long spin-ups that are necessary for these models to be in equilibrium (Séférian et al., 2016; Gupta et al., 2013). In some studies with fewer simulations than required for CMIP, higher-resolved ocean models have been shown to affect the ocean carbon sink and physics and their variability (Lachkar et al., 2007, 2009; Dufour et al., 2015; Griffies et al., 2015). While it remains impossible to evaluate the effect of higher resolution over a large ensemble of ESMs, as such an ensemble does not exist yet, I tested the relationship identified here with the highest-resolved Earth system model within CMIP6, GFDL-CM4 (Held et al., 2019), which has a horizontal resolution of 0.25° that allows resolving eddies in tropical and subtropical oceans but still has to parametrize some eddy activity in subpolar and polar oceans. GFDL-CM4 was not included in the overall analyses as it did not provide simulations under SSP1-2.6 and SSP3-7.0, presumably because of its large computational costs. The trends in the ocean carbon sink in GFDL-CM4 lie mostly within ±1σ of the relationship between changes in the atmospheric CO₂ growth rate and trends in the ocean carbon sink, with only a few decades being in the ±2σ range (Fig. 8). As for the other ESMs, the relationship in GFDL-CM4 only holds if the ocean carbon sink remains below 4.5 Pg C yr⁻¹. Although a potential change in this relationship at an even higher resolution cannot be excluded with certainty until simulations with higher resolution are performed, the robustness of the relationship even for higher-resolved ESMs such as GFDL-CM4 gives no indication that the relationship will not hold at even higher resolution.

Figure 8 The relationship between changes in the atmospheric CO₂ growth rate and decadal trends in the global ocean carbon sink in a high-resolution model. (a) Decadal trends in the multi-model mean ocean carbon sink compared to changes in decadal trends in atmospheric CO₂, which represent the decadal average growth rate of atmospheric CO₂. The dark-blue to yellow circles without a surrounding black line show multi-model averages for all years of the historical period from 1850 to 2014 and for all years from 2015 to 2100 for all four SSPs. The brown line shows a linear fit for all years when the global ocean carbon sink is smaller than 4.5 Pg C yr⁻¹, and the brown shading is the 1σ projection uncertainty. The red dots with black lines around them show values for the high-resolution ESM, GFDL-CM4 (Held et al., 2019), under the historical simulation, SSP2-4.5, and SSP5-8.5.

7 Discussion and conclusion

The analysis with ESM suggests that changes in the atmospheric growth rate of CO₂ can indeed explain most of the decadal trends in the ocean carbon sink, as previously proposed by McKinley et al. (2020). ESMs support the hypothesis by McKinley et al. (2020) that the weak decadal trend in the 1990s and the stronger trend in the 2000s is mainly driven by changes in the atmospheric CO₂ growth rate (Figs. 3 and 5). The importance of atmospheric CO₂ extends over all ocean basins, as also previously shown by McKinley et al. (2020). While McKinley et al. (2020) focused on the last decades and suggested that the trends in the ocean carbon sink depend on differences in the atmospheric growth rate of CO₂ compared to the long-term trend in the growth rate, I could generalize this idea here and show that it is the change in the growth rate compared to the previous decade that drives the trends in the ocean carbon sink over a wide range of timescales and SSPs. Moreover, this analysis extends the timeline of previous analyses and shows how atmospheric CO₂ drives trends in the ocean carbon sink in a range of different future scenarios, from high-mitigation low-emission scenarios to high-emission scenarios. In addition, the use of ESMs allowed the quantification of the link between changes in the atmospheric CO₂ growth rate and decadal trends in the ocean carbon, which allows for a better separation of the effects of internal and external forcing of past decadal trends in the ocean carbon sink. However, if atmospheric CO₂ rises too high and the impact of climate change on the ocean carbon sink increases, atmospheric CO₂ is no longer the dominant driver of trends in the ocean carbon sink due to changes in the buffer factor and ocean ventilation (Revelle and Suess, 1957; Heinze et al., 2015; Joos et al., 1999; McNeil and Matear, 2013; Frölicher et al., 2015).

Although atmospheric CO₂ is shown here to be the main driver of the decadal trends in the ocean carbon sink, climate variability also plays an important role in the decadal trends in the ocean carbon sink. With a standard deviation of ±0.16 Pg C yr⁻¹ across all ESMs, climate variability drives 17 % of the trend in the ocean carbon sink when the growth rate of atmospheric CO₂ is largest and drives all changes in carbon trends when the change in the growth rate of atmospheric CO₂ is zero. Known drivers of this internal climate variability are, for example, El Niño (Feely et al., 1999; Ishii et al., 2014; McKinley et al., 2004), the Atlantic Multi-decadal Oscillation (Breeden and McKinley, 2016; Keppler et al., 2023), changes in the overturning circulation (DeVries et al., 2017), and changes in the Southern Annular Mode associated with changes in Southern Ocean winds (Le Quéré et al., 2007; Keppler and Landschützer, 2019; Landschützer et al., 2015) and with stronger consequent upwelling of older waters (Lovenduski et al., 2008, 2007), as well as changes in the Southern Ocean stratification (Gruber et al., 2019b). Across the ESMs used here, the variability in decadal trends is highest in the Southern Ocean, followed by the tropical regions, and again followed by the northern subpolar gyres (Fig. 9); confirming that the decadal trends in the ocean carbon sink are indeed most variable due to internal climate variability in the regions where they are expected based on the previous studies mentioned above. ESMs simulate decadal trends in important climate modes that are even larger than their observation-based counterparts, which suggests that ESMs also capture the climate-driven decadal trends in the ocean carbon sink. As ESMs slightly overestimate decadal trends in important climate modes and still suggest that changes in the atmospheric growth rate are the dominant drivers of decadal trends in the ocean carbon sink, climate variability and associated changes in ocean circulation appear not to be the first-order drivers of decadal trends in the ocean carbon sink over the last decades as previously suggested (Landschützer et al., 2015; DeVries et al., 2017, 2019).

The results presented here have implications for previous estimates of the ocean carbon sink, especially those from pCO₂ products that suggested very strong decadal trends in the ocean carbon sink (Landschützer et al., 2015; Gruber et al., 2019b, 2023). The trend estimates of the ocean carbon sink by pCO₂ products are larger than the likely trends based on the relationship identified here between changes in the atmospheric CO₂ growth rate and the decadal trends in the ocean carbon sink. Thus, either pCO₂ products overestimate decadal trends in the ocean carbon sink or ESMs underestimate these trends. In each decade from 1990 to 2019, there are five out of eight pCO₂ products from the Global Carbon Budget (Friedlingstein et al., 2023) that estimate decadal trends that are outside of the 2σ range that is estimated here based on ESMs (Fig. 5), giving these results a likelihood of less than 5 % to occur if the ESM results are indeed robust. Some pCO₂ products estimate trends that are within the 5σ, 6σ, and 7σ ranges that correspond to events that occur once every 4776 years (5σ), once every 1.38 million years (6σ), and once every 1.07 billion years (7σ). While it is already extremely unlikely that decadal trends in all 3 decades from 1990 to 2020 lie outside the 2σ range, the estimates within the 5σ to 7σ range are virtually impossible based on the ESM-derived range. Only two pCO₂ products (NIES-ML3 from Zeng et al., 2022, and OS-ETHZ-Gracer from Gregor and Gruber, 2021) lie within the 1σ and 2σ ranges in the 1990s and 2000s and only very slightly above the 2σ range in the 2010s. The slightly higher trend in the 2010s in these products may very well be a consequence of the uneven sampling in space and time (Hauck et al., 2023). While the trends in these two pCO₂ products are closer to what is expected based on ESMs, only an in-depth analysis will eventually allow judging the performance of each pCO₂ product with certainty.

Figure 9 Variability in the decadal trends in the zonally integrated ocean carbon sink in Earth system models. The multi-model mean (thick blue line) and the 1σ standard deviation of the variability in the zonally integrated ocean carbon sink across the 251 years of the pre-industrial control simulation across all 12 ESMs. In addition, the maximum and minimum variability in the ESMs are shown at each latitude (thin blue lines).

Here, I have demonstrated that ESMs are capable of simulating the size of decadal trends in important climate modes that have a strong impact on the variability in the ocean carbon sink (Fig. 7) and that higher resolution does not alter the identified relationship (Fig. 8). While this analysis does not guarantee that ESMs do not underestimate the decadal trends in the ocean carbon sink, it suggests that ESMs can simulate the size of the variability in the ocean carbon sink. This conclusion also challenges earlier findings that GOBMs might underestimate decadal trends in the ocean carbon sink (DeVries et al., 2019). Other studies (Gloege et al., 2021; Hauck et al., 2023) support the hypothesis of an overestimation of decadal trends in the ocean carbon sink by pCO₂ products. Hauck et al. (2023) recently demonstrated with one GOBM and two pCO₂ products that sampling biases of pCO₂ have caused trends in the ocean carbon sink to be overestimated. A similar finding was made previously, when data from one GOBM, which was sampled in the same way as the real-world ocean was sampled, was extrapolated by one pCO₂ product to reconstruct the ocean carbon sink (Gloege et al., 2021). This reconstructed ocean carbon sink by the pCO₂ products also had a larger variability than the directly simulated ocean carbon sink by the GOBM (Gloege et al., 2021). Thus, it appears that most of the pCO₂ products overestimated and still overestimate decadal trends in the global ocean carbon sink. Therefore, estimates of the variability and size of the flux of natural carbon based on the difference in the total air–sea CO₂ flux from pCO₂ products and the change in interior ocean anthropogenic carbon (Müller et al., 2023; Gruber et al., 2019a), defined in this special case only as the additional carbon from increasing atmospheric CO₂ and not from climate change, are also likely too large.

The dependence found here of the ocean carbon sink on atmospheric CO₂ also has implications for studies that extrapolate present-day observation-based estimates of the ocean carbon sink back in time over the entire historical period to estimate a cumulative ocean carbon sink since the beginning of the Industrial Revolution using the difference in atmospheric CO₂ since pre-industrial times as a scaling factor (Gruber et al., 2009; Mikaloff Fletcher et al., 2006). While this scaling more or less works for most of the historical period, it breaks down from 1920 to 1960 and in the 1990s (Fig. 1). In addition, such estimates might be highly sensitive to the year for which the ocean carbon sink was estimated based on observations. If that year falls in one of these anomalous periods, as the year 1995 in Mikaloff Fletcher et al. (2006) does, the scaling might be biased low or high. Therefore, these extrapolations of present-day fluxes over the historical period should be used with caution or with a slightly more complex extrapolation method that takes the change in the atmospheric CO₂ growth rate into account.

The importance of changes in the atmospheric CO₂ growth rate for the trends in the global ocean carbon sink also affects our understanding of the uncertainty about the ocean carbon sink and the role of internal variability in the future. Previous studies used CMIP simulations with prescribed atmospheric CO₂ to quantify the importance of internal variability for the uncertainty about the projections of the ocean carbon sink in comparison to the importance of model and scenario uncertainty (Gooya et al., 2023; Lovenduski et al., 2016; Schlunegger et al., 2020). As these prescribed atmospheric CO₂ time series in CMIP simulations are much smoother than observed atmospheric CO₂ time series (Fig. 1), changes in the atmospheric CO₂ growth rate are also less variable. Thus, these concentration-driven CMIP SSPs suppress the internal variability in the atmospheric CO₂ growth rate caused by variabilities in atmospheric temperature, precipitation, El Niño, and volcanic eruptions (Keeling et al., 1995; Kuo et al., 1990; Raupach et al., 2008; Zeng et al., 2005; Bacastow, 1976; Yang and Wang, 2000). The suppressed variability in the atmospheric CO₂ growth rate in concentration-driven SSPs also suppresses the variability in the ocean carbon sink in the future, leading to an underestimation of the importance of the internal variability in ESMs for the overall uncertainty about ocean carbon sink projections over the 21st century (Gooya et al., 2023; Lovenduski et al., 2016; Schlunegger et al., 2020). This underestimation of the variability in the ocean carbon sink due to prescribed atmospheric CO₂ can be avoided if ESMs are run in an emission-driven mode that automatically introduces a strong variability in the atmospheric CO₂ growth rate, as in the model intercomparison project using the Adaptive Emission Reduction Approach (Terhaar et al., 2022a; Silvy et al., 2024).

While changes in the atmospheric CO₂ growth rate and ocean heat uptake might allow the estimation of changes in the decadal variability in the ocean carbon sink, it remains unknown how climate variability and individual modes can be used to predict inter-annual variability in the near-term ocean carbon sink (Lovenduski et al., 2019). In addition, other external forcings, such as volcanic eruptions, are an important factor in the inter-annual variability in the ocean carbon sink but also contribute to decadal trends (McKinley et al., 2017; Fay et al., 2023; Frölicher et al., 2011, 2013).

The influence of changes in the atmospheric CO₂ growth rate on the ocean carbon sink also has profound implications on the near-term future of the ocean carbon sink. With less strongly increasing or even peaking carbon emissions, the atmospheric CO₂ growth rate will also peak and potentially decline. The growth rate of atmospheric CO₂ on Mauna Loa has shown a robust negative trend since 2016, and the last time that the growth rate of atmospheric CO₂ was as small as in 2022 was the year 2008 (https://gml.noaa.gov/ccgg/trends/gr.html, last access: 10 June 2024). If this change from a rise in the atmospheric CO₂ growth rate towards a decline in the atmospheric CO₂ growth rate continues, the forced trend in the ocean carbon sink will also be negative. In addition, GOBMs and pCO₂ products suggest that the internal climate variability led to particularly positive trends in the ocean carbon sink in the 2000s and 2010s (Figs. 3 and 5; orange and pink dots lie above the derived relationship). This climate variability will eventually reverse at some point and lead to a larger decline in the ocean carbon sink. Furthermore, ocean heat uptake is projected to increase over the next decade or 2, independent of the chosen future pathway, also leading to stronger outgassing of carbon from the ocean in the near future and negative trends in the ocean carbon sink (Fig. 6). Although CMIP6 ESMs have tended to simulate a smaller ocean heat content over the last 2 decades (Lyu et al., 2021), this might not be an overestimation of the heat uptake by the ESMs but an especially low uptake due to climate variability. Indeed, the recent strong increase in ocean heat content and sea surface temperatures in 2023 (Cheng et al., 2024) might be the beginning of a shift from a period of low heat uptake due to climate variability to a period of high ocean heat uptake. Together, the decreasing atmospheric CO₂ growth rate, the potential change in internal climate variability, and the increased ocean heat uptake will likely cause a substantial negative trend in the ocean carbon sink over the next decade. If, however, emissions and atmospheric CO₂ rise, the continuous increase in the atmospheric CO₂ growth rate will cause the ocean carbon sink to increase as well.

Overall, this study demonstrates how ESMs can be used to better understand the past and future of the ocean carbon sink and drivers of its variability. They hence provide not only a valuable addition to pCO₂ products and GOBMs, but also a unique tool to statistically assess uncertainties and drivers of variability, also potentially in the interior ocean (Müller et al., 2023). The robustness in these results is further corroborated by their capability to simulate the size of decadal trends in important climate modes.