Decadal Trends in the Oceanic Storage of Anthropogenic Carbon From 1994 to 2014

The oceanic uptake and resulting storage of the anthropogenic CO2 (Cant) that humans have emitted into the atmosphere moderates climate change. Yet our knowledge about how this uptake and storage has progressed in time remained limited. Here, we determine decadal trends in the storage of Cant by applying the eMLR(C*) regression method to ocean interior observations collected repeatedly since the 1990s. We find that the global ocean storage of Cant grew from 1994 to 2004 by 29 ± 3 Pg C dec−1 and from 2004 to 2014 by 27 ± 3 Pg C dec−1 (±1σ). The storage change in the second decade is about 15 ± 11% lower than one would expect from the first decade and assuming proportional increase with atmospheric CO2. We attribute this reduction in sensitivity to a decrease of the ocean buffer capacity and changes in ocean circulation. In the Atlantic Ocean, the maximum storage rate shifted from the Northern to the Southern Hemisphere, plausibly caused by a weaker formation rate of North Atlantic Deep Waters and an intensified ventilation of mode and intermediate waters in the Southern Hemisphere. Our estimates of the Cant accumulation differ from cumulative net air‐sea flux estimates by several Pg C dec−1, suggesting a substantial and variable, but uncertain net loss of natural carbon from the ocean. Our findings indicate a considerable vulnerability of the ocean carbon sink to climate variability and change.

period to 4.4 ± 0.9 Pg C dec −1 from 2003 to 2014 (Wanninkhof et al., 2010;Woosley et al., 2016). While Woosley et al. (2016) reported a rather steady uptake behavior in the South Atlantic, Gao et al. (2022) found that the rates of C ant storage accelerated from the 1990s to the 2000s. In all regional studies, the temporal variability of the C ant storage was attributed to changing ventilation patterns of the upper ocean. However, differences in time periods and statistical methods applied in the regional studies limit their synoptic assessment and prevent combining them into a global reconstruction of the oceanic increase in C ant storage since 2007.
Such an extension of the reconstruction of the global increase in C ant storage beyond 2007 is the key aim of this study. Our work profits from including ∼100,000 additional observations of dissolved inorganic carbon (DIC) and related biogeochemical variables collected over the 2010s, which were compiled, quality-controlled, and made available by GLODAP, the Global Ocean Data Analysis Project (Lauvset et al., 2021). By consistently determining the storage increase between 1994 and 2004, and between 2004 and 2014, and benefitting from the reconstructed storage of C ant for 1994, we can investigate for the first time the temporal evolution of the global increase in the oceanic storage of C ant . This permits us to address whether the ocean has maintained its vital sink function in a changing climate. Our global-scale reconstruction of the oceanic storage of C ant further serves as an important independent reference point for the ocean carbon sink estimates established by other means, especially in the context of the GCB (Friedlingstein et al., 2022) and the Intergovernmental Panel on Climate Change (IPCC) (Canadell et al., 2021).

Overview of the Approach
Our global-scale analysis of the changes in the content of C ant (ΔC ant ) is based on measurements of the dissolved inorganic carbon (DIC) content and related hydrographical and biogeochemical properties gathered from 1989 to 2020 and synthesized in the data product GLODAPv2.2021 (Lauvset et al., 2021;Olsen et al., 2016). This data product includes high-quality measurements from reoccupied sections for the purposes of diagnosing long term climate signals such as the accumulation of C ant . The majority of the data used in this study stem from the JGOFS/ WOCE global CO 2 survey conducted in the 1980s and 1990s Wallace, 1995), the repeat hydrography program GO-SHIP that began in 2003 and is now completing its second cycle (Sloyan et al., 2019;Talley et al., 2016), as well as a number of additional programs, including INDIGO, SAVE, TTO, JOIS, and GEOSECS (Key et al., 2004, and references therein). In addition to DIC, our analysis requires observations of salinity (S), temperature (T), total alkalinity (TA), oxygen (O 2 ), the apparent oxygen utilisation (AOU), silicate (Si(OH) 4 ), nitrate (NO 3 − ), and phosphate ). To extract the ΔC ant signal from these data, we use the eMLR(C*) method (Clement & Gruber, 2018;Gruber et al., 2019) with a few modifications (see details below and in Text S2 in Supporting Information S1).
Our application of the eMLR(C*) method employs the following steps: 1. The semi-conservative tracer C* (Gruber et al., 1996;Gruber & Sarmiento, 2002) is calculated from DIC as C* = DIC − 117 × PO 4 3− − 0.5 × (TA + 16 × PO 4 3− ). This conversion from DIC to C* removes a substantial part of the inorganic carbon variability that is due to ocean interior redistributions and biogeochemical transformations of natural carbon. 2. The observations are clustered in neutral density slabs (Text S2.1 in Supporting Information S1) and ocean regions ( Figure S1 in Supporting Information S1), and assigned to one of the three sampling periods 1989-1999, 2000-2009, or 2010-2020 (Figure 1). 3. Within each sampling period, the observed C* is adjusted to the respective reference year (t ref ) 1994, 2004, or 2014 assuming a transient steady state increase of C ant (Gammon et al., 1982). . 5. The 10 best common MLR models for two compared sampling periods are selected within each density slab and ocean region based on the summed root mean squared error (RMSE), after excluding MLRs with strong multicollinearity between the predictors (Text S2.2 in Supporting Information S1).

10.1029/2023AV000875
4 of 28 6. The decadal change in anthropogenic carbon is computed as the difference between the average C* distribution for each sampling period, that is, ΔC ant = C*(t ref,n+1 )-C*(t ref,n ), where the C* distributions are predicted ("mapped") for each t ref by applying the selected MLRs to a common set of predictor climatologies. 7. In surface waters, ΔC ant is predicted based on an transient equilibrium approach (McNeil et al., 2003), which assumes that the increase of the surface ocean pCO 2 follows that of the atmosphere closely. This approach aims to avoid biases introduced by the seasonal and interannual variability of C* and the predictor variables.
The most important differences we introduced relative to the methods described by Clement and Gruber (2018) and Gruber et al. (2019) are detailed in Section 2.3 and include a more thorough selection of the predictor variables and MLR models, a more robust and standardised quantification of uncertainties, and an assessment of structural reconstruction uncertainties through tests with synthetic data, which were performed in parallel to the analysis of the real-world observations. Our reconstruction of ΔC ant with the eMLR(C*) method involves a number of choices regarding the configuration of the method. In Sections 2.2 and 2.3 we describe the standard configuration that we use to derive the results we report as our best informed estimates. In Section 2.4, we describe a number of reasonable alternative configurations. We use the offsets between results obtained with these plausible alternatives and the standard configuration as a basis for determining the uncertainty.  Figure S1 in Supporting Information S1. The shaded background indicates the assigned sampling periods 1989-1999 (t ref,1 = 1994), 2000-2009 (t ref,2 = 2004), and 2010-2020 (t ref,3 = 2014). (b) Map of observations as used in our standard case for the three sampling periods. Cruises that fulfilled all flagging criteria are displayed in black, while cruises for which one parameter did not fulfill all criteria but were still included (see also Table S1 in Supporting Information S1) are highlighted in color.

Data
This study relies on ocean interior observations of DIC, TA, S, T, O 2 , Si(OH) 4 , NO 3 − , and PO 4 3collected in the global ocean from 1989 through 2020 ( Figure 1) and provided through GLODAPv2.2021 (Lauvset et al., 2021;Olsen et al., 2016). The apparent oxygen utilisation was calculated from oxygen, salinity, and conservative temperature (Graham & McDougall, 2013) according to the solubility from Weiss (1970) and used as an additional predictor variable. Observations were filtered based on the GLODAP flagging scheme, using only highest-quality data identified as those with a f-flag value of 2, which indicates acceptable, measured data according to a simplified version of the WOCE flagging scheme, and a qc-flag value of 1, which indicates adjusted or unadjusted data that have undergone GLODAP's full secondary quality control. With a few exceptions, only samples were included for which all required variables fulfill the strictest quality criteria. We deviated from this only for cases where our tests of the eMLR(C*) method with synthetic data (see Text S5 in Supporting Information S1) revealed that omitting observations because of only one missing variable increased the biases of the ΔC ant reconstructions. In these cases either data that did not fulfill the strictest quality criteria were included (7,609 samples, 3.7% of all samples), or missing data were filled (7,305 samples, 3.5%) using CANYON-B (Bittig et al., 2018) predictions (see Text S1.3 in Supporting Information S1). In total, we used 206'836 samples for our standard case reconstruction.
Although the GLODAP data have already undergone a secondary quality control and are-if required-adjusted to improve their internal consistency, we applied a number of additional adjustments to the DIC, TA, and phosphate observations. These data adjustments are supported by multiple independent lines of evidence and quantified based either on our decade-by-decade reanalysis of deep water crossovers originally determined by GLODAP, or on previously unaccounted offsets in the measurements of certified reference materials (CRM) for DIC and TA (see Text S1.2 in Supporting Information S1 for details; Johnson et al., 1998Johnson et al., , 2002Millero et al., 1998). This affected the majority of the DIC and TA measurements from the Indian Ocean in the 1990s (12843 samples, 6.2% of all samples), and the DIC, TA, and phosphate measurements in the North Pacific from the 2010s (35395 samples, 17.2% of all samples). Our proposed adjustments of the Indian Ocean data were formally accepted by GLODAP and applied in the release of GLODAPv2.2022 (Lauvset et al., 2022). In contrast, the proposed adjustments for the North Pacific are still under discussion by the GLODAP community, and while well justified and likely attributable to measurement inconsistencies (Fong & Dickson, 2019;Olsen et al., 2019;Sharp & Byrne, 2020), currently need to be viewed as preliminary. We thus discuss the relevance of the North Pacific adjustments for our results and conclusions in detail. The magnitude of these additional adjustments are generally small (<2 μmol kg −1 for DIC, <4 μmol kg −1 for TA, and <1% for PO 4 3− ) and below the adjustment limits normally considered by GLODAP. Still, these adjustments proved to be critical in our work, since an offset of 1 μmol kg −1 in DIC integrated over 3,000 m amounts to a column inventory offset of ∼3 mol m −2 , which is of similar magnitude as some of the decadal changes we are aiming to detect.
For the purpose of predicting the C* distributions, we used the objectively analyzed climatology for the 1981-2010 period for salinity and temperature from the World Ocean Atlas 2018 Zweng et al., 2019), in combination with O 2 , Si(OH) 4 , NO 3 − , and PO 4 3− from the global interior ocean mapped climatology based on GLODAPv2 . The climatological distribution of AOU was calculated in accordance with the observational data. For atmospheric CO 2 (CO 2,atm ) data we use the globally averaged marine surface annual mean data provided by NOAA/GML (Lan et al., 2022).

