The ocean regulates the uptake, storage, and build up of carbon in the atmosphere on annual to millennial timescales, helping mitigate climate change by significantly modulating the long-term trend in the atmospheric CO₂ growth rate. The ocean took up on average, 3.0 ± 0.4 PgC yr⁻¹ or 29% of the total anthropogenic CO₂ emissions for the decade beginning in 2011, and cumulatively has sequestered over 1/3 of fossil emissions since the industrial era (Le Quéré et al., 2018; Friedlingstein et al., 2022). Although terrestrial ecosystems play the major role in modulating interannual fluctuations in the atmospheric CO₂ growth rate (Keppel-Aleks et al., 2014), when taking into account land-use change emissions to the atmosphere, which counteract the residual terrestrial sink, the ocean is the only long-term sink for anthropogenic carbon (Ballantyne et al., 2012; Crisp et al., 2022).

Interannual variability (IAV) in the atmospheric CO₂ growth rate is superimposed on the long-term positive trend. The atmospheric CO₂ surface network shows a decadal-average growth rate of 2.41 ± 0.26 ppm/y (mean ± 2 std dev) on an annual basis due to the impacts of internal climate variability during 2003–2016 (Buchwitz et al., 2018). This IAV has been shown to reflect, primarily, variations in the rate of terrestrial CO₂ uptake due to climate stressors (Luo et al., 2022) and disturbance (Keppel-Aleks et al., 2014), but also embodies IAV in net ocean carbon exchange and fossil fuel emissions (Doney et al., 2009; Crisp et al., 2022). The global ocean carbon flux IAV is estimated at about 0.3 PgC/yr (equivalent to about 0.15 ppm CO₂ mixed into the global atmosphere), with the tropical Pacific Ocean which is modulated by El Niño-Southern Oscillation (ENSO), accounting for a large fraction of the variability (Rödenbeck et al., 2014; Bennington et al., 2022). However, the influence of key regions of ocean flux IAV on spatial patterns of atmospheric CO₂ IAV may be obscured unless local scale variations in atmospheric CO₂ are analyzed.

The net flux of CO₂ between the atmosphere and the ocean is a function of (1) wind-driven gas exchange kinetic rates and (2) the difference of the partial pressure of CO₂ in the air and the ocean surface (ΔpCO₂), with the ocean component being far more regionally heterogeneous (Wanninkhof and McGillis, 1999; Garbe et al., 2014). Local variations in oceanic pCO₂ are related to physical or biogeochemical processes including sub-surface water upwelling, which can enrich the surface water in CO₂, temperature-driven solubility variations, alkalinity (which is closely related to salinity and controls the speciation of dissolved inorganic carbon), and phytoplankton photosynthesis and respiration (Doney et al., 2009; Crisp et al., 2022). These variables are all affected by interannual climate variability (McKinley et al., 2020). Because the response of the global carbon cycle to climate fluctuations may provide insight into the long-term response to climate change, understanding the global and regional characteristics of ocean-driven atmosphere CO₂ IAV can help to improve our understanding of the climate–carbon cycle processes and our ability to project the fate of the ocean CO₂ sink in the future.

Although studies have suggested that ocean flux IAV may impart an observable impact on atmospheric CO₂ IAV (Crisp et al., 2022), gaps remain in quantifying the influence from the ocean since most atmospheric CO₂ observations are made on land and coasts, with fewer island and ship-based measurements in the remote open ocean. Early research deduced that the CO₂ flux variation over the ocean, especially the equatorial Pacific Ocean, is one of the main causes of the atmospheric CO₂ IAV (Francey et al., 1995). Later studies based on inverse models, seawater system measurements, and air-sea CO₂ flux estimates (Feely et al., 2002; Rödenbeck et al., 2003) suggested air-sea CO₂ flux is not the primary driver for interannual to seasonal variations in atmospheric CO₂. Nevison et al. (2008) used an atmospheric transport model with an underlying mechanistic ocean flux model to show that the amplitude of atmospheric IAV owing to ocean fluxes was around 10% of the IAV amplitude at northern hemisphere surface stations, and up to 50% of the observed IAV in the Southern Hemisphere. In neither hemisphere however, was the IAV owing to ocean fluxes highly correlated with the observations (Nevison et al., 2008). Setting aside uncertainties within the atmosphere, there remain substantial uncertainties as to the magnitude of air-sea CO₂ flux variability itself based on ocean biogeochemical models (Hauck et al., 2020; Fay and McKinley, 2021), observation-based flux products (McKinley et al., 2020; Bennington et al., 2022; Hauck et al., 2023) and atmospheric CO₂ inverse models (Peylin et al., 2013).

