Data-based estimates of interannual sea–air CO₂ flux variations 1957–2020 and their relation to environmental drivers

This study considers year-to-year and decadal variations in as well as secular trends of the sea–air CO₂ flux over the 1957–2020 period, as constrained by the pCO₂ measurements from the SOCATv2021 database. In a first step, we relate interannual anomalies in ocean-internal carbon sources and sinks to local interannual anomalies in sea surface temperature (SST), the temporal changes in SST ($\mathrm{dSST}/\mathrm{d}t$), and squared wind speed (u²), employing a multi-linear regression. In the tropical Pacific, we find interannual variability to be dominated by $\mathrm{dSST}/\mathrm{d}t$, as arising from variations in the upwelling of colder and more carbon-rich waters into the mixed layer. In the eastern upwelling zones as well as in circumpolar bands in the high latitudes of both hemispheres, we find sensitivity to wind speed, compatible with the entrainment of carbon-rich water during wind-driven deepening of the mixed layer and wind-driven upwelling. In the Southern Ocean, the secular increase in wind speed leads to a secular increase in the carbon source into the mixed layer, with an estimated reduction in the sink trend in the range of 17 % to 42 %. In a second step, we combined the result of the multi-linear regression and an explicitly interannual pCO₂-based additive correction into a “hybrid” estimate of the sea–air CO₂ flux over the period 1957–2020. As a pCO₂ mapping method, it combines (a) the ability of a regression to bridge data gaps and extrapolate into the early decades almost void of pCO₂ data based on process-related observables and (b) the ability of an auto-regressive interpolation to follow signals even if not represented in the chosen set of explanatory variables. The “hybrid” estimate can be applied as an ocean flux prior for atmospheric CO₂ inversions covering the whole period of atmospheric CO₂ data since 1957.

1 Introduction

The atmospheric CO₂ content has risen during the recent decades, primarily due to anthropogenic emissions (IPCC, 2013). However, the actual rise has been co-determined by the exchange of CO₂ between the atmosphere and natural systems, notably the ocean and the land vegetation. The uptake of atmospheric CO₂ into the ocean is primarily driven by the solution disequilibrium across the sea–air interface. As the surface-ocean carbon content is lagging behind the atmospheric rise, the ocean uptake is, to first order, increasing in parallel with the atmospheric CO₂ rise. However, natural climate variability and anthropogenic climate change alter the uptake rate on year-to-year and decade-to-decade timescales as well as in its secular evolution. This leads to a feedback loop: atmospheric CO₂ influences the climate via the greenhouse effect, while the climate in turn influences the carbon-relevant natural systems in the ocean and on land. This feedback loop could dampen or accelerate climate change.

In order to understand the future climate trajectory, we therefore need to quantitatively understand the carbon response of the natural systems. For example, how will secular trends towards higher wind speeds in the Southern Ocean affect the sea–air CO₂ exchange in this region (Le Quéré et al., 2007; Hauck et al., 2013; and many others)? While the relevant timescale is secular (multi-decadal) trends, year-to-year or decade-to-decade variability in CO₂ fluxes can be used as “natural experiments” to understand the climatic controls of the land and ocean carbon cycle. This can be done by quantifying variations in carbon fluxes from suitable observations and statistically relating them to variations in quantities describing relevant environmental conditions. Even though the climate–carbon cycle feedback loop involves the global CO₂ fluxes only (because atmospheric CO₂ is mixed globally within about 1 year), the statistical analysis needs to be done on a spatial scale fine enough to accommodate the spatial inhomogeneity of the involved processes.

Suitable observational data therefore need to provide sufficient spatial and temporal detail and span several decades. Regarding ocean CO₂ fluxes, there are essentially two types of such data: (1) sustained atmospheric CO₂ measurements at various locations worldwide (Keeling, 1978; Conway et al., 1994; Francey et al., 2003; and many more) and (2) sustained and spatially extensive measurements of the CO₂ partial pressure (pCO₂) in the surface ocean (Bakker et al., 2016). As changes and gradients in atmospheric CO₂ reflect the sum of the regional CO₂ sources and sinks at the surface, atmospheric CO₂ data have been combined with simulations of atmospheric tracer transport and inverse techniques to estimate spatial and temporal variations in the CO₂ fluxes (“atmospheric inversion”; Newsam and Enting, 1988; Rayner et al., 1999; Bousquet et al., 2000; Rödenbeck et al., 2003; Baker et al., 2006; and many others). Even though most of the atmospheric inversions start in the 1990s or 2000s, when more and more stations became operational, the longest time series of atmospheric CO₂ measurements are available from 1957 (as used in Rödenbeck et al., 2018a). However, atmospheric inversions are known to have limited capability to correctly assign signals to land or ocean (Peylin et al., 2013). While the resulting error is relatively small for the land fluxes, it strongly affects the estimated ocean flux variability because the ocean variability is much smaller than the land variability.

Therefore, the surface-ocean pCO₂ data (Bakker et al., 2016) currently provide the most detailed information about the spatio-temporal variability in the sea–air CO₂ exchange. To cope with the very inhomogeneous distribution of these pCO₂ data in space and time, including substantial gaps, several methods have been developed to map (interpolate) the data into continuous spatio-temporal fields of pCO₂ (Takahashi et al., 2009; Watson et al., 2009; Valsala and Maksyutov, 2010; Landschützer et al., 2013; Nakaoka et al., 2013; Rödenbeck et al., 2013; Majkut et al., 2014; Iida et al., 2015; Jones et al., 2015; Zeng et al., 2015; Denvil-Sommer et al., 2019; Gregor et al., 2019; and several others). Most of these mappings employ either (i) an auto-regressive interpolation that fills unobserved areas or periods based on the neighbouring data within some prescribed correlation radii in space and time or (ii) a regression of pCO₂ against suitable explanatory variables that have been observed more densely and over the entire target period (using linear regression, neural networks, or machine learning). These two types of mappings offer complementary advantages, as regressions against explanatory variables possess predictive skill allowing longer data gaps to be filled (and potentially extrapolation into data-void periods), while auto-regressive mappings can reproduce all signals in the data even if they are not represented in the chosen explanatory variables (Rödenbeck et al., 2015). From the mapped pCO₂ fields, the sea–air CO₂ flux is then calculated via a gas exchange parameterization. In addition to studying the ocean carbon cycle, these flux estimates have also been used as an interannually varying ocean prior in atmospheric CO₂ inversions to potentially improve land CO₂ flux estimates (Rödenbeck et al., 2014).

With regard to the aim of understanding how the oceanic carbon cycle may respond to decadal and secular climatic changes as laid out above, however, the current pCO₂ mappings have two limitations. As a first limitation, the current pCO₂ mappings only provide spatio-temporal variations in the pCO₂ field and the sea–air CO₂ flux but do not explicitly quantify the relationships between these variations and underlying environmental drivers. This is true even for the regressions against explanatory variables: even though these relationships are implicitly contained in the synaptic weights of neural networks or similar parameters in machine learning algorithms, they are not accessible from these algorithms in interpretable form.

The second limitation arises from the fact that very few pCO₂ data exist before the mid-1980s (Bakker et al., 2016). In the equatorial Pacific, critical due to its large variability, sufficient coverage does not start before 1992. Despite their predictive skill, even the available pCO₂ regressions against explanatory variables only cover a time period not longer or even shorter than the pCO₂ data period, some for example because chlorophyll a data have only been available in the satellite era since 1997. Thus, none of the currently available pCO₂ mappings start before 1980. Consequently, they cannot be used as a data-based ocean prior in atmospheric CO₂ inversions over the full period of atmospheric data (1957–present). Further, the pCO₂ mappings do not cover the 1960–present period considered in ongoing synthesis projects like the annual carbon budget by the Global Carbon Project (GCP) (Friedlingstein et al., 2020), which so far exclusively relies on process model simulations during the first decades.

As a contribution to overcome these two limitations, this study has a 2-fold aim:

• First, extending the CarboScope pCO₂ mapping (Rödenbeck et al., 2013, 2014), we have developed a multi-linear regression explicitly estimating the sensitivities of the carbon sources and sinks in the oceanic mixed layer against the variations in relevant explanatory variables. This allows a data-based view of the processes plausibly underlying year-to-year variability in different parts of the ocean.

• Second, we have combined this multi-linear regression with an additive auto-regressive correction into a “hybrid” mapping, inheriting the complementary advantages of both auto-regressive and regression-based pCO₂ mappings. As the regression extrapolates the variability back to 1957 by only using explanatory variables available throughout the entire time frame, the hybrid mapping yields an observation-based estimate of the spatio-temporal variability in sea–air CO₂ fluxes since 1957.

After describing the mapping methods (Sect. 2), we present how the multi-linear regression traces the origin of interannual variations in the oceanic carbon system to the individual environmental quantities used as explanatory variables (Sect. 3.1). We present the spatial patterns in the regression coefficients (sensitivities) and discuss possible underlying mechanisms controlling the oceanic carbon system (Sect. 3.2). We evaluate the predictive skill of the multi-linear regression step as one of its most important requirements (Sect. 3.3). Finally, we present the interannual variations in sea–air CO₂ fluxes estimated by the hybrid mapping (Sect. 3.4) and compare it to the variations captured by the multi-linear regression (Sect. 3.5). In the discussion, we consider whether the presented multi-linear regression indeed meaningfully reflects biogeochemical processes (Sect. 4.1), which fraction of interannual variability it is able to capture (Sect. 4.2), to which extent the sensitivities depend on the timescale (Sect. 4.3), and how some uncertainties may affect the result (Sect. 4.4–4.6). In the Appendix, we focus on the global total sea–air CO₂ flux estimated by the hybrid mapping in terms of its mean (Sect. A1) and secular trend (Sect. A2), discussing its uncertainty and comparing it with literature values obtained by other methods.

2 Method

2.1 pCO₂ mapping

2.1.1 Overview

The pCO₂ mapping schemes used in this study are variants of the CarboScope pCO₂ mapping described in Rödenbeck et al. (2013). The estimates are based on the pCO₂ data (converted from the original fugacity data; see Table 1) in the SOCAT data collection version v2021 (Bakker et al., 2016, 2020). The elements common to all mapping variants are summarized in the following and illustrated in Fig. 1; for details we refer to Rödenbeck et al. (2013).