Standard Configuration of the eMLR(C*) Method
In the following, we describe the main changes of the eMLR(C*) method in comparison to the previous analysis by Gruber et al. (2019). Further minor configuration changes are presented in Text S2 in Supporting Information S1.

Temporal Clustering and C* Adjustment to Reference Year
For the temporal clustering of the data, we assigned each observation to one of the following three sampling periods (start and end years included): 1989-1999 ( ref,1 = 1994) 2000-2009 ( ref,2 = 2004) 2010-2020 ( ref,3 = 2014) with the assigned reference years (t ref ) given in parenthesis.
Based on these three sampling periods, we estimated the C ant storage changes between the reference years 1994-2004 and 2004-2014 (Figure 1). The ΔC ant estimates represent the changes over exactly 10 years, from mid-year of the first to mid-year of the second reference year. In addition, we determined ΔC ant directly for the 20 years period 1994-2014.

Spatial Clustering and Subsetting
For the fitting of the MLR models and mapping of ΔC ant in the standard configuration, we clustered the observations and predictor climatologies horizontally into the Atlantic, Pacific and Indian Ocean, according to mask "3" of our basin mask definitions ( Figure S1 in Supporting Information S1) derived from the World Ocean Atlas . To assess the contribution of the impact of this choice on the uncertainty of the reconstructed changes in C ant , we investigated five other basin configurations ( Figure S1 in Supporting Information S1). For clustering in the vertical dimension, we used the same neutral density (Jackett & McDougall, 1997) levels as employed by Gruber et al. (2019). Surface water samples collected shallower than 100 m were excluded from the MLR fitting to avoid seasonally biased observations.

Mapping C* and ΔC ant
The spatial distribution of C* was mapped by predicting the best MLR models of each reference year with climatological distributions of the predictor variables. For this purpose, the 10 best MLR models within each spatial cluster were selected as those with the lowest summed RMSE for the two paired sampling periods following Clement and Gruber (2018). The 10 individually mapped C* distributions were then averaged and subtracted to derive the mean ΔC ant distribution. In contrast to prior applications of the method, we included negative mapped ΔC ant values, since (a) ΔC ant can regionally be negative when water with low C ant displaces water with high C ant , (b) setting negative values to zero could lead to positively biased ΔC ant inventories, and (c) our tests with synthetic data revealed a tendency to lower biases when negative values were retained.

Surface Equilibrium ΔC ant
In the standard configuration, the equilibrium ΔC ant distribution at the sea surface was computed based on a rearranged definition of the Revelle factor, γ, as ΔC ant,eq (t ref,n -t ref,n+1 ) = 1/γ × DIC/pCO 2 × ΔCO 2,atm (t ref,n+1 -t ref,n ), where DIC, pCO 2 and γ are the climatological surface values  adjusted to the mean CO 2,atm of each analysis period. This adjustment of the climatological surface CO 2 -system parameters to the mean CO 2,atm was achieved by calculating as a first step the surface pCO 2 in 2002 based on the climatological values for temperature, salinity, DIC and TA, which are normalized to the same year . In a second step, the surface ocean pCO 2 was shifted according to the change in CO 2,atm , and DIC and γ were recalculated based on the new surface pCO 2 . Thus, our surface equilibrium approach takes changes in the surface ocean buffer capacity into consideration. All CO 2 -system calculations were done with the R-package seacarb (Gattuso et al., 2021) using the CO 2 dissociation constants from Lueker et al. (2000), the fluoride association constant from Perez and Fraga (1987) or Dickson and Riley (1979) at temperatures below 9°C and the acidity constant of hydrogen sulphide from Dickson (1990). To assess the uncertainties associated with this equilibrium ΔC ant estimate, we also used OceanSODA, an independent observation-based estimate of the increase in surface DIC .
While previous studies defined a distinct depth and neutral density threshold to separate water masses for which the surface equilibrium or eMLR(C*) reconstructions are used to determine ΔC ant , we blend both estimates smoothly over the top 200 m. For this purpose, the equilibrium ΔC ant is calculated at the sea surface only, while the eMLR-based ΔC ant is initially mapped across the entire water column. In a post processing step, the surfaceand eMLR-based ΔC ant estimates are averaged proportionally according to the water depth across the upper 200 m (e.g., 75% surface-based and 25% eMLR(C*)-based estimate at 50 m water depth).

Computation of Global ΔC ant Inventories and Sensitivities
Column inventories and inventories of ΔC ant in this study represent integrals across the upper 3,000 m of the water column. ΔC ant reconstructions below 3,000 m are not included in integrals to avoid the imprint of ΔC ant uncertainties that are small in terms of amount content but considerable in terms of integrated inventory changes. Instead, we follow previous studies and account for C ant storage changes below 3,000 m by adding 2% to our global ΔC ant inventories. This deep ocean scaling represents the fraction of the total C ant inventory in 1994 beneath 3,000 m according to Sabine et al. (2004). We further scale our global inventories for the storage of C ant in unmapped regions according to previously determined fractions of the global C ant storage that occurs in these regions, namely 2% in the Arctic Ocean (Tanhua et al., 2009), 1.5% in the Mediterranean Sea (Palmiéri et al., 2015), 1% in the Nordic Seas (Olsen et al., 2010), and 0.3% in the Sea of Japan (Park et al., 2006). In sum, the upscaling amounts to 7% of our directly mapped global ΔC ant inventory. Regional inventories refer to the integral of directly mapped ΔC ant distributions and no areal scaling was applied, for example, the regional inventory of the Atlantic Ocean does not account for storage in the Mediterranean Sea. Thus, our global inventory differs from the sum of the regional inventories by 7%.
To relate the global decadal change in the oceanic C ant inventory (Inv(∆C ant )) to the primary driver of the C ant uptake, that is, the increase in atmospheric CO 2 over the same decade (∆CO 2,atm ), we compute the storage sensitivity β = Inv(∆C ant )/∆CO 2,atm . Globally, the sensitivity β measures how much additional anthropogenic carbon the ocean takes up and stores Ríos et al., 2012) in response to the increase in atmospheric CO 2 . As long as atmospheric CO 2 increases close to exponentially, as was the case for the past 50 years, the sensitivity β is expected to remain constant, provided that ocean circulation remains in steady-state (Keeling, 1979) and ocean chemistry does not change. This implies that under these circumstances, changes in the oceanic storage of C ant scale linearly with changes in atmospheric CO 2 , that is, Inv(∆C ant ) = β * ∆CO 2,atm . Conversely, deviations from this linear scaling indicate changes in ocean circulation or other changes affecting β, such as the reduction in the oceanic buffering capacity owing to ongoing ocean acidification. Note, however, that this metric will be less useful when the atmospheric CO 2 growth rates will start to fall (see Gruber et al. (2023) for discussion). For our regional and basin-scale analyses, we will use the regional expression of this sensitivity, that is, the area-normalized storage sensitivity β area , to achieve direct comparability of the storage sensitivity across different regions. It is defined as β area = β/A, where A is the surface area of the region under consideration.