Quantifying ocean-driven atmospheric CO₂ IAV remains challenging since direct observations of pCO₂ in surface waters and the corresponding difference with the atmosphere are sparse, and leading to large spatiotemporal gaps. CO₂ fluxes can be estimated/calculated? from observation-based pCO₂ products. There are many such products available, each of which estimate near global 1°x1°, monthly pCO₂ fields from the sparse pCO₂ measurements available, using various techniques including statistical interpolation, linear and non-linear regressions, and machine learning-based methodologies (Rödenbeck et al., 2015). For all products, the CO₂ flux is then calculated from the interpolated pCO₂ fields (Fay et al., 2021). Here, we focus on only three representative data-based flux products (Landschützer et al., 2014; Rödenbeck et al., 2014, 2016, Denvil-Sommer et al., 2019) which are described in more detail below. Robust estimates of mean and seasonality can be derived from observation-based products (Fay and McKinley, 2021; Gloege et al., 2021), while significant uncertainties in terms of higher and lower frequency variability remain (Rödenbeck et al., 2015; Hauck et al., 2020; Bennington et al., 2022).

New space-based observations of the column-integrated CO₂ mole fraction, XCO₂, from NASA’s Orbiting Carbon Observatory 2 (OCO-2) mission may provide a unique vantage to observe the atmospheric imprint of IAV in air-sea CO₂ fluxes. OCO-2 is NASA’s first dedicated Earth remote sensing satellite to study atmospheric carbon dioxide from space (Eldering et al., 2017). It was designed to collect space-based measurements of atmospheric CO₂ with high precision and near-global coverage (Crisp et al., 2012, 2017). Compared to surface in-situ CO₂ observations, such as from the NOAA greenhouse gas network, which is mostly sited in coastal, island, and inland locations, OCO-2 can observe directly over the open ocean, which may improve the spatiotemporal attribution of ocean fluxes.

The research presented here leverages these two advances in carbon cycle data products to answer the following scientific questions: (1) What are the fingerprints of ocean carbon fluxes on IAV of atmospheric CO₂, and from what regions are these imprints most prominent? (2) What differences emerge from different observation-based products, and how large are these differences compared to the IAV owing from atmospheric transport? (3) Is the IAV signal from ocean fluxes detectible in column CO₂ observed from the state-of-the-art OCO-2 satellite? To answer these questions, we run atmospheric transport simulations with underlying air-sea fluxes from several observation-based products. We use the output from these simulations to quantify the imprint of regional and global sea-air CO₂ fluxes on atmospheric CO₂ IAV. As a final step, we use OCO-2 observations of XCO₂ to contextualize the simulated ocean-driven XCO₂ variations.

