# Climate-driven variability of the Southern Ocean CO 2 sink

Nicolas Mayot, Corinne Corinne Le Quéré, Christian Rödenbeck, Raffaele Bernardello, Laurent Bopp, Laique M. Djeutchouang, Marion Gehlen, Luke Gregor, Nicolas Gruber, Judith Hauck, et al.

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

To cite this version:

Nicolas Mayot, Corinne Corinne Le Quéré, Christian Rödenbeck, Raffaele Bernardello, Laurent Bopp, et al. Climate-driven variability of the Southern Ocean CO 2 sink. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2023, 381 (2249), 10.1098/rsta.2022.0055. hal-04106303

HAL Id: hal-04106303

https://hal.science/hal-04106303

Submitted on 27 May 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Research

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

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

Check for
updates

**Cite this article:** Mayot N *et al*. 2023
Climate-driven variability of the Southern
Ocean CO₂ sink. *Phil. Trans. R. Soc. A 381*:
20220055.

https://doi.org/10.1098/rsta.2022.0055

Received: 24 October 2022

Accepted: 3 April 2023

One contribution of 13 to a discussion meeting
issue 'Heat and carbon uptake in the Southern
Ocean: the state of the art and future
priorities'.

**Subject Areas:**

oceanography, biogeochemistry

**Keywords:**

Southern Ocean, carbon sink, climate, oxygen,
interannual, decadal

**Author for correspondence:**

C. Le Quéré

e-mail: c.lequere@uea.ac.uk

Climate-driven variability of
the Southern Ocean CO₂ sink

N. Mayot<fnref n="1" />, C. Le Quéré<fnref n="1" />,⁺, C. Rödenbeck<fnref n="2" />,

R. Bernardello<fnref n="3" />, L. Bopp<fnref n="4" />, L. M. Djeutchouang<fnref n="5" />,<fnref n="6" />,

M. Gehlen<fnref n="7" />, L. Gregor<fnref n="8" />, N. Gruber<fnref n="8" />, J. Hauck<fnref n="9" />,

Y. Iida¹⁰, T. Ilyina¹¹, R. F. Keeling¹², P. Landschützer¹¹,¹³,

A. C. Manning<fnref n="1" />, L. Patara¹⁴, L. Resplandy¹⁵,

J. Schwinger¹⁶,¹⁷, R. Séférian¹⁸, A. J. Watson¹⁹,

R. M. Wright<fnref n="1" /> and J. Zeng²⁰

<footnote>†Present address: School of Environmental Sciences, University of East Anglia, Norwich NR4 7TJ, UK.</footnote>

<footnote n="1">¹Centre for Ocean and Atmospheric Sciences, School of Environmental Sciences, University of East Anglia, Norwich NR4 7TJ, UK</footnote>

<footnote n="2">²Max Planck Institute for Biogeochemistry, PO Box 600164,
Hans-Knöll-Str. 10, 07745 Jena, Germany</footnote>

<footnote n="3">³Department of Earth Sciences, Barcelona Supercomputing Center,
Barcelona, Catalonia, Spain</footnote>

<footnote n="4">⁴Laboratoire de Météorologie Dynamique/Institut Pierre-Simon
Laplace, CNRS, Ecole Normale Supérieure/Université PSL, Sorbonne
Université, Ecole Polytechnique, Paris, France</footnote>

<footnote n="5">⁵Department of Oceanography, University of Cape Town, Cape Town
7701, South Africa</footnote>

<footnote n="6">⁶SOCCO, Council for Scientific and Industrial Research, Cape Town
7700, South Africa</footnote>

<footnote n="7">⁷Laboratoire des Sciences du Climat et de l'Environnement,
LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, 91191
Gif-sur-Yvette, France</footnote>

<footnote n="8">⁸Environmental Physics, ETH Zürich, Institute of Biogeochemistry
and Pollutant Dynamics and Center for Climate Systems Modeling
(C2SM), Zurich, Switzerland</footnote>

<footnote n="9">⁹Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und
Meeresforschung, Postfach 120161, 27515 Bremerhaven, Germany</footnote>
¹⁰Atmosphere and Ocean Department, Japan Meteorological Agency, 1-3-4 Otemachi, Chiyoda-Ku, Tokyo 100-8122, Japan

¹¹Max Planck Institute for Meteorology, Hamburg, Germany

¹²Scripps Institution of Oceanography, University of California, San Diego, CA, USA

¹³Flanders Marine Institute (VLIZ), Jacobsenstraat 1, 8400 Ostend, Belgium

¹⁴GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany

¹⁵Department of Geosciences and High Meadows Environmental Institute, Princeton University, Princeton, NJ, USA

¹⁶Bjerknes Centre for Climate Research, Bergen, Norway

¹⁷NORCE Norwegian Research Centre, Jahnebakken 5, 5007 Bergen, Norway

¹⁸CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France

¹⁹College of Life and Environmental Sciences, University of Exeter, Exeter EX4 4RJ, UK

²⁰Earth System Division, National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba, Ibaraki 305-8506, Japan

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

NM, 0000-0003-4255-5939; CR, 0000-0001-6011-6249; LB, 0000-0003-4732-4953; JH, 0000-0003-4723-9652; PL, 0000-0002-7398-3293; ACM, 0000-0001-6952-7773; AJW, 0000-0002-9654-8147

The Southern Ocean is a major sink of atmospheric CO₂, but the nature and magnitude of its variability remains uncertain and debated. Estimates based on observations suggest substantial variability that is not reproduced by process-based ocean models, with increasingly divergent estimates over the past decade. We examine potential constraints on the nature and magnitude of climate-driven variability of the Southern Ocean CO₂ sink from observation-based air-sea O₂ fluxes. On interannual time scales, the variability in the air-sea fluxes of CO₂ and O₂ estimated from observations is consistent across the two species and positively correlated with the variability simulated by ocean models. Our analysis suggests that variations in ocean ventilation related to the Southern Annular Mode are responsible for this interannual variability. On decadal time scales, the existence of significant variability in the air-sea CO₂ flux estimated from observations also tends to be supported by observation-based estimates of O₂ flux variability. However, the large decadal variability in air-sea CO₂ flux is absent from ocean models. Our analysis suggests that issues in representing the balance between the thermal and non-thermal components of the CO₂ sink and/or insufficient variability in mode water formation might contribute to the lack of decadal variability in the current generation of ocean models.

This article is part of a discussion meeting issue 'Heat and carbon uptake in the Southern Ocean: the state of the art and future priorities'.

# 1. Introduction

The Southern Ocean CO₂ sink represents about 40% of the global oceanic CO₂ sink. Large decadal variations have been evidenced from observations [1–6], which are not related directly to changes in emissions of CO₂ from human activities but rather to variable climate conditions and/or variable external forcings and their influence on the growth rate of atmospheric CO₂ [7–9]. It is essential to better characterize and simulate correctly the variability of the Southern Ocean CO₂ sink in process-based ocean models to improve our understanding of the global carbon-cycle and its future evolution [9–13].

Two types of approaches widely used to estimate the oceanic CO₂ sink and its spatio-temporal variability include ‘pCO₂ products’, which are based on observations compiled in the Surface Ocean CO₂ Atlas (SOCAT) [14], and Global Ocean Biogeochemistry Models (GOBMs), which
are based on simulating the carbon cycle and its response to observed climate variability and changes in atmospheric CO₂ [15]. Results from both methods suggest a relative stagnation of the Southern Ocean CO₂ sink in the 1990s [2] and its reinvigoration in the 2000s [3]. However, they strongly disagree on the magnitude of these temporal variations, which are below 0.04 PgC yr⁻¹ in GOBMs and around 0.08–0.18 PgC yr⁻¹ in pCO₂ products [11,15]. While pCO₂ products might overestimate decadal variability due to sparse and unevenly distributed data [16,17], independent constraints from atmospheric CO₂ inversions tend to support larger variability compared with GOBMs [15].

Here, we make use of a temporal decomposition methodology and of observation-based air–sea O₂ fluxes to gain further insights on the variability of the Southern Ocean CO₂ sink. First, we isolate the climate-driven variability of the Southern Ocean CO₂ sink, i.e. the part of the CO₂ sink caused only by fluctuations in climate, and decompose it into its short-term interannual component (i.e. year-to-year) and its longer-term decadal/sub-decadal component. This temporal decomposition aims to better identify the issues and underlying processes [16,18,19]. Second, we compare the climate-driven variability of air–sea fluxes of CO₂ and O₂, from both observational products and GOBMs, recognizing that CO₂ and O₂ are affected by the same processes, but in different proportions. Indeed, both O₂ and CO₂ in the ocean are influenced by thermal processes (e.g. warming of the surface ocean) in similar ways, while they are both also influenced by non-thermal processes (e.g. biological photosynthesis and respiration, and ocean circulation) largely in opposite ways [20,21].

The overall objective of this study is thus to provide potential constraints on the nature and magnitude of the climate-driven variability of the Southern Ocean CO₂ sink by examining the coherence between CO₂ and O₂ flux variability estimated by data products and ocean models. For this, we evaluate the ability of 10 GOBMs to simulate interannual and decadal variability in air–sea fluxes of CO₂ and O₂ inferred from observations, and examine the overall coherence between CO₂ and O₂ variability. We hypothesize that if observation-based air–sea fluxes of CO₂ and O₂ from completely independent methods suggest similar interannual and/or decadal variabilities, then these are true signals of climate-driven variability of the Southern Ocean CO₂ sink that GOBMs should simulate.

## 2. Data and methods

This study focuses on the period 1985 to 2018 and is based on monthly gridded data of air–sea fluxes of CO₂ and O₂ in the Southern Ocean from GOBMs and observation-based products. The Southern Ocean is defined here as the ocean area south of 30° S. The total oceanic CO₂ sink (Total Flux) can be described as:

$$ \text{Total Flux} = \text{Flux}_{\text{ant}}^{\text{ss}} + \text{Flux}_{\text{ant}}^{\text{ns}} + \text{Flux}_{\text{nat}}^{\text{ss}} + \text{Flux}_{\text{nat}}^{\text{ns}} \quad (2.1) $$

where Flux<sub>ant</sub> and Flux<sub>nat</sub> are the air–sea fluxes of anthropogenic and natural CO₂, respectively. The superscript ss (steady state) denotes fluxes under unchanging climate conditions (on time scales longer than a year), whereas ns (non-steady state) denotes fluxes that are solely affected by changing climate conditions. Therefore, Flux<sub>ant</sub><sup>ss</sup> captures the effect of rising atmospheric CO₂ alone on the ocean CO₂ sink, Flux<sub>nat</sub><sup>ss</sup> captures the flux of natural CO₂ in a constant climate and Flux<sub>nat</sub><sup>ns</sup> + Flux<sub>ant</sub><sup>ns</sup> captures the climate-driven variability of the ocean CO₂ sink.

### (a) Air–sea CO₂ fluxes

Each of the 10 GOBMs used here comprises an ocean physical model coupled with a marine biochemistry module (table 1). Models are forced with observed atmospheric CO₂ mole fraction, and winds and other weather conditions from atmospheric reanalysis datasets (called ‘atmospheric forcing’). The 10 GOBMs differ through their use of different ocean physical models, representation of biogeochemistry, forcing products, spin-up strategies and spatial resolutions [15], all of which influence model representation of CO₂ and O₂ fluxes.
<table><caption>**Table 1.** List of the GOBMs and pCO₂ products used in this study with some of their characteristics.</caption><thead><tr><th colspan="4">GOBMs</th></tr><tr><th>name</th><th>physical ocean model</th><th>biogeochemistry model</th><th>atmospheric forcing</th></tr></thead><tbody><tr><td>CESM-ETHZ</td><td>CESMv1.3</td><td>BEC</td><td>JRA55</td></tr><tr><td>CNRM</td><td>NEMOv3.6</td><td>PISCESv2</td><td>JRA55</td></tr><tr><td>EC-Earth3-CC</td><td>NEMOv3.6</td><td>PISCESv2</td><td>JRA55</td></tr><tr><td>FESOM-RECOm</td><td>FESOM-1.4</td><td>RECOm-2</td><td>JRA55</td></tr><tr><td>MPIOM-HAMOCC</td><td>MPIOM</td><td>HAMOCC6</td><td>NCEP</td></tr><tr><td>IPSL</td><td>NEMOv3.6</td><td>PISCESv2</td><td>JRA55</td></tr><tr><td>NorESM1-OCv1.2</td><td>MICOM</td><td>HAMOCC</td><td>NCEP</td></tr><tr><td>MOM6-Princeton</td><td>MOM6-SIS2</td><td>COBALTv2</td><td>JRA55</td></tr><tr><td>NEMO-PlankTOM12</td><td>NEMOv3.6</td><td>PlankTOM12</td><td>NCEP</td></tr><tr><td>ORCA025-GEOMAR</td><td>NEMO-ORCA025</td><td>MOPS</td><td>JRA55</td></tr></tbody><thead><tr><th colspan="4">pCO₂ products</th></tr><tr><th>name</th><th>gas exchange parameterization</th><th>wind product</th><th>atmospheric CO₂ fields</th></tr></thead><tbody><tr><td>CMEMS-LSCE-FFNN</td><td>quadratic exchange formulation [22]</td><td>ERA 5</td><td>[23]</td></tr><tr><td>CSIR-ML6</td><td>quadratic exchange formulation [22]</td><td>ERA 5</td><td>NOAA</td></tr><tr><td>Jena-MLS</td><td>quadratic exchange formulation [22]</td><td>JMA55-do reanalysis</td><td>Jena CarboScope</td></tr><tr><td>JMA-MLR</td><td>quadratic exchange formulation [22]</td><td>JRA55</td><td>JMA-GSAM</td></tr><tr><td>MPI-SOMFFN</td><td>quadratic exchange formulation [22]</td><td>ERA 5</td><td>NOAA</td></tr><tr><td>NIES-NN</td><td>quadratic exchange formulation [22]</td><td>ERA 5</td><td>NOAA</td></tr><tr><td>OS-ETHZ-GRaCER</td><td>quadratic exchange formulation [24]</td><td>JRA55, ERAS, NCEP1</td><td>NOAA</td></tr><tr><td>Watson2020</td><td>Nightingale formulation [25]</td><td>CCMP</td><td>NOAA</td></tr></tbody></table>

All 10 GOBMs performed standardized simulations as part of the RECCAP2 project (https://reccap2-ocean.github.io/), following a common protocol (electronic supplementary material). Models conducted a simulation 'A' designed to capture the Total Flux term (equation (2.1)). For this, models were all forced with similar observed atmospheric CO₂ mole fraction, and with variable weather conditions (electronic supplementary material, table S1). Models performed another simulation 'C' to capture the $Flux_{ant}^{ss} + Flux_{nat}^{ss}$ terms, where models were again forced with observed atmospheric CO₂ mole fraction, but this time with weather conditions reflecting a climatological year (e.g. looping over the same year). The climate-driven variability $Flux_{nat}^{ns} + Flux_{ant}^{ns}$ was obtained by subtracting simulation C from simulation A.

All eight pCO₂ products used here were part of the RECCAP2 project. Each product estimates the oceanic CO₂ sink and its variability based on the observations of sea surface fugacity of CO₂ (fCO₂) from the SOCAT database [14]. First, fCO₂ observations are interpolated and extrapolated in time and space using statistical or machine learning methods. Second, the air-sea CO₂ fluxes are calculated by subtracting the corresponding atmospheric CO₂ mole fraction from the monthly (or daily) gridded ocean fCO₂ estimates, and multiplying the difference by a gas-exchange coefficient, which is a function of wind speed. The eight pCO₂ products differ through their use of different methods to produce the gridded maps of fCO₂, and in the variety of gas-exchange parametrizations and ancillary datasets required in the calculation of air-sea CO₂ fluxes (table 1). The pCO₂ products estimate the 'Total Flux' component of equation (2.1). In order to isolate the climate-driven variability component $Flux_{nat}^{ns} + Flux_{ant}^{ns}$, the model estimate of the air-sea CO₂
fluxes driven by atmpsheric CO₂ alone as captured by the multi-model average of simulation C was subtracted from the pCO₂ products estimates.

For the remaining of the manuscript, air-sea CO₂ fluxes refer to the climate-driven variability of the CO₂ fluxes for simplicity, unless specified otherwise.

## (b) Air-sea O₂ fluxes

All GOBMs provided monthly gridded air-sea O₂ fluxes from their simulation 'A' except ORCA025-GEOMAR. Observation-based estimates are taken from an atmospheric inversion method that optimized air-sea fluxes of Atmospheric Potential Oxygen (APO) to track the observed changes in APO concentration [26]. APO combines atmospheric observations of O₂ and CO₂ (APO = O₂ + 1.1 CO₂) in a manner that is conservative with regard to land biosphere exchanges (O₂:CO₂ = 1.1 mol mol⁻¹) [27]. The combustion of fossil fuels also influences APO. This contribution has been removed as part of the inversion method using fuel-specific O₂:CO₂ ratios (globally weighted average O₂:CO₂ ≈ 1.4 mol mol⁻¹). Therefore, APO data adjusted for the relatively well-known fossil fuel combustion are records of air-sea fluxes of O₂ and CO₂. It has been demonstrated that variability in APO air-sea flux is approximately equal to the variability in air-sea O₂ fluxes [26,28]. Here, we assume that the variability in APO air-sea flux estimated by an atmospheric inversion provides an observation-based estimate of the variability in O₂ air-sea fluxes.

The O₂ air-sea fluxes are not influenced by a strong anthropogenic signal (i.e. Fluxant) as is the case for CO₂ [26]. Following the principles of equation (2.1), the climate-driven variability in O₂ air-sea fluxes is captured by the term Fluxnat. Therefore, this observation-based estimate can be used as an independent constraint to evaluate the ability of GOBMs to correctly simulate the oceanic processes that influence air-sea exchanges of O₂, which impact the climate-driven variability of the CO₂ sink as well.

Here, we use the atmospheric CarboScope APO inversion (https://www.bgc-jena.mpg.de/CarboScope), which is based on the TM3 atmospheric tracer transport model. The inversion optimizes APO observations at five stations (period 1994–2018) or nine stations (period 1999–2018). Observations are from the Scripps O₂ program (https://scrippso2.ucsd.edu/). The temporal variations in air-sea O₂ flux estimated with the APO inversion method depend mainly on the inversion configuration, the quality of the APO data and the number and location of stations [26]. As the inversion is based on atmospheric observations, it is totally independent of the pCO₂ products, which are based on oceanic observations. Note that the spatial patterns of the O₂ air-sea fluxes cannot be studied with the atmospheric inversion because of an insufficient number of sampling stations in the Southern Hemisphere to obtain robust estimates of finer-scale spatial features [26].

## (c) Time series decomposition

The interannual and decadal/sub-decadal (hereafter named decadal) components of the climate-driven CO₂ and O₂ air-sea fluxes time series were isolated using the following signal decomposition methodology: (i) the long-term mean was removed to focus on the temporal variability, (ii) the seasonal cycle was removed by applying a 12-month moving average, (iii) the decadal component of this de-seasonalized time series was obtained by filtering this time series with a 48-month Hanning window, (iv) the interannual component was extracted by removing the decadal component from the original de-seasonalized time series, and (v) the interannual component was smoothed with a 5-month Hanning window to eliminate the small month-to-month variability. The Hanning window is a filtering function with a 'bell-shaped' curve used to smooth the signal by emphasizing the feature near the centre of the window. A 48-month wide window eliminates most year-to-year variability, and therefore isolates the decadal component. The significance of the correlation coefficients between two filtered time series takes into account the corresponding degree of freedom by considering the e-folding decay time of autocorrelation
<figure><img src="image_6.png" /></figure>

**Figure 1.** Southern Ocean CO₂ sink from 1985 to 2019. Positive values denote a sink for CO₂. The GOBM mean estimate (solid black line) is the average of 10 GOBMs with ±1 standard deviation of the model ensemble (grey shading). The dashed line represents the effect of increasing atmospheric CO₂ only estimated from the GOBM mean. The pCO₂ product mean estimate (blue line) is the average of eight pCO₂ products with ±1 standard deviation of the product ensemble (blue shading). For the pCO₂ product estimates, a river flux adjustment term of 31.6 Tmol yr⁻¹ was added to be comparable with the GOBM estimates (see reference [15] for more details). The seasonal cycle was removed from all estimates with a 12-month moving average (see Methods). (Online version in colour.)

[29]. The standard deviation of each time series is used as a measure of the magnitude of the variability.

We apply similarly steps (i–v) above to time series of the Southern Annular Mode (SAM) index [30] and to the atmospheric forcing time series of wind speed and Sea Surface Temperature (SST) from the National Centers for Environmental Prediction (NCEP), after removing their long-term trends using a linear fit. However, for the SAM index, step (ii) was not applied because there is no clear seasonal cycle and step (v) used a longer 18-month Hanning window because the SAM index signal is noisier than the climate-driven CO₂ and O₂ air–sea fluxes time series.

# 3. Results

## (a) Overview of the recent changes in the Southern Ocean CO₂ sink

On average over the year, the Southern Ocean is a sink for CO₂. The CO₂ uptake mostly occurs in the Subtropical Zone (between 30°S and the Subtropical Front) and the Subantarctic Zone (between the Subtropical and the Subantarctic Front). This is mainly driven by natural CO₂ uptake due to the cooling of subtropical waters, which are transported southwards, and the anthropogenic CO₂ uptake by recently upwelled waters with low anthropogenic CO₂ concentrations, which are transported northwards [1,31]. In the Subpolar Zone (north of the seasonally ice-covered zone), there is a net outgassing of CO₂ due to the upwelling of deep waters with a high concentration of dissolved inorganic carbon [1,31].

On average, the Southern Ocean CO₂ sink increased by 30.6 Tmol yr⁻¹ between the first and last decades (1985–1994 and 2009–2018, respectively) according to the pCO₂ products (figure 1). GOBMs estimated a slightly lower increase of 24.5 Tmol yr⁻¹, almost entirely caused by the response to the increasing atmospheric CO₂ mole fraction (26.5 Tmol yr⁻¹, dashed line in figure 1). Climate variability induced fluctuations in the Southern Ocean CO₂ sink of approximately 6.9 Tmol yr⁻¹ according to the pCO₂ products and 3.3 Tmol yr⁻¹ according to GOBMs, significantly larger than the fluctuations of about 1.9 Tmol yr⁻¹ also estimated by
<figure><img src="image_7.png" /></figure>

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

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

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

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

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

**Figure 2.** Comparison of the climate-driven air–sea CO₂ (left) and O₂ (right) fluxes estimated by GOBMs and observation-based products. (a,b) Climate-driven air–sea CO₂ and O₂ fluxes from GOBMs (black) and observation-based products (blue and magenta, see legend), decomposed into their (c,d) decadal and (e,f) interannual components, showing the mean across estimates and the ±1 standard deviation across the ensembles (shading). Fluxes are defined as positive from the atmosphere into the ocean. (Online version in colour.)

GOBMs and caused by variability in atmospheric CO₂ mole fraction only (estimated as the standard deviations of the detrended time series).

## (b) Simulated versus observed temporal variability of CO₂ and O₂ air–sea fluxes

Temporal variations in the CO₂ sink from the means of the GOBMs and pCO₂ products are clearly different, both in phase and in magnitude (figure 2a). Most of the discrepancies are related to the amplitude of the decadal variability (figure 2c) which is uncorrelated and three times lower in GOBMs (2.0 Tmol yr⁻¹) than in pCO₂ products (6.3 Tmol yr⁻¹). None of the 10 GOBMs simulates decadal variability that significantly correlates with that of the pCO₂ product mean, and all GOBMs underestimate the magnitude of the pCO₂ product mean by 24% to 77%. Moreover, all pCO₂ products except one have a higher decadal amplitude than any other GOBM (table 2). However, only half of the individual pCO₂ products are significantly and positively correlated
with the mean pCO₂ product signal, although the degree of freedom to test the significance of the correlation is small and some moderate or high positive correlation values (r ≥ 0.66) are considered non-significant. Longer time series are needed to increase the degree of freedom.

The decadal variability in O₂ fluxes is also lower in the GOBM mean (8.5 Tmol yr⁻¹) than in the atmospheric inversion (22.2 Tmol yr⁻¹), but they are significantly correlated (r = 0.63, figure 2d). However, the correlation between the GOBM mean and the inversion is not improved when using the version of the atmospheric inversion with nine stations. Only three of the nine GOBMs have a significant correlation with the atmospheric inversion using five stations (0.62 ≥ r ≥ 0.74), and two with the inversion using nine stations (0.74 ≥ r ≥ 0.77). All GOBMs underestimate the magnitude compared with the inversion estimate by 25% to 74%.

The interannual variabilities in CO₂ and O₂ air-sea fluxes are more similar between the GOBM mean and the observation-based estimates (figure 2e,f). For CO₂, the interannual variability is significantly correlated (r = 0.64) and similar in magnitude (2.1 Tmol yr⁻¹ for GOBMs and 1.7 Tmol yr⁻¹ for pCO₂ products). The median values of the interannual amplitude from all individual GOBMs and from all individual pCO₂ products are similar (p-value = 0.0634; Wilcoxon Rank Sum Test). Eight GOBMs have a positive correlation with the pCO₂ product mean (0.44 ≥ r ≥ 0.63). All individual pCO₂ products are positively correlated with the pCO₂ product mean (0.55 ≥ r ≥ 0.91; table 2). For O₂, the interannual variability is also significant correlated (r = 0.46 and r = 0.55 for the inversions with the five and nine stations, respectively). Six of the nine GOBMs are significantly correlated with the atmospheric inversion (0.36 ≥ r ≥ 0.56). However, the magnitude of the mean GOBM interannual variability (8.3 Tmol yr⁻¹) is three times lower than that estimated by the atmospheric inversion (26.1 Tmol yr⁻¹). Within the observation-based estimates, the winter season is more correlated with the interannual variations (r = 0.97 for CO₂, and r = 0.96 for O₂) than the summer season (r = 0.87 for CO₂, and r = 0.71 for O₂).

### (c) Simulated versus observed spatial variability of CO₂ and O₂ air-sea fluxes

For the decadal component of the CO₂ fluxes, there are clear differences between GOBMs and pCO₂ products (figure 3a,b). First, the values are in general two times lower in GOBMs than in the pCO₂ products. Second, in GOBMs, most of the highest magnitudes are south of or along the Subantarctic Front, apart from coastal regions south of Australia. In pCO₂ products, areas south of the Subantarctic Front are also of importance, but they extend well to the north of this front and in all basins. For the decadal component of O₂ air-sea fluxes in GOBMs, the areas south of the Subantarctic Front are of importance according to GOBMs (figure 3c).

For the interannual component of the CO₂ fluxes, the areas of highest variability occur south of the Subantarctic Front in both GOBMs and pCO₂ products (figure 3d,e). GOBMs highlight the importance of the Indian and Pacific sectors, a distinction less clear in the pCO₂ products. For the interannual component of O₂ fluxes (figure 3f), the area south of the Subantarctic Front is also of importance according to GOBMs, with also a larger influence of the Indian and Pacific Oceans, as for the decadal variability (figures 3c,f).

### (d) Relationship between observed variability of CO₂ and O₂ air-sea fluxes

The observed decadal variability in air-sea CO₂ flux from the pCO₂ product mean is significantly correlated with the decadal variability in air-sea O₂ fluxes inferred with the APO atmospheric inversion based on nine stations only (r = -0.81, figure 4a and table 3), but not with the five-station inversion. This significant correlation is mostly induced by three individual pCO₂ products that have correlation values r ≤ -0.72 (table 2).

The observed interannual variability in air-sea CO₂ flux from the pCO₂ product mean is significantly correlated with the interannual variability in air-sea O₂ flux estimated with the APO atmospheric inversion with both five stations (r = -0.40) and nine stations (r = -0.64) (figure 4b and table 3). The significant correlation with the five-station inversion is supported by three out
<table><caption>**Table 2.** Correlation (r) and amplitude of the decadal and short-term interannual variability in CO₂ and O₂ air–sea fluxes. Estimates inferred from the pCO₂ product mean and the APO inversion with five stations (the first two lines), are compared with the fluxes inferred from each individual GOBM, the GOBM mean, each individual pCO₂ product and the APO inversion with nine stations. The correlation coefficients are calculated with a Pearson correlation and take into account the degree of freedom for each time scale. Values in bold are statistically significant.</caption><thead><tr><th rowspan="2"></th><th colspan="2">CO₂ (Tmol yr<sup>-1</sup>)</th><th colspan="2">interannual</th><th colspan="2">O₂ (Tmol yr<sup>-1</sup>)</th><th rowspan="2">interannual<br>r (APO inversion<br>(5 stn.) versus)</th><th rowspan="2">amplitude</th></tr><tr><th>decadal<br>r (pCO₂ product<br>mean versus)</th><th>amplitude</th><th>r (pCO₂ product<br>mean versus)</th><th>amplitude</th><th>decadal<br>r (APO inversion<br>(5 stn.) versus)</th><th>amplitude</th></tr></thead><tbody><tr><td>pCO₂ product mean</td><td>—</td><td>6.31</td><td>—</td><td>1.68</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>APO inversion (5 stn.)</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>22.2</td><td>—</td><td>26.1</td></tr><tr><td colspan="9"><strong>GOBMs</strong></td></tr><tr><td>CESM-ETHZ</td><td>0.08</td><td>2.8</td><td><strong>0.44</strong></td><td>2.4</td><td><strong>0.72</strong></td><td>8.5</td><td>0.28</td><td>9.9</td></tr><tr><td>CNRM-electronic<br>supplementary<br>material</td><td>-0.12</td><td>4.8</td><td><strong>0.5</strong></td><td>3.4</td><td>0.53</td><td>14.9</td><td><strong>0.36</strong></td><td>9.4</td></tr><tr><td>EC-Earth3</td><td>0.35</td><td>2.9</td><td><strong>0.62</strong></td><td>2.7</td><td>0.45</td><td>11.9</td><td><strong>0.41</strong></td><td>9.8</td></tr><tr><td>FESOM_RECOM</td><td>-0.02</td><td>1.9</td><td><strong>0.63</strong></td><td>3.0</td><td>0.55</td><td>5.8</td><td><strong>0.42</strong></td><td>9.5</td></tr><tr><td>MPIOM-HAMOCC</td><td>0.23</td><td>4.1</td><td><strong>0.58</strong></td><td>5.6</td><td>-0.07</td><td>15.4</td><td><strong>0.39</strong></td><td>16.0</td></tr><tr><td>NEMO-PISCES</td><td>0.08</td><td>3.1</td><td><strong>0.48</strong></td><td>2.3</td><td>0.55</td><td>10.6</td><td><strong>0.43</strong></td><td>9.6</td></tr><tr><td>NorFSM</td><td>0.07</td><td>1.5</td><td>0.07</td><td>1.7</td><td>0.5</td><td>9.1</td><td>0.30</td><td>8.9</td></tr><tr><td>MOM6-Princeton</td><td>-0.05</td><td>3.6</td><td><strong>0.54</strong></td><td>2.7</td><td><strong>0.74</strong></td><td>8.2</td><td><strong>0.56</strong></td><td>10.0</td></tr><tr><td>NEMO-PlankTOM12</td><td>0.2</td><td>2.4</td><td>-0.32</td><td>2.1</td><td><strong>0.62</strong></td><td>16.6</td><td>0.20</td><td>11.3</td></tr><tr><td>GEOMAR</td><td>0.06</td><td>2.9</td><td><strong>0.57</strong></td><td>3.1</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>GOBM mean</td><td>0.12</td><td>2.1</td><td><strong>0.64</strong></td><td>2.1</td><td><strong>0.63</strong></td><td>8.5</td><td><strong>0.46</strong></td><td>8.3</td></tr></tbody><tfoot><tr><td colspan="9">(Continued)</td></tr></tfoot></table>
<table><caption>Table 2. (Continued.)</caption><thead><tr><th rowspan="2"></th><th colspan="2">CO<sub>2</sub> (Tmol yr<sup>-1</sup>)</th><th colspan="2">interannual</th><th colspan="2">O<sub>2</sub> (Tmol yr<sup>-1</sup>)</th><th rowspan="2">interannual<br>r (APO inversion<br>(5 stn.) versus)</th><th rowspan="2">amplitude</th></tr><tr><th>decadal<br>r (pCO<sub>2</sub> product<br>mean versus)</th><th>amplitude</th><th>r (pCO<sub>2</sub> product<br>mean versus)</th><th>amplitude</th><th>decadal<br>r (APO inversion<br>(5 stn.) versus)</th><th>amplitude</th></tr></thead><tbody><tr><td colspan="9"><strong>pCO<sub>2</sub> products</strong></td></tr><tr><td>CMEMS-LSCEFFNN</td><td>0.70</td><td>4.9</td><td><strong>0.72</strong></td><td>2.0</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>CSIR-ML6</td><td><strong>0.97</strong></td><td>6.7</td><td><strong>0.91</strong></td><td>2.1</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Jena-MLS</td><td>0.46</td><td>8.0</td><td><strong>0.78</strong></td><td>4.9</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>JMA-MLR</td><td>0.66</td><td>5.8</td><td><strong>0.61</strong></td><td>1.6</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>MPI-SOMFFN</td><td><strong>0.95</strong></td><td>13.7</td><td><strong>0.55</strong></td><td>2.3</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>NIES-NN</td><td>0.63</td><td>2.7</td><td><strong>0.68</strong></td><td>1.2</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>OS-ETHZ-GRaCER</td><td><strong>0.94</strong></td><td>5.8</td><td><strong>0.79</strong></td><td>2.0</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Watson2020</td><td><strong>0.96</strong></td><td>13.5</td><td><strong>0.66</strong></td><td>2.6</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td colspan="9"><strong>APO inversion</strong></td></tr><tr><td>APO inversion (9 stn.)</td><td>—</td><td>—</td><td>—</td><td>—</td><td><strong>0.79</strong></td><td>24.7</td><td><strong>0.68</strong></td><td>28.1</td></tr></tbody></table>
<figure><img src="image_13.png" /></figure>

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

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

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

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

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

**Figure 3.** Standard deviation of the decadal (left) and interannual (right) variability of CO₂ and O₂ air–sea fluxes. For CO₂ (a,b and d,e), the averaged maps from GOBMs and from pCO₂ products are shown, while for O₂ (c and f) only GOBMs could be used to study the spatial patterns. Note the differences in units and colour scales between O₂ and CO₂ maps. The black lines represent the average location of the Subantarctic Front (northern line) and of the September extent of sea ice (southern line). Maps from individual GOBM and pCO₂ product are available in the electronic supplementary material (figures S1–S6). (Online version in colour.)

of the eight pCO₂ products, while all pCO₂ products except one support the significant correlation with the nine-station inversion (table 3).

Significant correlations between CO₂ and O₂ fluxes are all negative, suggesting that non-thermal processes, which influence CO₂ and O₂ in opposite direction, are the dominant source of variability on average over the Southern Ocean.
<figure><img src="image_19.png" /><figcaption>**Figure 4.** Comparison of the variability in CO₂ and O₂ air–sea fluxes. The (a) decadal and (b) interannual components of CO₂ and O₂ air–sea fluxes are from the observation-based products (i.e. pCO₂ products and atmospheric inversion with five or nine stations). The y-axes on the right for O₂ air–sea fluxes are inverted. Fluxes are defined as positive from the atmosphere into the ocean. (Online version in colour.)</figcaption></figure>

### (e) Relationships with the SAM index, wind speed and SST

On decadal time scale, significant correlations are found between the SAM and the CO₂ and O₂ air–sea fluxes from GOBMs (r = -0.52 for CO₂ and r = 0.73 for O₂), but not between SAM and the observation-based air–sea fluxes (i.e. pCO₂ products and atmospheric inversion, figure 5d,f). Therefore, the simulated influence of the SAM on the decadal variability could be either spurious or due to missing or poorly represented processes in most GOBMs.

On interannual time scale, significant correlations are found between the SAM and the CO₂ and O₂ air–sea fluxes (figure 5e,g) in both models and observation-based fluxes (for GOBM fluxes, r = -0.71 for CO₂ and r = 0.70 for O₂; for observation-based fluxes, r = -0.56 for CO₂ and r = 0.44 for O₂). These results tend to support that the observed negative correlation between the interannual component of CO₂ and O₂ air–sea fluxes (see previous section) is also related to the SAM index.

Negative correlations between CO₂ and O₂ fluxes on interannual time scale are associated with wind speed and occur in the Subpolar Zone (that extents from the Subantarctic Front to the September extent of sea ice; 15% of sea ice concentration, figure 6). In detail, in this region, correlations are mostly negative between CO₂ fluxes and wind speed, and mostly positive between O₂ fluxes and wind speed. Finally, and still in the Subpolar Zone, correlations are generally positive between CO₂ fluxes and SST, and negative between O₂ fluxes and SST. This suggests that interannual variations in the Southern Ocean CO₂ sink are induced by processes
<table><caption>**Table 3.** Correlation between observation-based variability of CO₂ and O₂ air–sea fluxes inferred from the pCO₂ product mean and the APO inversion, respectively (with five or nine stations). The analysis was done with the decadal and interannual components (see also figure 4). The correlation coefficients are calculated with a Pearson correlation and take into account the degree of freedom for each time scale. Values in bold are statistically significant.</caption><thead><tr><th rowspan="2"></th><th colspan="2">decadal</th><th colspan="2">interannual</th></tr><tr><th>APO inversion (5 stn.)</th><th>APO inversion (9 stn.)</th><th>APO inversion (5 stn.)</th><th>APO inversion (9 stn.)</th></tr></thead><tbody><tr><td><strong>pCO<sub>2</sub> products</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>CMEMS-LSCEFFNN</td><td>-0.03</td><td>-0.07</td><td>-0.1</td><td><strong>-0.37</strong></td></tr><tr><td>CSIR-ML6</td><td>-0.45</td><td>-0.69</td><td><strong>-0.41</strong></td><td><strong>-0.64</strong></td></tr><tr><td>Jena-MLS</td><td>-0.41</td><td>-0.62</td><td><strong>-0.38</strong></td><td><strong>-0.46</strong></td></tr><tr><td>JMA-MLR</td><td>-0.26</td><td>-0.08</td><td>-0.17</td><td>-0.34</td></tr><tr><td>MPI-SOMFFN</td><td>-0.58</td><td><strong>-0.72</strong></td><td>-0.25</td><td><strong>-0.55</strong></td></tr><tr><td>NIES-NN</td><td>-0.23</td><td>-0.37</td><td>-0.26</td><td><strong>-0.43</strong></td></tr><tr><td>OS-ETHZ-GRaCER</td><td>-0.31</td><td><strong>-0.76</strong></td><td><strong>-0.5</strong></td><td><strong>-0.62</strong></td></tr><tr><td>Watson2020</td><td>-0.41</td><td><strong>-0.75</strong></td><td>-0.25</td><td><strong>-0.43</strong></td></tr><tr><td>pCO<sub>2</sub> product mean</td><td>-0.47</td><td><strong>-0.81</strong></td><td><strong>-0.4</strong></td><td><strong>-0.64</strong></td></tr></tbody></table>

occurring in the Subpolar Zone, with stronger CO₂ outgassing events related to stronger winds, colder SST and stronger O₂ ingassing events (and vice versa).

# 4. Discussion

The temporal variability of the Southern Ocean CO₂ sink driven by climate variations can be decomposed into two components: a short-term interannual component and a decadal component (figure 2).

## (a) Interannual variability of the Southern Ocean CO₂ sink

Our analysis suggests that pCO₂ products accurately represent the climate-driven interannual variations of the Southern Ocean CO₂ sink. This is supported by the strong similarities between the interannual variability of CO₂ and O₂ fluxes provided by two completely independent observation-based products (figure 4b and table 3). Our analysis further suggests that this interannual variability is predominately regulated by wintertime deep-water ventilation south of the Subantarctic Front. This is supported (i) by the negative sign of the correlation between CO₂ and O₂ air–sea fluxes, which indicates the dominance of non-thermal processes [18], (ii) by the highest variability rates south of the Subantarctic Front, which corresponds to the location of ocean ventilation and outgassing of CO₂ (ingassing of O₂) [31,32], and (iii) by the higher correlation during the winter season pointing at a physical rather than a biological process. Therefore, the climatic processes regulating the wintertime ventilation of the Southern Ocean might also control the interannual variations of the Southern Ocean CO₂ sink [33].

Our results indicate that the SAM could exert a strong control on the interannual variations of Southern Ocean CO₂ sink (figure 5e,g). A positive (negative) SAM is related to an intensification (weakening) of the winds, which drive the upwelling intensity and carbon storage through the Ekman transport and other processes regulating the Mixed Layer Depth [2,34–36]. Years with stronger CO₂ outgassing events in the Subpolar Zone are related to years experiencing stronger wind, colder SST (because of the upwelling of cold deep waters) and stronger O₂ ingassing events (figure 6). Colder SST reinforces the ingassing of O₂ in the Subpolar Zone (i.e. the thermal and non-thermal components reinforce each other), but dampens the outgassing of CO₂.
<figure><img src="image_20.png" /></figure>

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

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

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

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

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

<figure><img src="image_26.png" /><figcaption>**Figure 5.** Comparison of the SAM index and the climate-driven variability of CO₂ and O₂ air–sea fluxes from the GOBM mean, pCO₂ product mean and APO inversion. (a) The detrended time series of the SAM index, which has been decomposed into (b) a decadal component and (c) an interannual component, and compared with the decadal and interannual components from CO₂ (d and e, the right-hand axes are inverted) and O₂ (f and g) air–sea fluxes. Fluxes are defined as positive from the atmosphere into the ocean. (Online version in colour.)</figcaption></figure>

Our results show that GOBMs simulate the interannual CO₂ variability well, both in phase and in magnitude (figure 5e), including the SAM control of this variability. The correlation between the SAM index and the simulated and observed CO₂ sink are similar for the interannual components (r = -0.64 and r = -0.51, respectively). This may be because the SAM is an atmospheric mode of variability and GOBMs are forced with atmospheric reanalysis dataset for wind forcing and surface heating whose variations are partly induced by the SAM.

The SAM is known to develop non-zonal atmospheric forcing in the Southern Ocean [37,38], inducing regional variations of the Mixed Layer Depth that influence the CO₂ and O₂ air–sea fluxes [18,19,32,39], but such regional features are currently missing in pCO₂ products (figure 3b), possibly because of insufficient observational coverage (electronic supplementary material, figure S7). GOBMs simulated some regional variabilities, with a strong influence of the Indo-Pacific sectors on the interannual component of the CO₂ sink (figure 3a). Recently, the Indo-Pacific sectors were pointed out as the main regions for CO₂ outgassing when considering new indirect estimates of air–sea CO₂ fluxes derived from Biogeochemical-Argo float observations [33,40]. New observations maintained over several years will be needed to confirm these regional influences [33].

Our findings are slightly conflicting with previous results that argued for the influence of the summer season on the short-term interannual variability of the CO₂ sink [19]. However, this later study [19] only compares the year-to-year variability of monthly pCO₂ during the
<figure><img src="image_27.png" /><figcaption>**Figure 6.** Correlation maps between the interannual components of the CO₂ and O₂ air–sea fluxes, and with the interannual components of wind speed and SST. The black lines represent the averaged location of the Subantarctic Front and of the September extent of sea ice. Fluxes are defined as positive from the atmosphere into the ocean. (Online version in colour.)</figcaption></figure>

months associated with the annual maximum and minimum of pCO₂ and not the flux itself. In winter, small changes of pCO₂ can be amplified by the strong wintertime winds. Our suggestion that wintertime deep water ventilation events have a key role also implies that variability in gas transfer velocity has a second-order effect on interannual variations of O₂ air–sea fluxes, contrary to [41] who estimated that it was the most important process controlling the variability of O₂ air–sea fluxes. However, their sensitivity analysis was done with model outputs instead of performing separate model simulations with non-varying piston velocity or deep-water ventilation. Nonetheless, they highlighted that changes in the ventilation of O₂-depleted deep water strongly influenced the temporal variations of O₂ air–sea fluxes south of the Subantarctic Front. Moreover, they also found a relationship between air–sea O₂ fluxes and the SAM index and argued that it is driven by the influence of the SAM index on the upwelling rate of O₂-depleted deep waters.

Although the processes controlling the interannual variability of CO₂ and O₂ air–sea fluxes are emerging from our analysis, uncertainties remain on their magnitude. Whereas the variability in the CO₂ sink is comparable between GOBMs and pCO₂ products, GOBMs underestimate the interannual variability in O₂ air–sea fluxes by a factor of two to three, as pointed out elsewhere [26,41–43]. Further work is needed to examine the plausible causes of this discrepancy.

## (b) Decadal variability of the Southern Ocean CO₂ sink

The decadal component of the air–sea O₂ fluxes (figure 4a) tends to support the existence of a decadal climate-driven variability of the Southern Ocean CO₂ sink. Longer time series of atmospheric APO will be needed before this can be firmly confirmed, but the existence
of significant large decadal variability is consistent with several studies based on different methods (e.g. [5], see review in [1,9]). Decadal variations in the Southern Ocean CO₂ sink are mostly marked by a saturation period in the 1990s [2] and a reinvigoration period in the 2000s [3]. According to the pCO₂ products, the climate-driven decadal variations are three times larger than the interannual variations (figure 2).

Although our analysis does not fully resolve the magnitude of the decadal variability, it does suggest it is real and larger than GOBMs simulate, even though it may not be as large as estimated by pCO₂ products due to uneven sampling [16,17]. The relative performance of the GOBMs in reproducing O₂ decadal variability is better than that of CO₂ decadal variability, suggesting that the representation of the balance between thermal and non-thermal processes might be partly at fault in models. The balance between thermal and non-thermal components of the CO₂ sink has been shown to play an important role in driving seasonal [1,44] and interannual [6] variability. In contrast to CO₂ fluxes, for which the thermal and non-thermal components oppose each other, the thermal and non-thermal components involved in the O₂ air-sea fluxes reinforce each other [26,45], which might explain why some GOBMs were able to reproduce the observed decadal variability in O₂ air-sea fluxes and less the decadal changes in the CO₂ sink.

The spatial extent of the decadal variability in the air-sea CO₂ flux north of the Subantarctic Front suggests that Subantarctic Mode Water formation might also influence the decadal variability in the Southern Ocean CO₂ sink. Areas just north of the Subantarctic Front are where the upper cell of the Southern Ocean overturning circulation contributes to the subduction of mode waters, which could have substantial influence on the decadal variability of the CO₂ sink. The subduction of surface water into intermediate layers was already mentioned to explain some of the regionally enhanced CO₂ sink in the eastern Pacific between 2012 and 2016 [6]. Such subduction phenomena occur at specific locations in the Southern Ocean [46,47], which could explain the asymmetrical spatial pattern observed in pCO₂ products. Subantarctic Mode Water is particularly important for the transport and recirculation of absorbed anthropogenic CO₂ [48], but less so on its temporal variability, which is lower than the variability associated with natural CO₂ [49]. However, variations in the uptake flux of anthropogenic CO₂ could enhance the climate-driven variability of the natural CO₂ flux [1]. For instance, an enhanced upwelling rate south of the Subantarctic Front could either slightly increase the CO₂ sink north of this front by increasing the subsequent subduction of surface water that absorbed anthropogenic CO₂, or slightly decrease this CO₂ sink by shortening the residence time of the surface water [46,49,50].

The current generation of GOBMs have coarse resolutions and most likely do not correctly simulate the formation of Subantarctic Mode Water [51], which is sensitive to oceanic currents, wind speed, Mixed Layer Depth, sea ice and eddies [46,48,52]. Other studies highlight the importance of eddies after a positive SAM event that compensate for the intensified upwelling [53,54]. Therefore, improvements of the physical and ice ocean models might be needed to correctly simulate the decadal variability of the CO₂ sink, as well as the future evolution in the Southern Ocean [55,56].

Other possible sources of variability missing in GOBMs include the coupling with atmospheric dynamics, internal tracer anomalies (memory) and ecosystem variability. Atmospheric coupling is unlikely to be driving large decadal variability of the Southern Ocean CO₂ sink since Earth System Models also do not generate a ratio between the decadal and interannual variability similar to the one suggested by the pCO₂ products ensemble mean [18]. Internal anomalies in dissolved inorganic carbon and/or O₂ concentration could in theory trigger variability in air-sea fluxes when those anomalies become in contact with the atmosphere, but more work would be needed to verify if such anomalies in the ocean interior exist and how they relate to the patterns of variability identified here. Finally, changes in marine ecosystems in response to variability in ocean properties could act to enhance or dampen the thermal and/or non-thermal components and therefore amplify the total signals [5,34]. Current generation GOBMs represent ecosystems that are largely driven by upwelled nutrients and do not yet include the more complex ecosystem responses, such as vertical migrations and salps/krill dipole that characterize the Southern Ocean
(e.g. [57]). The importance of these processes for the variability in CO₂ and O₂ fluxes has not yet been examined.

## 5. Conclusion

The degree of concordance between observed and simulated variability in air-sea fluxes of CO₂ and O₂, at different time scales, was used to gain insights on climate-driven changes in the Southern Ocean CO₂ sink. The interannual variations of the Southern Ocean CO₂ sink derived from pCO₂ products and models are consistent with the interannual variations in air-sea O₂ flux derived from observations. The current generation of GOBMs can simulate the influence of stronger (weaker) winds during years of positive (negative) SAM that induce, in the Subpolar Zone, stronger (weaker) upwelling of deep waters and drive the short-term interannual variations of the Southern Ocean CO₂ sink. The decadal variations of the Southern Ocean CO₂ sink, suggested by several pCO₂ products, tends to be supported by the observed decadal variations of the air-sea O₂ flux. However, GOBMs do not reproduce these decadal CO₂ variations. Although the climate-driven processes associated with these decadal variations remain unclear, pCO₂ products suggest an influence from regions associated with the formation of Subantarctic Mode Water, a physical process that might be poorly represented in GOBMs, while the relative performance of models in reproducing the decadal variability of O₂ compared with CO₂ suggests issues in representing the balance between thermal and non-thermal processes. More in-situ pCO₂ data are required to confirm the influence of different Southern Ocean regions while more atmospheric APO data could help constrain the size of the decadal variability.

**Data accessibility.** All MATLAB scripts and data discussed in the results are publicly available. This GitHub repository contains instructions on how to access them: https://github.com/nmayot/PTA_SouthernOcean. The data are provided in the electronic supplementary material [58].

**Authors' contributions.** N.M.: conceptualization, formal analysis, investigation, methodology, visualization, writing—original draft; C.L.Q.: conceptualization, funding acquisition, investigation, methodology, project administration, supervision, writing—review and editing; C.R.: data curation, writing—review and editing; R.B.: data curation, writing—review and editing; L.B.: data curation, writing—review and editing; L.M.D.: data curation, writing—review and editing; M.G.: data curation, writing—review and editing; L.G.: data curation, writing—review and editing; N.G.: data curation, writing—review and editing; J.H.: data curation, writing—review and editing; Y.I.: data curation, writing—review and editing; T.I.: data curation, writing—review and editing; R.F.K.: data curation, writing—review and editing; P.L.: data curation, writing—review and editing; A.C.M.: data curation, writing—review and editing; L.P.: data curation, writing—review and editing; L.R.: data curation, writing—review and editing; J.S.: data curation, writing—review and editing; R.S.: data curation, writing—review and editing; A.W.: data curation, writing—review and editing; R.M.W.: data curation, writing—review and editing; J.Z.: data curation, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** N.M., C.L.Q., R.B., N.G., L.G. and A.C.M. acknowledge the funding from the European Commission through the H2020 project 4C (grant no. 821003). J.H., L.G., M.G. and N.G. acknowledge the funding from the European Commission through the H2020 project COMFORT project (grant no. 820989). C.L.Q. was funded by the UK Royal Society (grant no. RP\R1\191063). N.M. and R.M.W. were funded by UK's Natural Environment Research Council (SONATA: grant no. NE/P021417/1). J.S. received funding from the Research Council of Norway through project INES (grant no. 270061) and HPC resources provided by the National Infrastructure for HPC and Data Storage in Norway, UNINETT Sigma2 (grant no. nn/ns2980k). L.R. acknowledges the Princeton University Catalysis Initiative. N.G. and L.G. acknowledge funding from ETH Zürich. M.G. and R.S. acknowledge the ESM2025 project under the grant agreement number 101003536. M.G. also acknowledges funding from the European Union's Horizon 2020 Blue Growth research and innovation programme under grant agreement number 862923 (project AtlantECO). J.H. acknowledges support by the Initiative and Networking Fund of the Helmholtz Association (Helmholtz Young Investigator Group Marine Carbon and Ecosystem Feedbacks in the Earth System [MarESys], grant number VH-NG-1301). The integration of the ORCA025-GEOMAR experiment was performed at the North German Supercomputing Alliance (HLRN) and was financially supported by the German Research Foundation (project PA 3075/2-1).
The APO measurements were supported by a series of grants to the Scripps Institution of Oceanography from the US NSF and NOAA, most recently OPP-1922922 and NA20OAR4320278.

**Acknowledgements.** Thanks to Erik Buitenhuis and David Willis for their work on the NEMO-PlankTOM12 model development. The research presented in this paper was carried out on the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia.

## References

1. Gruber N, Landschützer P, Lovenduski NS. 2019 The variable Southern Ocean carbon sink. *Annu. Rev. Mar. Sci.* **11**, 159–186. (doi:10.1146/annurev-marine-121916-063407)

2. Le Quéré C et al. 2007 Saturation of the Southern Ocean CO₂ sink due to recent climate change. *Science* **316**, 1735–1738. (doi:10.1126/science.1136188)

3. Landschützer P, Gruber N, Bakker DCE. 2016 Decadal variations and trends of the global ocean carbon sink. *Glob. Biogeochem. Cycles* **30**, 1396–1417. (doi:10.1002/2015GB005359)

4. Rödenbeck C et al. 2015 Data-based estimates of the ocean carbon sink variability – first results of the Surface Ocean pCO₂ mapping intercomparison (SOCOM). *Biogeosciences* **12**, 7251–7278. (doi:10.5194/bg-12-7251-2015)

5. DeVries T et al. 2019 Decadal trends in the ocean carbon sink. *Proc. Natl Acad. Sci. USA* **116**, 11 646–11 651. (doi:10.1073/pnas.1900371116)

6. Keppler L, Landschützer P. 2019 Regional wind variability modulates the Southern Ocean carbon sink. *Sci. Rep.* **9**, 7384. (doi:10.1038/s41598-019-43826-y)

7. Canadell JG et al. 2021 Global Carbon and Other Biogeochemical Cycles and Feedbacks. In *Climate change 2021: The physical science basis contribution of working group I to the sixth assessment report of the intergovernmental panel on climate change* (eds V Masson-Delmotte et al.), pp. 673–816. Cambridge, United Kingdom and New York, NY: Cambridge University Press.

8. McKinley GA, Fay AR, Eddebbar YA, Gloege L, Lovenduski NS. 2020 External forcing explains recent decadal variability of the ocean carbon sink. *AGU Adv.* **1**, e2019AV000149. (doi:10.1029/2019AV000149)

9. Gruber N, Bakker DCE, DeVries T, Gregor L, Hauck J, Landschützer P, McKinley GA, Müller JD. 2023 Trends and variability in the ocean carbon sink. *Nat. Rev. Earth Environ.* **4**, 119–134. (doi:10.1038/s43017-022-00381-x)

10. DeVries T. 2022 The ocean carbon cycle. *Annu. Rev. Environ. Resour.* **47**, 317–341. (doi:10.1146/annurev-environ-120920-111307)

11. Hauck J et al. 2020 Consistency and challenges in the ocean carbon sink estimate for the global carbon budget. *Front. Mar. Sci.* **7**, 571720. (doi:10.3389/fmars.2020.571720)

12. Li H, Ilyina T. 2018 Current and future decadal trends in the oceanic carbon uptake are dominated by internal variability. *Geophys. Res. Lett.* **45**, 916–925. (doi:10.1002/2017GL075370)

13. McKinley GA, Fay AR, Lovenduski NS, Pilcher DJ. 2017 Natural variability and anthropogenic trends in the ocean carbon sink. *Annu. Rev. Mar. Sci.* **9**, 125–150. (doi:10.1146/annurev-marine-010816-060529)

14. Bakker DCE et al. 2022 Surface Ocean CO₂ Atlas Database Version 2022 (SOCATv2022) (NCEI Accession 0253659) [cited 2022 Aug 22]. (doi:10.25921/1h9f-nb73)

15. Friedlingstein P et al. 2022 Global carbon budget 2022. *Earth Syst. Sci. Data* **14**, 4811–4900. (doi:10.5194/essd-14-4811-2022)

16. Gloege L et al. 2021 Quantifying errors in observationally based estimates of ocean carbon sink variability. *Glob. Biogeochem. Cycles* **35**, e2020G.B006788. (doi:10.1029/2020GB006788)

17. Hauck J, Nissen C, Landschützer P, Rödenbeck C, Bushinsky SM, Olsen A. 2023 Sparse and unevenly distributed observations at high-latitudes dominate biases in air-sea CO₂ flux estimates. *Phil. Trans. R. Soc. A* **381**, 20220068. (doi:10.1098/rsta.2022.0068)

18. Resplandy L, Séférían R, Bopp L. 2015 Natural variability of CO₂ and O₂ fluxes: what can we learn from centuries-long climate models simulations? *J. Geophys. Res. Oceans* **120**, 384–404. (doi:10.1002/2014JC010463)

19. Gregor L, Kok S, Monteiro PMS. 2018 Interannual drivers of the seasonal cycle of CO₂ in the Southern Ocean. *Biogeosciences* **15**, 2361–2378. (doi:10.5194/bg-15-2361-2018)

20. Resplandy L, Keeling RF, Stephens BB, Bent JD, Jacobson A, Rödenbeck C, Khatiwala S. 2016 Constraints on oceanic meridional heat transport from combined measurements of oxygen and carbon. *Clim. Dyn.* **47**, 3335–3357. (doi:10.1007/s00382-016-3029-3)
21. Sarmiento JL, Gruber N. 2006 Ocean Biogeochemical Dynamics.

22. Wanninkhof R. 2014 Relationship between wind speed and gas exchange over the ocean revisited. *Limnol. Oceanogr. Methods* **12**, 351–362. (doi:10.4319/lom.2014.12.351)

23. Chevallier F. 2013 On the parallelization of atmospheric inversions of CO₂ surface fluxes within a variational framework. *Geosci. Model Dev.* **6**, 783–790. (doi:10.5194/gmd-6-783-2013)

24. Wanninkhof R. 1992 Relationship between wind speed and gas exchange over the ocean. *J. Geophys. Res. Oceans* **97**(C5), 7373–7382. (doi:10.1029/92JC00188)

25. Nightingale PD, Malin G, Law CS, Watson AJ, Liss PS, Liddicoat MI, Boutin J, Upstill-Goddard RC. 2000 In situ evaluation of air-sea gas exchange parameterizations using novel conservative and volatile tracers. *Glob. Biogeochem. Cycles* **14**, 373–387. (doi:10.1029/1999GB00091)

26. Rödenbeck C, Le Quéré C, Heimann M, Keeling RF. 2008 Interannual variability in oceanic biogeochemical processes inferred by inversion of atmospheric O₂/N₂ and CO₂ data. *Tellus B Chem. Phys. Meteorol.* **60**, 685–705. (doi:10.1111/j.1600-0889.2008.00375.x)

27. Stephens BB, Keeling RF, Heimann M, Six KD, Murnane R, Caldeira K. 1998 Testing global ocean carbon cycle models using measurements of atmospheric O₂ and CO₂ concentration. *Glob. Biogeochem. Cycles*. **12**, 213–230. (doi:10.1029/97GB03500)

28. Manning A, Keeling RF. 2006 Global oceanic and land biotic carbon sinks from the Scripps atmospheric oxygen flask sampling network. *Tellus B Chem. Phys. Meteorol.* **58**, 95–116. (doi:10.1111/j.1600-0889.2006.00175.x)

29. Leith CE. 1973 The standard error of time-average estimates of climatic means. *J. Appl. Meteorol. Climatol.* **12**, 1066–1069. (doi:10.1175/1520-0450(1973)012<1066:TSEOTA>2.0.CO;2)

30. Marshall GJ. 2003 Trends in the southern annular mode from observations and reanalyses. *J. Clim.* **16**, 4134–4143. (doi:10.1175/1520-0442(2003)016<4134:TITSAM>2.0.CO;2)

31. Gray AR, Johnson KS, Bushinsky SM, Riser SC, Russell JL, Talley LD, Wanninkhof R, Williams NL, Sarmiento JL. 2018 Autonomous biogeochemical floats detect significant carbon dioxide outgassing in the high-latitude southern ocean. *Geophys. Res. Lett.* **45**, 9049–9057. (doi:10.1029/2018GL078013)

32. Bushinsky SM, Gray AR, Johnson KS, Sarmiento JL. 2017 Oxygen in the southern ocean from argo floats: determination of processes driving air-sea fluxes. *J. Geophys. Res. Oceans* **122**, 8661–8682. (doi:10.1002/2017JC012923)

33. Prend CJ et al. 2022 Indo-pacific sector dominates southern ocean carbon outgassing. *Glob. Biogeochem. Cycles* **36**, e2021GB007226.

34. Hauck J, Völker C, Wang T, Hoppema M, Losch M, Wolf-Gladrow DA. 2013 Seasonally different carbon flux changes in the Southern Ocean in response to the southern annular mode. *Glob. Biogeochem. Cycles* **27**, 1236–1245. (doi:10.1002/2013GB004600)

35. Lenton A, Matear RJ. 2007 Role of the Southern Annular Mode (SAM) in Southern Ocean CO₂ uptake. *Glob. Biogeochem. Cycles* **21**, GB2016. (doi:10.1029/2006GB002714)

36. Lovenduski NS, Gruber N, Doney SC, Lima ID. 2007 Enhanced CO₂ outgassing in the Southern Ocean from a positive phase of the Southern Annular Mode. *Glob. Biogeochem. Cycles* **21**, GB2026. (doi:10.1029/2006GB002900)

37. Fogt RL, Marshall GJ. 2020 The southern annular mode: variability, trends, and climate impacts across the southern hemisphere. *WIREs Clim. Change* **11**, e652. (doi:10.1002/wcc.652)

38. Fogt RL, Jones JM, Renwick J. 2012 Seasonal zonal asymmetries in the southern annular mode and their impact on regional temperature anomalies. *J. Clim.* **25**, 6253–6270. (doi:10.1175/JCLI-D-11-00474.1)

39. Sallée JB, Speer KG, Rintoul SR. 2010 Zonally asymmetric response of the Southern Ocean mixed-layer depth to the Southern Annular Mode. *Nat. Geosci.* **3**, 273–279. (doi:10.1038/ngeo812)

40. Chen H, Haumann FA, Talley LD, Johnson KS, Sarmiento JL. 2022 The deep ocean's carbon exhaust. *Glob. Biogeochem. Cycles* **36**, e2021GB007156. (doi:10.1029/2021GB007156)

41. Nevison CD, Munro DR, Lovenduski NS, Keeling RF, Manizza M, Morgan EJ, Rödenbeck C. 2020 Southern annular mode influence on wintertime ventilation of the southern ocean detected in atmospheric O₂ and CO₂ measurements. *Geophys. Res. Lett.* **47**, e2019GL085667. (doi:10.1029/2019GL085667)

42. Andrews O, Buitenhuis E, Le Quéré C, Suntharalingam P. 2017 Biogeochemical modelling of dissolved oxygen in a changing ocean. *Phil. Trans. R. Soc. A* **375**, 20160328. (doi:10.1098/rsta.2016.0328)
43. McKinley GA, Follows MJ, Marshall J. 2000 Interannual variability of the air-sea flux of oxygen in the North Atlantic. *Geophys. Res. Lett.* **27**, 2933–2936. (doi:10.1029/2000GL011492)

44. Mongwe NP, Vichi M, Monteiro PMS. 2018 The seasonal cycle of pCO₂ and CO₂ fluxes in the Southern Ocean: diagnosing anomalies in CMIP5 Earth system models. *Biogeosciences* **15**, 2851–2872. (doi:10.5194/bg-15-2851-2018)

45. Gruber N, Gloor M, Fan SM, Sarmiento JL. 2001 Air-sea flux of oxygen estimated from bulk data: implications For the marine and atmospheric oxygen cycles. *Glob. Biogeochem. Cycles* **15**, 783–803. (doi:10.1029/2000GB001302)

46. Morrison AK, Waugh DW, Hogg A, Jones DC, Abernathey RP. 2022 Ventilation of the Southern Ocean Pycnocline. *Annu. Rev. Mar. Sci.* **14**, 405–430. (doi:10.1146/annurev-marine-010419-011012)

47. Sallée JB, Speer K, Rintoul S, Wijffels S. 2010 Southern ocean thermocline ventilation. *J. Phys. Oceanogr.* **40**, 509–529. (doi:10.1175/2009JPO4291.1)

48. Sallée JB, Matear RJ, Rintoul SR, Lenton A. 2012 Localized subduction of anthropogenic carbon dioxide in the Southern Hemisphere oceans. *Nat. Geosci.* **5**, 579–584. (doi:10.1038/ngeo1523)

49. DeVries T, Holzer M, Primeau F. 2017 Recent increase in oceanic carbon uptake driven by weaker upper-ocean overturning. *Nature* **542**, 215–218. (doi:10.1038/nature21068)

50. Gruber N et al. 2019 The oceanic sink for anthropogenic CO₂ from 1994 to 2007. *Science* **363**, 1193–1199. (doi:10.1126/science.aau5153)

51. Sallée JB, Shuckburgh E, Bruneau N, Meijers AJS, Bracegirdle TJ, Wang Z, Roy T. 2013 Assessment of Southern Ocean water mass circulation and characteristics in CMIP5 models: historical bias and forcing response. *J. Geophys. Res. Oceans* **118**, 1830–1844. (doi:10.1002/jgrc.20135)

52. Cerovečki I, Meijers AJS. 2021 Strong quasi-stationary wintertime atmospheric surface pressure anomalies drive a dipole pattern in the subantarctic mode water formation. *J. Clim.* **34**, 6989–7004. (doi:10.5194/egusphere-egu2020-12128)

53. Dufour CO, Sommer JL, Gehlen M, Orr JC, Molines JM, Simeon J, Barnier B. 2013 Eddy compensation and controls of the enhanced sea-to-air CO₂ flux during positive phases of the Southern Annular Mode. *Glob. Biogeochem. Cycles* **27**, 950–961. (doi:10.1002/gbc.20090)

54. Meredith MP, Hogg AM. 2006 Circumpolar response of Southern Ocean eddy activity to a change in the Southern Annular Mode. *Geophys. Res. Lett.* **33**, L16608. (doi:10.1029/2006GL026499)

55. Bourgeois T, Goris N, Schwinger J, Tjiputra JF. 2022 Stratification constrains future heat and carbon uptake in the Southern Ocean between 30°S and 55°S. *Nat. Commun.* **13**, 340. (doi:10.1038/s41467-022-27979-5)

56. Terhaar J, Frölicher TL, Joos F. 2021 Southern Ocean anthropogenic carbon sink constrained by sea surface salinity. *Sci. Adv.* **7**, eabd5964. (doi:10.1126/sciadv.abd5964)

57. Yang G, Atkinson A, Pakhomov EA, Hill SL, Racault MF. 2022 Massive circumpolar biomass of Southern Ocean zooplankton: implications for food web structure, carbon export, and marine spatial planning. *Limnol. Oceanogr.* **67**, 2516–2530. (doi:10.1002/lno.12219)

58. Mayot N et al. 2023 Climate-driven variability of the Southern Ocean CO₂ sink. Figshare. (doi:10.6084/m9.fishshare.c.6597304)