Determination of Uncertainty and Method Testing
The most important sources of uncertainty associated with the eMLR(C*) method are structural in nature and involve choices associated with the configuration of the method (Clement & Gruber, 2018;Gruber et al., 2019). We identified the following six configuration choices as critical (a) the regional clustering of the data, (b) the approach to perform data adjustments, (c) the nutrient used to compute the target variable C*, (d) the approach to estimate surface ocean ΔC ant , (e) the gap-filling of flagged data, and (f) the choice of predictor climatologies. In addition, we consider (g) an uncertainty contribution in our global ΔC ant inventories arising from the scaling to account for C ant storage changes in unmapped waters.
To assess the impact of these choices and to obtain an estimate of uncertainty of our reconstructions, we have reconstructed a set of total 10 alternative estimates of ΔC ant using modified choices for each of the six configurations listed above. These modifications of the eMLR(C*) configuration are described in detail in Text S4 in Supporting Information S1. We base our estimate of uncertainty on the ΔC ant offsets between the standard case and the reconstructions obtained with the configuration changes. The individual offsets are considered as independent uncertainty contributions and combined as the square root of the sum of the squares (RSS) to derive the standard uncertainty (±1 ) of our reconstructions, which we consider as a 68% confidence interval. Throughout the results and discussion, we report all results together with this ±1 uncertainty. In figures, we display also the expanded ±2 uncertainty representing a confidence interval of 95%. For each variable that is derived from a primary ΔC ant estimate (e.g., the decadal difference between two ΔC ant estimates, the ocean-borne fraction, etc.) we combine individual uncertainties through standard error propagation.
In order to independently assess the quality of our ΔC ant reconstructions against a known truth, we used synthetic data generated from the GOBM CESM-ETHZ (Doney et al., 2009;Hauck et al., 2020;Yang et al., 2017) following the approach of Clement and Gruber (2018). The GOBM is a hindcast model forced with reanalyzed atmospheric data (Tsujino et al., 2018) and the observed atmospheric CO 2 trajectory. The synthetic data set was generated by subsetting the model output in space and time according to availability of real-world observations. The eMLR(C*) approach was then applied to the synthetic data set to reconstruct ΔC ant . The comparison of the reconstructed ΔC ant to the known model truth allows us to determine the biases of the reconstruction. Details and results of this assessment are given in the Text S5 in Supporting Information S1.
In order to assess the sensitivity of our reconstructions to other sources of uncertainty, such as the limited sampling of the spatio-temporal variability of the changes in DIC, we performed additional eMLR(C*) analyses, wherein we pushed the configurations beyond the limits of what we consider a reasonable modification. For example, we limited the observations to those from repeatedly occupied sections, used other DIC variants than C* as target variable and omitted the data adjustments (see Section 4.2). In contrast, other choices associated with the application of the eMLR(C*) method have shown to be of minor relevance for the uncertainty of the reconstructed ΔC ant and were thus not considered in the uncertainty assessment (Text S2 in Supporting Information S1). This applies, for example, to the choice of the a priori C ant estimate used for the adjustment of C* to the reference year (Clement & Gruber, 2018).

Column Inventories
The storage changes in anthropogenic carbon (ΔC ant ) integrated over the upper 3,000 m (referred to as column inventories) reveal strong similarities during both decades of our analysis, that is, between 1994between and 2004between , and between 2004between and 2014. The ΔC ant column inventories vary markedly with latitude. The highest mean decadal C ant storage changes of about 10 mol m −2 dec −1 located near the center of the subtropical gyres (30-40°N/S) are about twice as high as those in the equatorial regions (10°N/S) and in the Southern Ocean south of 60°S. In the Southern Hemisphere, the consistently high ΔC ant column inventories in all subtropical gyres form a circumpolar band, while in the Northern Hemisphere the storage changes per unit area in the Atlantic exceed those in the Pacific roughly by a factor of two. The general patterns of our ΔC ant reconstructions are reminiscent of those reconstructed for the pre-industrial to 1994 period  and for the 1994 to 2007 period (Gruber et al., 2019), supporting in first approximation the expectation of a steady-state increase in the oceanic storage of C ant . We also note that our ΔC ant column inventory reconstructions for the 1994-2007 period ( Figure  S6 in Supporting Information S1) confirm those reported by Gruber et al. (2019) for the same period, supporting the consistency of our approaches and estimates despite the modifications of the eMLR(C*) method and our use of an updated database.
While the general pattern of the increases in the C ant column inventories are similar in the two decades, there are distinct differences ( Figure 2c). These decadal differences are most pronounced in the Atlantic Ocean, where we find a shift of the highest ΔC ant column inventories from the Northern to the Southern Hemisphere. In the subtropical latitudes (20-50°N) of the North Atlantic, the mean, area-weighted ΔC ant column inventory in the second decade is 3.7 ± 0.8 mol m −2 dec −1 lower compared to the first one (1994-2004: 12.1 ± 0.5; 2004-2014: 8.4 ± 0.6 mol m −2 dec −1 ; ±1σ uncertainty). In contrast, the mean ΔC ant column inventory in the subtropical latitudes of the South Atlantic (20-50°S) is 4.7 ± 2.2 mol m −2 dec −1 higher in the second decade (1994-2004: 8.5 ± 0.9; 2004-2014: 13.2 ± 1.9 mol m −2 dec −1 ). In other regions of the ocean, the decadal differences are within or close to the bounds of the uncertainty of our estimates inferred from a set of alternative reconstructions. We thus refrain from further analysis. Overall, the spatial patterns in the decadal difference of our ΔC ant column inventories are reminiscent of the anomaly structure that Gruber et al. (2019) derived from a comparison of their ΔC ant column inventories to the steady-state projection of the total C ant inventory in 1994 from Sabine et al. (2004).

Vertical Distribution
The high column inventories of ΔC ant in the centers of the subtropical gyres seen in Figures 2a and 2b are due to a deeper penetration of ΔC ant in these regions. This is illustrated by plotting ΔC ant along a global section (Figures 3a and 3b) that connects the zonal mean sections of the Atlantic and Pacific Ocean with a meridional mean section crossing the Indian Ocean sector of the Southern Ocean. Within the subtropical gyres, ΔC ant exceeding 5 μmol kg −1 dec −1 reaches at least 300 m deeper than in the equatorial regions. This is primarily a consequence of the passive tracer transport of C ant along isopycnal surfaces (depicted in Figure 3c), which brings C ant more rapidly into the ocean's interior in regions of downward sloping isopycnals (Bopp et al., 2015;DeVries & Primeau, 2011). In the South Pacific, the ΔC ant signal reaches deeper into the water column compared to the North Pacific, which is a persistent pattern for both decades and attributed to the formation of Subantarctic Mode and Antarctic Intermediate Water (SAMW and AAIW). In contrast, in the Atlantic Ocean we identify a shift of the deepest penetration of the ΔC ant level of 5 μmol kg −1 dec −1 from the Northern to Southern Hemisphere ( Figure 3c).
These changes in the downward extensions of ΔC ant cause the decadal differences identified in the column inventories of ΔC ant in the Atlantic Ocean ( Figure 2c). The decrease of the mean ΔC ant column inventory in the Stippling on the maps with small/large dots indicates regions in which the column inventories are lower than the respective 2σ-/1σ-uncertainty. Likewise, shading around the zonal integrals represents the 1σ-and 2σ-uncertainty ranges. The selected contours represent every second density slab used to cluster the data in the vertical dimension. Stippling with small or large dots indicates regions in which the ΔC ant and ΔΔC ant are lower than the 2σ-or 1σ-uncertainty, respectively. The North-South sections through the Atlantic and Pacific Ocean show zonal mean values across the entire basins, while the Southern Ocean sector is represented by a meridional mean section ranging from 55°S to 65°S that is compressed by a factor of 2.5 in distance. The location of the global section and areas for the computation of the zonal and meridional mean sections are displayed in Figure S1 in Supporting Information S1.
North Atlantic is a consequence of the weaker C ant storage increase in the North Atlantic Deep Water (NADW), evident in the zonal mean section as negative decadal differences (−2 to −5 μmol kg −1 dec −1 ) at neutral densities >27.5 kg m −3 (Figure 3c). In contrast, the decadal increase of the ΔC ant column inventory in the South Atlantic can be attributed to an intensified rate of C ant storage in the Subantarctic Mode Waters (SAMW) and Antarctic Intermediate Waters (AAIW). Here, positive ΔC ant differences in the zonal mean sections are well confined to the neutral density slabs ranging from 26.5 to 27.5 kg m −3 (Figure 3c). The decadal ΔC ant differences in the NADW and the AAIW are larger than the uncertainty of our reconstructions, while in most other water masses, the decadal differences are not significant (stippling in Figure 3c).
Examining the vertical distribution of the reconstructed ΔC ant in terms of area-weighted mean content profiles (Figure 4a), we find that globally C ant in the upper 50 m of the ocean increased on average by 10-11 μmol kg −1 dec −1 in both decades. As the storage change near the surface is determined based on an assumed equilibrium with the atmospheric CO 2 , it shows latitudinal patterns that reflect gradients in the surface ocean buffer capacity (Figures 3a and 3b). However, regional differences of the surface ΔC ant averaged over hemispheric ocean basins are small (Figure 4a). Furthermore, the almost identical surface ΔC ant for both decades is due to a compensation of the higher atmospheric CO 2 growth rate and the reduced surface ocean buffer capacity in the second decade of our analysis. In contrast to the rather steady accumulation of C ant in the surface ocean, the decadal differences in the mean ΔC ant profiles tend to be more pronounced at depth in some regions. In the North Atlantic, ΔC ant below 1000 m is about 3 μmol kg −1 dec −1 lower during the second decade (Figure 4a), while the mean ΔC ant signal in the South Atlantic was significantly higher by about 3 μmol kg −1 dec −1 between 500 and 1500 m during the second decade.