Figure 2 shows the global distribution of air-sea CO₂ flux IAV amplitude (calculated as the standard deviation of IAV anomalies) over 40 years from 1982 to 2021 in the observation-based products. The three observation-based products show broadly similar magnitudes and spatial patterns in air-sea flux IAV in the Northern Hemisphere and Southern Hemisphere low to mid latitudes ( Supplementary Figure S3 ). The flux IAV amplitude is highest over the northern midlatitude ocean, Eastern Pacific, and Southern Ocean, with the largest inter-product differences in the northern hemisphere. In the Southern Hemisphere, the zonally integrated air-sea flux IAV is largest at 60°, around 0.6-0.7 mol m^-2 y^-1. The flux IAV decreases by almost a factor of three in the subtropical regions in both hemispheres to around 0.2 mol m^-2 s^-1. The large zonal IAV and the small zonal standard deviation around a latitude circle in the Southern Ocean, especially for SOMFFN and CMEMS (shading in Figures 2D, F ), is due to the fact that in the Northern Hemisphere, though there is larger flux IAV amplitude (solid blue line in Figures 2D–F ), there is also more variability from gridcell to gridcell (shaded blue area in Figures 2D–F ). Thus, the gridcell-level variability partially cancels when integrated around a latitude band, leading to a smaller northern ocean impact on atmospheric XCO_2. Consequently, the Southern Hemisphere has a bigger imprint of atmospheric XCO₂ IAV, as we will show below.

The XCO₂ timeseries simulated from global air-sea fluxes (i.e., the sum of all tagged regions) shows IAV that ranges from +/-0.2 ppm for JENA ( Figure 3 ) with smaller variations from +/-0.1 ppm for SOM-FFN and CMEMS. The IAV from the products is small compared with the OCO-2 XCO₂ IAV (Guan et al., 2023). Both spatially averaged timeseries ( Figure 3 ) and maps of IAV ( Supplementary Figure S4 ) are notably different across the three observation-based flux products when propagated through an atmospheric transport model. The simulated XCO₂ from JENA is different compared to that derived from the other two observation-based flux products, due to differences especially in the Southern Hemisphere fluxes. We note, however, that XCO₂ simulated from any one flux product is generally coherent across latitude bands. The corresponding globally averaged XCO₂ IAV amplitude is about 0.11 ppm for JENA, but less than 0.07 ppm for SOMFFN and CMEMS ( Figures 4A–C ). Each data product shows the ocean-driven XCO₂ IAV amplitude largest over the tropical Pacific and the Southern Ocean. When the IAV amplitude is calculated for the period from September 2014 to 2021, the overlapping period between OCO-2 observations and observation-based products, the amplitude is similar for SOMFFN and JENA (around 0.12 ppm globally averaged), and similar for CMEMS (around 0.08 ppm). The difference in simulated IAV for a given observation-based product calculated across two different time periods suggests that the ocean may leave a larger imprint on XCO₂ IAV in recent years compared to historical averages.

The Southern Hemisphere oceanic regions are the dominant contributors to XCO₂ IAV across the three models. We calculated the XCO₂ IAV amplitude arising separately from each regional tracer, then determined the region that contributes the largest and second largest IAV amplitude to each global atmospheric gridcell (left-hand column of Figure 5 and Supplementary Figures S5–7 ). Across all three models, the Southern Hemisphere low- and mid-latitude regions are the primary contributors to IAV over 70-95% of the global area ( Table 1 ), including large swaths of the Northern Hemisphere for SOMFFN and CMEMS. In fact, the SOMFFN simulation shows that Northern Hemisphere ocean regions are not a dominant contributor to XCO₂ IAV even locally. In contrast, the XCO₂ IAV from the simulation with JENA and CMEMS fluxes show a modest contribution from the Northern Hemisphere, with the mid latitude region emerging as the dominant contributor to about 15% of the global area. The dominant region in each gridcell accounts for up to 60% of the atmospheric XCO₂ IAV amplitude for CMEMS ( Figure 5F ), and roughly 35-50% of the IAV amplitude for SOMFFN and JENA ( Figures 5B , D ), suggesting more equable contributions from different regions. The second most-dominant contributors to atmospheric CO₂ IAV are similarly the Southern Hemisphere, or tropics or Northern Hemisphere mid-latitude across all three ocean-flux-GEOS-Chem simulations ( Table 2 ).