Figure 1 Illustration of the quantities involved in the mixed-layer scheme (time series panels) and the calculations done to connect them (thick-framed boxes). At the arrows on the right of each calculation box, we give its most important environmental input fields (see Table 1). The time series represent the example pixel enclosing the TAO140W mooring location (2^∘ N, 140^∘ W) in the tropical Pacific; they are taken from the results of this study but shown here for illustration only. Left: quantities on the original daily time steps, plotted for five example years. Right: the same quantities displayed as smoothed yearly averages, which is the way all results are shown in this paper. The background shading indicates the El Niño–Southern Oscillation (ENSO) phase (multivariate El Niño index (MEI) by Wolter and Timlin, 1993).

Bakker et al. (2016, 2020)Hirt and Rexer (2015)Good et al. (2013)Titchner and Rayner (2014)de Boyer Montégut et al. (2004)Tsujino et al. (2018)Tsujino et al. (2018)Rödenbeck et al. (2018b)Lee et al. (2006)Conkright et al. (2002)Lee et al. (2006)Garcia et al. (2006)DeVries (2022)

Table 1 Input data sets.

SST: sea surface temperature; MLD: mixed-layer depth; LOCEAN: Laboratoire d'océanographie et du climat: expérimentations et approches numériques; NCEP: National Centers for Environmental Prediction; SOCAT: Surface Ocean CO₂ Atlas; WOCE: World Ocean Circulation Experiment; WOA: World Ocean Atlas; param.: parameterizations; expl. var.: explanatory variable.

Parameterizations of sea–air gas exchange (quadratic wind speed dependence as in Wanninkhof, 1992) and solubility (Weiss, 1974), a calculation of the chemical equilibrium of the carbonate chemistry in seawater (Orr and Epitalon, 2015) as well as a mixed-layer budget of dissolved inorganic carbon (DIC) (Rödenbeck et al., 2013), are used to express the pCO₂ field and the sea–air CO₂ flux field as a function of the ocean-internal flux of DIC, f_int (Fig. 1). The ocean-internal DIC flux f_int is meant to comprise all sources and sinks of DIC into or out of the oceanic mixed layer, through biological conversion within the mixed layer or through mixing-in of waters with different DIC concentration. It is expressed as the sum of a fixed (a priori) flux field and a set of predefined spatio-temporal patterns of adjustment each scaled by an adjustable parameter (the sets of patterns are detailed for each variant of the mapping below). Then, the mismatch between the calculated pCO₂ field (at the respective pixels and time steps containing the SOCAT pCO₂ samplings) and the corresponding measured pCO₂ values (black dots in the pCO₂ panel of Fig. 1) is gauged by a quadratic cost function. The (a posteriori) estimates of the mapping are calculated from those values of the adjustable parameters that minimize this cost function. In the example of Fig. 1, the two estimates (coloured) follow the data points (black dots) more closely than the prior (grey).

Table 2 Mapping runs used in this study. The main results are given in bold; the other runs are used to assess uncertainty (“uncertainty cases”; Sect. 2.2), to illustrate specific points of discussion (“test cases”; Sect. 2.2), or to assess predictive skill (“cross-validation”; Sect. 2.3).

^* Regression run only over 1998–2019, 1985–2009, or 1988–2018, respectively.

Spatial and temporal interpolation between the very inhomogeneously sampled data is implemented in the following way. By choosing a set of spatial patterns of adjustment that are centred at all the individual ocean pixels but simultaneously affect the respective neighbouring pixels within some correlation radius (to be detailed below), in conjunction with additional Bayesian terms in the cost function that penalize large adjustments to the adjustable parameters, the parameter fields (the ocean-internal DIC flux field or the fields of sensitivities, respectively; see below) are forced to be smooth. These smoothness constraints spread the information from data-covered pixels to neighbouring unconstrained pixels (see Fig. 5 of Rödenbeck et al., 2013), thereby interpolating spatial data gaps. (The set of patterns of adjustment indirectly defines the Bayesian a priori covariance matrix; see Rödenbeck, 2005, for background.) Interpolation in time is achieved analogously by temporal smoothness constraints (even though, for practical reasons, a mathematically equivalent Fourier formulation is used).

The four mapping variants used here (Table 2) differ in the choices of the prior for f_int and the set of spatio-temporal patterns of adjustment. Our development started from a variant (Sect. 2.1.2) essentially identical to Rödenbeck et al. (2013) used as the CarboScope pCO₂ mapping before version v2020, except for some technical changes described later (Sect. 2.1.6–2.1.7). As an intermediate modification, we introduced a prior stabilizing the secular trend (Sect. 2.1.3); the result of this variant will be used to help discuss specific aspects. The main results of this study come from the multi-linear regression (Sect. 2.1.4) and the hybrid mapping (Sect. 2.1.5). Figure 2 summarizes the differences and the flow of information between the four variants.

Figure 2 Information flow between the presented mapping runs. Thick-framed boxes denote calculations, with arrows denoting their input and output data sets (mostly spatio-temporal fields). The violet, green, orange, and blue arrows represent the spatio-temporal sea–air CO₂ fluxes estimated by our main mapping runs M_ZE, M_E, M_R, and M_H (these same four colours are also used in all line plots in this paper). All mapping runs use the SOCAT database of pCO₂ measurements (point data, black) as their primary information source. The runs are algorithmically identical, except for the representation of the ocean-internal DIC sources and sinks (f_int) and the corresponding set of adjustable unknowns; this part of the mapping algorithm has therefore been represented explicitly as separate parts of the boxes at their left-hand sides, labelled by the respective equation number. Runs M_ZE and M_E are using the representations with explicit interannual degrees of freedom, either without f_int prior (Eq. 1) or using the decadally smoothed OCIM result as f_int prior (cyan; Eq. 2). Run M_R is using the representation involving regression terms (Eq. 3), which requires the explanatory variables (magenta) as further input fields; it uses the same f_int prior as the run M_E. The representation of f_int in the “hybrid” run M_H again has explicit interannual degrees of freedom but no long-term and seasonal degrees of freedom (Eq. 4). Importantly, it uses the ocean-internal DIC sources and sinks estimated by the “multi-linear regression” M_R as its f_int prior (cyan again). All mapping calculations use the various input fields shown in Fig. 1, which are however omitted here for clarity. More technically, the mappings also need a linearization of the non-linear dependence of pCO₂ on DIC, consisting of three fields (the derivatives dpCO₂ $/$ dDIC as well as reference fields for pCO₂ and DIC) together depicted by the grey arrows. Box L denotes the calculation of dpCO₂ $/$ dDIC and the reference DIC field from a reference pCO₂ field (light blue) as described in Sect. 2.1.6 (the further input fields required by this calculation are not depicted). For the main mappings M_ZE, M_E, M_R, and M_H, the reference pCO₂ field comes from the data-based estimate of the pre-mapping M_P. The linearization for M_P, in turn, uses the atmospheric pCO₂ field as pCO₂ reference (light blue again).

2.1.2 The “zero-prior explicitly interannual” pCO₂ mapping (ZE)

The starting variant has a general set of (many) patterns of adjustment, allowing an arbitrary smooth spatio-temporal internal DIC flux field $f_{\mathrm{int}}^{\mathrm{ZE}}$ (Rödenbeck et al., 2013). This field $f_{\mathrm{int}}^{\mathrm{ZE}}$ is implemented as the sum of a constant term (subscript “LT” for “long-term”) and terms for seasonal (subscript “Seas”) and interannual anomalies (non-seasonal, subscript “IAV”):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{f}_{\mathrm{int}}^{\mathrm{ZE}}(x,y,t)& =f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}(x,y,t)\\ & +f_{\mathrm{int},\mathrm{LT}}^{\mathrm{ADJ}}(x,y)+f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{ADJ}}(x,y,s).\end{array}\end{array}$$

As indicated by the superscript “adj” or “ADJ” (difference explained below), all these terms involve degrees of freedom being adjusted in the cost function minimization sketched above. A priori, all adjustable terms are zero, such that the prior of $f_{\mathrm{int}}^{\mathrm{ZE}}$ is zero as well.

The interannual term $f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}(x,y,t)$ can represent non-seasonal anomalies on all month-to-month, year-to-year, or decadal timescales, including secular trends. The level of its temporal smoothness corresponds to a priori correlation length scales of about 4 weeks, implemented through a mathematically equivalent Fourier series with dampened higher-frequency components (where Fourier terms dampened to less than 2 % are discarded entirely). This amounts to 722 scalable Fourier terms for our 71-year calculation period 1951–2021. The seasonal term $f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{ADJ}}$ only contains seasonal Fourier components; thus it only depends on the time s within the year and repeats itself every year. Along the seasonal cycle, it has the same temporal correlation length as the interannual term of about 4 weeks, amounting to 10 scalable Fourier terms. The constant term $f_{\mathrm{int},\mathrm{LT}}^{\mathrm{ADJ}}$ is not time-dependent by definition (1 temporal degree of freedom).

Spatially, the level of smoothness in all three terms corresponds to a priori correlation length scales of about 640 km in longitude and latitude.

As symbolized by the capitalized superscript “ADJ”, the a priori uncertainties in the seasonal Fourier terms of $f_{\mathrm{int},\mathrm{LT}}^{\mathrm{ADJ}}$ and $f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{ADJ}}$ are chosen to be enlarged relative to the non-seasonal Fourier terms of $f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}$, corresponding to larger expected amplitudes of seasonal variations in f_int compared to non-seasonal ones. In terms of the implied a priori autocorrelation function, these enhanced a priori uncertainties in seasonal variations are equivalent to non-zero temporal correlations between the flux at any given time of year and the same time of year in all other years (in addition to the 4-week decaying correlations mentioned above). Due to these periodic correlations, f_int in time periods without data does not fall back to the prior (here zero) but to the mean seasonal cycle $f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{ADJ}}$ as constrained by the data-covered periods.

2.1.3 The “explicitly interannual” pCO₂ mapping (E)

In order to stabilize the secular trend in the early decades (as discussed in Sect. A2 below), we now add a fixed (i.e. non-adjustable) term (superscript “fix”):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{f}_{\mathrm{int}}^{\mathrm{E}}(x,y,t)& =f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}(x,y,t)\\ & +f_{\mathrm{int},\mathrm{LT}}^{\mathrm{ADJ}}(x,y)+f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{ADJ}}(x,y,s)\\ & +f_{\mathrm{int},\mathrm{Decad}}^{\mathrm{fix}=\mathrm{OCIM}}(x,y,t).\end{array}\end{array}$$