Regional and Global Inventories
Reflecting the rapid decrease of ΔC ant with depth (Figure 4a), almost 50% of the global C ant storage change occurs in the upper 500 m of the water column ( Figure 4b). Over the upper 1,000 m, this share increases to around 75%, except in the North Pacific, where the entire inventory increase is fully confined to the top 1,000 m. In all other ocean regions depicted in Figure 4b, a significant fraction of the C ant storage change occurs below 1,000 m (∼25% on a global basis) as a result of the more rapid water mass transport and mixing between the surface and ocean interior in these regions. The decadal differences in the ΔC ant inventories for 500 m depth layers reflect the profiles of the amount content and were found to be significant only in the three depth layers below 1500 m of the North Atlantic, as well as between 500 and 1,500 m in the South Atlantic ( Figure 4b).
Once integrated over the top 3,000 m of entire ocean basins ( Figure 5 and Table 1), the South Pacific stands out as the region with the highest increase of C ant storage for both decades (1994-2004: 8.6 ± 1.2 Pg C dec −1 and 2004-2014: 7.4 ± 1.0 Pg C dec −1 ), followed by the Indian Ocean (7.2 ± 0.9 and 5.7 ± 0.6 Pg C dec −1 ). In contrast, the North Pacific accounts for the smallest contribution to the increase in the oceanic storage of C ant in both decades (2.9 ± 0.8 and 3.2 ± 1.8 Pg C dec −1 ). The decadal differences in the increases in C ant storage are not significant at the 2 -uncertainty level in any of these three regions. This is different for the North Atlantic ( Figure 5, Table 1), which represented the third largest sink region during the 1994-2004 decade (4.8 ± 0.2 Pg C dec −1 ), but experienced a C ant accumulation rate that was significantly reduced by 0.9 ± 0.4 Pg C dec −1 during the 2004-2014 period (3.9 ± 0.4 Pg C dec −1 ). In contrast, the storage rate in the South Atlantic increased by 1.5 ± 0.8 Pg C dec −1 from the first (3.9 ± 0.5 PgC dec −1 ) to the second decade (5.4 ± 0.6 Pg C dec −1 ) of our analysis. As a consequence, the South Atlantic represents the third largest contributor to the global accumulation of C ant from 2004 to 2014. The ranking of the ocean basins in terms of their contribution to the global ocean C ant sink reflects primarily the differences in surface area of the basins.
The ranking of the regions changes markedly in terms of the area-normalized storage sensitivity β area (Table 1 and Figures S4 and S5 in Supporting Information S1), which is a measure of how much additional C ant the ocean is storing within any given region in response to the increase in atmospheric CO 2 (see Section 2.4). The North Atlantic has the highest regional storage sensitivity (25%-40% larger than the global mean), and the North Pacific the lowest one (about 45%-55% lower than the global mean). This largely reflects the differences in intermediate, mode, and deep water formations between the different basins, and will be further discussed in Section 4.1. When integrating ΔC ant globally and scaling it for unmapped regions and deep water storage, we determine an ocean sink for anthropogenic CO 2 that amounts to 29.3 ± 2.5 and 27.3 ± 2.5 Pg C dec −1 for the 1994-2004 and the 2004-2014 decade, respectively (Figure 5, Table 1). These inventory changes of C ant are indistinguishable for both decades with a difference of −1.9 ± 3.6 Pg C dec −1 , indicating that the global ocean continued to act as a strong sink for anthropogenic CO 2 in the recent past. The global area-normalized storage sensitivity β area decreased markedly and significantly, however, from 0.37 ± 0.03 mol m −2 ppm −1 for the decade 1994-2004 to 0.31 ± 0.03 mol m −2 ppm −1 during the second decade 2004-2014 (Table 1), suggesting a slowdown of the global ocean C ant uptake relative to what one would expect on the basis of the growth in atmospheric CO 2 . For the 20-year period from 1994 to 2014, the directly estimated global ΔC ant inventory is identical to the sum of the estimates from the individual decades (56.7 ± 3.5 Pg C dec −1 ).

Uncertainty Assessment
Nearly all identified main configuration choices of the eMRL(C*) method matter for the reconstruction of the decadal increases in C ant and contribute to the uncertainties of the global inventory changes of about ±10% (Table 1, Figure 6; Figures S9 and S10 in Supporting Information S1). The smallest uncertainty contributions stem from alternative definitions of C* and from uncertainties associated with the predictor climatologies used for mapping ( Figure 6). All other configuration choices matter more, although differently for the two decades, largely reflecting differences in the data distribution and data consistency.
The single largest contribution to the global inventory uncertainty comes from the way we applied the adjustments to the DIC, TA, and PO 4 3measurements in order to ensure the highest level of data consistency. If the adjustments were determined and applied separately for each cruise instead of in a bulk manner as done in the standard configuration, the global ΔC ant inventories change by ∼−0.5 Pg C dec −1 during the 2004-2014 period,   and assuming proportional growth with atmospheric CO 2 (see text for details). and by ∼+1.5 Pg C dec −1 for the 2004-2014 period. Nearly all of the changes during the earlier decade originate at the Indian Ocean, while those for the second period originate at the North Pacific ( Figure 6; Figure S9 in Supporting Information S1), as expected given that these were the regions where we applied these adjustments.
The second most important uncertainty contribution comes from the regional clustering of the observations with uncertainty contributions to the regional inventories ranging from 0.2 to 0.8 Pg C dec −1 (Figure 6; Figure S9 in Supporting Information S1). However, for the global ΔC ant inventories, the uncertainties arising from the regional clustering partially cancel out, such that the global uncertainty contribution is lower than the sum of the uncertainties in the individual basins (<1.5 Pg C dec −1 for both decades, Figure 6).
Among the other choices, the approach for the surface ΔC ant reconstruction contributes about 1 Pg C dec −1 to the uncertainty of the global ΔC ant inventory, determined by comparing our estimate based on the assumption of surface ocean equilibrium of C ant with a reconstruction based on observation-based surface changes of DIC from OceanSODA . Choices associated with the gap filling also contribute about 1 Pg C dec −1 to the global uncertainty, but this error source is limited largely to the first decade and the South Pacific ( Figure  S9 in Supporting Information S1), where we had to include a substantial number of cruises that took place in the 1990s and provided only calculated TA data (Figure 1). Finally, our scaling of the ΔC ant inventories for unmapped regions and the deep ocean introduces an additional uncertainty contribution of about 1 Pg C dec −1 , however, confined to the global ΔC ant inventories and of very similar magnitude for both decades.
In addition to investigating the contribution of the configuration choices to the uncertainty of our standard configuration, we also determined the sensitivity of our results to a set of additional decisions we had to take for our reconstructions (see Text S4.2 in Supporting Information S1). We did not add these results to our formal uncertainty estimate since we consider our choices as well justified, and the alternatives as clearly inferior choices. Among all these additional decisions, we find that the biggest sensitivity of the results is associated with the For consistency with the regional estimates, the global β area estimate is based on the directly mapped ΔC ant inventory, whereas the global β estimate corresponds to the global ΔC ant inventory as given in this table, obtained by scaling the mapped inventory to global coverage. All values refer to the standard cases of our ΔC ant reconstruction and the corresponding 1 -uncertainty ranges. Decadal differences are tagged with ** or * when they exceed the combined 2 -or 1 -uncertainty of both decades, respectively.