Not all IAV in XCO₂ arises from IAV in local surface fluxes; IAV in patterns of atmospheric transport is also an important contributor to XCO₂ IAV. We calculate the relative contributions of ocean flux variability and IAV in atmospheric transport from GEOS-Chem simulations with climatological monthly-mean cyclostationary ocean CO₂ fluxes as surface boundary conditions. IAV in these simulations arises only from IAV in atmospheric transport acting on the spatial gradients in XCO₂ set up by the climatological annual cycle of ocean air-sea CO₂ exchange. All three products show nearly identical spatial patterns in the XCO₂ amplitude arising from cyclostationary air-sea CO₂ fluxes ( Figures 6A–C ), suggesting that seasonal flux patterns result in similar atmospheric gradients (Landschtzer et al., 2023) across the three models. The western Pacific warm pool is a hotspot for dynamics-driven variability, likely due to ENSO-driven changes in atmospheric transport, and this region has an IAV amplitude of 0.07-0.08 ppm across all three simulations. Given the large variability in the magnitude of IAV among SOMFFN, JENA, and CMEMS ( Figure 4 ), the ratio of dynamics-driven variations to total IAV diverged among the ocean-flux-GEOS-Chem simulations ( Figures 6D–F ). For SOMFFN and CMEMS, which showed smaller IAV when driven with interannually varying fluxes, transport contributes about 50% of IAV in the Southern Hemisphere subtropics and subpolar regions, and 50-80% of IAV in the tropics and high latitudes of both hemispheres. For JENA, which had high IAV in the full simulations ( Figure 4B ), transport contributes less than 20% of total IAV except in the tropical Pacific, where it contributes 50%.

In the atmospheric transport simulations with cyclostationary fluxes, the Southern Hemisphere mid-latitude region is the dominant contributor to IAV for about 60% to 96% of the global area. Consistent across the three observation-based products, the second dominant contributor is the Southern Hemisphere low-latitude region. These results suggest that IAV in transport interacts most strongly with atmospheric gradients set up by Southern Hemisphere fluxes. JENA and CMEMS show that the Northern Hemisphere mid latitude regions are also important contributors when run with cyclostationary fluxes, and are the dominant contributor to roughly 15% to 25% of the global area. Figure 7 shows that in the cyclostationary simulations, the influence of Northern and Southern Hemisphere regions are more localized to the originating hemisphere. Taken together with the results from the time-varying fluxes, these results suggest that IAV in the Southern Hemisphere fluxes is synergistic with IAV from atmospheric transport, whereas, in the Northern hemisphere, IAV in the fluxes counteracts the IAV due to transport alone.

The ocean-driven XCO₂ IAV simulated by the three observation-based products is largely unrelated to observed patterns of IAV as observed by OCO-2, which includes other sources of IAV. The OCO-2 IAV amplitude, which is of order 1 ppm for the 2014-2021 period, was an order of magnitude larger than that simulated from any of the three-ocean observation-based products which is of order 0.1ppm ( Figures 4 , 8 ). A possible explanation could be that the OCO-2 observations contain the imprint not only of ocean fluxes but also land and fossil emissions. The OCO-2 XCO₂ IAV amplitude, calculated as the standard deviation of the IAV timeseries, suggests that XCO₂ interannual variability over ocean basins is smaller than that over continents (around 0.4 ppm vs 1.2 ppm; Figure 8 ), although this may be due to larger error variance due to complex topography and land surface albedo variations (Guan et al., 2023; Mitchell et al., 2023). We calculate the correlation coefficient, slope, and fractional ratio between simulated XCO₂ IAV and OCO-2 XCO₂ IAV ( Figure 9 , Supplementary Figure S8 ) at the gridscale to explore where the imprint of the ocean might be detectible. Over most of the planet, the observed XCO₂ IAV is only weakly correlated with the XCO₂ IAV simulated from ocean fluxes, with many regions showing a slight negative correlation. These correlations are generally not statistically significant, with p>0.05. Simulations from the three different flux products show different meridional patterns of correlation, with JENA simulations showing slight positive correlations within the tropics and CMEMS showing slight positive correlations in the Southern Hemisphere, although we note these relationships are not statistically significant. Across all regions, the slope between the observed and simulated XCO₂ IAV is less than 0.2 units? ( Figures 9D–F ), consistent with the magnitude of the observed IAV being dominated by terrestrial signals. The maps of correlations and slopes suggest that detecting ocean IAV directly from space-based measurements is challenging given other sources of variability. Disentangling these multiple impacts to differentiate between plausible ocean flux products, such as the three analyzed here, will require a combination of observations and modeling.