Consequently, the prior of $f_{\mathrm{int}}^{\mathrm{E}}$ is given by this fixed term. It is obtained from the sea–air flux product by DeVries (2022), which is based on an abiotic carbon cycle model that captures the rising atmospheric CO₂ boundary condition and is embedded in a data-driven model of time-mean ocean circulation (OCIM). The OCIM fluxes have been decadally smoothed (indicated by subscript “Decad”) because the OCIM result originally represents sea–air fluxes including SST-related interannual variations, which are created by our parameterizations already (Fig. 1).

2.1.4 The “multi-linear regression” (R)

In the third variant, the ocean-internal DIC flux is represented as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{f}_{\mathrm{int}}^{\mathrm{R}}(x,y,t)& =\sum _i{\mathit{\gamma }}_i^{\mathrm{adj}}(x,y)\cdot V_i(x,y,t)\\ & +f_{\mathrm{int},\mathrm{LT}}^{\mathrm{ADJ}}(x,y)+f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{ADJ}}(x,y,s)\\ & +f_{\mathrm{int},\mathrm{Decad}}^{\mathrm{fix}=\mathrm{OCIM}}(x,y,t).\end{array}\end{array}$$

Compared to Eq. (2), the degrees of freedom representing interannual variations (f_int,IAV) are replaced here by a multi-linear function involving three explanatory fields (V_i):

• sea surface temperature (SST),

• its temporal change ($\mathrm{dSST}/\mathrm{d}t$), and

• squared wind speed (u²).

There is a 2-fold motivation behind this choice of explanatory variables: (1) variations in carbon-relevant processes (e.g. carbon and nutrient input into the mixed layer, stratification, mixing, entrainment, wind-driven deepening of the mixed layer) are expected to also be related to these variables, and (2) observation-based data sets for SST and u are available over our entire calculation period 1951–2021 (u at least from reanalysis). The specific input data sets used in our base case are given in Table 1.

The simultaneous use of SST and $\mathrm{dSST}/\mathrm{d}t$ is motivated as it is changes in SST that are related to DIC fluxes (i.e. changes in DIC). Moreover, the sum of SST and $\mathrm{dSST}/\mathrm{d}t$ mathematically allows a temporal shift between SST and f_int for a dominant Fourier mode (similar to sine and cosine terms).

To prevent confusion, we point out that the multi-linear regression as introduced here is set up in terms of the ocean-internal DIC flux f_int (see Sect. 2.1.1), not in terms of pCO₂ or sea–air flux as done in various other studies in the literature. This also means that important processes (SST dependence of solubility and carbonate chemistry, wind speed dependence of gas exchange) are not included into the regression Eq. (3) but are already taken care of by the parameterizations listed in Sect. 2.1.1 and Fig. 1.

All the explanatory fields V_i are implemented on a monthly timescale, smoothly transformed onto our daily time steps. The scaling factors ${\mathit{\gamma }}_i^{\mathrm{adj}}$ between the internal DIC flux and these explanatory fields V_i are taken as the adjustable degrees of freedom in the cost function minimization (very analoguous to the “NEE-T inversion” of Rödenbeck et al., 2018b). These unknown scaling factors are allowed to vary spatially (with correlation length of about 2000 km in longitude and 1000 km in latitude, thus more smoothly than the direct adjustments of f_int in the explicitly interannual mapping of Sect. 2.1.3), but are constant in time (1 temporal degree of freedom per explanatory field per pixel). All three regression terms are normalized such that the a priori uncertainty in their global integral on 1 July (averaged over the 1 July time steps of all years within the analysis period 1957–2020) is the same as that of f_int,IAV in Eq. (2) (1 July is an arbitrary choice, in line with the normalization with respect to the flux in the middle of the final year used in CarboScope so far.)

In order to avoid influences of the spin-up transient on the regression coefficients (estimated sensitivities), the regression terms (first line of Eq. 3) only cover the analysis period 1957–2020, while the remaining years before and after are filled by explicitly interannual degrees of freedom just as $f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}$ in Eq. (2). For clarity, this detail has been omitted from Eq. (3).

2.1.5 The “hybrid” pCO₂ mapping (H)

The final variant aims to combine the temporal extrapolation capability of the multi-linear regression (Sect. 2.1.4) and the flexibility to reproduce observed signals of the explicitly interannual mapping (Sect. 2.1.3). Technically being an explicitly interannual mapping itself, its representation of the ocean-internal DIC flux,

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{f}_{\mathrm{int}}^{\mathrm{H}}(x,y,t)& =f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}(x,y,t)\\ & +f_{\mathrm{int},\mathrm{LT}}^{\mathrm{adj}}(x,y)+f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{adj}}(x,y,s)\\ & +f_{\mathrm{int}}^{\mathrm{fix}=\mathrm{R}}(x,y,t),\end{array}\end{array}$$

is similar to Eq. (2), but with the following two changes:

• As the essential change, the interannually varying result of the multi-linear regression (Sect. 2.1.4) is used as a prior for the internal DIC flux ($f_{\mathrm{int}}^{\mathrm{fix}=\mathrm{R}}(x,y,t)$) instead of the decadally smoothed OCIM result only containing decadal variations and the secular trend.

• As a merely technical change, the a priori uncertainties in the mean flux $f_{\mathrm{int},\mathrm{LT}}^{\mathrm{adj}}(x,y)$ and the seasonality $f_{\mathrm{int},\mathrm{Seas}}^{\mathrm{adj}}(x,y,s)$ are not enhanced with respect to non-seasonal variability $f_{\mathrm{int},\mathrm{IAV}}^{\mathrm{adj}}(x,y,t)$ any more (indicated by the lower-case superscript “adj” in all three terms) because the prior $f_{\mathrm{int}}^{\mathrm{fix}=R}(x,y,t)$ already contains a long-term mean and a mean seasonal cycle.

In essence, the hybrid mapping thus adds an interannually varying correction to the multi-linear regression. Due to this construction, the hybrid result will fall back to the multi-linear regression during periods without data, but it is nevertheless able to fit pCO₂ signals on month-to-month, year-to-year, and decadal timescales that have not yet been reproduced via the explanatory variables of the multi-linear regression.

Methodological note. Mathematically, the hybrid run is equivalent to estimating the additive correction to the multi-linear regression from the pCO₂ residuals of the multi-linear regression. That is, the signals being used by the hybrid run are those that could not yet be explained by the multi-linear regression. The hybrid run is thus similar to a hypothetical joint run simultaneously having regression degrees of freedom (like the multi-linear regression) and explicitly interannual degrees of freedom (like the explicitly interannual estimate). We abandoned the concept of such a joint run, however, because it would face two problems: (1) its result would depend on the relative a priori weighting between the two groups of degrees of freedom, for which there is no clear information, and (2) the explicitly interannual degrees of freedom would necessarily also absorb part of the signals actually proportional to the explanatory variables. Running the multi-linear regression and the hybrid step sequentially, as done here, reduces both problems.

2.1.6 The pre-mapping (P): determining the linearization of the carbonate chemistry

In contrast to Rödenbeck et al. (2013), we now allow for the secular trend in the Revelle factor. We deem this necessary due to our longer period of interest 1957–2020, during which the mixed-layer carbon content notably increased, leading to shifts in the relation between variations in the ocean-internal DIC flux (f_int) and the sea–air CO₂ flux. As our scheme extrapolates the seasonality (and in the “multi-linear regression” also the interannual variations) from the data-constrained recent decades to the almost unconstrained earlier decades through correlations in f_int (see the last paragraph of Sect. 2.1.2), the shifting relation has the potential to alter the amplitude of flux variations in the earlier decades.

As in Rödenbeck et al. (2013), the non-linear dependence of pCO₂ on DIC is linearized around reference fields $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}$ and DIC_Ref:

$$\begin{array}{}\text{(5)}& p{\mathrm{CO}}_{\mathrm{2}}=p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}+\left({\displaystyle \frac{\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{d}\mathrm{DIC}}}\right)(\mathrm{DIC}-{\mathrm{DIC}}_{\mathrm{Ref}}).\end{array}$$

The linearization is needed to be able to use the fast minimization algorithm in the CarboScope software. Previously in Rödenbeck et al. (2013), the reference fields $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}$ and DIC_Ref were temporally constant and had been taken from observation-based data sets not guaranteed to be mutually consistent, and the derivative $(\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}\mathrm{DIC})$ had been calculated from these via approximation formulas. In order to now include the secular trend in Revelle factor (and simultaneously to remove the mentioned approximations), we employ the mocsy package (Orr and Epitalon, 2015), which provides routines to accurately calculate pCO₂ and $(\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}\mathrm{DIC})$ from a given field of DIC (and from fields of alkalinity, SST, salinity, silicate, phosphate, and air pressure, which we take from external sources; Table 1). Using an adjusted Newton algorithm calling mocsy iteratively, we obtain an algorithm to calculate (reference) DIC and $(\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}\mathrm{DIC})$ from a given (reference) pCO₂ value at each location and time (box L in Fig. 2). The $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}$ field is obtained as the posterior pCO₂ field of a “pre-mapping” run (P, the leftmost one in Fig. 2). The $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}$ and $(\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}\mathrm{DIC})$ fields used in this pre-mapping run, in turn, are calculated from a preliminary reference identical to atmospheric pCO₂. This yields a reasonable starting point because the atmospheric pCO₂ field does already contain the secular CO₂ rise, which is the most important feature in this context.

Potentially, we might expect to need a loop with further pre-mappings, each getting its $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}$ field from the posterior pCO₂ field of the respective previous one. However, we confirmed by explicit testing that the fields are not appreciably altered any more after the first pre-mapping; thus a single pre-mapping is sufficient.

All other mapping runs of this study use the same spatio-temporal linearization fields $p{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{Ref}}$, DIC_Ref, and $(\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}\mathrm{DIC})$ as calculated by the pre-mapping.

2.1.7 Technical details common to all variants

As in Rödenbeck et al. (2013), the pCO₂ data comprise the individual observations from file https://www.ncei.noaa.gov/data/oceans/ncei/ocads/data/0235360/SOCATv2021.tsv, last access: 1 June 2021, including all observations flagged A–D. The additional file flagged E was not used.

In contrast to Rödenbeck et al. (2013), the analysis period now starts in 1957 (chosen in light of the potential use of the results as a prior in atmospheric inversions). The actual calculation period of all runs starts in 1951. According to explicit tests, this allows the initial transient of the mixed-layer DIC budget equation to decay by 1957.

As in Rödenbeck et al. (2013), the calculation period includes 1 more year (“spin-down”, here 2021) after the valid period constrained by the data (until end of 2020), in order to avoid numerical edge effects.