Table 1
Inventories of the Change in Anthropogenic Carbon Storage (ΔC ant Inventory) for the First (1994First ( -2004 and Second (2004Second ( -2014  data adjustments in the North Pacific. Our reconstructions based on the unadjusted data reveal that the ΔC ant inventory in the North Pacific for the 2004-2014 period would be about 4 Pg C dec −1 higher than the standard case reconstruction (Figures S11 and S12 in Supporting Information S1) that includes the adjustments. Hence, the global ΔC ant inventory in the last decade would be higher than that obtained for the 1994-2004 period without North Pacific data adjustment. As a consequence, the growth of the C ant inventory in the last decade would be roughly proportional to growth in atmospheric CO 2 , such that we would interpret the storage sensitivity as unchanged over both decades and the industrial period. While multiple lines of evidence suggest that our data adjustments are justified and reflect inter-decadal measurement inconsistencies, the relevance of these adjustments for the interpretation of trends in the ocean carbon sink underpin the importance of long-term consistent ocean observations. Furthermore, we tested the sensitivity of our results to data coverage by reconstructing ΔC ant with observations only from cruise sections that were reoccupied during both sampling periods. This reoccupation filter has generally a low impact on the reconstructed column inventories and basin inventories (Figures S11 and S12 in Supporting Information S1), suggesting an overall sufficient data coverage. An exception to this are the ΔC ant reconstructions in the Indian Ocean, which are more sensitive to the reoccupation filter ( Figure S11 in Supporting Information S1) due to the lack of data from the Arabian Sea during the central sampling period, that is, the 2000s (Figure 1b).
Our observation-based uncertainty and sensitivity findings are corroborated by our tests with synthetic data generated from an ocean hindcast model (Text S5 in Supporting Information S1), following the approach developed by Clement and Gruber (2018). Comparing the biases of our reconstructed model inventories to the estimated uncertainty based on the configuration choices, we find that the bias of 7 (11) out of 12 reconstructed Figure 6. Inventory uncertainty contributions at the 1σ-uncertainty level determined as offsets between our standard case reconstruction and six configuration choices of the eMLR(C*) method (colors) for each ocean region and both decades (panels). All offsets are shown as absolute values. Where applicable, + or − symbols on top of each bar indicate that higher or lower inventory changes, respectively, were obtained when the configuration change was applied. The uncertainty contribution that accounts for our upscaling for unmapped water masses is shown for the global inventory.
ΔC ant inventories is within the 1σ-(2σ-) uncertainty range ( Figure S18 in Supporting Information S1), which meets the expectation of a 68% (95%) confidence interval. Our uncertainty ranges and confidence intervals are thus considered suitable criteria to evaluate the significance of our ΔC ant reconstructions. Furthermore, the eMLR(C*) method proves capable of reconstructing the global ΔC ant patterns ( Figure S13 in Supporting Information S1). The reconstruction biases of the ΔC ant column inventories are below 2 mol m −2 dec −1 for 87% of the total ocean surface area, within 2-4 mol m −2 dec −1 for 13%, and only exceed the latter threshold for <0.5% of the surface area. In contrast, the true ΔC ant column inventories in our model are larger than these thresholds of 2 and 4 mol m −2 dec −1 over more than 90% and 55% of the surface area of the ocean, respectively (Figure S13a in Supporting Information S1). In agreement with this assessment based on synthetic data, the observation-based ΔC ant column inventories exceed the 1σ-and 2σ-uncertainty level over more than 90% and 75% of the total ocean surface area (Figure 2a). ΔC ant column inventories that are lower than the local uncertainty are confined to regions with low ΔC ant column inventories, primarily in the North Pacific (Figure 2a). We conclude that the global distribution patterns of ΔC ant are robustly reconstructed, which is in line with previous assessments of the method (Clement & Gruber, 2018) despite the shorter sampling periods applied in this study.
The tests with synthetic data subsampled from a GOBM are useful to demonstrate the skill of the eMLR(C*) method to reconstruct the inventory and spatial patterns of ΔC ant . However, the tendency of the GOBMs to simulate a lower decadal variability of the ocean carbon sink compared to observation-based surface CO 2 flux products (Hauck et al., 2020) suggest that this testing approach may not be well suited to investigate the method's ability to detect decadal differences in the C ant storage rates (ΔΔC ant ). As a consequence of the model's low decadal variability, the spatial patterns in the column inventory biases show a strong correlation with the reconstructed ΔΔC ant (Figures S14 and S15 in Supporting Information S1). However, the observation-based ΔΔC ant column inventories in the Atlantic Ocean (Figure 2b) are higher than the ΔΔC ant biases in the tests with synthetic data by about a factor of two ( Figure S15 in Supporting Information S1), suggesting that the observation-based ΔΔC ant patterns carry, at least in part, a true signal.

Comparison With Regional Estimates
As a final component to assess the robustness of our estimates, we compare our ΔC ant reconstructions to previous regional estimates that-same as our study-resolve changes for at least two periods and apply an MLR approach to ocean interior observations (see Text S6 in Supporting Information S1 for details). Regional studies that fulfill these criteria are available for the North and South Atlantic (Gao et al., 2022;Wanninkhof et al., 2010;Woosley et al., 2016), as well as the North and South Pacific , but not for the Indian Ocean (Murata et al., 2010;Sabine et al., 1999). We conclude from this comparison that the magnitude, patterns and trends in our ΔC ant reconstructions agree with those determined in regional studies, and that differences can-where they exist-be attributed to differences in the chosen integration depth, differences in the definition of the target variable C*, and sometimes also to the uncertainty associated with the computation of a whole basin inventory from a single reoccupied transect (Gao et al., 2022;Woosley et al., 2016).
The most pronounced difference to a regional estimate exists in the South Pacific, where Carter et al. (2019) determined a change of the C ant inventory of 5.4 ± 0.6 Pg C dec −1 from 1995 to 2005, whereas we determine a substantially higher inventory change of 8.6 ± 1.2 Pg C dec −1 for almost the same period (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004). A main difference between these studies is the calculation of C* without  or with (this study) a TA contribution. Our sensitivity reconstruction of ΔC ant in the South Pacific without considering the contribution of TA for the calculation of C* indeed reveals an inventory change that is about 3 Pg C dec −1 lower than our standard case reconstruction. Carter et al. (2019) relied on a synthetic data sensitivity test that included synthetic measurement uncertainties to conclude that the eMLR results are more robust for their implementation when the TA adjustment is omitted from the C* calculation. However, methodological differences between the methods used in that and our study limit the applicability of their sensitivity tests for our approach. We further contend that the additional attention paid to TA quality control results in better TA consistency between cruises than is assumed by Carter et al. (2019) and that the possibility of DIC changes driven by the calcium carbonate cycle (that are neither well represented in the synthetic data nor sufficiently correlated with nutrient and oxygen changes that they would be removed by following the eMLR approach) should not be neglected.
The second largest difference from a regional ΔC ant inventory exists in the North Atlantic Ocean, where Wanninkhof et al. (2010) determined a C ant storage rate of only 1.9 ± 0.4 Pg C dec −1 for the 1993-2003 period based on a single reoccupied cruise section. Their estimate is drastically lower than ours for the 1994-2004 decade (4.8 ± 0.2 Pg C dec −1 ). Our tests of the eMLR(C*) method reveal that the uncertainty of our estimates can only explain a minor part of this offset. We conclude that the offset is primarily due to differences in the integration depth, structural uncertainties in the regional estimate (Wanninkhof et al., 2010) and extrapolation errors from a single reoccupied cruise section to a whole basin inventory. In fact, Woosley et al. (2016) found that the whole basin inventories of the North Atlantic differ by ∼30% when comparing estimates obtained from a single section to those obtained from three sections.
In general, the individual differences between our regional ΔC ant inventories and those obtained in the previous regional studies are within the uncertainty range of our global inventory. Therefore, we do not assume that offsets at the regional scale challenge the robustness and interpretation of our global inventories.

Decadal Trends or the Ocean at a Time of Change
Our ocean-interior observation-based reconstruction of the changes in the global oceanic storage of anthropogenic carbon reveals for the two consecutive decades since 1994 indistinguishable accumulation rates, that is, 29.3 ± 2.5 Pg C dec −1 and 27.3 ± 2.5 Pg C dec −1 for 1994-2004 and 2004-2014, respectively ( Figure 5, Table 1).
Putting this continued accumulation of anthropogenic carbon into the context of the total oceanic C ant storage in 1994 (118 ± 19 Pg C), we confirm the expected linear correlation with the increase in atmospheric CO 2 (Figure 7, see also Gruber et al., 2023). But for a more detailed and quantitative discussion of the recent decadal trends in the  (Table 1) to the total C ant inventory in 1994 (118 ± 19 Pg C) according to Sabine et al. (2004). The red line in (a) shows the time history of atmospheric CO 2 . The cumulative uncertainty in (b) for the 1994-2014 period (dark blue ribbon) assumes zero uncertainty in 1994.
ocean carbon sink, it is informative to look at potential changes in the slope of this relationship, that is, changes in the global sensitivity β and the associated area-normalized sensitivity β area (see Section 3.3). For the global sensitivity β, we compute values of 1.6 ± 0.1 Pg C ppm −1 and 1.3 ± 0.1 Pg C ppm −1 for the two decades, respectively (Table 1, Figure 7). Their average confirms the long-term mean value of 1.4 ± 0.1 Pg C ppm −1 diagnosed by Gruber et al. (2023), but the significant decrease of about 15 ± 11% between the two decades indicates a weakening of the ocean sink for C ant . The reason that this difference in β is significant, while the 7 ± 12% reduction (−1.9 ± 3.6 Pg C dec −1 ) in the global ΔC ant inventory is not, is due to the ∼10% higher growth rate in atmospheric CO 2 from 2004 to 2014 compared to that during the previous decade (Table 1).
As the ocean acidifies in response to taking up CO 2 , a decrease of the ocean sink efficiency is expected due to the decrease of the ocean buffer capacity (Jiang et al., 2019). Over the past 40 years, seawater that followed the same pCO 2 increase as the atmosphere would have experienced a ∼6% reduction of the DIC increase per change in pCO 2 roughly every 10 years. This estimate is based on fundamental marine CO 2 -system considerations and reflects an increase in the Revelle factor (Sarmiento & Gruber, 2006). The 6% decadal weakening of the ability of the surface ocean carbonate chemistry to buffer the increase in pCO 2 can explain about half of the observed decrease in the sink sensitivity β. The other half is most likely attributable to changes in the ocean's circulation and upper ocean stratification (Sallée et al., 2021) that appears to have led to a less efficient downward transport of C ant , which we discuss further in the following.
Roughly half of the decrease of the global ocean carbon sink stems from the reduced decadal storage changes in the North Atlantic (−0.9 ± 0.4 Pg C dec −1 ). Here, we find a significant weakening of the area-normalized sink sensitivity β area (−0.14 ± 0.04 mol m −2 ppm −1 ) when comparing the first (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) to the second decade (2004-2014) of our analysis (Table 1, Figures S4 and S5 in Supporting Information S1). Furthermore, our β area estimates for both decades are well below that obtained for the 1800-1994 period , indicating a progressive weakening of the sink efficiency in the North Atlantic. The most plausible explanation for this progressive weakening is a tendency of the Atlantic Meridional Overturning Circulation (AMOC) to weaken since the 1980s (Jackson et al., 2019(Jackson et al., , 2022Latif et al., 2022). Attributing the decadal ΔC ant differences to changes in AMOC strength is supported by the localization of the negative ΔΔC ant signal in the North Atlantic Deep Water (Figure 3c). This view is also in line with two previous sets of regional studies: (a) Pérez et al. (2010Pérez et al. ( , 2013 found that the C ant storage rates in the North Atlantic subpolar gyre during the phase of a low North Atlantic Oscillation (NAO) from 1997 to 2006 were ∼48% lower than those during the first half of the 1990s, when a high NAO phase was dominant, although the mechanistic processes linking the NAO, subpolar convection strength, gyre circulation and the AMOC are not yet fully understood. However, an important caveat regarding our finding of a progressively weakening North Atlantic C ant sink is the fact that the available ocean interior observations in the North Atlantic stem mostly from the first half of our last sampling period, the 2010s. However, after 2013 the NAO switched to a strong positive phase (Holliday et al., 2020) and in line with this, the C ant column inventories in the Central Labrador Sea rapidly increased (Raimondi et al., 2021). Likewise, a deep convection event in the Irminger Sea in winter 2014/15 injected anthropogenic carbon into the ocean interior and almost tripled the storage rates compared to those determined from previous hydrographic sections (Fröb et al., 2016). It is thus possible that our reconstructions do not capture a very recent reinvigoration of the North Atlantic C ant sink, due to the temporal distribution of the available observations.
The decadal difference of the South Atlantic C ant sink is significant at the 1σ-uncertainty level (+1.5 ± 0.8 Pg C dec −1 ) and slightly exceeds the increase expected from the growth in atmospheric CO 2 alone, expressed in an increase of the sink sensitivity (+0.1 ± 0.08 mol m −2 ppm −1 ). In contrast to the North Atlantic, the decadal change in the South Atlantic is not of progressive nature when putting it into context of the total C ant storage until 1994 ( Figure 5), that is, only the second decade reveals a tendency toward an elevated storage sensitivity. Due to the strong spatial coherence between the positive ΔΔC ant signal and the Subantarctic Mode and Antarctic Intermediate Waters (Figure 3c), we attribute the decadal differences found in the South Atlantic to increased ventilation rates of these water masses (DeVries et al., 2017;Panassa et al., 2018;Patara et al., 2021;Shi et al., 2021).
Although the decadal inventory changes in the Indian Ocean and South Pacific are much less robust than those in the Atlantic, they represent in sum a contribution of about 2.6 ± 1.9 Pg C dec −1 to the decline of the global inventory. As the negative ΔΔC ant signals in these regions are associated primarily with Antarctic Bottom Water and Lower Circumpolar Deep Waters ( Figure S7 in Supporting Information S1), the decadal differences in the C ant storage changes are likely a consequence of circulation changes as well, albeit determined with lower uncertainty than in other regions.
An additional, and globally perhaps more uniform, contribution to the decrease may stem from the observed increase in upper ocean stratification (Sallée et al., 2021). Sallée et al. found that the density contrast across the base of the mixed layer had increased by about 9% per decade between 1970 and 2018. Although their estimate pertains only to the summer, such an increase in stratification is bound to decrease the transport of C ant from the surface to depth, that is, the most important bottleneck for the uptake of C ant from the atmosphere.
We conclude that in addition to the decrease of the ocean buffer capacity, ocean circulation changes and the increase in stratification are likely to play an important role in driving the decrease in the sink sensitivity, that is, the strength of the oceanic sink for anthropogenic carbon per unit change in atmospheric CO 2 .

Comparison With Observation-Based Surface Flux Estimates: Implications for Changes in the Natural Carbon Storage
In the following, we compare our estimates of the global oceanic sink of anthropogenic carbon to an ensemble of observation-based surface CO 2 flux products assembled by the Global Carbon Budget (Friedlingstein et al., 2022). During the 1994-2004 decade, the ocean interior accumulation of anthropogenic carbon (29.3 ± 2.5 Pg C dec −1 ) exceeds the time-integrated net air-sea flux of CO 2 of 21.4 ± 2.8 Pg C dec −1 (Figure 8 and Table 2). For this comparison, we followed GCB procedures and adjusted the observation-based air-sea fluxes for a preindustrial steady-state outgassing of riverine CO 2 of 6.1 Pg C dec −1 (Jacobson et al., 2007;Resplandy et al., 2018) without considering the uncertainty contribution from this adjustment. We further excluded the air-sea flux estimates provided by Watson et al. (2020) when calculating the ensemble mean and standard deviation of the flux products (see separate analysis below). Analogous to Gruber et al. (2019), the difference between the ocean interior estimates and the surface fluxes of 7.9 ± 3.8 Pg C dec −1 can plausibly be interpreted as a loss of natural carbon from the ocean to the atmosphere (Table 2). Such natural carbon fluxes are captured by the surface flux products, but are not included in the eMLR(C*) based estimates of the accumulation of C ant in the ocean's interior (Clement & Gruber, 2018). The 1σ-uncertainty of this residual term is almost half as large as the signal itself, suggesting that while the determined flux is significant, its magnitude is not well constrained. However, postulating a loss of natural carbon for the first decade of our analysis is qualitatively in line with previous studies, which concluded that the stagnation of the ocean carbon sink during the 1990s is due to an anomalously strong outgassing of natural carbon primarily in the Southern Ocean (Landschützer et al., 2015;Le Quéré et al., 2007;Lovenduski et al., 2008). It should be noted that our approach to determine the natural carbon flux as a residual term between the net surface flux and the storage change of anthropogenic carbon is only possible at global scale, because storage changes at regional scales can also occur due to ocean interior transport of anthropogenic carbon.
For the 2004-2014 decade, there is no significant difference between the cumulative surface CO 2 fluxes (26.5 ± 1.3 Pg C dec −1 ) and our ocean interior ΔC ant inventory (27.3 ± 2.5 Pg C dec −1 ), suggesting only a minor global net flux of natural CO 2 across the air-sea interface (0.9 ± 2.9 Pg C dec −1 ). The fact that during the second decade the gain of anthropogenic and the loss of natural carbon weakened simultaneously can plausibly be explained by a weakened ventilation of the ocean interior, induced by changes in circulation and/or stratification (Sallée et al., 2021), which results in lower upward transport rates of natural carbon to the surface and, vice versa, reduced downward transport of C ant into the ocean interior. This coupling was already hypothesized in previous studies that identified a synchronised reduction of both carbon flux components during periods of a weak upper-ocean overturning circulation (DeVries et al., 2017;Lovenduski et al., 2008). Specifically, DeVries et al. (2017) suggested a more vigorous global overturning in the 1990s that drove an increased outgassing of natural CO 2 and uptake of anthropogenic CO 2 , whereas a weaker overturning in the 2000s was found to have the opposite effect. Although the periods of our and their study do not fully overlap, the tendencies toward a weaker anthropogenic carbon uptake agree.
The synchronisation of the uptake of anthropogenic and the outgassing of natural carbon was also observed at a regional scale in the Irminger Sea, where Fröb et al. (2018) detected a sharp increase of the C ant inventory from 2012 to 2015, accompanied by a decline in the natural carbon inventory. This was attributed to a deep convection event during 2015 (Fröb et al., 2016), and underlines that variability at regional scale can superimpose upon the postulated global trend toward a declining anthropogenic carbon uptake.
Our determination of the natural carbon flux as a residual component relies sensitively on the underlying estimates of the net surface CO 2 flux. Watson et al. (2020) recently suggested that the current pCO 2 -based flux estimates are all biassed low due to their lack of consideration of near-surface temperature gradients. When taking these temperature gradients into consideration, they estimated an ocean carbon sink that is up to 50% larger than the unadjusted estimate. If we used their adjusted flux estimate instead of the ensemble mean of the other flux products as a basis for our calculation of the natural carbon flux (Figure 8), we would determine a release of natural carbon of 5.5 ± 3.8 Pg C dec −1 for the period 1994-2004 and an uptake of −4.0 ± 2.9 Pg C dec −1 for the period 2004-2014. Hence, we would find a weaker loss of natural carbon during the first decade and an uptake of natural carbon rather than a near-zero flux during the second decade. The integral of the natural carbon flux over the two decades would be near zero rather than a net loss. Thus, if we used the temperature adjusted CO 2 flux estimate of Watson et al. (2020), the absolute values of our globally-integrated natural carbon flux estimates would change, but this substitution would not change the conclusion that these fluxes are subject to substantial decadal variability. We note that a recent re-evaluation of the temperature adjustments by Dong et al. (2022) concluded that the effects are likely smaller than suggested by Watson et al. (2020). With the nature and magnitude of these  ) and (c) represent a CO 2 flux product that is not included in the ensemble and accounts for near-surface temperature gradients (Watson et al., 2020). Note: The eMLR(C*) estimates represent storage changes of anthropogenic carbon only, while the GCB estimates include fluxes of natural and anthropogenic CO 2 (see detailed discussion in the main text).
temperature adjustments not yet settled by the community, we stay with the ensemble of unadjusted CO 2 fluxes, but emphasize that the absolute magnitude of the natural CO 2 fluxes inferred from our work may be subject to future revisions. Another need for later revision might arise in the future if recent indications that some ensemble members of the surface flux products overestimate the true decadal variability proved to be correct (Gloege et al., 2021). A final need for revision might arise if our suggested adjustments of the North Pacific data turn out to be erroneous. This would cause no change in the C nat outgassing estimated for the first decade, but it would make the C nat flux larger (more outgassing) during the second decade. In conclusion, there is strong indication that the natural CO 2 flux has been responding sensitively to recent climate change and variability, but the current estimates underpinning the estimation of this component come with too many caveats to provide a hard constraint.

Implications for the Global Carbon Budget and Climate Change
During the first decade assessed in our study (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004), the flux estimates of the GOBMs (19.5 ± 2.7 Pg C dec −1 ) and the observation-based flux products (21.4 ± 2.8 Pg C dec −1 ) are significantly lower than the eMLR(C*)-based estimate (29.3 ± 2.5 Pg C dec −1 ). In the previous section, we proposed climate-driven outgassing of natural CO 2 from the ocean as a plausible explanation for the discrepancy to the surface flux products. The same argument would apply to the GOBMs, as their flux estimate includes the natural CO 2 flux component. For the 2004-2014 period, the difference between the cumulative observation-based surface fluxes and the ocean interior storage change of C ant disappears, suggesting only a minor global net outgassing of natural CO 2 . However, the GOBMs (22.8 ± 2.7 Pg C dec −1 , Figure 8) still diagnose a 17% weaker ocean carbon sink during the second decade compared to our ocean interior estimates. The most likely explanation for this discrepancy lies in a combination of three challenges that the majority of the current generation of GOBMs are facing: (a) the surface to deep ocean transport of C ant is rather sluggish, as these ocean models tends to underestimate the ventilation rates of the ocean interior, which was illustrated in structurally similar Earth System Models (Fu et al., 2022;Terhaar et al., 2022), (b) the GOBMs reveal a generally lower decadal variability of the CO 2 fluxes compared to observation-based estimates (Gloege et al., 2021;Hauck et al., 2020), and (c) the non steady-state fluxes of natural CO 2 (i.e., the fluxes that are caused by climate variability) appear to be low in the GOBMs. For example, in the CESM-ETHZ model used in this study, the cumulative flux of natural CO 2 over the 1994-2004 period amounts only to a net uptake of 0.1 Pg C dec −1 , and ranges in terms of annual fluxes from an outgassing of 0.2 Pg C yr −1 to an uptake of 0.4 Pg C yr −1 , which is substantially lower than the observation-based estimate determined in this study.
Considering these three challenges of the GOBMs, the decadal offsets compared to both observation-based estimates (surface and interior) could be explained as follows: During the first decade, the GOBMs simulate a low anthropogenic CO 2 uptake and a low outgassing of natural CO 2 . The compensation of these two biases leads to an Note. Estimates of emissions, and the atmospheric, land and net ocean sink are based on the Global Carbon Budget 2021 (Friedlingstein et al., 2022). The net ocean sink estimates represent the cumulative surface fluxes based on surface-pCO 2 observations, adjusted for the outgassing of riverine carbon. The oceanic sink estimates of C ant are from this study, and the oceanic outgassing of natural carbon was determined as the residual between the C ant and the net ocean sink. The uncertainties of the oceanic sink estimates follow the approach of this study, and for all other estimates apply the relative uncertainties for the 2000s according to Table 6 in the GCB (Friedlingstein et al., 2022). Numbers in parentheses indicate the airborne, land-borne, and ocean-borne fractions in % of the total emissions, with propagated uncertainties from the total emissions and the sink terms.

Table 2
Main Sources and Sinks of CO 2 for the Periods 1994-2014 apparent agreement with the surface-flux and ocean interior estimates. For the second decade, the anthropogenic CO 2 uptake in the GOBMs remains low, but this bias is no longer compensated for by the bias in the outgassing of natural CO 2 . While this interpretation of the decadal differences between the three groups of estimates is internally consistent, we emphasize that most of the evaluated offsets are of similar magnitude as their uncertainties. A more bottom-up assessment requires the comparison of the different carbon flux components for the ensemble of GOBMs reported in the GCB, an effort that is currently underway in the framework of Phase 2 of the REgional Carbon Cycle Assessment and Processes (RECCAP2) project (Poulter et al., 2022).
To put our ocean carbon sink estimates over the last two decades further into the context of the global carbon cycle, we compare them with the evolution of anthropogenic carbon emissions, as well as with the land and the atmospheric carbon sink extracted from the Global Carbon Budget 2021 (Friedlingstein et al., 2022) for the same periods (Table 2). From this comparison we derive the airborne, ocean-borne and land-borne fraction of the total emissions ( Table 2). The ocean-borne fraction of anthropogenic carbon is computed by dividing, for the respective periods, the oceanic sink with the total anthropogenic emissions following the procedures in the GCB (Friedlingstein et al., 2022). In analogy, we report the airborne and land-borne fractions by relating the increase in the atmospheric CO 2 inventory or the land carbon sink to the anthropogenic emissions of CO 2 . Note that the global sensitivity β is directly related to these fractions, since it is the ratio between the ocean-and airborne fraction. For the total anthropogenic emissions we consider the emissions due to the combustion of fossil fuels (including the cement carbonation sink) and land use change.
An important aspect for the contextualization of our ocean sink estimates is the increase of the total emissions by about 25% from the first to the second decade of our analysis. In contrast, the growth of the atmospheric sink for CO 2 was only about 10% higher during the second decade, which is reflected in a decrease of the airborne fraction from 48 ± 4% to 44 ± 4%. Because the emissions grew more rapidly than the atmospheric CO 2 , the ocean-borne fraction of C ant (Table 2) decreases even more pronouncedly than the global ocean's uptake sensitivity β (Table 1), namely from 36 ± 4% to 27 ± 4% for the 1994-2004 and 2004-2014 decade, respectively, which corresponds to a reduction of the uptake fraction of −9 ± 6% (or ∼25% in relative terms). In contrast, the land sink evolved very consistently with the total emissions over the two decades, such that the land-borne fraction remained at a stable level (31 ± 6% and 30 ± 6%). Due to the identified oceanic outgassing of natural carbon during the first but not the second decade of our analysis, the net ocean sink for anthropogenic and natural carbon increases in a remarkably stable manner with the total emissions. Accordingly, the ocean-borne fraction of the total emissions in terms of the net oceanic CO 2 uptake remained unchanged at a level of 26 ± 4% and 26 ± 3%. According to our assessment, the sum of all three sink estimates exceeds the total emissions by about 5% during the 1994-2004 period, whereas the sources and sinks of CO 2 during the 2004-2014 period match almost perfectly.
Although we do not identify a change in the net ocean-borne fraction of the total emissions between 1994 and 2014, it is not for granted that the ocean carbon sink will remain constant for the decades to come (Ridge & McKinley, 2021). DeVries et al. (2017) hypothesized that a trend toward a more stratified ocean is likely to strengthen the CO 2 sink in the near future by trapping natural CO 2 in the deep ocean, but further concluded that this process may ultimately limit the net oceanic carbon sink, when the reduced uptake of anthropogenic CO 2 that continues to accumulate in the atmosphere outweighs the reduced outgassing of natural carbon. Our findings demonstrate that these compensating processes are in progress and we deem it of utmost importance to continue the monitoring of the ocean interior accumulation of carbon to keep track of them.

Caveats and Recommendations
Building on the quantitative uncertainty assessment of our ΔC ant reconstructions (Section 3.4, Text S4 and S5 in Supporting Information S1), we highlight in the following some caveats of our study and provide recommendations on how to overcome them in future studies.
While all sampling periods assigned for this study are relatively well covered with observations, the large changes we reconstruct in the North Atlantic between the first and second period need to be viewed with caution since the number of data records that provide all required variables for the eMLR(C*) analysis is very limited after ∼2015. During this period, a reinvigoration of the anthropogenic carbon accumulation has been reported in regional studies (Fröb et al., 2018;Raimondi et al., 2021), albeit only for small subregions of the whole North Atlantic. The inclusion of North Atlantic observations collected since 2015 into an eMLR(C*)-based ΔC ant reconstruction will contribute to further improve our understanding of the basin-wide C ant storage changes in this highly dynamic region.
We further expect substantially new and improved insights from the completion of another cycle of the repeat hydrography program over the 2020s. In contrast to our two decadal ΔC ant reconstructions, which both build on the same data for the central sampling period (2000s) and are thus not fully independent, reconstructing ocean interior trends with yet another decade of observations would resolve this issue. Furthermore, future investigations based on more recent data will profit from the improved data quality, in particular when becoming independent from the observations of the 1990s, which tend to be less consistent than the more recent measurements (Lauvset et al., 2021).
The continued tracking of the oceanic C ant storage, for example, by providing a global ΔC ant reconstruction for the 2014-2024 period, would also shed light on the very recent divergence of GOBMs and surface-flux products which increased to more than 1 Pg C yr −1 in 2020 (Friedlingstein et al., 2022;Hauck et al., 2020). A burning question in this regard is whether the high uptake determined by the surface flux products around 2020 can be confirmed by ocean interior estimates. Scaling our C ant accumulation estimates from the 2004-2014 period to the 2010s according to the atmospheric CO 2 increase, we would indeed project an uptake of ∼32 Pg C dec −1 , which is very similar to the mean observation-based net surface flux over the same period.
Another recommendation emerges from the high sensitivity of our results to the adjustments we applied to a subset of the observations provided through GLODAPv2.2021. This pronounced sensitivity highlights the importance of data quality and consistency for the ocean interior observing system. Continued efforts to maintain and improve the quality of seawater biogeochemical measurements, such as through the continued use of reference materials and undertaking inter-laboratory comparisons (Bockmon & Dickson, 2015), are indispensable. Furthermore, the timely submission, compilation, and harmonization of data through GLODAP appears crucial. The release of version 3 of GLODAP including a complete revision of the data adjustments is anticipated in 2024. Based on our findings, we suggest a critical revision of the observation from the Pacific with a particular focus on the TA measurements.
Tightly linked to the observational data consistency is the accuracy of deep ocean ΔC ant reconstructions. Small biases in ΔC ant can indeed exert a strong impact on the basin inventory changes due to the large volume of the deep ocean. Below 1,000 m, the mean ΔC ant reconstructed in this study is lower than 5 μmol kg −1 dec −1 across all ocean basins (Figure 4a). Despite the low ΔC ant rates compared to surface waters, the ocean below 1,000 m represents a potentially significant contribution to the global C ant inventory as it accounts for roughly 75% of the total ocean volume. On a global average the content and inventory changes between 1,000 and 3,000 m carry a significant positive signal and contribute about 25% to the total inventory integrated over the top 3,000 m. To derive our global inventories, we have chosen to account for the storage change below 3,000 m (∼30% of the total ocean volume) by scaling the inventory with +2% according to the total C ant accumulation at depth in 1994 . This approach is consistent with previous studies (Gruber et al., 2019) and represents a compromise between neglecting deep water storage changes and potentially introducing biases from integrating small and highly uncertain ΔC ant below 3,000 m. It is important to note that the general decadal trends reported in this study for the regional inventories are maintained when integrating the reconstructed ΔC ant across the full water column, that is, without the deep ocean scaling (data not shown). Nevertheless, we deem it important that future observation-based studies explicitly include also the accumulation of anthropogenic carbon in the deep ocean below 3,000 m water depth, taking advantage, for example, of measurements of transient tracers, such as SF 6 and CFCs.
Finally, our study revealed that fluxes of natural carbon are key to understanding the oceanic response to a changing climate. The comparison of our estimates of the ocean interior accumulation of anthropogenic carbon with the net surface fluxes of CO 2 allowed us to distinguish a decade with presumably strong net outgassing of natural carbon (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)) from a decade with low net fluxes of natural carbon (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). However, our quantification of the natural carbon flux as a residual quantity between two entirely independent estimates remains prone to uncertainties that are in the same order of magnitude as the flux itself (McNeil & Matear, 2013). This uncertainty is partly due to remaining uncertainties in the estimates of the surface flux products, some of which might tend to overestimate decadal variability and could be prone to systematic biases such as a missing consideration of near-surface temperature gradients. Furthermore, our testing of the eMLR(C*) method based on synthetic data from a single GOBM with low decadal variability in the anthropogenic carbon storage (see Section 3.4) and low changes in the natural carbon pool (see Section 4.3) impedes a comprehensive assessment of the method's ability to quantify decadal variability and separate changes in the inventories of natural and anthropogenic carbon. Tests of the method over a larger ensemble of GOBMs, with artificially perturbed data, or time-variant predictor variables may allow further investigation of this aspect. As natural carbon fluxes are expected to vary substantially at sub-decadal time scales, a reoccupation of selected repeat hydrography sections with increased frequency or the extension of the BGC Argo program to global coverage could provide an important observational basis for future studies of ocean interior carbon dynamics. Furthermore, progress in the development of statistical methods to separate storage changes of natural and anthropogenic carbon based on a consistent interpretation of ocean interior observations alone could provide new valuable insight. The application of neural networks to reconstruct ocean interior dynamics of DIC are a meaningful first step in this direction (Broullón et al., 2019(Broullón et al., , 2020Keppler et al., 2020Keppler et al., , 2023, albeit being based on even fewer data than surface pCO 2 products.

Conclusion and Outlook
This study provides the first global reconstruction of the decadal evolution of the ocean interior storage changes of anthropogenic carbon covering the decades 1994 to 2004 and 2004 to 2014. We provide uncertainty estimates for all reported estimates, including regional inventories and spatial distributions of ΔC ant , and decompose the uncertainties into contributions from various configuration choices associated with the eMLR(C*) method.
We find that the oceanic sink for anthropogenic carbon remained strong during both decades (29 ± 3 PgC dec −1 and 27 ± 3 PgC dec −1 , respectively). But the sink strength relative to the growth in atmospheric CO 2 (sensitivity β) and the uptake fraction of anthropogenic emissions weakened from the first to the second decade by about 15% and 25%, respectively. We attribute these changes to a decrease of the ocean buffer capacity and a reduction in the surface ocean to deep transport, induced by changes in ocean circulation (most apparent in the Atlantic) and an increase in upper ocean stratification. In contrast to our findings for the accumulation of C ant , observation-based estimates of the surface fluxes of CO 2 indicate that the net ocean sink for anthropogenic and natural carbon increased proportionally with the anthropogenic emissions. This implies that the net ocean uptake fraction remained stable throughout both decades. We attribute the difference between the anthropogenic and the net carbon sink to substantial fluxes of natural carbon that vary between decades, but remain hard to quantify.
Our results can serve as new reference points for the annual ocean sink estimates published in the Global Carbon Budget and provide guidance to further develop and assess GOBMs, which most likely underestimate the anthropogenic carbon sink. Furthermore, our reconstructions of the continuing accumulation of C ant can be used to infer acidification trends in the ocean interior at global scale.
Future studies of the ocean interior storage of C ant may allow us to address questions arising from our analysis, including the drivers for the very recent increase in the net uptake flux as determined based on surface pCO 2 -observations and the question whether compensating processes of the ocean carbon cycle remain effective, such as the regional shift of the anthropogenic carbon storage from the North to the South Atlantic and the apparent coupling between the fluxes of natural and anthropogenic carbon. Mandatory requirements to address these topics are (a) the continued and extended collection of biogeochemical ocean interior observations, that is, the completion of a fourth cycle of the repeat hydrography program and the expansion of the biogeochemical Argo program to global coverage, (b) the continued compilation of the observations into a harmonized and quality-controlled data product, and (c) the continued improvement and further development of statistical methods, for example, to separate the storage changes of anthropogenic and natural carbon. the collection and harmonisation of high-quality ocean interior observations made available through GLODAP. We are grateful for the formal reviewer comments provided by Andrew Watson and two anonymous reviewers, as well as the informal comments provided on a preprint of this manuscript by Judith Hauck and Jens Terhaar. JDM and NG acknowledge support from the European Union's Horizon 2020 research and innovation programme under grant agreements no. 821003 (project 4C) and no. 821001 (SO-CHIC). FFP was supported by the BOCATS2 (PID2019-104279GB-C21) project funded by MCIN/AEI/10.13039/501100011033 and contributed to WATER:iOS CSIC PTI. AO and SKL were supported by the project N-ICOS-2 (Research Council of Norway Grant 296012). SKL also acknowledges internal funding support from NORCE. MI was supported by JPMEERF21S20810. RW, RAF, and BC were supported by the Office of Ocean and Atmospheric Research (OAR) of NOAA, including the Global Observation and Monitoring Program (GOMO), FundRef 100018302. BC and RAF contributions are PMEL contribution 5454 and CICOES contribution 2022-1244. TT acknowledges support by EU Horizon 2020 through the EuroSea action (grant agreement 862626).