Previous analysis of OCO-2 observations over the tropical Pacific suggested that space-based observations can detect an ocean-driven decrease in XCO₂, potentially as large as 0.5 ppm, during the initial phases of the strong 2015 El Niño (Chatterjee et al., 2017). Our simulations show that atmospheric XCO₂ decreases over the Nino 3.4 region during the strong 1997-1998 El Niño (by about 0.3 ppm), but show a weaker decrease (about 0.1ppm) during 2015-2016 El Niño. The decrease in the atmospheric XCO₂ IAV timeseries corresponds to the reduced ocean outgassing over the Niño 3.4 region (up to 1.5 mol/m²/y) ( Figure 10 ), which is roughly similar between the two El Niño periods and leads to similar decreases in the tropical XCO₂ tracer ( Supplementary Figure S9 ). During the strong 1997-98 El Niño, however, a temporary reduction in Northern Hemisphere CO₂ outgassing also reduced XCO₂ over the Niño.4 region, amplifying the apparent response ( Supplementary Figure S9 ). As expected, XCO₂ anomalies during the two El Niño events are mainly driven by changes in the air-sea CO₂ flux rather than changes in the atmospheric circulation ( Supplementary Figure S10 ). However, these flux-driven anomalies are small (less than +/- 0.03 ppm). Our simulations suggest that even during strong El Niño events, detecting the direct atmospheric imprint of changes in ocean outgassing within the tropics will be challenging given the small direct signal predicted from all three observational flux products, compounding variability from other ocean and terrestrial regions, and retrieval uncertainties in XCO₂.

We simulate ocean-driven IAV in atmospheric CO₂ based on three estimates of air-sea CO₂ fluxes from interpolated observation-based products. Although these observation-based products all use the same ocean pCO₂ data to estimate gridded fluxes and show largely similar spatial patterns of flux IAV, they result in very different estimates of atmospheric XCO₂ IAV over the 40 year period from 1982-2021 ( Figures 3 , 4A–C ). The inter-model spread was much reduced when looking at the four years from 2014-2021 overlapping with OCO-2 observations ( Figures 4D–F ). Our results largely corroborate the results from Nevison et al. (2008), who showed a small imprint for the ocean on atmospheric CO₂ IAV. Our results show that even in the remote Southern Ocean, which shows large IAV in the fluxes, the atmospheric IAV signature imparted by the ocean is small.

Further, our simulations show that up to 80 percent of the IAV of the IAV in atmospheric XCO₂ results from IAV in atmospheric transport acting on atmospheric gradients derived from cyclostationary ocean fluxes ( Figure 6 ), not from IAV in the ocean fluxes themselves. The three observation-based products showed a very similar pattern of transport-induced IAV, with the dominant spatial contributors (SH low-lat and SH mid-lat) showing high consistency among the three data products ( Table 1 ). This is somewhat expected, since the MERRA/GEOS-Chem atmospheric transport patterns were common among all three simulations, but also requires that the three observation-based products provide similar mean spatial seasonal patterns in oceanic CO₂ fluxes ( Supplementary Figure S11 ).

Our study shows that atmospheric XCO₂ IAV is affected most strongly by air-sea CO₂ fluxes in the Southern Hemisphere ( Figure 7 ), although the observation-based products show tradeoffs between the Southern Hemisphere low- and mid-latitude regions in terms of which region dominates. In the simulations with time-varying air-sea fluxes, the Southern Hemisphere, and to a lesser extent the tropics, dominate ocean-driven XCO₂ IAV even in the Northern Hemisphere. This pattern generally holds in the simulations with cyclostationary fluxes, although the Northern Hemisphere regions have a relatively larger contribution.