In order to cover the entire calculation period since 1951, we now use SST from Hadley EN.4.2.1 (Good et al., 2013) and sea ice concentration from HadISST 2.2.0.0. (Titchner and Rayner, 2014, https://www.metoffice.gov.uk/hadobs/hadisst2/data/HadISST.2.2.0.0_sea_ice_concentration.nc.gz, last access: 5 June 2020).

Compared to Rödenbeck et al. (2013), the spatial resolution of all the mapping calculations has been increased to 2.5^∘ longitude × 2^∘ latitude (previously on the grid of the TM3 atmospheric transport model, 5^∘ × 4^∘). Moreover, the adjustments are now done over the entire ocean (i.e. we do not fix part of the temporally ice-covered regions any more).

2.2 Uncertainty and test cases

In order to explore how robust the results of the multi-linear regression (Sect. 2.1.4) are, we also perform uncertainty cases where certain set-up parameters are modified within ranges deemed as plausible as the base case (Table 2):

RegrSSTNOAA

– using SST from NOAA_ERSST v5 (Huang et al., 2017) as an alternative data set for the explanatory variables SST and $\mathrm{dSST}/\mathrm{d}t$ (but no change to any other SST-dependent items such as solubility);

RegrU2NCEP

– using wind speeds from NCEP reanalysis (Kalnay et al., 1996) as an alternative data set for the explanatory variable u² (but no change to wind-dependent gas exchange);

RegrAdddSSTdt2

– additional regression term based on ($\mathrm{dSST}/\mathrm{d}t$)²;

RegrAddU4

– additional regression term based on u⁴;

RegrAddpaCO2

– additional regression term based on decadally smoothed p_aCO₂;

RegrNoDecad

– removing any decadal variability and secular trends from the explanatory fields V_i, such that the multi-linear regression term only represents interannual variability on a timescale of a few years;

RegrShort

– shorter spatial correlation lengths for the sensitivities ${\mathit{\gamma }}_i^{\mathrm{adj}}$ (Supplement Fig. S5);

RegrLoose

– a priori uncertainty in the sensitivities increased by a factor of 4 (i.e. the strength of the mathematical regularization is reduced);

MLDq2

– dividing mixed-layer depth by 2;

MLDx2

– multiplying mixed-layer depth by 2 (lacking a clear uncertainty range of mixed-layer depth, MLDq2 and MLDx2 represent a rather strong change, maybe already outside the actual uncertainty);

GasexLow

– weaker gas exchange by scaling the gas transfer velocity field such that its global mean matches the lower limit of the range 16.5±3.2 cm h⁻¹ (Naegler, 2009) rather than the central value;

GasexHigh

– stronger gas exchange (analogously, using upper limit);

GasexU1

– replacing the u² dependence of gas exchange by a $\left|u\right|$ dependence (while keeping the global mean gas transfer velocity the same).

GasexU3

– replacing the u² dependence of gas exchange by a $|u{|}^{\mathrm{3}}$ dependence (while keeping the global mean gas transfer velocity the same).

To help in the discussion of specific aspects, we performed further test cases (not necessarily as plausible as the base case):

RegrOnlySST, RegrOnlydSSTdt, RegrOnlyU2

– the explanatory variables used individually (i.e. the regression terms of the remaining two were omitted);

RegrAddChl_98r19

– addition of Chl a as a further explanatory variable (Fig. S7; chlorophyll concentration has been taken from the GlobColour project (Maritorena et al., 2010), which combined retrievals from the SeaWiFS (NASA), MODIS (NASA), MERIS (ESA), OLCI (ESA), and VIIRS (NOAA/NASA) satellites into a harmonized data set; as the Chl a data are only available for the years 1998–2019, the regression is restricted to this period, plus spin-up and spin-down periods);

RegrHeat_85r09

– replacing $\mathrm{dSST}/\mathrm{d}t$ with the net sea–air heat flux taken from the OAFlux project (Yu and Weller, 2007); regression period restricted to 1985–2009 according to the availability of the heat flux data set;

RegrCurl_88r18

– replacing u² with wind stress curl calculated from Cross-Calibrated Multi-Platform (CCMP) v2.0 wind speeds (Atlas et al., 2011); regression period restricted to 1988–2018 according to the availability of CCMP;

98r19, 85r09, 88r18

– using the same regression terms as in the base case but restricting the time period of regression to the same years as used for RegrAddChl_98r19, RegrHeat_85r09, and RegrCurl_88r18, respectively.

Uncertainties in the hybrid mapping (Sect. 2.1.5) were explored analogously by re-running the hybrid step with several of the uncertainty cases of the regression listed above (Table 2). Part of the involved set-up changes (mixed-layer depth, gas exchange) also affect the hybrid calculation itself.

2.3 Gauging the predictive skill of the multi-linear regression

In order to test whether the multi-linear regression against explanatory variables (Sect. 2.1.4) is actually meaningful, we determine its predictive skill. For this, the multi-linear regression is re-run six times, each time omitting the pCO₂ data from one of the 5-year periods 1985–1989, 1990–1994, 1995–1999, 2000–2004, 2005–2009, or 2010–2014. That is, each of the six test runs possesses an artificial data gap of 5 years, a duration chosen to be longer than typical features of year-to-year variability like El Niño. We can then compare the predictions during the data gaps with the results of the completely constrained run.

3 Results

The main results of this study are of two different types:

• From the multi-linear regression, we obtain spatial maps of the sensitivities γ_i (Eq. 3) relating the variations in the surface-ocean carbon system to variations in SST, $\mathrm{dSST}/\mathrm{d}t$, and u² (Sect. 3.2).

• The hybrid mapping yields a spatio-temporal estimate of the sea–air CO₂ flux over 1957–2020, in particular its evolution from year to year (Sect. 3.5).

Further results are presented for illustration and to elucidate the robustness of the main results.

3.1 Origin of interannual variations as estimated by the multi-linear regression

The multi-linear regression attempts to trace the interannual variations in the surface-ocean carbon system (and hence in the sea–air CO₂ flux) to the interannual variations in the chosen explanatory variables SST, $\mathrm{dSST}/\mathrm{d}t$, and u². In Fig. 3 (left panels), the estimated contributions of the three explanatory variables to the ocean-internal DIC flux are depicted for a subdivision of the ocean into five latitudinal bands. In the centre panels, the resulting contributions to the sea–air CO₂ flux are shown, as calculated by the parameterizations and the budget equation in our mapping scheme (Sect. 2.1.1). These contributions and the prior sum up to the total sea–air CO₂ flux, shown in the right panels together with our set of uncertainty results (Sect. 2.2).

Figure 3 Left: estimated contributions of the three explanatory variables in the multi-linear regression (as well as the prior, plotted here without its mean) to the ocean-internal DIC flux in five latitudinal bands (top to bottom). Centre: corresponding contributions to the sea–air CO₂ flux. Right: total sea–air CO₂ flux estimated by the multi-linear regression (base case, orange) together with the uncertainty cases listed in Sect. 2.2 (the cases with the largest impact on interannual variability – RegrSSTNOAA, RegrU2NCEP, RegrAddpaCO2, RegrNoDecad – are plotted explicitly in different colours; since the cases related to gas exchange – GasexLow, GasexHigh, GasexU1, GasexU3 – shift the long-term mean of the flux, the range of this shift has been indicated by the length of the vertical orange bars just to the right of each panel for clarity; the remaining uncertainty cases having rather small impact – RegrAdddSSTdt2, RegrAddU4, RegrLoose, RegrShort, MLDq2, MLDx2 – have been subsumed into the pale orange band depicting their envelope). All curves show interannual variations. The background shading indicates the El Niño phase according to the multivariate El Niño index (MEI) by Wolter and Timlin (1993). In the left and centre panels, fluxes are given in per-area units to emphasize the local process perspective, while fluxes in the right panels are given as regional integrals to emphasize their share in the total ocean flux.

When disregarding the secular increase in the ocean carbon sink, the largest year-to-year variations in the regionally integrated sea–air carbon flux are found in the tropics (Fig. 3, middle right), in particular the tropical Pacific (Fig. S4). Correspondingly, the year-to-year variations in the ocean-internal carbon flux (f_int) from the three terms in the multi-linear regression (Eq. 3) are largest in the tropics as well (Fig. 3, middle left). Of the three explanatory variables, the contribution of year-to-year variability in temporal SST changes ($\mathrm{dSST}/\mathrm{d}t$, black) is the largest. Concurrent with the warming (dSST $/$ dt>0) at the onset of each El Niño event (grey background stripes), we find a negative carbon flux anomaly (reduction in the carbon source in this region) because smaller amounts of cold, carbon-rich water are upwelling. At the end of each El Niño event, we find an analogous coupling of the cooling (dSST $/$ dt<0) and an additional carbon source to the mixed layer. The contribution of year-to-year variability in SST itself (red) is second-largest in the tropics, causing anomalous carbon sinks during El Niño events and anomalous carbon sources during La Niña conditions afterwards. This could be interpreted as a small correction to the $\mathrm{dSST}/\mathrm{d}t$ contribution: the sum of the $\mathrm{dSST}/\mathrm{d}t$ and SST contributions (not shown) is similar to the $\mathrm{dSST}/\mathrm{d}t$ contribution alone but slightly shifted in time by a few months. The smallest contribution to the year-to-year variability in the tropics is estimated for squared wind speed (u², light blue), with a temporal pattern relatively similar to that of the SST contribution. Due to the co-variation between SST and u² on a year-to-year timescale, these two explanatory variables could be partly confounded by the regression, though the detailed locations where their respective sensitivities are high do not actually overlap much (see Sect. 3.2 below).

In the high-latitude bands (top and bottom left panels of Fig. 3), the wind speed contribution is estimated to be larger than in the tropics, now on the same order of magnitude as the SST and $\mathrm{dSST}/\mathrm{d}t$ contributions or even larger. As a notable feature in the Southern Ocean (bottom left), the secular increase in wind speed leads to a secular increase in the carbon source into the mixed layer. Across our set of uncertainty cases (Sect. 2.2), the linear trend of the wind speed contribution over the 1960–2019 period in the ocean south of 45^∘ S is estimated in the range 0.002 to $\mathrm{0.005}\left(\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}\right)\mathrm{yr}{}^{-\mathrm{1}}$ (see Supplement Fig. S3, bottom, light-blue bars). As a secular trend in f_int (bottom left in Fig. 3) causes a secular trend in sea–air flux of the same size (bottom centre), it represents a reduction by 17 % to 42 % of the trend towards an increasing Southern Ocean sink strength (relative to the trend of $-\mathrm{0.012}\left(\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}\right)\mathrm{yr}{}^{-\mathrm{1}}$ estimated by OCIM over 1960–2019). A slowing-down of the Southern Ocean sink increase (compared to the increase expected from rising atmospheric CO₂) has also been found in model simulations and attributed to an increase in upwelling of old carbon by the accelerating winds (Le Quéré et al., 2007; Hauck et al., 2013; and many others). We need to note, however, that our multi-linear regression estimates the wind-speed-related trend only indirectly: as the sensitivities ${\mathit{\gamma }}_{u^{\mathrm{2}}}$ are presumably largely constrained by year-to-year variations (because they do not change much if the linear trend of the explanatory variables is removed; see sensitivity case RegrNoDecad, Sect. 4.3), the slope of the secular trend can only be correct to the extent that the sensitivity ${\mathit{\gamma }}_{u^{\mathrm{2}}}$ is identical for year-to-year and secular variations.

The year-to-year anomalies from the f_int contributions (Fig. 3, left) carry through to the sea–air CO₂ flux (centre) in a delayed and dampened fashion due to the buffer effect of carbonate chemistry in combination with the limited gas exchange. We also note again that the sea–air CO₂ flux contains additional year-to-year variability from solubility and gas exchange anomalies as represented by the involved parameterizations (also see Fig. 1 and Sect. 3.2.4 below).

3.2 Patterns of the sensitivity of ocean-internal DIC sources and sinks to interannual variations in SST, $\mathrm{dSST}/\mathrm{d}t$, and u² estimated by the multi-linear regression – which underlying processes do they suggest?

The estimated sensitivities of the ocean-internal DIC flux (f_int) against interannual variations in the chosen explanatory variables of the multi-linear regression (sea surface temperature SST, temporal changes in sea surface temperature $\mathrm{dSST}/\mathrm{d}t$, and squared wind speed u²) are shown in Fig. 4. Here we consider the most prominent features in these sensitivity patterns and mention oceanic processes that are compatible with these and may thus control surface-ocean biogeochemistry. Even though regression analysis cannot prove causation, we argue later (Sect. 4.1) why such a tentative attribution may be meaningful here. Also see Sect. 4.3–4.6 for further discussion on uncertainties.

Figure 4 Estimated sensitivities of the ocean-internal DIC flux f_int against interannual variations in the temporal changes in sea surface temperature (a), in the sea surface temperature itself (b), and in squared wind speed (c). Positive (negative) sensitivities mean that increases in the respective explanatory variable are associated with a stronger source (stronger sink) of DIC in the mixed layer.

3.2.1 Sensitivity of f_int to $\mathrm{dSST}/\mathrm{d}t$

We start with $\mathrm{dSST}/\mathrm{d}t$ (Fig. 4, top) as the explanatory variable contributing the largest year-to-year variability (Sect. 3.1). Events of decreasing SST are estimated to be associated with more positive ocean-internal DIC fluxes in the tropical Pacific (within a tilted band located around the Equator in the western tropical Pacific and around about 15^∘ S in the eastern tropical Pacific) and in most parts of the higher latitudes in both hemispheres (blue and cyan areas in Fig. 4a). Such a correlation would arise from variations in the upwelling of waters that are both colder and more carbon-rich than the mixed layer.

In the rest of the ocean, the absolute value of the sensitivity ${\mathit{\gamma }}_{\mathrm{dSST}/\mathrm{d}t}$ is small (light blue or light red). We assume that these sensitivities mainly reflect insignificant correlations, especially due to the higher uncertainty in regions of sparse data coverage or in regions where $\mathrm{dSST}/\mathrm{d}t$ is mainly driven by atmospheric heating or cooling. In particular, positive sensitivities are not compatible with any known oceanic mechanism.

3.2.2 Sensitivity of f_int to SST

The estimated sensitivity γ_SST between the interannual variations in the ocean-internal DIC flux and SST itself is rather patchy, with both positive and negative areas (Fig. 4, middle). This may reflect the fact that various biological processes contribute to f_int, depending on temperature in different ways and thus potentially cancelling each other. For example, carbon fixation (net primary productivity, NPP) will invigorate with increasing temperature (until a threshold is reached); as NPP represents a sink (i.e. a negative contribution to f_int), it would thus cause negative γ_SST sensitivities. Carbon export (or export ratio at least) is generally anticorrelated with temperature (Laws et al., 2000), thus causing positive γ_SST sensitivities, though also the opposite behaviour seems possible.

Positive interannual sensitivity to SST would also be compatible with a nutrient effect. Upwelling and mixing-in from below both decreases SST and increases the availability of nutrients. Thus, negative anomalies in SST tend to be associated with higher biological production and thus enhanced removal of carbon (negative anomalies in f_int). However, upwelling also brings up carbon, which is usually assumed to dominate the carbon signal. For example, Hauck et al. (2013) showed that – in the model – in the Southern Ocean south of 55^∘ S, there would be more biological export per increase in the Southern Annular Mode (SAM), which goes along with more upwelling. Yet, whether the carbon effect or the opposing nutrient effect dominates the upwelling signal is still controversially discussed.

As the statistical inference by our regression can only respond to the sum of all contributing processes, we therefore cannot draw specific conclusions from the estimated γ_SST pattern. In addition, the regression may adjust γ_SST to effectively shift the $\mathrm{dSST}/\mathrm{d}t$ contribution in time (Sect. 3.1).

3.2.3 Sensitivity of f_int to u²

Higher wind speeds are estimated to be associated with more positive ocean-internal DIC fluxes (stronger sources into or weaker sinks out of the mixed layer) along the Equator in the Pacific; in the eastern upwelling zones of the North Pacific, South Pacific, and South Atlantic; and in circumpolar bands in the high latitudes of both hemispheres (red and yellow areas in Fig. 4c). Such a positive sensitivity is compatible with wind-driven deepening of the mixed layer, Ekman pumping, or speeding-up of the wind-driven upwelling, such that more carbon-rich waters are mixed in from below during stronger winds.

In contrast, higher wind speeds tend to be associated with more negative ocean-internal DIC fluxes (i.e. weaker sources or stronger sinks) at the western extratropical fringes of all ocean basins (blue areas). In these regions of mode water formation, higher wind speeds lead to more subduction of anthropogenic CO₂ away from the surface into the ocean interior.

3.2.4 Additional variability in the sea–air flux

We note again that the sensitivities discussed here are those of the ocean-internal DIC sources and sinks f_int (Fig. 1 bottom or Fig. 3 left). The sea–air CO₂ fluxes (Fig. 1 top or Fig. 3 centre) contain additional variability also driven by interannual variations in SST (e.g. via the changes in CO₂ solubility and chemical equilibrium) or in wind speed (via the gas transfer velocity of gas exchange). As this additional variability is already generated by the parameterizations contained in our algorithm (Sect. 2.1.1), these processes are, within uncertainties, not reflected in the sensitivities against SST or u² again.

Even though the sea–air CO₂ flux is the quantity most directly relevant to the atmospheric CO₂ budget and its consequences for global climate, this additional variability partly disguises the variability caused by ocean-internal processes as those discussed above. This also means that the ocean-internal DIC sources and sinks f_int are potentially easier to be related to environmental variables than the sea–air CO₂ flux or the pCO₂ field traditionally chosen as a target variable of linear or non-linear regressions because it is a directly observed quantity.

3.3 How much predictive skill does the multi-linear regression have?

The results of the multi-linear regression are only meaningful if the regression actually possesses some predictive skill to bridge unconstrained periods. Only then can they be considered to represent generalizing relationships.

In order to test this, we performed runs with artificial data gaps of 5 years length (Sect. 2.3). Figure 5 illustrates this using runs discarding all pCO₂ data during 1995–1999. For context, we first consider the explicitly interannual mapping (E), which draws all information about year-to-year variations from the data and therefore does not have any predictive skill. Indeed, it essentially defaults to the prior (having upside-down El Niño response as it misses any variations related to the ocean-internal sources and sinks) during the data gap (Fig. 5a), except for a shift in long-term mean (see Sect. 2.1.2, last paragraph, for explanation). In contrast, the multi-linear regression (Fig. 5b) almost completely reconstructs the 1995–1999 flux variations based on the relationships between the ocean-internal DIC flux and the driving variables learned on the basis of the remaining data outside 1995–1999.

Figure 5 Interannual variations in the sea–air CO₂ flux in the tropical Pacific estimated by the explicitly interannual mapping (a), the multi-linear regression (b), and the hybrid mapping (c), using either all pCO₂ data (base cases, colour) or all data but the ones during 1995–1999 (black).

As demonstrated by Fig. S1 in the Supplement, this predictive skill generally holds for all parts of the ocean and other 5-year data gaps. This means that no particular pCO₂ data point is causing features in the variability and the estimated sensitivities (Sect. 3.2 above) by its own.

3.4 Sea–air CO₂ flux variations estimated by the hybrid mapping

After presenting the interannual sensitivities from the multi-linear regression, we now turn to interannual flux variations as estimated by the hybrid mapping involving an additional interannually varying correction (Sect. 2.1.5). Figure 6 (blue) shows its estimated interannual (i.e. slower-than-seasonal) variations in the sea–air CO₂ flux, subdividing the ocean into basins and latitude bands.

Figure 6 Yearly sea–air CO₂ flux as estimated from pCO₂ data by the explicitly interannual mapping (green; discarded before 1985, when the data constraint is very weak), the multi-linear regression (orange), and the hybrid mapping (blue). Fluxes have been integrated over a set of regions subdividing the ocean into basins (left to right) and latitude bands (top to bottom). For better temporal orientation across panels, the grey vertical background stripes indicate the positive phases of El Niño–Southern Oscillation according to the multivariate El Niño index (MEI) by Wolter and Timlin (1993).

The most prominent feature of interannual variability is the secular trend towards more CO₂ uptake in all ocean regions. Considering variations around this secular trend, the tropical Pacific is the region providing the largest contribution to total ocean variability (compare Le Quéré et al., 2000) on both a decadal timescale and a year-to-year timescale. The year-to-year variations are strongly tied to El Niño as indicated by the background stripes (Feely et al., 1999). When considering trends within individual decades, the decadal increase in the CO₂ sink slowed down in the 1990s and early 2000s and accelerated again afterwards (Landschützer et al., 2016; DeVries et al., 2019), even though it may be questioned whether such trends over chosen 10-year periods truly represent decadal variations rather than apparent trends arising from high-amplitude anomalies on the faster year-to-year timescale.

3.5 How do the year-to-year sea–air CO₂ flux variations estimated by regression and hybrid mapping compare with each other?

In addition to the variations in the sea–air CO₂ flux estimated by the hybrid mapping (blue), Fig. 6 also shows those estimated by the multi-linear regression (Sect. 2.1.4, orange) and the explicitly interannual pCO₂ mapping (Sect. 2.1.3, green). From the late 1980s onwards, when progressively more pCO₂ data are available to constrain interannual variations explicitly, the hybrid mapping (blue) shows some corrections over the multi-linear regression (orange). For the large El Niño-related variability in the tropical Pacific, these corrections are generally small compared to the estimated variations themselves. This indicates that the multi-linear regression already captures a notable fraction of the year-to-year flux variations in this region, even though it underestimates the size of most of these anomalies (the interannual standard deviation between 1985 and 2019 from the multi-linear regression is only about 82 % of that from the hybrid mapping in the tropical Pacific). Figure 7 (dots) confirms that the hybrid mapping fits the pCO₂ data closely (the blue dots are located right under the black dots), while the multi-linear regression (orange dots) also follows the variability in the data (black dots) but does not match them as closely as the hybrid mapping.

Figure 7 Estimated pCO₂ in the tropical Pacific, averaged spatially and over calendar years. The coloured lines give full regional averages from the explicitly interannual mapping (green), the multi-linear regression (orange), and the hybrid mapping (blue). The coloured dots are from the same estimates but averaged only over the pixels and time steps covered by pCO₂ data in the respective year. The smaller black dots give the corresponding averages over the data. We note that the green, blue, and black dots are not visible individually because they are almost exactly located on top of each other, indicating that the model–data residuals of the explicitly interannual and hybrid mappings are very small. The differences between dots and lines reflect the bias of the incompletely sampled average compared to the full regional average, which the mapping algorithm is trying to address.

In the intermediate and high latitudes (top and bottom panels of Fig. 6), in contrast, the multi-linear regression (orange) does not pick up most of the year-to-year anomalies. This may indicate that the set of explanatory variables used in the regression misses essential modes of variability there. However, some of the variations estimated with explicitly interannual degrees of freedom (green and blue) may also be spurious effects from the temporally very uneven data coverage.

Although the hybrid mapping (blue) has the same interannual degrees of freedom (i.e. the same flexibility) as the explicitly interannual mapping, it does not always bring the fluxes back to the explicitly interannual result (green), especially in the region south of the tropical Pacific (Fig. 6). Since the two estimates are actually very close to each other where data exist (as illustrated in Fig. 7; the green dots are essentially invisible under the co-located blue and black dots, despite the differences between the green and blue lines), the differences in areal averages as in Fig. 6 reflect differences in data-void areas and periods being filled by the mappings. However, while the explicitly interannual mapping falls back to the prior not constrained by pCO₂ data, the hybrid mapping falls back to the multi-linear regression, which is at least indirectly constrained via the statistical relationships between the ocean-internal DIC flux and the chosen explanatory variables (Sect. 3.3 above). This may also prevent some undue spatial extrapolation from the tropical Pacific into unconstrained areas by the explicitly interannual scheme. Thus, we expect the hybrid mapping (blue) to be more realistic than the multi-linear regression and the explicitly interannual mapping in terms of their detailed interannual anomalies.

In view of applying the multi-linear regression as a prior of the hybrid mapping, its predictive skill (Sect. 3.3) is only meaningful to the extent that it is actually able to explain all signals in the data. For example, since the regression underestimates the year-to-year anomalies in the tropical Pacific compared to the explicitly interannual estimate as discussed above, it will fill data gaps with somewhat too small an amplitude (Fig. 5c). This indicates that the variability extrapolated into the earlier decades without data will likely be underestimated, too, even though this is still a clear qualitative improvement compared to the explicitly interannual mapping (Fig. 5a).

3.6 What can the pCO₂ mappings say about the secular flux trend?

In light of climate change, quantitative information about the secular flux trend is relevant. Unfortunately, as discussed in more detail in the appendix (Sect. A2), the secular trend in our mapping results is mostly determined through the prior derived from the OCIM estimate based on ocean interior data (DeVries, 2022). Due to the lack of pCO₂ data in the early decades, the mapping does not add credible information about the secular trend. Due to some slight inconsistency in the use of the prior, the secular trend is even slightly overestimated (Sect. A2); this remains to be addressed in future versions of our flux product.

4 Discussion: robustness of the multi-linear regression

4.1 How meaningful is the multi-linear regression in terms of biogeochemical processes?

Statistical inferences by multi-linear regression are at the risk of overfitting, i.e. adjustment of coefficients to follow minor signals in the data, or even noise. If that was the case, the estimated sensitivities would not reflect underlying biogeochemical processes. Various findings indicate however that the results of the presented multi-linear regression do reflect actual signals:

• The patterns of the sensitivities (Fig. 4), at least those with respect to $\mathrm{dSST}/\mathrm{d}t$ and u², are quite systematic spatially and in many respects interpretable (Sect. 3.2). This is especially true in the tropical Pacific, but also throughout the entire ocean.

• Test regression runs only using one of the explanatory variables (RegrOnlySST, RegrOnlydSSTdt, RegrOnlyU2) yield sensitivities very similar to the base case using all explanatory variables (Supplement Fig. S2). This indicates that the regression terms are essentially mutually independent, such that each explanatory variable picks up a more or less unique portion of the signals contained in the data.

• The regression possesses predictive skill (Sect. 3.3 above), which also means that it does not depend on any particular portion of the pCO₂ data or the explanatory fields alone.

• The estimates are relatively robust against alternative data sets for the explanatory variables (cases RegrSSTNOAA and RegrU2NCEP in Supplement Fig. S4). This corroborates the notion that the regression is likely not dominated by any particular feature in these fields.

• The regression hardly responds to a 4-fold change in the regularization strength (case RegrLoose in Supplement Fig. S4). In the overfitting regime, one would expect a substantial dampening effect when the regularization is stronger.

• The regression results are quite robust against further changes in the set-up (Supplement Figs. S4, S5, and S6).

The relatively small set of explanatory variables used here (Sect. 2.1.4) is certainly helpful to avoid overfitting as it means a relatively small number of degrees of freedom (cf. Thacker, 2012). Also the use of temporally constant sensitivity coefficients helps to keep the number of degrees of freedom sufficiently low. For example, in test runs with seasonally resolved sensitivity coefficients, the data could be fitted more closely but the predictive skill deteriorated (not shown).

Clearly, for any given spatial area, the presence of pCO₂ data over a sufficient variety of environmental conditions is a prerequisite to estimate meaningful sensitivity coefficients. The “reduction in uncertainty” diagnostic for interannual variations given in Rödenbeck et al. (2014) provides at least a rough indication. A reduction-in-uncertainty diagnostic could also be performed for the sensitivity coefficients directly, which however remains for follow-on work.

4.2 Which fraction of the year-to-year variability can be captured by the multi-linear regression?

As seen in Fig. 6 and quantified explicitly in Fig. 8 (top), the amplitude of year-to-year variability in the global sea–air CO₂ flux estimated by our multi-linear regression (orange) is lower than that estimated by the hybrid mapping possessing the degrees of freedom to follow any interannual signals (blue). This indicates that the pCO₂ data also contain signals of year-to-year variability that cannot be represented in terms of the variations contained in the set of explanatory variables used in the regression. Possibly, however, the hybrid mapping may also exaggerate the amplitude of signals by spreading them over too large an area in data-poor parts of the ocean.

Figure 8 Amplitudes of variability in the sea–air CO₂ flux on year-to-year timescales around its secular trend, from the multi-linear regression (orange group of bars; solid: base case; hashed: uncertainty cases) and the hybrid mapping (blue), compared to other pCO₂ mapping methods (salmon; CMEMS v2020, Denvil-Sommer et al., 2019; CSIR-ML6 v2020, Gregor et al., 2019; JMA-MLR v2020, Iida et al., 2020; and MPI-SOMFFN v2020, Landschützer et al., 2013) as well as the ocean biogeochemical process models collated in Friedlingstein et al. (2020) (mint green). The amplitudes are represented by temporal standard deviations of detrended yearly fluxes over the 1990–2019 period. The top panel gives the global flux, the bottom panel the tropical Pacific.

The situation is different in the tropical Pacific (Fig. 8, bottom). Here, the multi-linear regression (orange) already captures a large part of the variability found in the hybrid mapping (blue). This indicates that our explanatory variables are reasonably suited to represent the ENSO-related variability dominating in this region.

To elucidate the ability of the multi-linear regression to capture year-to-year anomalies, we compare it with other pCO₂ mappings based on linear or non-linear regressions of pCO₂ (itself) against various sets of explanatory variables (Landschützer et al., 2013; Iida et al., 2020; Denvil-Sommer et al., 2019; Gregor et al., 2019). Globally (Fig. 8, top), the variability obtained by the other pCO₂ mappings (salmon) is larger than that from our multi-linear regression (orange). Closer inspection (not shown) reveals that these larger amplitudes mostly reflect variability on multi-year (decade-to-decade) timescales occurring coherently in both northern and southern extratropics, while the multi-linear regression does not involve such globally correlated contributions. Accordingly, when splitting up the global flux into regional contributions, the amplitudes from the other pCO₂ mappings and our multi-linear regression are quite comparable. For example, in the tropical Pacific (Fig. 8, bottom) our regression yields year-to-year variability larger than any of the other pCO₂ mappings considered. Based on reconstructions of model-based pseudo-data, Gloege et al. (2021) found for one of the other methods included in Fig. 8 that the amplitude of Southern Ocean decadal variability was overestimated by 15 % to 58 %.

Could alternative or additional explanatory variables help to capture a larger fraction of variability by the multi-linear regression?

• As the explanatory variables of the base case are all physical variables, we tested using chlorophyll a concentration as a biological variable (run RegrAddChl_98r19; Supplement Fig. S7). A practical problem with chlorophyll a is that data sets are only available for the most recent years (from 1998); therefore it is not used in our base case. The test suggests, however, that chlorophyll is not actually adding much information about the year-to-year variations in the sea–air CO₂ flux beyond what is already provided by the explanatory variables of the base case (SST, $\mathrm{dSST}/\mathrm{d}t$, and u²). A reason may be that chlorophyll variability is already covered in the other variables as nutrients are also a function of upwelling, stratification, etc. It is also important to keep in mind that chlorophyll concentration is not directly observed but only indirectly inferred from optical properties of the seawater. Due to that, part of the variability in the chlorophyll data may originate from processes unrelated to the carbonate system, which makes it less helpful as a predictor in the regression considered here.

• Conceivably, more general non-linear relationships between pCO₂ and the explanatory variables may allow the capturing of signals not represented in linear relationships as used in our base case. Uncertainty cases involving additional regression terms proportional to ($\mathrm{dSST}/\mathrm{d}t$)² (run RegrAdddSSTdt2) and (u²)² (run RegrAddU4), respectively, only marginally increase year-to-year variability (within the narrow band in Fig. S4). Also from the set of other pCO₂ mappings (salmon in Fig. 8), there is no indication that the non-linear regressions (CMEMS-FFNN, CSIR-ML6, and MPI-SOMFFN) would generally capture more variability than the linear ones (JMA-MLR and ours). We conclude that non-linearities in the pCO₂ relationships are not essential for explaining year-to-year anomalies in the pCO₂ field on a regional scale.

• Using heat flux as an explanatory variable instead of $\mathrm{dSST}/\mathrm{d}t$ (RegrHeat_85r09) deteriorates the ability of the multi-linear regression to reproduce ENSO-related variability (Supplement Fig. S8).

• Replacing u² by the wind stress curl (RegrCurl_88r18) does not change the flux IAV much (Supplement Fig. S10). A further alternative explanatory variable may be “Ekman pumping”, which however diverges at the Equator and was not tested.

A common methodological feature of all present-day regression-based pCO₂ mappings including ours is that the carbon variables are only related to the concurrent values of the explanatory variables, disregarding any dependence on past values of the explanatory variables possible due to memory effects. This might be a serious limitation, but allowing for memory effects is not straightforward. For example, regression terms with lagged explanatory variables would only allow discrete lag times, and using an extensive spectrum of lag times would possibly exceed the number of well-determined degrees of freedom. Theoretically, fitting comprehensive process models to the pCO₂ data would include emerging memory effects, but this faces various conceptual and computational challenges (see a recent application of a low-dimensional Green's function approach by Carroll et al., 2020). (Note that the amplitudes simulated by the hindcast ocean biogeochemical models included in Fig. 8 are roughly similar to those from our multi-linear regression and smaller than those from the hybrid scheme.) We notice that our algorithm involves some elements that do represent history effects (the budget equation Eq. A18 in Rödenbeck et al., 2013, accumulating past f_int contributions; the seasonal “history flux” Eq. A20 in Rödenbeck et al., 2013; and the use of both SST and $\mathrm{dSST}/\mathrm{d}t$ as explanatory variables; see Sect. 3.1 above). However, if memory effects are important, they are evidently not yet adequately captured by those elements.

4.3 To which extent do the sensitivities γ_i depend on the timescale?

In our formulation of the regression (Eq. 3), the sensitivities γ_i are applied to the fields V_i of the explanatory variables including all their variations on year-to-year, decadal, and secular timescales. Conceivably, however, the relationships between f_int and the explanatory variables may differ for year-to-year, decadal, or secular variations. In ocean areas where the data period is long enough to possibly constrain decadal timescales directly, the estimates may therefore reflect some mixture of timescales, which would be hard to interpret.

We assessed this by the uncertainty case RegrNoDecad, where any decadal variability (including any secular trend) has been removed from the three explanatory variables. As this case can only pick up year-to-year signals to constrain the sensitivities, any changes compared to the base case may indicate such potential timescale conflicts. In most regions, this is not evident (Supplement Fig. S6). Exceptions are the southern Pacific and the tropical Indian (for the wind-speed sensitivity ${\mathit{\gamma }}_{u^{\mathrm{2}}}$) and the western tropical Pacific (for the SST sensitivity γ_SST). As the explanatory variable $\mathrm{dSST}/\mathrm{d}t$, which dominates the large tropical variability, does not have much secular trend, it is not prone to timescale dependence anyway.

An alternative way to assess the impact of secular trends in the explanatory variables is the uncertainty case RegrAddpaCO2 having an additional regression term proportional to decadally smoothed atmospheric CO₂ (p_aCO₂). As p_aCO₂ is rising steadily over the calculation period, this run is able to adjust the secular trend independently of the trends in SST, $\mathrm{dSST}/\mathrm{d}t$, or u², thus breaking any potential timescale conflicts. Indeed, the sensitivities estimated by RegrAddpaCO2 (not shown) are similar to those from the base case as well, and any differences from the base case are similar to those of RegrNoDecad.

We note that in ocean areas with data periods of a few years only, a possible timescale dependence will not affect the sensitivities themselves, but it may still affect secular trends in the fluxes if sensitivities estimated for year-to-year variations are applied to secular trends in the explanatory variable. We do not have a means to detect whether this is the case.

4.4 Spurious effects from uncertainties in the parameterizations

Errors in the sea–air CO₂ flux resulting from deficiencies in our chosen parameterizations of solubility and gas exchange lead to compensating spurious contributions to f_int because it is the sum of both fluxes which changes the mixed-layer carbon content in our budget equation (see Fig. 1 or  Rödenbeck et al., 2013). This will then also lead to spurious contributions to the estimated sensitivities γ_i. For example, spurious u² sensitivity may arise if the wind speed dependence of our gas exchange parameterization is not strong enough such that it is reinforced by additional changes in the ocean-internal carbon flux (or vice versa).

Luckily, the interannual variability in the sea–air CO₂ flux is much smaller than that of f_int due to the buffer effect (see Fig. 1). Therefore, in relative terms, the error in the sea–air CO₂ flux translates into a much smaller error in f_int and in the sensitivities γ_i.

4.5 Spurious effects from missing interannual alkalinity variations

The estimated ocean-internal DIC flux f_int – and thus the estimated sensitivities γ_i in the regression – contains some spurious contributions to compensate any errors in our representation of carbonate chemistry because the pCO₂ data constrain the pCO₂ field rather than the DIC field (Fig. 1). Even though we represent the carbonate chemistry – up to the linearization – by exact equations (Sect. 2.1.6), some error arises because we only use a seasonal alkalinity climatology, while alkalinity also varies interannually due to (1) changing degrees of dilution due to freshwater fluxes (evaporation, precipitation, ice formation, and ice melt) as well as (2) mixing-in of alkalinity-rich deep waters and possibly biological influences.

1. Freshwater fluxes dilute not only alkalinity but also DIC, in equal proportions. At the same time, the sensitivities of pCO₂ to changes in alkalinity and DIC are almost equal in absolute value but of opposite sign (Sarmiento and Gruber, 2006). Therefore, the total effect of freshwater fluxes on pCO₂ is small compared to that on alkalinity and DIC, respectively. Therefore, as we neglect both the freshwater contributions to f_int and the freshwater-related alkalinity variations, the combined error in pCO₂ should be small.

2. Alkalinity variations related to mixing from below are linked to DIC variations as well because deep waters are rich in both DIC and alkalinity, compared to the mixed layer. In contrast to the freshwater effects, however, the regression terms γ_iV_i in Eq. (3) do contain mixing contributions to f_int, such that the absence of the corresponding alkalinity variations does affect our pCO₂ field being matched to the data. On the seasonal timescale (where there is no problem anyway as we are using a monthly alkalinity climatology), alkalinity variations in the tropical and subtropical oceans are dominated by freshwater effects; only at higher latitudes are alkalinity variations increasingly affected by mixing (Lee et al., 2006). For the interannual timescales relevant here, the relative role of mixing is unclear. A better understanding – and hopefully solution – of this problem remains for further work.

We note that the spurious compensatory contributions to f_int do not affect the pCO₂ field being constrained by the observations. Thus, they essentially do not affect the estimated sea–air CO₂ fluxes either.

4.6 Further sources of uncertainty

The interannual variations estimated before the pCO₂ data period (i.e. before about 1990) represent extrapolations based on the estimated sensitivities γ_i and the variations in the explanatory variables. As the data sets used for the explanatory variables are generally based on fewer and more uncertain observations in the earlier decades, the uncertainty in our results is expected to be larger in the earlier decades as well. A meaningful quantification of this uncertainty is deemed impossible.

5 Conclusions

In this study, we considered the interannual variability in the sea–air CO₂ flux over the 1957–2020 period, constrained by the pCO₂ measurements from the SOCATv2021 database (Bakker et al., 2016). Extending the pCO₂ mapping scheme of Rödenbeck et al. (2013, 2014), we employed (1) a multi-linear regression against interannual anomalies of sea surface temperature (SST), the temporal changes in SST ($\mathrm{dSST}/\mathrm{d}t$), and squared wind speed (u²), yielding maps of interannual sensitivities, and (2) a subsequent explicitly interannual additive correction, yielding a “hybrid” estimate of spatio-temporal variations in the contemporary sea–air CO₂ flux (formal resolution 2.5^∘ longitude × 2^∘ latitude × 1 d).

• According to our multi-linear regression, interannual variability in the tropical Pacific is dominated by a positive correlation of ocean-internal DIC fluxes to $\mathrm{dSST}/\mathrm{d}t$, as arising from variations in the upwelling of colder and more carbon-rich waters into the mixed layer.

• In the eastern upwelling zones as well as in circumpolar bands in the high latitudes of both hemispheres, we find a positive sensitivity to wind speed, compatible with the entrainment of carbon-rich water during wind-driven deepening of the mixed layer. To the extent that this sensitivity inferred from year-to-year variations also applies to secular trends, the wind trend in the Southern Ocean (south of 45^∘ S) implies a wind-related reduction in the flux trend by about 17 % to 42 % (weaker increase in sink).

• As a pCO₂ mapping method, the hybrid mapping combines (a) the ability of regression to bridge data gaps and extrapolate into the early decades without much pCO₂ data constraint and (b) the ability of an auto-regressive interpolation to follow signals even if not represented in the chosen set of explanatory variables. This way, at least the large contributions of the tropical Pacific to the global year-to-year variability in the oceanic CO₂ exchange can be extrapolated over the entire 1957–2020 period, even though the extrapolated variability prior to about 1985 is probably underestimated.

Appendix A: The global ocean carbon sink estimated by the hybrid mapping

Here we discuss the global total of the sea–air CO₂ flux as estimated by the hybrid mapping and compare it to various literature estimates. In order to allow a quantitative comparison, we focus on specific features, namely the mean flux (Sect. A1) and the secular flux trend (Sect. A2).

A1 The mean sink (1994–2007)

Figure A1 shows the contemporary global sea–air CO₂ flux estimated by the hybrid mapping (solid blue bar) averaged over the 1994–2007 period. According to the set of uncertainty cases shown (hashed blue bars), the uncertainty in the mean flux from the hybrid mapping is dominated by the uncertainty in gas exchange (cases GasexLow, GasexHigh, GasexU1, and GasexU3; diagonally hashed bars), while all other uncertainty cases do not affect the mean sink estimate very much.

Figure A1 Mean global sea–air CO₂ flux over 1994–2007 from the hybrid mapping (blue group of bars; solid: base case; hashed: uncertainty cases), compared to other pCO₂ mapping methods (salmon; CMEMS v2020, Denvil-Sommer et al., 2019; CSIR-ML6 v2020, Gregor et al., 2019; JMA-MLR v2020, Iida et al., 2020; and MPI-SOMFFN v2020, Landschützer et al., 2013), and the ocean biogeochemical process models collated in Friedlingstein et al. (2020) (not including a river-induced sea–air flux (see Table A1); mint green). To the left of the hybrid mapping, we also give OCIM (DeVries, 2022; total contemporary flux) as well as intermediate results from this study. The long-dashed horizontal line indicates the estimate of $-\mathrm{2.6}\pm \mathrm{0.4}\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}$ from ocean interior data by Gruber et al. (2019) (anthropogenic carbon only) and the dotted line the same estimate shifted by 0.62 PgC yr⁻¹ (average of Jacobson et al., 2007 and Resplandy et al., 2018) as an assumed contribution from outgassing of terrestrial carbon transported to the ocean by rivers. Positive fluxes denote oceanic CO₂ outgassing into the atmosphere; negative fluxes denote CO₂ sinks into the ocean.

The spread between the flux estimates from other pCO₂ mapping methods (group of salmon bars) together with the base case of our hybrid mapping (solid blue bar) only indicates uncertainties due to the mapping algorithms as all the estimates use the same global scaling of the gas transfer velocity from Naegler (2009). Notably, this spread does not exceed the differences due to scaling sea–air gas exchange within the uncertainty range of Naegler (2009) (cases GasexLow, GasexHigh).

The comparison between the results of the hybrid mapping and further literature values is hampered as pCO₂ mappings are estimating the total contemporary net CO₂ flux (F_net) through the sea–air interface, while other methods may only include certain components of it. Adopting the notation by Hauck et al. (2020), Table A1 gives the six components of F_net and their respective inclusion in the literature estimates considered here (note that the terms “anthropogenic” or “contemporary” are also defined differently in part of the literature).

Gruber et al. (2019)

Table A1 The components of the contemporary net sea–air CO₂ flux ($F_{\mathrm{net}}=F_{\mathrm{ant},\mathrm{ss}}+F_{\mathrm{ant},\mathrm{ns}}+F_{\mathrm{nat},\mathrm{ss}}+F_{\mathrm{nat},\mathrm{ns}}+F_{\mathrm{riv},\mathrm{ss}}+F_{\mathrm{riv},\mathrm{ns}}$) according to Hauck et al. (2020) and whether or not they are included in the individual estimates shown in Figs. 8, A1, and A2.

¹ For a mean ocean circulation over the industrial era; i.e. the variations are due to temperature and gas transfer velocity variations only. ² Mean 1994–2007 flux only.

From the increase in the anthropogenic carbon inventory in the ocean between the extensive ocean surveys in 1994 and 2007, Gruber et al. (2019) estimate an anthropogenic CO₂ uptake of $F_{\mathrm{ant},\mathrm{ss}}+F_{\mathrm{ant},\mathrm{ns}}=-\mathrm{2.6}\pm \mathrm{0.3}\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}$ over the interjacent period, shown in Fig. A1 as a long-dashed line. This estimate conceptually differs from the hybrid mapping by the river-induced flux $F_{\mathrm{riv},\mathrm{ss}}+F_{\mathrm{riv},\mathrm{ns}}$ and the non-steady-state modifications F_nat,ns to the natural sea–air fluxes, while F_nat,ss is zero at the global scale (Table A1). The river-induced flux is very uncertain, with literature estimates ranging between $\mathrm{0.45}\pm \mathrm{0.18}\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}$ (Jacobson et al., 2007) and $\mathrm{0.78}\pm \mathrm{0.41}\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}$ (Resplandy et al., 2018), though the real uncertainty may be even larger. If the Gruber et al. (2019) estimate is shifted by a mid-range river-induced value of 0.62 PgC yr⁻¹ (resulting in the dotted line), the base case value from the hybrid estimate is matched more closely. Nevertheless, given the uncertainty ranges of gas exchange, river-induced outgassing, and the Gruber et al. (2019) estimate, we cannot draw any conclusions from the remaining difference.

The CO₂ flux difference between the hybrid estimate and the dotted line in Fig. A1 may also contain a contribution from systematic differences between pCO₂ in the bulk ocean water (as typically measured at a few metres depth) and pCO₂ at the diffusive surface layer (as relevant for gas exchange), arising due to systematic differences in water temperature and salinity (Woolf et al., 2016). Further, the cooler ocean skin temperature translates the atmospheric pCO₂ to a different equilibrium DIC concentration than that implicitly calculated based on bulk temperature (Robertson and Watson, 1992). Watson et al. (2020) estimated that the sum of these two effects would shift pCO₂-based estimates of the mean global CO₂ flux by −0.8 to $-\mathrm{0.9}\mathrm{PgC}\mathrm{yr}{}^{-\mathrm{1}}$ (stronger sink). So far, however, it is unclear how well the water temperature at the relevant vertical positions can actually be determined (an important source of uncertainty not included in the range of Watson et al., 2020) and how it varies in space and time. In any case, we note that our study mainly considers the variability in the flux, for which the effect of a time-constant correction as in Watson et al. (2020) would cancel out.

Figure A1 further shows the global fluxes simulated by a set of global ocean biogeochemical models (GOBMs) collated in the annual global carbon budget (Friedlingstein et al., 2020; mint green). Like OCIM or Gruber et al. (2019), the GOBMs' results do not include the river-induced flux component, but they do conceptually include the non-steady-state modification F_nat,ns of carbon uptake and natural carbon cycling (Table A1). The range of results covered by the GOBMs slightly exceeds the range of the hybrid estimates due to the gas exchange uncertainty. The medians of the GOBM ensemble and the pCO₂ mapping ensemble differ by less than the mid-range river-induced value of 0.62 PgC yr⁻¹.

A2 The secular sink trend (1960–2019)

Regarding the 1960–2019 secular sink trend, our estimate from the hybrid mapping (1) is not able to add much independent information and (2) even slightly overestimates this trend relative to OCIM used in the prior:

1. According to Fig. A2a, the 1960–2019 trend from the base case (solid blue bar) is quite similar to that of the base case prior (open grey bar). Among the uncertainty cases (hashed blue bars), the largest deviations are seen when mixed-layer depth is changed (MLDq2 and MLDx2); these deviations are in fact mostly inherited from their respective priors as well (not shown).

2. Figure A2a further reveals that the prior (open grey bar) has a slightly steeper trend than the OCIM estimate (magenta) even though the prior has been derived from OCIM (Sect. 2.1.3). This discrepancy arises because we are using OCIM's sea–air fluxes as a prior of the ocean-internal flux f_int even though these two quantities differ by the carbon accumulation in the mixed layer. Since the carbon accumulation accelerates (following the accelerating increase in atmospheric pCO₂), this leads to a difference not only in mean flux (Fig. A1) but also in trend. Due to the lack of information to correct the 1960–2019 secular trend from the pCO₂ data as discussed under (1), this issue leads to an overestimation of the trend in the hybrid estimate compared to OCIM. Most GOBMs (mint green) simulate an even flatter 1960–2019 trend than OCIM.

Figure A2 Secular linear trend of the global sea–air CO₂ flux over 1960–2019 (a) and over the more recent 1990–2019 period better constrained by pCO₂ data (b), from the same estimates as in Fig. A1. The error bars give the formal error in the slope calculated from 5-year flux averages; i.e. they reflect the uncertainty due to the pentadal variability around the linear trend, thereby roughly taking into account the serial correlations of the flux, e.g. on the El Niño timescale.

Looking at the linear trend over the better-constrained, more recent period 1990–2019 (Fig. A2c), the estimate from the hybrid mapping becomes more independent from the prior. The pCO₂-based hybrid estimates tend to show steeper trends than both OCIM and the GOBM simulations. Most other pCO₂ mappings (salmon) estimate the trend to be even more negative than the hybrid mapping. However, given the substantial pentad-to-pentad variations in the global flux (as reflected in the error bars), it is not fully clear how well defined the trend over the 1990–2019 period actually is.

The level of constraint in the trend over the different periods is corroborated by the “zero-prior” mapping not using the secular trend from OCIM as a prior (Fig. A3). Even though the zero-prior explicitly interannual mapping (violet) and the explicitly interannual mapping (green) start from priors with very different secular trends (shown in dark and light grey, respectively), their estimated multi-decadal trends during the recent decades are still very close. In well-constrained regions like the tropical Pacific (bottom panel) they are practically identical, while some deviations occur in poorly constrained regions such as the Indian Ocean, adding up to the slight deviations in the global total flux (top). Only in the early decades where there are hardly any pCO₂ data to constrain the estimates do the two mappings stick to the differing multi-decadal trends (and also to the year-to-year variations) of their respective priors.

Figure A3 Interannually filtered sea–air fluxes as in Fig. 6 (two example panels only) estimated by the zero-prior explicitly interannual pCO₂ mapping (violet) and the explicitly interannual pCO₂ mapping (green) as well as their respective priors (dark and light grey, respectively; note that the designation “zero-prior” refers to $f_{\mathrm{int}}^{\mathrm{ZE},\mathrm{pri}}=\mathrm{0}$, while the a priori sea–air fluxes shown here are in fact non-zero owing to the rise in atmospheric CO₂ and the variations in SST). The line width roughly distinguishes the early period with insufficient data constraint (thin) and the more recent period with better constraint (thick). Vertical background stripes as in Fig. 6.

As the better-constrained trend over the recent decades (after about 1992) is essentially the same as that in the prior of the explicitly interannual mapping, the flat multi-decadal trend of the zero-prior mapping in the early decades is very unlikely to be true. This illustrates that a prior with the correct secular trend (such as the OCIM result used here) is indeed needed to extrapolate the ocean CO₂ sink into the data-poor first decades of our extended period of interest 1957–2020.