Given the small signature of ocean fluxes on atmospheric XCO₂ IAV, attributing ocean-induced IAV based on space-based observations will be challenging. Mitchell et al. (2023) provide a detailed assessment of time and space scales of variation in OCO-2 XCO₂ over North American land and coastal ocean using a geostatistical approach. They identify synoptic scale variations as contributing up to 2 ppm² variance, mesoscale transport, and correlated error as contributing up to 1 ppm² variance, and random noise as up to 1 ppm. Given that our simulations suggest the imprint of ocean IAV is less than 0.1 ppm, these results suggest that directly observing ocean IAV is not possible with the current technology for space-based remote sensing of atmospheric CO₂. Rather, XCO₂ observations will require precision better than 0.1 ppm, or 0.025% on a ~400 ppm background to detect and attribute ocean-driven variation.

Furthermore, because transport is an important contributor to the patterns of atmospheric CO₂ IAV, using methods such as atmospheric inversions to back out ocean fluxes in an optimal estimation framework requires fidelity in atmospheric transport modeling (Schuh et al., 2019). Here, patterns of atmospheric transport-induced IAV from cyclostationary air-sea CO₂ fluxes ( Figure 6 ) are similar because all flux products were transported through the same GEOS-Chem transport model. The choice of a different transport model would likely result in different spatial patterns, albeit with a similarly small magnitude compared to the IAV amplitude of the observations.

During the four years from 2014-2021 overlapping with OCO-2 observations, all three observation-based flux products show a multi-year decrease in the net ocean flux. There is low correspondence between the simulated XCO₂ IAV driven by ocean fluxes with total IAV detected from OCO-2. We hypothesize that multi-decadal variability in air-sea CO₂ fluxes may be more detectable in atmospheric XCO₂ than the year-to-year variability presented here, but testing this hypothesis requires ongoing space-based observations of XCO₂.

Our results temper the optimistic results from Chatterjee et al. (2017), who showed that a large ocean flux signal could be discerned from space-based XCO₂ data. Our study shows that the Southern Hemisphere oceans were the largest ocean-driven XCO₂ IAV in our simulations ( Figure 4 ), and that tropical Pacific that was most sensitive to IAV in atmospheric transport ( Figure 6 ).

Nevertheless, all three ocean products we analyze suggest a smaller XCO₂ IAV than what was observed for the 2015-16 El Niño. These results suggest that flux anomalies associated with changing modes in ocean oscillations in other basins may impart smaller variations in the atmosphere that are difficult to detect using OCO-2 or a similar satellite.

We evaluate the imprint that air-sea CO₂ fluxes from the whole ocean and different oceanic subregions leave on the atmospheric XCO₂ interannual variation. We simulate ocean-driven XCO₂ using the GEOS-Chem atmospheric transport model and air-sea CO₂ fluxes estimated from observation-based surface pCO₂ products and compare against the observed total XCO₂ IAV based on OCO-2 column-mean CO₂ observation from late 2014 to 2022. The ocean-driven IAV amplitude (standard deviation) in atmospheric XCO₂ caused by air-sea CO₂ exchange and IAV in atmospheric transport is generally between 0.08 and 0.12 ppm. While this magnitude is up to 40% of the total IAV in OCO-2 XCO₂ over the tropical and subtropical ocean basins, it is well below random noise in individual OCO-2 soundings and systematic errors in the satellite observations themselves. These results indicate that direct observation of air-sea CO₂ flux variations from total column XCO₂ would be very challenging with current space-based sensors.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

YG: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. GM: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing. AF: Conceptualization, Data curation, Methodology, Resources, Software, Writing – review & editing. SD: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing, Validation. GK: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing.