Constraining the trend in the ocean CO_2 sink during 2000–2022

Domain: earth_science
Published: 2024-09-28
Editorial outcome: minor_revision

Abstract
The ocean will ultimately store most of the CO2 emitted to the atmosphere by human activities. Despite its importance, estimates of the 2000−2022 trend in the ocean CO2 sink differ by a factor of two between observation-based products and process-based models. Here we address this discrepancy using a hybrid approach that preserves the consistency of known processes but constrains the outcome using observations. We show that the hybrid approach reproduces the stagnation of the ocean CO2 sink in the 1990s and its reinvigoration in the 2000s suggested by observation-based products and matches their amplitude. It suggests that process-based models underestimate the amplitude of the decadal variability in the ocean CO2 sink, but that observation-based products on average overestimate the decadal trend in the 2010s. The hybrid approach constrains the 2000−2022 trend in the ocean CO2 sink to 0.42 ± 0.06 Pg C yr−1 decade−1, and by inference the total land CO2 sink to 0.28 ± 0.13 Pg C yr−1 decade−1.

Introduction
The ocean plays a crucial role in the global carbon cycle and climate. Every year, it absorbs, on average, 25% of the anthropogenic CO2 emitted to the atmosphere from burning fossil fuels and land-use change1,2. The ocean will be the main storage for most of the CO2 emitted to the atmosphere from human activities, absorbing more than half of the cumulative emissions on a timescale of 1000 years, and between 60% and 85% on timescales of 10,000 years or longer3. The ocean CO2 sink responds mainly to the increase in atmospheric CO2 and is the result of processes taking place on different timescales: the dissolution of anthropogenic CO2 at the surface, its buffering by carbonate chemistry, and its transport to the intermediate and deep ocean by ocean circulation4. This uptake is modulated by the variable growth rate in atmospheric CO25. Furthermore, the ocean carbon reservoir is also sensitive to climate variability and climate change6,7, which modulates the response of the ocean CO2 sink to anthropogenic CO2 emissions, generates variability, and can alter both the decadal rate of change and long-term storage of carbon in the ocean. An accurate estimate of the trends and variability in the ocean CO2 sink is needed to better understand how the ocean carbon cycle responds to the various drivers of change and to predict its evolution and its long-term fate. An accurate assessment of the ocean CO2 sink and its variability also provides an independent constraint on the other terms of the global carbon budget, in particular the trend and variability of the CO2 sink by the terrestrial biosphere8.

Over the last two decades, considerable progress has been made in improving our quantitative assessment of the ocean CO2 sink and its trend and variability. The IPCC Fourth Assessment Report9 assessed a mean ocean CO2 sink for the 1990s based on observations of 2.2 ± 0.4 Pg C yr−1, which has stood the test of time1,3. The first estimates of the ocean CO2 sink were based on ocean general circulation models10, followed by estimates inferred from O2/N2 observations11,12 and broader geochemical data10,13. Later, the development of global ocean biogeochemical models (GOBMs) provided the first estimates of the ocean CO2 sink variability in response to climate variations and long-term increases in atmospheric CO2 concentration14,15. These models simulate the processes that regulate the full carbon cycle in the ocean (both natural and anthropogenic) and respond to changes in climate, weather, and variations in CO2 levels in the atmosphere. Results from such process-based ocean models have suggested that the ocean CO2 sink is sensitive to climate change and variability but that this variability is much smaller than the trend induced by the rising atmospheric CO216.

More recently, observation-based estimates of the ocean CO2 sink have become available from multiple data products. This has been made possible by the annual release of quality-controlled observations of CO2 fugacity (fCO2) at the sea surface—analogous to the partial pressure of CO2—compiled within the Surface Ocean CO2 Atlas (SOCAT) since 201517,18, and the emergence of machine learning and other advanced numerical approaches to extrapolate these relatively sparse space-time observations of fCO2. Estimates from these observation-based products (fCO2-products), which all used the SOCAT database as a starting point, confirmed some aspects of the ocean CO2 sink simulated by the GOBMs, including the mean ocean CO2 sink within observational uncertainties, the presence of a hiatus in its growth rate in the 1990s19,20, and of variability much smaller than its long-term trend2. It should be noted that these fCO2-products only assess the air-sea CO2 flux, and not where anthropogenic CO2 is ultimately stored in the ocean, which would require additional measurements of carbon in the water column, as well as more assumptions21–24.

Despite recent progress, the two types of estimates used within the global carbon budget annual update by the carbon research community diverge over the last two decades (2000–2022)22,23,25, with the fCO2-product ensemble suggesting a decadal rate of growth of the ocean CO2 sink almost twice as high as that simulated by the GOBMs ensemble (Table 1)2,25. The systematic nature of this discrepancy suggests a structural error in one or both methodologies used. In addition, the fCO2-products suggest a higher amplitude of decadal variability on average than that simulated by the GOBMs over the period 1990–2022 (Table 1, Supplementary Fig. S1)2. Results from models using data assimilation26,27 also suggest an underestimated decadal variability in GOBMs air-sea CO2 flux28,29. This lack of consistency between the fCO2-products and the GOBMs ensemble originates in the mid- and high-latitude regions of both hemispheres (poleward of 30°N and 30°S)2. We focus here on understanding the factor-of-two inconsistency in estimates of the 2000–2022 trend in the ocean CO2 sink between fCO2-products and GOBMs, which occurs despite the increasing number of fCO2 observations (i.e., from ~4500 gridded observational data points a year in the 1990s, to 10,000 in the 2000s, and 15,000 in the 2010s). We introduce and use a hybrid approach that uses, as a starting point, the NEMO-PlankTOM12.1 GOBM forced with meteorological and climate reanalysis data from the National Centres for Environmental Prediction (NCEP30) and observed atmospheric CO2 concentration. The hybrid approach then constrains the model outcome on a yearly basis using SOCAT observations to provide an annual estimate of the ocean CO2 sink (see methods). NEMO-PlankTOM12.1 is the latest update of an established GOBM that was used from the onset in the annual updates of the global carbon budget analysis and which was designed for the study of the ocean CO2 sink variability31. Its overall performance in simulating ocean physics and biogeochemistry is comparable to that of other GOBMs in the global carbon budget analysis (see methods, Table 1, Supplementary Figs. S2 and S3)2,19. The hybrid approach preserves the coherence of the physical and biogeochemical processes represented in NEMO-PlankTOM12.1 and goes beyond the traditional model evaluation by constraining the model output fields of fCO2 against the observed fCO2 data provided by SOCAT. A similar hybrid approach was recently published32, but with a machine learning algorithm used to derive the factors influencing the fCO2 variability. Here, the mechanism as represented in the NEMO-PlankTOM12.1 model, including the mixed-layer dynamics and the large-scale circulation, the carbonate chemistry, and the organic carbon transfer to depth resulting from biological processes (see methods for a description of the model) remained unchanged and thus also constrained the results. The hybrid approach used here thus provides an estimate of the ocean CO2 sink that is based on both observations and on current theoretical understanding of the response of the ocean CO2 sink to climate change and variability. Note that the hybrid approach used here, which optimises a target variable using multiple model simulations and a cost function, has been used in previous studies to constrain global ocean primary production33 and air-sea fluxes of N2O34 and CCl435.

Results
Constraints on the decadal variability of the global ocean CO2 sink

In order to use observations to constrain the annual global ocean CO2 sink simulated by NEMO-PlankTOM12.1 standard model simulation, perturbed simulations were deliberately produced. This was obtained by perturbing model parameters. Perturbed simulations provided a range of possible values for the ocean CO2 sink around the estimate from the standard model simulation, and covered the expected range suggested by the global carbon budget analysis. We show here results obtained with the perturbation of the half-saturation constant of bacterial remineralization, which is more homogeneous and, therefore, more robust (see methods and the Supplementary information for details of the sensitivity analyses). Then, for each year, the optimal CO2 sink was found within this range of possibilities by optimising the calculated mean square error (MSE) between the simulated fCO2 and the SOCAT observations. The hybrid approach also provides a quantitative estimate of uncertainty (see methods).

The hybrid approach increases the decadal variability of the ocean CO2 sink simulated by the NEMO-PlankTOM12.1 process-based ocean model (see methods for the definition of decadal variability). Originally, over the period 1990–2022, NEMO-PlankTOM12.1 simulated amplitudes of decadal variability for the ocean CO2 sink of 0.11 Pg C yr−1. This value is at the high end of the decadal variability simulated by the other GOBMs used in the global carbon budget analysis (Table 1). The hybrid approach further increases this simulated decadal variability by 18% to 0.13 ± 0.02 Pg C yr−1, to a value close to the decadal variability estimated by the fCO2-products (0.14 ± 0.06 Pg C yr−1).

We tested the robustness of the decadal variability produced by the hybrid approach with respect to (i) the choice in the selected model’s configuration and parameter perturbed, and (ii) the annual availability and distribution of the SOCAT data. To do this, we first applied the hybrid approach to a total of six different model setups (see methods). This first sensitivity analysis suggested a comparable increase in decadal variability (0.16 ± 0.05 Pg C yr−1, Supplementary Table 3). Secondly, the hybrid approach was applied by considering observations from three consecutive years (see methods). This second sensitivity analysis also suggested a comparable increase in the decadal variability (to 0.14 ± 0.02 Pg C yr−1, Supplementary Fig. S1). Overall, these two sensitivity analyses confirmed the robustness of the amplitude of the decadal variability suggested by the hybrid approach (Table 1). In contrast, the amplitude of the year-to-year variability was less robust because of insufficient data to constrain the hybrid approach on a yearly basis, especially in the 1990s (see Supplementary information).

Constraints on the regional ocean CO2 sink

The regional ocean CO2 sinks were constrained by applying the same hybrid approach but separately for three latitude bands, using only the observations in the North (>30°N), Tropics (30°S to 30°N), or South (<30°S) to constrain the regional ocean CO2 sink simulated by NEMO-PlankTOM12.1 (Fig. 1b–d). The hybrid approach substantially modified the simulated ocean CO2 sink in the mid- and high-latitude regions, particularly in the South, but with the Tropics remaining similar to the original NEMO-PlankTOM12.1 results. The hybrid approach increased the decadal variability simulated by NEMO-PlankTOM12.1 in the North and South regions (Supplementary Table 2).

In the North, where there are more SOCAT observations, the hybrid approach barely modified the mean ocean CO2 sink simulated by NEMO-PlankTOM12.1. However, in the South, where there are fewer observations, the hybrid approach increased the mean ocean CO2 sink simulated by NEMO-PlankTOM12.1 by 44% in the period 2000–2022. The hybrid approach applied at every 5° of latitude to constrain the mean climatological ocean CO2 sink during 2000–2022 (rather than its annual values) shows that the mismatch between the model and the fCO2 observations in the 40–60°S band could be the main cause of the underestimation of the ocean CO2 sink in the South (Fig. 2). Note that it is the Southern region that has the most influence on the global ocean CO2 sink.

Decadal trends of the global ocean CO2 sink

In the 1990s and 2000s, the hybrid approach enhanced the decadal trends simulated by NEMO-PlankTOM12.1, bringing the trends closer to those suggested by the fCO2-products. In the 1990s, the hybrid approach decreased the NEMO-PlankTOM12.1 simulated trend from 0.02 to −0.19 ± 0.17 Pg C yr−1 decade−1, to a value with the same sign and within the range of the fCO2-product ensemble trend of −0.12 ± 0.35 Pg C yr−1 decade−1. In the 2000s, the hybrid approach increased the simulated trend from 0.27 to 0.80 ± 0.21 Pg C yr−1 decade−1, to a value also close to and within the range of the fCO2-product ensemble of 0.71 ± 0.38 Pg C yr−1 decade−1 (Fig. 3, Table 1). In the 2010s, the hybrid approach decreased the simulated trend by NEMO-PlankTOM12.1 from 0.53 to 0.44 ± 0.15 Pg C yr−1 decade−1, which is below the strong decadal trend of 0.68 ± 0.19 Pg C yr−1 decade−1 suggested by the fCO2-product ensemble, but above the decadal trend of 0.34 ± 0.08 Pg C yr−1 decade−1 suggested by the GOBMs ensemble. As a consequence, the trend between 2000 and 2022 from the hybrid approach (0.42 ± 0.06 Pg C yr−1 decade−1) lies between the trend simulated by NEMO-PlankTOM12.1 and the GOBM ensemble (0.28 ± 0.05 Pg C yr−1 decade−1), and that suggested by the fCO2-product ensemble (0.54 ± 0.13 Pg C yr−1 decade−1). However, unlike the GOBM and fCO2-product ensembles where the growth is similar between the two decades, the trend in the hybrid approach is made of a decade of strong growth (in the 2000s) followed by a decade of weak growth (the 2010s). The distinct decadal trend variations between the 2000s and the 2010s suggested by the hybrid approach are robust to different configurations of the original model and to the choice of perturbed parameters (see methods). In comparison, only two out of seven fCO2-products suggest that the trend in the 2000s should be higher than the trend in the 2010s over the global ocean.

Origin of the discrepancy among estimates

The differences in the decadal trends for the 2010s among the hybrid approach, the NEMO-PlankTOM12.1 and the fCO2-products were mostly associated with the mid- and high-latitude regions (Fig. 3). Further investigation to determine the origin of the discrepancy was conducted, first by scrutinising estimates in the North where the density of SOCAT observations was the highest. This region encompasses 17% of the global ocean area. Because of the higher density of observations in the North compared with other latitudes, it is possible in this region to compare the trends in ocean CO2 flux in areas that are generally well-sampled to those in areas that are poorly sampled. A similar analysis cannot be done for the South because of the insufficient data coverage. In the North, the hybrid approach and average GOBMs produce lower trends compared to the average of the fCO2-products (Fig. 3b).

The SOCAT observations, by themselves, do not confirm the existence of a strong decadal trend in the 2010s in the North. The strong decadal trend in the ocean CO2 flux estimated by the fCO2-products is primarily driven by diverging trends between the CO2 fugacity at the surface of the ocean compared to that in the atmosphere (ΔfCO2)36. The SOCAT observations converted into ΔfCO2 in the North region between 2000 and 2019 show a positive and higher trend in the 2000s compared to the 2010s (Fig. 4). Similar temporal patterns were visible in the ΔfCO2 data from the fCO2-products subsampled to SOCAT sampling points, with a decadal trend in ΔfCO2 in the 2000s significantly higher than in the 2010s (Kruskal–Wallis test, p value < 0.01), as expected. However, when not subsampled, the fCO2-products suggested a decadal trend in ΔfCO2 in the 2000s that is not significantly higher than in the 2010s (Kruskal–Wallis test, p value = 0.14), with an overlap in the estimated uncertainties in the two decades, explaining the small differences in the CO2 sink trend between the 2000s and 2010s. This is induced by the fact that three of the seven fCO2-products suggested a greater trend in the 2010s compared to the 2000s, and a fourth fCO2-product suggested a strong trend in both decades. In comparison, when subsampled or not subsampled, the GOBMs suggested a decadal trend in ΔfCO2 in the 2000s significantly higher than in the 2010s (Kruskal–Wallis tests, p value < 0.001). In addition, the decadal trend in ΔfCO2 in the North during the 2010s is significantly higher in the fCO2-product ensemble than in the GOBMs ensemble (Kruskal–Wallis test, p value < 0.01). The differences between the subsampled and not subsampled results suggest that different extrapolation methods outside of data-rich regions could account for the higher decadal trend in ocean CO2 sink in the North over the 2010s in the fCO2-products ensemble compared to the GOBMs ensemble and the hybrid approach (Fig. 3b, Supplementary Table 2).

Within the four fCO2-products that suggested a strong positive decadal trend in ΔfCO2 in the 2010s in the North, the ocean area associated with the highest trend values overlaps the northern Pacific Ocean (Fig. 5a) which was undersampled in the 2010s (white areas in Fig. 5b). In quantitative terms, these four fCO2-products suggested that the Pacific Ocean contributed 0.04 Pg C yr−1 decade−1 of the 2010s trend in the North, while the other three fCO2-products (below the fCO2-product average) suggested that the Pacific Ocean had a negative trend in the 2010s (−0.01 Pg C yr−1 decade−1), as simulated by the GOBM average.

In addition to the extrapolation problems in the undersampled northern regions mentioned above, on a global scale, the estimate of the positive trend in the 2010s from the fCO2-product ensemble has been revised downwards in successive publications of the global carbon budget analysis between 2021 and 2023, while their trends for the 1990s and 2000s have remained relatively similar (Fig. 6). For each of the global carbon budget analyses published between 2021 and 2023, the fCO2-product ensemble average has always been produced from seven estimates. However, two of the seven fCO2-products were introduced, replacing previously submitted products that were not updated, and five were slightly updated. Among these five fCO2-products, on average, the 2010 trend between the 2021 and 2023 publications decreased by −0.05 Pg C yr−1 decade−1. Thus, it was mainly the turnover in the last two fCO2 products between 2021 and 2023 that led to a visible decrease in the ensemble average of −0.25 Pg C yr−1 decade−1 for the 2010 trend. Consequently, the downward revision observed for the observation-based estimate was mainly due to a change in two fCO2-product methodologies and, to a lesser extent, to the annual updates of the SOCAT database and fCO2-product methods, suggesting that the trend of the 2010s estimated with the fCO2-products is not robust at this stage.

Methods
Amplitude of the decadal variability and decadal trend estimates

The amplitude of the decadal variability was estimated as the standard deviation of the decadal component of the annual ocean CO2 sink time series. The decadal components were extracted from the time series using the following signal decomposition methodology52: (i) the linear trend and long-term mean were removed to isolate the temporal variability, (ii) the decadal component of this detrended time series was obtained by filtering this time series with a 10-year Hanning window. The Hanning window is a filtering function with a ‘bell-shaped’ curve used to smooth the signal by emphasizing the feature near the center of the window. All uncertainty ranges are reported here with a ±1σ (68%) confidence interval, as in global carbon budget analysis2.

The decadal trends of ocean CO2 sink and ΔfCO2 were estimated with the Theil-Sen slope estimator61,62. To calculate ΔfCO2, the atmospheric CO2 mole fraction (xCO2) used in the global carbon budget analysis (data from the U.S. National Oceanic and Atmospheric Administration, Global Monitoring Laboratory63) is first converted to pCO2 taking into account the atmospheric surface pressure corrected for the vapour pressure of seawater (Patm), and then to fCO2 with a fugacity coefficient estimated as follows64,

fCO2 = pCO2 * exp(Patm * (B + 2 * (1 - xCO2)2 * δ) / (R * T)) (1)

where T is the sea surface temperature (in Kelvin, from OISST 1.265), B and δ are the virial coefficients for carbon dioxide66 and R is the gas constant64. The necessary sea surface salinity data come from EN4 (EN.4.2.2.g1067), and the surface atmospheric pressure data come from ERA568.

For the GOBM and fCO2-product ensembles (from the global carbon budget analysis), the amplitude of the decadal variability, and the decadal trends were calculated for each member of the ensemble. The mean ensemble values are reported with their standard deviation. For the hybrid approach, the amplitude of the decadal variability, and decadal trends were calculated with the annual constrained values. An estimate of the standard deviation around these values was obtained by re-calculating 10,000 times these estimates with annual ocean CO2 sink values that had been randomly selected (from a uniform distribution) within the estimated confidence intervals of the annual constrained value.

Description and evaluation of the NEMO-PlankTOM12.1 model

The NEMO-PlankTOM12.1 model consists of a global ocean general circulation model, NEMO v3.6, with an embedded biogeochemical model, PlankTOM12.1, forced by atmospheric meteorological data from the NCEP reanalysis product30.

NEMO-PlankTOM12.1 used the NEMO model69 in its global configuration on the ORCA tripolar grid, with a longitudinal resolution of 2° and an average latitudinal resolution of 1.5°, the latter being enhanced up to 0.3° in the tropics and at high latitudes, and a temporal resolution of 96 min. This physical ocean model comprises a total of 31 vertical z levels with a vertical resolution of 10 m for the first 100 m, decreasing progressively to a resolution of 500 m at a depth of 5 km. The NEMO model is based on the Navier-Stokes equations and a non-linear equation of state. It explicitly calculates vertical mixing using a turbulent closure model. Subgrid-scale eddy-induced mixing is represented with a parameterisation70. NEMO is coupled to the Louvain-La-Neuve sea ice model (LIM71).

The PlankTOM12.1 biogeochemical model simulates the full marine cycles of carbon, oxygen, phosphorus and silicon, and simplified cycles for iron and nitrogen. This biogeochemical model was obtained by merging two versions of the PlankTOM model series, which had been developed in parallel, one focused on the role of jellyfish72, and one focused on the role of pteropods73. This version has been used in the global carbon budget analysis 202239 and 20232. Its ecosystem component is based on the representation of 12 Plankton Functional Types (PFTs), including six phytoplankton, five zooplankton and one bacteria. Spatiotemporal variations in PFT concentrations are induced by the simulated response of each PFT to environmental conditions, including temperature, nutrient availability, light, and interactions between PFTs. PlankTOM12.1 explicitly represents dissolved organic carbon and two size classes of particulate organic carbon, one small and one large. These components are influenced by the particle aggregation process, and the large particles are also influenced by the effect of mineral ballasting. Simulated dissolved inorganic carbon and alkalinity are influenced by air-sea exchanges of CO2, calcification (production and dissolution), primary production, and remineralisation of organic matter (grazing by zooplankton and remineralisation by bacteria). The alkalinity is also influenced by denitrification. A full description of PlankTOM12.1 biogeochemistry has been published34. Model simulations are too short to fully represent the input of river fluxes and their subsequent outgassing of CO2 in the open ocean. Instead, constant river fluxes of dissolved and organic carbon and nutrients are prescribed as input at the location of river mouth, and corresponding fluxes are removed from the bottom sediments to conserve mass. The version of the NEMO-PlankTOM12.1 code used here is the same as that used in the latest global carbon budget analysis, forced with NCEP reanalysis. For this, the model was spun up first from 1750 to 1947 with a 30 years (1948−1977) climatological annual cycle of atmospheric forcing from the NCEP reanalysis product, followed by the use of annual forcing from 1948 until 2022.

The validation of this model version was first carried out by ensuring that the simulated surface chlorophyll-a concentration, primary production, and nutrient distributions were reasonably simulated, as in previous model versions74. Second, we examined the RMSE relative to the SOCAT gridded fCO2 observations and the temporal variability of the ocean CO2 sink between 1990–2022. These two variables are used to evaluate GOBMs in the global carbon budget analysis. The RMSE value associated with NEMO-PlankTOM12.1 (38.5 μatm) is within the range of GOBMs used in the global carbon budget analysis (31.3 μatm – 45.0 μatm). The simulated temporal variabilities of the ocean CO2 sink from NEMO-PlankTOM12.1 are also comparable to the other GOBMs of the global carbon budget (Supplementary Table 1), and the simulated mean ocean CO2 sink in the 1990s (1.91 Pg C yr−1) falls within the observational range (1.5 to 2.9 Pg C yr−1)3.

Finally, the performance of NEMO-PlankTOM12.1 was evaluated with the metrics adopted by the global carbon budget in 2023: the simulated Atlantic Meridional Overturning Circulation, Southern Ocean sea surface salinity, the Southern Ocean stratification index, and surface ocean Revelle factor. The values simulated by NEMO-PlankTOM12.1 are within the range of the values simulated by the other GOBMs and are close to the observed values, with the exception of the Southern Ocean stratification index for which NEMO-PlankTOM12.1 has the lowest value but remains within comparable range (Supplementary Fig. S3).

Hybrid approach

A hybrid approach is developed to constrain the annual ocean CO2 sink simulated by the NEMO-PlankTOM12.1 model on the basis of the model-observation mismatch for surface fCO2. This approach is not implemented to significantly improve the model-observation mismatch, but to correct for annual biases in the simulated ocean CO2 sink, after the standard simulation is done, thus providing an adjusted annual estimate with uncertainty. This methodology has been previously used to constrain the climatological ocean primary production33 and air-sea fluxes of N2O34 and CCl435. Because there are more observations available for surface fCO2 observations than for N2O or CCl4, this hybrid approach can be performed annually.

The surface fCO2 observations used here are the ones compiled within the SOCAT v2023 database75. This database is a gridded product (1°×1°) with a monthly temporal resolution. All monthly model outputs used here were regridded to the same spatial resolution as SOCAT.

First, to perform this hybrid approach, four perturbed simulations of NEMO-PlankTOM12.1 are produced with higher average MSEs (between the model and SOCAT observations, with equal weight given to each gridded observational data) over the simulated period, and lower or higher annual ocean CO2 sink. These perturbed simulations range from 1.3 Pg C yr−1 to 3.2 Pg C yr−1 on average during 2000–2022, and span the expected range suggested by the global carbon budget analysis (i.e., 2.6 ± 0.4 Pg C yr−1 on average, with individual years ranging from 1.8 Pg C yr−1 to 3.0 Pg C yr−1). They are obtained by perturbing model parameters, changing the half-saturation constant for bacteria remineralisation of organic carbon from 5×10−6 mol L−1 to 18×10−6 mol L−1. This parameter was chosen because of its strong influence on the ocean CO2 sink that is relatively uniform over the entire ocean. Second, for each year, a plot of the annual MSE values (on the y axis) against the annual ocean CO2 sink (x axis) associated with the four model simulations (the optimal simulation and the four perturbed simulations) is produced and a cubic function is fitted through these data points. Third, the constrained annual ocean CO2 sink is estimated by finding the local minimum (turning point) associated with the concave upward section of the fitted cubic function. This local minimum corresponds to a theoretical model simulation with an annual ocean CO2 sink that presents the smallest MSE (MSEmin). Note that if the fit did not have a local minimum the ocean CO2 sink from this year is not constrained. Years with a constrained ocean CO2 sink not within the range of the ocean CO2 sink from the perturbed simulations are kept but considered uncertain values.

Finally, the ± 1σ (68%) confidence interval associated with the determined MSEmin value is estimated with this formula33,

MSE68% = 0.468 × (n - 2) × MSEmin / (n - 4) + (n - 2) (2)

where n is the number of gridded observational data points, and MSE68% corresponds to the MSE value for theoretical model simulations located at the borders of the confidence interval, and their associated annual ocean CO2 sink is estimated by using the cubic function (i.e., determined where f(x) = MSE68% with x being the annual ocean CO2 sink). See Supplementary Fig. S4 for a graphical interpretation of this hybrid approach. The global performance of the hybrid approach is evaluated with the two standard metrics used by the global carbon budget annual analysis2—i.e., the RMSE relative to the SOCAT fCO2 observations and the estimated temporal variability of the ocean CO2 sink (results in Supplementary Table 1).

The hybrid approach reduced the RMSE between NEMO-PlankTOM12.1 and SOCAT observations (the RMSE values associated with the perturbed simulations, PlankTOM12.1 and the hybrid approach are 39.9 μatm, 38.5 μatm and 38.0 μatm, respectively). For comparison, the GOBMs and fCO2-products listed in the global carbon budget analysis in 2023 had mean RMSE values of 39.0 μatm and 20.3 μatm, respectively. Note that fCO2-products were based on SOCAT observations and are therefore not independent, which explains their lower RMSE values.

This hybrid approach is also used to constrain regional ocean CO2 sink values at three latitude bands, and at every 5° of latitude. For this regional analysis, the number of available SOCAT observations is lower, therefore only the quadratic function is used to provide wider confidence intervals and reflect the lower confidence in the constrained values.

We performed sensitivity analyses to test the robustness of our results to the choice of perturbed model parameters and model configurations. The perturbed simulations were repeated with parameters of phytoplankton respiration, and with a combination of both bacterial half-saturation and phytoplankton respiration. The model configuration was changed by using ERA5 reanalysis as weather-forcing data. In total, we thus have applied the hybrid approach to six different setups, with three choices of perturbation parameters and two choices of forcing configurations (Supplementary Fig. S5). Regardless of the parameter and configuration used, the results consistently produced the lowest trend in the 1990s, and a higher trend in the 2000s than in the 2010s, although the exact trends within each decade varied (Supplementary Fig. S6 and Supplementary Table 3). We show here results of the model forced with NCEP, which has a lower RMSE (38.5 μatm) compared to the configuration forced with ERA5 (40.0 μatm), and the perturbation of the half-saturation constant of bacterial remineralisation, which produces changes in fCO2 that are more uniform across the ocean.

Finally, we also carried out an analysis to assess the influence of the annual application of the hybrid approach by considering, instead, three consecutive years. Therefore, rather than using a series of plots of the annual MSE value against the annual ocean CO2 sink (e.g., Supplementary Fig. S4), we used a series of plots of the 3-year MSE value against the 3-year average ocean CO2 sink.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Authors Nicolas Mayot, Erik T. Buitenhuis, Rebecca M. Wright, Judith Hauck, Dorothee C. E. Bakker,\nand Corinne Le Quéré\nTitle Constraining the trend in the ocean CO2 1 sink during 2000-2019\nManuscript NCOMMS-23-48601\nKey results:\nMayot et al. present a hybrid method fusing processed-based models and observations to\nconstrain the trend in oceanic CO2 uptake for the period 2000-2019 to 0.30 ± 0.05 Pg C yr-1 decade-\n1. Critically, this revised trend in the ocean CO2 sink enables them to infer the terrestrial CO2 sink\nto 0.33 ± 0.12 Pg C yr-1 decade-1 over the same 2000-2019 period.\nThe hybrid model forces the state-of-the-art NEMO-PlankTOM12.1 GOBM with NCEP reanalysis\ndata and observed atmospheric CO2 concentrations. It provides an annual estimate of the ocean\nCO2 sink by constraining the NEMO-PlankTOM12.1 yearly outcome with SOCAT observations\nValidity:\nThe overall analysis appears sound, but there are several discrepancies that require better\nexplanation and justification. The hybrid results for the 1990s and 2000s react significantly to the\nfCO2 data, adjusting the NEMO-PlankTOM12.1 results into closer agreement with the fCO2 data,\nas expected. However, the 2010s trends not only appears to disregard the data, but drives the NEMO-\nPlankTOM12.1 trend 50% lower and in the opposite direction of an adjustment towards the\nobservations. This result dominates the 2000-2019 trend analysis as well, where the hybrid model\nagain adjusts the final result away from the data trend. The explanations given for this in the text is\nnot compelling – if the hybrid approach did disregard the fCO2 observations, then it should have\ndelivered a result very close to the NEMO-PlankTOM12.1 result.\nThe discussion surrounding the North and South regions and how the hybrid model weighted the\nobservations in these regions is also inconsistent. The authors clearly state that there are more\nobservations in the North, yet the hybrid model does not appear to weight them heavily as the\nhybrid results in Figs 1b and 3b largely follow the NEMO-PlankTOM12.1 lines. Yet the hybrid results\ndo track the fCO2 data and not the NEMO-PlankTOM12.1 results in the Southern region (Figs 1d, 3d)\ndespite the fact that these data are sparser and potentially biased due to undersampling (Section\n4). This needs to be explained with greater clarity since the high latitude regions dominate the\nglobal results (Figs 1a, 3a)\nSignificance:\nThe results of this study will have a major impact on the understanding and quantification of the\nglobal ocean CO2 sink and its implications for quantifying the terrestrial CO2 sink. It should be\npublished after revisions.\nData and Methodology/Analytical Approach:\nThe approach and data use appear reasonable. Previous reviews have in principle vetted the hybrid\napproach used here (Refs 23-25).\nSuggested improvements\nImprovements are suggested below in the detailed line by line comments. In particular, there are\nseveral important results and conclusions that the authors should consider adding to the Abstract.\nDetailed Comments:\nLine 10 The study period is 2000-2019, but the first decade mentioned in the abstract is the 1990s.\nThe report should include the entire period analyzed, unless strong justification is provided for\nfocusing on a particular sub-period.\nLines 25-26, Ref 5 “However, the ocean carbon reservoir is also sensitive to climate variability and\nclimate change”\nThe outcome reported by McKinley et al. [2020] - that there will be an immediate reduction in ocean\ncarbon uptake as atmospheric pCO2 responds to cuts in anthropogenic emissions - seems more\nrelevant to the theme of this paragraph than the review paper by Gruber et al since the McKinley\nstudy deals explicitly with the timescale and magnitude of the ocean response to changes in\nanthropogenic CO2 emissions\nMcKinley, G.A., Fay, A.R., Eddebbar, Y.A., Gloege, L. and Lovenduski, N.S., 2020. External forcing\nexplains recent decadal variability of the ocean carbon sink. Agu Advances, 1(2),\np.e2019AV000149.\nLine 56 “…last two decades…”\nLines 56-59, Table 1 “the two types of estimates diverge over the last two decades (2000-2019)19–\n21, with the fCO2-product ensemble suggesting a decadal rate of growth of the ocean CO2 sink of\n0.68 ± 0.13 Pg C yr-1 decade-1, double the 0.35 ± 0.05 Pg C yr-1 decade-1 simulated by the GOBMs\nensemble (Table 1)”\nNote: the hybrid model derived trend, 0.30 ± 0.05 Pg C yr-1 decade-1, is very close to the GOBM\nensemble value of 0.35 ± 0.05 Pg C yr-1 decade-1, which makes one question how much influence\nthe observations have on the hybrid model result, particularly since one would expect the\nobservation-based trend of 0.68 ± 0.13 Pg C yr-1 decade-1 to drive the hybrid model result higher\nthan the GOBM ensemble result, not lower.\nTable 1 The first 2 columns provide information on the 1990-2019 trends, but these are largely\nignored in the analysis\nTable 1 The -50% adjustment that the Hybrid approach makes to the decadal trend in the 2010s\n(0.26 → 0.12 ± 0.12) makes no sense given the fCO2 value of 0.79 ± 0.41. One would expect a\npositive adjustment of ~0.4 in line with the adjustment that the method made in the 2000s (0.31 →\n0.70 ± 0.12 with fCO2 value of 0.81 ± 0.40). This must be justified as the 2010s trend also appears\nto distort the 2000-2019 trend significantly lower than the NEMO-Plank estimate and the fCO2\ntrend\nLines 185-194 do not justify or even explain this result, they only note it.\nLine 92-93 “… by constraining the model output fields of fCO2 against the observed fields”\nThe authors should quantify the strength of the fCO2 constraint given the sparse nature of the\nobserved fCO2 fields (Lines 50-51, Ref 16) and the fact that machine learning has been used to\nextrapolate the observations. How much of the hybrid result is due to the actual observations vs the\nextrapolation method?\nLines 93-96 As above, it is unclear how much the ocean CO2 sink estimate is constrained by the\nactual (sparse) observations vs the machine learning method used to extrapolate fCO2 fields\nLines 145-47, Fig 1 “The hybrid approach substantially modified the simulated ocean CO2 sink in\nthe high latitude regions, with the tropics remaining similar to the original NEMO-PlankTOM12.1\nresults”\nFig 1d shows that the Southern oceans (Hybrid – NEMO-Plank TOM ~ 0.5 PgC yr-1) drive most of the\nglobal difference between the NEMO-Plank TOM and the hybrid results (Fig 1a); the Northern\noceans contribute only ~0.1 PgC yr-1 to the model-hybrid difference (See lines 157-158)\nFig1 Suggest plotting Figs 1b-1d with the same y-axis scale (-0.5 to 2.5) to illustrate the magnitude\nof the fluxes in each domain and their relative importance to the global result (Fig 1a). This will also\nemphasize the degree of agreement between NEMO-Plank TOM and the Hybrid model in each\nregion\nLine 152-53 It looks like the Hybrid approach increased the Northern ocean CO2 sink by ~0.1 Pg C\nyr-1, not 0.01 PgC yr-1 (Fig 1b)\nLines 157-158 Suggest adding “the mismatch between the model and the fCO2 observations in the\n35-55°S band is the main cause of the underestimation of the ocean CO2 sink in the South and\nglobally” to the abstract\nFig 3 As with Fig 1, suggest plotting Fig 1b-d with the same y-axis\nFig 3b,d Why does the hybrid product track the fCO2 values in the Southern region (Fig 3d) but not\nthe Northern region (Fig 3b), where it tracks the NEMO-Plank TOM results instead? This appears to\nbe a breakdown of the hybrid approach, especially given the high density of observations in the\nNorthern region (Lines 200-202)\nFig 5, Line 208 Arctic Ocean not included in analysis\nLines 247-250 Deficiencies of the GOBMs in the high latitudes are a major result of this analysis\nand should be emphasized more. These deficiencies also demand increased sampling in the high\nlatitudes to confirm the hybrid results presented here and correct GOBM deficiencies.\nLines 273-74 “Our results suggest that fCO2-products overestimate the decadal trend of the 2010s\nin regions where there are few observations”\nThis suggests that the hybrid method defaulted to the NEMO-Planck TOM model in undersampled\nhigh latitudes in both hemispheres; however, there should have been some adjustment made for\nthe actual observations even if they are sparse.\nLines 281-84 The authors appear to assume that the observations in the Southern Ocean are\nbiased based on Refs 31 and 33, but the hybrid model locks onto the measurement data (Fig 3d). If\nthe observations are biased, then why trust the results of the hybrid method in the Southern region?\nA similar argument about biased observations in the North carries less weight given there are more\nobservations in the 30-90N latitude range than the 30-90 S range. This is a significant inconsistency\nthat is not clearly resolved in the text.\nLines 293-95 “GOBMs represent the correct processes, but either they do not respond sufficiently\nto changes in external forcing, or the balance among thermal and non-thermal processes in\nresponse to that forcing is imperfect.”\nThis is a general statement, but should address NEMO-Planck TOM specifically given its central\nimportance in the hybrid method\nLines 320-322 Suggest adding “the stronger trend in the 2000s and the weaker trend in the 2010s\nestimated with the hybrid approach help to reduce the trend associated with the carbon budget\nimbalance between 2000 and 2019 from -0.24 Pg C yr-1 decade-1 to -0.003 Pg C yr-1 decade-1” to\nthe Abstract\nLine 323 How sensitive is the estimated trend in the terrestrial CO2 sink to the value of the trend\nassigned to the ocean CO2 sink?\nLine 331 The authors miss the opportunity to advocate for more and more extensive fCO2\nobservations, particularly at the high latitudes. Strongly urge that a final paragraph be added to\nSection 4 that addresses research needs unresolved by the current study.\nLines 389-90 “Model simulations are too short to fully represent the input of river fluxes and their\nsubsequent outgassing of CO2 in the open ocean”\nThis is a potentially serious model deficiency as riverine inputs will significantly alter coastal CO2\nuptake and will have a measurable impact on the global ocean CO2 sink", "reviewer_2": "The manuscript addresses an important point regarding the disagreement between different\nestimates of ocean carbon and shows that data availability may be one of the culprits for the\ndiscrepancy.\nHowever, in the current state, the manuscript is focused too much on technical details and lacks\nbroader context, so it may be more suitable for a more methodology-oriented journals, such as\nGlobal Biogeochemical Cycles. Below are some of points that would need to be addressed to\nbroaden the context of the manuscript before I can recommend it for publication in Nature\nCommunications.\n1. Why was the comparison only done with NEMO-PLANK model? There are other (higher-\nresolution) ocean circulation models for physics-biogeochemistry coupling, for example, ECCO-\nDARWIN (Carroll et al, 2020, JAMES). This model is higher resolution than NEMO-PLANK and is\ndata-assimilating. I suggest the authors compare how the global estimates from their models\ncompare with ECCO-DARWIN and whether their hybrid approach performs better than ECCO-\nDARWIN.\n2. As Nature Communications is for a broader audience, the authors should discuss the meaning of\nfCO2 and the limitation of this metric in relation to ocean carbon storage in the introduction.\nSpecifically, fCO2 is an ocean-surface measurement and does not represent the changes in ocean\ncarbon content at depth, importantly below the permanent pycnocline. In this context, the authors\nshould discuss recent studies done for the ocean carbon storage in 3D space, e.g., Gruber et al\n(2019, Science), Zemskova et al (2022, Nat Comm).\n3. The manuscript contains many acronyms and numbers, which deteriorates the readability of the\nmanuscript. For example, are the acronyms A-IV and A-DV really necessary? The authors could just\nsay \"interannual\" and \"decadal\" variability in words. I also suggest that the authors comb through\nthe manuscript to limit the number of numbers that they present to those that are crucial to the\ncentral points of the paper. These points need to be emphasized, as they are currently lost in the\nplethora of values presented in text. This trimming would make the paper more suitable for a\nbroader higher-impact journal like Nature Communications rather than a methodology-oriented\njournal.", "reviewer_3": "Mayot et al. review, November 2023\nThe manuscript by Mayot and coauthors evaluates discrepancies between the surface fCO2 and\nair-sea flux of carbon in earth system models and observation-based data products. The authors\nuse an optimization approach to identify estimates that are as consistent as possible with both\nobservational and model constraints.\nI think this manuscript will be of interest to the readership and appropriate in scope for Nature\nCommunications. It is also timely and well written. I do have concerns about how particular\nmethodological choices and resulting caveats are communicated. In general I think my concerns\ncan be resolved by in text clarifications, though I do recommend inclusion of some diagnostics\nwhich may be more time consuming to produce.\nMy major comments are as follows (in order of my perception of priority):\nA. The authors chose to exclude the coastal ocean from their analysis--everything within 3 degrees\nof the coast and shallower than 1000 m. Given the outsize role of these regions for marine carbon\nfluxes, I believe the quantitative consequences of that choice should be clarified and visualized for\nreaders.\nB. In my opinion, this paper would be more impactful and robust if the authors included the results\nof sensitivity analyses used to inform their approach. Specifically:\ni. Several possible mechanisms are proposed to drive discrepancies between models and data-\nproducts. Varying subsampling and model parameters would help identify the most likely of these\nmechanisms.\nii. Presenting the results of the sensitivity analysis used to select and tune key parameters for\noptimization and perturbation--within the constraints identified in the literature--would increase\nconfidence in the results.\nC. I think the model manipulations described here constitute a sufficiently different approach\nand parameter set that validation should be presented. e.g., an evaluation of how major variables like\ntemperature and salinity, as well as fCO2, differ between this and prior PlankTOM configurations\nand climatological hydrography.\nD. Various results are compared, but statistical tests are not presented. Including these would help\nsupport interpretations of whether or not values are similar or meaningfully different.\nE. I think the nudging/data assimilation approach could be better referenced to illustrate how it\ndraws on or differs from prior literature.\nF. The authors appear to use different sources and calculations for water and air fCO2, instead of\nthe concurrently and self-consistently calculated values from SOCAT. I recommend the authors\nevaluate whether it is possible to consolidate their data sources in this case.\n___\nDetailed comments by line number (because this is a long list, I don't expect responses to all of\nthese).\n52. It would be helpful to further distinguish what is meant by fCO2-products here in contrast to\nobservational data only. Are the authors only referring to gridded, interpolated, and climatological\ndatasets that use different methods but all fundamentally rely on SOCAT? Are there additional\nsubstantive differences readers should be aware of?\n57. What are the specific fCO2 products being considered as part of the ensemble? And the\nspecific GOBMs included in the ensemble? I think each product considered should be referenced\nand included in a table associated with this manuscript. If both ensembles are identical to those\nused in the Global Carbon Budget 2022 publication, the I think that this point should further be\nclearly stated in the inline text with an associated reference number.\n65.\na. The result in the 2010s trend column--in which the result is significantly lower than the model\nensemble, preferred model, or fCO2-product ensemble and thus seems to contradict either\nobservational or model evidence--strikes me as a red flag that requires clearer explanation in the\ntext. I'm not sure section 3.4 is sufficient to convince me. Any proposed resolution would be more\nconvincing with a model sensitivity study demonstrating a plausible mechanistic basis.\nb. The larger interannual variability in the presented approach than the other models and data\nproducts is also concerning, though less so given the reported uncertainties.\nc. I suggest a concise table title, rather than expanding this into a caption.\n85. Given that this approach is so central to the analysis, I think a brief overview here (a couple\nmore sentences), prior to the trailing methods section, would be appropriate.\n90. The trends in the preferred model are 20-60% lower than the ensemble in each decade. Is it fair\nto call these similar trends? Or perhaps it makes sense to restate this as a similar multidecadal\ntrend, since the decadal trends diverge. I also suggest double checking the overall 2000-2019\nNEMO trend--based on Fig. S1 it appears that either the overall trend or detrended decadal\ncomponent trend should be a bit smaller than that of the ensemble (but I may be mistaken).\n96. A caveat to this comparison to is that these prior studies performed sensitivity analyses across\na wide range of parameters and included observational constraints (of parameters other than the\nobjective variable) on the the tuned parameters. If this study does either, it was not obvious to me.\nSo I think it is appropriate to include some comment on how those prior studies guided the authors\napproach, and any limitations.\n103. Similar to my comment about line 85, I think the authors should introduce a concise\nsummary of the perturbation and tuning approaches here, beyond stating that they occurred. This\ncontext seems important for helping readers understand why state estimates sometimes arise outside\nthe bounds established by the perturbations.\n113. This is a matter of personal preference, but now that the authors have defined this variability in\nterms of amplitude, could the acronyms be shortened to IV and DV? Or, see comment for line 133.\n115. \"Very close\" = indistinguishable? Considering these as a one sided test of difference\n(with table of tests results in supplement?) might clarify the meaning of statements about\nsimilarity ('very close', 'close to', 'marginally overlapping', 'in line with', 'within the range of',\netc.). For consistency of interpretation by readers, the authors might consider choose a\nsingle phrase to indicate not significantly different comparisons.\n124. I think it would be useful to label the perturbed simulation lines with a descriptive\nlabel of how they were perturbed (e.g., microbial K1/2 = ...). Perhaps this is semantic, but in\nthe final sentence of the caption I think \"constrained but outside the bounds of the\nperturbation experiments\" might be more descriptive than 'uncertain constrained'.\n133. Wouldn't a moving-average as described here often result in lower variability even without\nchanges in the data distribution? I suggest that smoothing the data in this way is no longer a\ncomparison of interannual variability in the same sense as it was previously introduced. To be very\nclear, the authors might decide to replace A-DV and A-IV, with a notation that clearly displays the\ntimescale considered. e.g., V_{10y}, V_{3y}, V_{1y}.\n156.\na. A mismatch may help diagnose mechanisms, but it is not itself a cause of anything. Nor does the\npresence of a mismatch--absent a clearly aligned expectation grounded in a particular\nmechanism--conclusively indicate that the southern mid-latitude CO2 sink is underestimated in\nmodels as opposed to overestimated in sparse observations. I think this statement needs to be\nrevised--perhaps the authors meant instead that the overall regional and global disagreements\nbetween data sources reflect these latitudes in particular rather than elsewhere.\nb. At face value, this finding seems to contradict that of Fay et al. 2021 (cited by the authors\nas reference 27) with respect to the geographic specificity of model-data mismatches.\n160. Are the uncertainty bounds smaller than the symbols on most of these bands? Does it make\nsense that they are this small given the large interannual and interdecadal variability previously\ndiscussed?\n194. See comment for line 65a.\n204. If relying on significance for argument, I suggest including test result(s) (null, p,\ntest-statistic).\n216. I think that this statement implicitly assumes a constant gas transfer velocity. I'm not\nfamiliar enough with global windspeed or momentum flux trends to know whether that is\nlikely but there are observed and modeled changes in these parameters over the last few\ndecades in at least some regions within this latitude band (e.g., Eastern boundary current\nregions). The assumptions here probably need some reinforcing citations.\n219. a. Please include trends and relevant statistical tests to reinforce this conclusion.\nb. To me, this distinction is a bit arbitrary (or at least reflect a choice of emphasis). Visually,\nany given decade seems unlikely to diverge greatly from the overall 1990-2020 trend. i.e., if\ncomparing each decade to null hypothesis of no difference from the secular trend, are any of these\ndecadal variations actually significant?\n225. a. I think this point about spatial bias in the interpolated products is valuable and it would be\ngood to start a new paragraph here to visually highlight this idea.\nb. This point (bifurcation of the interpolation products into high and low clusters) would be\nstrengthened by a supplemental figure.\n242. My perception is that statistical approaches to 'nudging' earth, ocean, and climate model\nresults to incorporate observational constraints has a long history, which should probably be\nreferenced here. Examples and context include: http://dx.doi.org/10.5065/D60Z716B,\nhttps://doi.org/10.1175/1520-0493(2000)128>3664:ASNTFD<2.0.CO;2,\nhttps://www.myroms.org/forum/viewtopic.php?t=320, https://doi.org/10.5194/gmd-16-1857-2023.\nThe author's will note that these sources point out that nudging to observations is sometimes\nconsidered a suboptimal data assimilation strategy. I recognize that the model optimization used\nhere is more complicated, and may not be well summarized as \"nudging\".\nOther works have used more data assimilation techniques to incorporate marine carbonate system\nvariables in concert with physical parameters: https://doi.org/10.1016/j.jmarsys.2016.02.011,\nhttps://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JC012650,\nhttps://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2009GB003531,\nhttps://www.jstor.org/stable/24860162, https://www.jstor.org/stable/24861004.\nGiven this context, I think the authors need to not just reference related concepts and works of this\nnature, but more clearly delineate how this work is novel in approach and focus. I acknowledge that\nthe length-limitations and format of this journal may have obscured significant advances in the\napproach used by the authors--if so and I missed these, I take that as further evidence that these\nneed to be clearly and concisely delineated in the text.\n258. Please see comment regarding line 156a. Without testing one of the mechanisms proposed, I\nthink that this statement is not fully supported as written.\n324. This is unclear to me as written. I suggest including an equation and calculated fluxes (and\ntheir data sources) in the supplement.\n348. SOCAT reports both the wet and dry atmospheric fCO2 measured concurrently with the\nseawater fCO2 selected by the authors, eliminating the need for this calculation and the potential\nerrors in introduces because of discrepancies in choices of equation and virial coefficients used by\ndifferent groups, differences in the spatial gridding, etc. I strongly recommend the authors revisit\ntheir data source and this calculation considering the alternative SOCAT availability--if this is not\npossible, please explain why (e.g., are there no comparable interpolated products for the SOCAT air\nfCO2?).\n363. Randomly selected with what distribution? I ask because the authors previously chose to use\ntests that are suitable for non-normal data. If the data is not well approximated by a normal\ndistribution, that should be noted here and an appropriate distribution chosen for this Monte Carlo\napproach.\n385. Some additional details about the treatment of major inorganic carbon components would be\nappropriate here, even if covered in the reference.\n398. a. If re-optimized as described, this is a new model variant that really requires independent\nvalidation from prior published works. Did this process knock any parameters back out of\nreasonable ranges? Are basic properties like T & S well described in at least a climatological sense?\nI acknowledge model validations are a lot of work, but as a reviewer it is difficult for me to have\nconfidence in the results presented here without validation of parameters beyond fCO2 of this new\nproduct presented in the supplement.\nb. This optimization scheme is empirical but it isn't clear to me that the optimized parameters are\nrealistically constrained. For example, if I am interpreting the temperature sensitivities\ncorrectly (e.g., based on parameter description by Le Quere 2016), that's equivalent to a Q10 values that are\nlower than typically measured or modeled on either the organismal or ecosystem level (usually 2-3,\nthe values here are equivalent to 1-2). Another example is the use of a constant e-folding scale for\nlight, even though this varies widely in the ocean. Things like this make me question whether\nappropriate constants are being used and if the model optimization is able to distinguish between\ncompensatory drivers.\nc. Also, Table S1 does not include the 47 parameters that the text says were optimized. I think that if\nother variables were unchanged (I think only changed parameters have been tabulated?), then the\ntext should be updated and the table clear on what is included and why.\n403. Please elaborate on why these parameters were chosen, at least in a supplemental text. That\nmakes these choices more transparent, and helps other modeling teams understand which\ndifferences are important.\n410. I would prefer representative results to back up this claim were included.\n422. But Table S1 appears to show changed physical parameters as well...please resolve this\n(appearance of) inconsistency.\n425. Is this different from the A-IV discussed earlier? If so please consider removing the new\nacronym, as it is only used a couple times in the methods.\n456. What is the order of magnitude of observations required to do this robustly?\n460. a. This is a choice that I think should be clearly communicated in the main text, and\nas a mask or outline in the maps in Fig. 5 and Fig. S3. It is unclear from this whether the\nmodel outputs and data products are similarly masked out in the same coastal band, so\nplease clarify that these analyses are comparing like to like.\nb. Areas of <1000 m depth and within 3 degrees of shore are disproportionately important for\nmarine carbon cycling. Ignoring them in order to characterize more pelagic regions may make it\neasier to identify trends and present more constrained appearing values, but at the cost of\nconsidering real, evnironmentaly relevant variability and significant contributions to the global\ncarbon flux. Reading this, I have diminished confidence in claims about the changed magnitude of\nthe climatological flux (and inferred land flux), trends (which may be compensated or exacerbated\nby coastal processes), and the relative agreement of variability estimates. I think that this is too\nimportant a caveat to be buried in the Methods and not discussed in the Discussion.\nc. How do the results change if coastal areas are considered as well, or separately? I would like to\nsee even a quick discussion of this and followup plots in the supplement. It may be that focusing on\nthe pelagic ocean is defensible and interesting, but this would be a more compelling manuscript if\nreaders were able to consider the coasts as well--and if the authors were to point out some\npotential open questions that arise from those considerations.\n467. Is this the most relevant variable to change?\nCould the authors present sensitivity analysis supporting the selection of this parameter instead\nof/without others?\nWhat are the actual values of the perturbations used? This doesn't seem to be reported\nanywhere, but needs to be.\nIs it reasonable to assume a constant value (as shown in Table S1) or are biogeographic variations\nin this (and other constants) expected?\nWhat are the observational constraints on this parameter and how were the magnitudes of changes\nselected within those constraints?\nI think this choice needs significantly more explanation. There many fit results (in time and space)\nplotted outside the range of these perturbed models, which leads me to ask if the perturbations\nreally are providing useful bounds or not--and if they are, whether the hybrid approach would\nbenefit from additional constraints or flexibility in parameters. I understand that this is a poorly\nconstrained parameter, but there does seem to be a finite range that is consistent with measured\nocean DOC values, which could help justify bounds on the perturbations.\n470. Please provide a representative graphical example of this approach in the supplement (or a\ncouple, a good fit and a bad fit).\n477. Related to comment on line 467, I'd like to see some more quantitative assessment of the\nstatement here.\n486. There is no Fig. S4 in the supplement provided for review, but this would be welcome and\nwould address my request for line 470 as well.\n501. What quadratic function? Is there a mathematical justification for a different calculation? i.e.,\nis this an arbitrarily larger confidence bound?\nFig. S1. The letters above the last two panels are repeated.\nFig. S3. I suggest making these panels larger.\nIn the data link, the hybrid results have a mislabeled column 1 in the tab for results by latitude.\nThe code looks plausible (I did not not attempt to run it), and additional annotation would be\nappreciated to help users looking to translate this to an open source programming language for\ntesting.", "reviewer_author_response": "Dear Editor,\nThank you for inviting a revision of our work for further consideration in your journal. We\nsincerely thank the reviewers for their invaluable comments on our manuscript. Our revision\naddresses all the points raised by the reviewers and reduces the use of acronyms and\ntechnical language to enhance accessibility.\nIn particular, following your recommendation and the reviewer’s comments, we have made\nthe following important improvements:\n● We conducted and detailed a sensitivity analysis of the hybrid approach, which\nshows that results on decadal trends are robust to the choices of parameters and\nmodel set up.\n● We have included a validation of our model based on the model evaluation process\nused in the Global Carbon Budget 2023, which demonstrates that the\nNEMO-PlankTOM12.1 model is appropriate for estimating the evolution of the\noceanic CO sink, and that it is largely in line with the results of other models.\n2\n● We have improved the description of the hybrid approach, both in the Methods and in\nthe main text.\n● We have used all available fCO observations including those from the coastal zones\n2\n(which were excluded in the original submission).\n● We have clarified the presentation of the results about the North and South regions.\n● We have included a new figure that illustrates the low robustness of the 2010s trend\nin the ocean CO sink estimated with fCO products, giving even more weight to our\n2 2\nown conclusions using our hybrid approach.\nFurthermore, since the submission of our manuscript, the Global Carbon Budget 2023 (GCB\n2023) has been published. We thus updated our analysis to use the latest available\nestimates for the NEMO-PlankTOM12.1, the Global Ocean Biogeochemistry Models\n(GOBMs) and the fCO -products published in GCB 2023, and expanded the period of\n2\nanalysis to year 2022. Our main focus is now on constraining the ocean CO sink trend\n2\nduring 2000-2022 (instead of 2000-2019 in the original submission). Finally, we limited the\nuse of the hybrid approach to start in 1980 as the number of fCO observations prior to that\n2\nwere extremely low and discontinued.\nAll these modifications strengthen our results and improve the readability of our manuscript.\nSee below our point-by-point response to each reviewer, with the reviewer's comments in\nblue, the new text in red and our response in black.\nResponse to Reviewer #1:\nThe overall analysis appears sound, but there are several discrepancies that require better\nexplanation and justification. The hybrid results for the 1990s and 2000s react significantly to\nthe fCO2 data, adjusting the NEMO-PlankTOM12.1 results into closer agreement with the\nfCO2 data, as expected. However, the 2010s trends not only appears to disregard the data,\nbut drives the NEMO-PlankTOM12.1 trend 50% lower and in the opposite direction of an\n1\nadjustment towards the observations. This result dominates the 2000-2019 trend analysis as\nwell, where the hybrid model again adjusts the final result away from the data trend. The\nexplanations given for this in the text is not compelling – if the hybrid approach did disregard\nthe fCO2 observations, then it should have delivered a result very close to the\nNEMO-PlankTOM12.1 result.\nFor clarification, the hybrid approach uses directly the fCO observations from the SOCAT\n2\ndatabase, not the interpolated and extrapolated version of these observations proposed by\nthe fCO -products. Therefore, there is no assumption that the obtained estimates from the\n2\nhybrid approach for the ocean CO sink should follow the estimates from the fCO -products\n2 2\nat global or regional scale. A difference between the hybrid approach estimate and the\nfCO -product estimate could, for example, be linked to imperfect interpolation and\n2\nextrapolation of fCO observations by fCO -products. We also would like to mention that the\n2 2\nhybrid approach is applied similarly for each year, meaning that it cannot disregard the fCO\n2\nobservations for a specific period. A graphical interpretation of the hybrid approach method\nis presented in Fig S4 (this figure was regrettably missing in our first submission).\nThe discussion surrounding the North and South regions and how the hybrid model weighted\nthe observations in these regions is also inconsistent. The authors clearly state that there\nare more observations in the North, yet the hybrid model does not appear to weight them\nheavily as the hybrid results in Figs 1b and 3b largely follow the NEMO-PlankTOM12.1 lines.\nYet the hybrid results do track the fCO2 data and not the NEMO-PlankTOM12.1 results in\nthe Southern region (Figs 1d, 3d) despite the fact that these data are sparser and potentially\nbiased due to undersampling (Section 4). This needs to be explained with greater clarity\nsince the high latitude regions dominate the global results (Figs 1a, 3a)\nAs mentioned above, we are not expecting that the hybrid approach results should\nexclusively follow the results from the fCO -products. In section 3.4, which has been revised\n2\nby modifying figure 4, we identified discrepancies in the North between fCO -products,\ninterpolation and extrapolation of fCO observations by fCO -products in under-sampled\n2 2\nnorthern regions. For the South, in section 3.2, we demonstrate the deficiency of\nNEMO-PlankTOM12.1 in simulating a correct ocean CO sink in this region.\nbeen revised. For the North, our suggested explanations are now clearly stated and we have\nadded information on the potential influence of coastal regions. For the South, a new\nseparate paragraph has been added, in which we clearly recognise the uncertainties\nassociated with the lack of fCO observations (new text in red): “In the northern latitudes,\n2\nwhere the availability of measurements is highest, the fCO2-product ensemble gives a\ndecadal trend in the 2000s not significantly different from that of the 2010s (Fig. 4). Four\nfCO2-products suggest a growing or strong trend during the 2010s, contrary to the fCO2\nobservations alone, which is explained here by their strong trends in areas that were\nundersampled during the 2010s (Fig. 5). Hence, we hypothesise that methodological issues\nin some fCO2-products could lead to an unrealistic amplification of the ocean CO2 sink trend\nin the 2010s (Bennington et al., 2022a). In addition, the ocean CO2 sink in the northern\nregion is also more influenced by coastal phenomena than the southern region, which\ndespite their importance remain largely uncertain (Dai et al., 2022; Regnier et al., 2022;\nregion, the way in which coastal fCO2 observations are taken into account by the various\nfCO2-products could induce some of the discrepancies among fCO2-products. This would\npartly explain the lack of coherence between the GOBMs and the fCO2-products over this\nrecent decade (Friedlingstein et al., 2023; Gruber et al., 2023; Keppler et al., 2023; Müller et\nal., 2023).\nIn the Southern Ocean, our hybrid approach suggests that existing fCO2 measurements\ncould corroborate a strong and positive decadal trend in this region in the 2010s, and more\ngenerally between 2000-2022. However, the paucity of fCO2 measurements impede our\nability to evaluate the decadal trend using observations only (Gloege et al., 2021), but recent\nstudies showed that undersampling could be responsible for strong biases in fCO2-products\nin that region (Hauck et al., 2023a).”\nAnd in the second to last paragraph: “Moreover, it [our analysis] confirms the importance of\nhigh-density fCO observations, which are notably lacking in the Southern Ocean, for\n2\ninforming the fCO -products and our hybrid approach.”\n2\n1) Line 10 The study period is 2000-2019, but the first decade mentioned in the abstract is\nthe 1990s. The report should include the entire period analyzed, unless strong justification is\nprovided for focusing on a particular sub-period.\nOur main study period is 2000-2022 (expanded to 2022 to take account of the latest models\nand fCO -products made available in the global carbon budget 2023 update). The\n2\njustification of this choice is that it is during this period that ocean CO sink estimates differ\n2\nthe most between GOBMs and fCO -products: as reported in the global carbon budget\n2\nanalysis 2023 (Friedlingstein et al., 2023 - doi.org/10.5194/essd-15-5301-2023), ocean CO\n2\nsink trends from GOBMs and fCO -products have diverged by a factor of 2 since 2002, and\n2\nby a factor of 2.5 since 2010. Furthermore, the number of fCO observations considerably\n2\nincreased from the early 2000 (the number of monthly grid cells observed per year more\nthan doubled between 2000 and 2005), making the discrepancy between various estimates\nmore puzzling. Therefore our analysis helps to resolve this divergence.\nThe justification sentence of our introduction has been expanded to briefly explain the\nrationale for the foci of the paper, as follows: “Here, we investigate the discrepancies\nbetween fCO -products and GOBMs, focusing on the inconsistency, by a factor of two, of the\n2\n2000-2022 trend, which occurs despite the growing number of fCO observations.”\n2\n2) Lines 26, Ref 5 “However, the ocean carbon reservoir is also sensitive to climate\nvariability and climate change” The outcome reported by McKinley et al. [2020] - that there\nwill be an immediate reduction in ocean carbon uptake as atmospheric pCO2 responds to\ncuts in anthropogenic emissions - seems more relevant to the theme of this paragraph than\nthe review paper by Gruber et al since the McKinley study deals explicitly with the timescale\nand magnitude of the ocean response to changes in anthropogenic CO2 emissions\nMcKinley, G.A., Fay, A.R., Eddebbar, Y.A., Gloege, L. and Lovenduski, N.S., 2020. External\nforcing explains recent decadal variability of the ocean carbon sink. Agu Advances, 1(2),\np.e2019AV000149.\natmospheric growth rate, as follows: “This uptake is modulated by the variable growth rate in\natmospheric CO (McKinley et al., 2020). Furthermore, the ocean carbon reservoir is also\n2\nsensitive to climate variability and climate change (Gruber et al., 2023)...” The Gruber et al.\n(2023) paper is a review paper that summarises the factors responsible for trends and\nvariability in the ocean CO sink, which includes also those of Mckinley et al. (2020).\n2\n3) Line 62 “…last two decades…” Done\n4) Lines 61, Table 1 “the two types of estimates diverge over the last two decades\n(2000-2019)19–21, with the fCO2-product ensemble suggesting a decadal rate of growth of\nthe ocean CO2 sink of 0.68 ± 0.13 Pg C yr-1 decade-1, double the 0.35 ± 0.05 Pg C yr-1\ndecade-1 simulated by the GOBMs ensemble (Table 1)”\nNote: the hybrid model derived trend, 0.30 ± 0.05 Pg C yr-1 decade-1, is very close to the\nGOBM ensemble value of 0.35 ± 0.05 Pg C yr-1 decade-1, which makes one question how\nmuch influence the observations have on the hybrid model result, particularly since one\nwould expect the observation-based trend of 0.68 ± 0.13 Pg C yr-1 decade-1 to drive the\nhybrid model result higher than the GOBM ensemble result, not lower.\nThe fCO -products are not used by the hybrid approach, which instead uses only the fCO\n2 2\nobservations published by SOCAT. Therefore, we are not expecting that the hybrid approach\nresults should exclusively follow the results from the fCO -products. The hybrid approach is\n2\nconstrained by both the SOCAT observations and by the simulated processes affecting the\ncarbon cycle in the ocean (as represented in the NEMO-PlankTOM12.1 model). In contrast,\nthe fCO -products are closely linked to fCO observations and are also sensitive to\n2 2\nuncertainties in gas exchange parameterization and the scarcity of data. The extrapolation in\ndata-poor regions performed by the fCO -products depends on the numerical methods\n2\nchosen and does not need to respect physical processes.\nLikewise, the GOBMs are not constrained by SOCAT observations, contrary to the hybrid\napproach. The fact that the hybrid approach deviates from the trends of both the\nNEMO-PlankTOM12.1 model and the GOBM ensemble shows that the hybrid approach can\nindeed be influenced by the observations and demonstrates that the method works and does\nnot simply reproduce the results of the original NEMO-PlankTOM12.1 model.\nMoreover, in this revised submission, values in Table 1 have been updated to take into\naccount the latest information available. The value of the ocean CO sink trend, between\n2\n2000 and 2022, obtained from the hybrid approach now lies between the GOBM and\nfCO -product estimates, while still slightly closer to the GOBMs. We explain the slightly\n2\nlarger updated trend by an additional 3 years as well as additional SOCAT data covering the\nlater years, and by the use of an updated version of the NEMO-PlankTOM12.1 model which\nbetter fits the SOCAT observations. The closer proximity of the hybrid approach results to\nthe fCO -products results stems also from the update of these products within the last three\n2\nupdates of the Global Carbon Budget analysis (see the end of section 3.3 and Figure 6).\n5) Table 1 The first 2 columns provide information on the 1990-2019 trends, but these are\nlargely ignored in the analysis\nparagraphs of the discussion. Because the first two columns give information about the\namplitudes of the interannual and decadal variations, those results are more reliable when\nusing longer time series. We use these results in two ways. First to show that GOBMs likely\nunderestimate variability (mentioned in the discussion). Second as an additional element to\nsupport the validity of the hybrid approach, because the hybrid approach reproduces the\npatterns of variability produced by the NEMO-PlankTOM12.1 model (r = 0.5, p = 0.004,\nPearson's correlation coefficient) but enhances its variability (also in the discussion). We\ntherefore left two columns of the table unchanged, but have reordered the information to be\neasier to follow.\n6) Table 1. The -50% adjustment that the Hybrid approach makes to the decadal trend in the\n2010s (0.26 0.12 ± 0.12) makes no sense given the fCO2 value of 0.79 ± 0.41. One would\nexpect a positive adjustment of ~0.4 in line with the adjustment that the method made in the\n2000s (0.31 0.70 ± 0.12 with fCO2 value of 0.81 ± 0.40). This must be justified as the 2010s\ntrend also appears to distort the 2000-2019 trend significantly lower than the NEMO-Plank\nestimate and the fCO2 trend\nWe do not expect the hybrid approach to make the same adjustments between different\nyears or decades, nor to follow the results of fCO -products for the reasons mentioned\n2\nabove (the hybrid approach is constrained by both SOCAT observations and physical\nprocesses, while the fCO -products are tightly linked to SOCAT observations and\n2\nnumerical extrapolations). While the NEMO-PlankTOM12.1 model may wrongly estimate a\ntrend for one period, it could correctly estimate a trend for another period. This explains why\nthe adjustments are different between different years and decades. It is also important to\nnote that the estimate of the ocean CO sink obtained for one year by the hybrid approach is\n2\nnot used by the hybrid approach to derive the estimate for the following or previous year,\nwhich again induces different adjustments for different periods.\n7) Lines 182-185 do not justify or even explain this result, they only note it.\nWe report the most important results within the results section (including those of lines\n182-185), and further explain and discuss them in the discussion section. For example, the\ntrends in the 2010s are reported for the various methods in results section 3.3, the\ndiscrepancies are reported in detail in section 3.4 and are discussed in the second\nparagraph of the discussion. Note, we have added additional results concerning the\nevolution of the fCO -product estimate over the last three Global Carbon Budget analyses for\n2\nthe 2010s trend.\n8) Line 87-90 “… by constraining the model output fields of fCO2 against the observed\nfields” The authors should quantify the strength of the fCO2 constraint given the sparse nature of the\nobserved fCO2 fields (Lines 50-51, Ref 16) and the fact that machine learning has been used to\nextrapolate the observations. How much of the hybrid result is due to the actual observations vs the\nextrapolation method?\nLines 91-95 As above, it is unclear how much the ocean CO2 sink estimate is constrained by the\nactual (sparse) observations vs the machine learning method used to extrapolate fCO2 fields\n2\nare not influenced by the use of an interpolation and extrapolation method. The fit minimises\nthe RMSE between the NEMO-PlankTOM12.1 model output and SOCAT observations (that\nare provided on a 1x1 grid at monthly resolution). The same weight is given to all\nobservations. We have reformulated the sentence above to make this clear. The sentence\nnow reads: “...goes beyond the traditional model evaluation by constraining the model output\nfields of fCO against the observed fCO data provided by SOCAT“.\n2 2\nWe also added a precision on the treatment of weight in the Methods section 5.3: “...average\nMSEs (between the model and SOCAT observations, with equal weight given to each\ngridded observational data)”\n10) Lines 143-45, Fig 1 “The hybrid approach substantially modified the simulated ocean\nCO2 sink in the high latitude regions, with the tropics remaining similar to the original\nNEMO-PlankTOM12.1 results” Fig 1d shows that the Southern oceans (Hybrid –\nNEMO-Plank TOM ~ 0.5 PgC yr-1) drive most of the global difference between the\nNEMO-Plank TOM and the hybrid results (Fig 1a); the Northern oceans contribute only ~0.1\nPgC yr-1 to the model-hybrid difference (See lines 157-158)\nWe added a mention of the larger importance of the South in this sentence. Although the\ndifferences are not on the same order of magnitude between the different regions, it is also\nimportant to mention that the variabilities are mostly modified in the high latitudes regions\nrather than in the tropics. The new sentence now reads: “The hybrid approach substantially\nmodified the simulated ocean CO sink in the high latitude regions, particularly in the South,\n2\nbut with the tropics remaining similar to the original NEMO-PlankTOM12.1 results”\n12) Fig1 Suggest plotting Figs 1b-1d with the same y-axis scale (-0.5 to 2.5) to illustrate the\nmagnitude of the fluxes in each domain and their relative importance to the global result (Fig\n1a). This will also emphasize the degree of agreement between NEMO-Plank TOM and the\nHybrid model in each region.\nFig 3 As with Fig 1, suggest plotting Fig 1b-d with the same y-axis.\nWe modified Fig 3 as suggested. That figure now clearly shows the variability as well as the\nrelative magnitude of the three regions. For Fig 1, we modified panels 1b-c to have the\nsame y-axis scale, but kept Fig 1d unchanged, as the mean flux value is very different in the\nSouth and the presentation of error bars (also larger in the South) required a different and\nmore extended y-axis. The point about the relative magnitude is now made by Fig 3. We\nalso added a mention of the larger influence of the South region on the global results in the\ntext: “...the mismatch between the model and the fCO observations in the 40-60°S band\n2\ncould be the main cause of the underestimation of the ocean CO sink in the South (Fig. 2).\n2\n13) Line 152-53 It looks like the Hybrid approach increased the Northern ocean CO2 sink by\n~0.1 Pg C yr-1, not 0.01 PgC yr-1 (Fig 1b)\ncomments.\n14) Lines 154-156 Suggest adding “the mismatch between the model and the fCO2\nobservations in the 35-55°S band is the main cause of the underestimation of the ocean\nCO2 sink in the South and globally” to the abstract\nOur manuscript focuses on explaining the trends in the ocean CO sink. The suggested\n2\nsentence, although interesting to explain the origin of a model bias, is peripheral to the main\nmessage of the manuscript. Because of the length restriction for the abstract for Nature\nCommunication, we did not add the suggested information to the abstract.\nLines 373-375 Suggest adding “the stronger trend in the 2000s and the weaker trend in the\n2010s estimated with the hybrid approach help to reduce the trend associated with the\ncarbon budget imbalance between 2000 and 2019 from -0.24 Pg C yr-1 decade-1 to -0.003\nPg C yr-1 decade-1” to the Abstract\nThe reduction in the global carbon budget imbalance is subject to additional uncertainties\nthat are not related to our analysis, but to the estimations of the other components of the\nglobal carbon budget (e.g., land sink). We therefore did not add the suggested information to\nthe abstract. Furthermore, we would like to mention that a manuscript targeting specifically\nthis question is in preparation within our research group. It will be submitted for publication\nsoon.\n16) Fig 3b,d Why does the hybrid product track the fCO2 values in the Southern region (Fig\n3d) but not the Northern region (Fig 3b), where it tracks the NEMO-PlankTOM results\ninstead? This appears to be a breakdown of the hybrid approach, especially given the high\ndensity of observations in the Northern region (Lines 200-202)\nAs detailed in our response above, we are not expecting that the hybrid approach results\nshould exclusively follow the results from the fCO -products. For the South, in section 3.2,\n2\nwe demonstrate the deficiency of NEMO-PlankTOM12.1 in simulating a correct ocean CO\n2\nsink in this region. In section 3.4, which has been revised by modifying figure 4, we identified\ndiscrepancies in the North between fCO -products, GOBM and the hybrid approach, and\n2\nrelate these discrepancies to potential imperfect interpolation and extrapolation of fCO\n2\nobservations by fCO -products in under-sampled northern regions.\nbeen revised. For the North, our suggested explanations are now clearly stated and we have\nadded information on the potential influence of coastal regions. For the South, a new\nseparate paragraph has been added, in which we clearly recognise the uncertainties\nassociated with the lack of fCO observations.\n2\n17) Fig 5, Line 225 Arctic Ocean not included in analysis\nThe figure caption was clarified and now reads: “Decadal trends in ΔfCO in the North\n2\n(>30°N, excluding the Arctic Ocean).”\n7\n18) Lines 337-340 Deficiencies of the GOBMs in the high latitudes are a major result of this\nanalysis and should be emphasized more. These deficiencies also demand increased\nsampling in the high latitudes to confirm the hybrid results presented here and correct\nGOBM deficiencies.\nIndeed, as mentioned in this paragraph of the discussion, the lack of variability in GOBMs in\nhigh latitude regions is a major issue and for reasons that remain partially unknown.\nHowever it is already well-known, with several publications have specifically addressed this\ndeficiency of the GOBMs (for example, DeVries et al., 2019, 2023; Hauck et al. 2020; Li et\nal., 2019; Friedlingstein et al. 2023). They are mentioned and discussed, but we added in the\ndiscussion that (lines 330-331) “Additional fCO sampling at high latitudes could help resolve\n2\nsome of these [high-latitude model] issues”.\nIn our revised discussion, we also now reaffirm the lack of observations in the high latitudes\nregions: “In the Southern Ocean, our hybrid approach suggests that existing fCO\n2\nmeasurements could corroborate a strong and positive decadal trend in this region in the\n2010s, and more generally between 2000-2022. However, the paucity of fCO\n2\nmeasurements impede our ability to evaluate the decadal trend using observations only\n(Gloege et al., 2021), but recent studies showed that undersampling could be responsible for\nstrong biases in fCO -products in that region (Hauck et al., 2023a).”\n2\n“Moreover, it [our analysis] confirms the importance of high-density fCO observations, which\n2\nare notably lacking in the Southern Ocean, for informing the fCO -products and our hybrid\n2\napproach.”\nHowever, we consider that the main objective and major results of our manuscript concern\nrecent trends in ocean CO sink. To better highlight the focus on our main results, we have\n2\nnow rearranged the order of the paragraphs in the discussion to start with our results on\ntrends in ocean CO sink (the first 5 paragraphs and the last 2 paragraphs). The results\n2\nrelating to the variability of the ocean CO sink are discussed afterwards (in 2 paragraphs).\n2\n19) Lines 299-300 “Our results suggest that fCO2-products overestimate the decadal trend\nof the 2010s in regions where there are few observations” This suggests that the hybrid\nmethod defaulted to the NEMO-PlankTOM model in undersampled high latitudes in both\nhemispheres; however, there should have been some adjustment made for the actual\nobservations even if they are sparse.\nActual observations are fully considered in the hybrid approach. Suggested errors in\nfCO -products had no influence on the hybrid approach, because we did not use them to\n2\nconstrain the ocean CO sink simulated by NEMO-PlankTOM12.1. The hybrid approach\n2\nconstrains the simulated ocean CO sink based on the SOCAT observations, even if\n2\nobservations are sparse (e.g., in the Southern Ocean). Moreover, this sentence was related\nto our analysis conducted in the North. We clarified our text in the sentence to avoid\nconfusion: “Our results suggest that some fCO -products overestimate the decadal trend of\n2\nthe 2010s in northern regions where there are few observations”\n20) Lines 315-317 The authors appear to assume that the observations in the Southern\nOcean are biased based on Refs 31 and 33, but the hybrid model locks onto the\n8\nmeasurement data (Fig 3d). If the observations are biased, then why trust the results of the\nhybrid method in the Southern region? A similar argument about biased observations in the\nNorth carries less weight given there are more observations in the 30-90N latitude range\nthan the 30-90 S range. This is a significant inconsistency that is not clearly resolved in the\ntext.\nThe reviewer is pointing out that the hybrid approach is following the results from the\nfCO -products in the Southern Ocean, but not in the North. As mentioned earlier, the\n2\nfCO -products results are not used by the hybrid approach to constrain the ocean CO2 sink\n2\nsimulated by NEMO-PlankTOM12.1, only direct observations are used. Therefore, we do not\nexpect the hybrid approach to produce estimates similar (or dissimilar) to those of\nfCO -products.\nconfirmed due to the lack of measurements in this region. In our revised discussion, we have\nadded a text arguing in favour of greater data collection, mainly in the Southern Ocean: “In\nthe Southern Ocean, our hybrid approach suggests that existing fCO2 measurements could\ncorroborate a strong and positive decadal trend in this region in the 2010s, and more\ngenerally between 2000-2022. However, the paucity of fCO2 measurements impede our\nability to evaluate the decadal trend using observations only (Gloege et al., 2021), but recent\nstudies showed that undersampling could be responsible for strong biases in fCO2-products\nin that region (Hauck et al., 2023a).”\nIn contrast, we have more confidence in the hybrid approach estimates for the North, as it\nfollows the temporal patterns visible in the SOCAT data, which are more abundant in this\nregion. In addition, we have suggested deficiencies for some fCO -products in their\n2\nestimates of the trend in the ocean CO sink in the 2010s in the North (as presented in\n2\nsection 3.4, and discussed in section 4)\n21) Lines 325-328 “GOBMs represent the correct processes, but either they do not respond\nsufficiently to changes in external forcing, or the balance among thermal and non-thermal\nprocesses in response to that forcing is imperfect.” This is a general statement, but should\naddress NEMO-PlanckTOM specifically given its central importance in the hybrid method\nWe added the mention of NEMO-PlankTOM12.1 in the sentence, which now reads as\nfollows “...suggest that NEMO-PlankTOM12.1 and other GOBMs represent the correct\nprocesses, but…”.\nNEMO-PlankTOM12.1 exhibits deficiencies in the simulation of the Southern Ocean CO\n2\nsink similar to those of the majority of GOBMs used in the analysis of the global carbon\nbudget, and we argue that the mention of other GOBMs in this sentence is justified. This was\nshown in particular by the analysis carried out by Hauck et al. (2023,\nhttps://doi.org/10.1029/2023GB007848). In this published analysis, NEMO-PlankTOM12.1\nhas been classified within the group of the DIC-weak GOBMs (as most GOBMs), where the\nstrong underestimation of non-thermal processes causes these models to be too strongly\ntemperature driven across the year. We have added this information in the manuscript to\nstrengthen and clarify our justification: “For example, in the Southern Ocean, ocean surface\n9\nfCO variations over the year in NEMO-PlankTOM12.1, and in most GOBMs, tend to be too\n2\nstrongly influenced by temperature changes (Hauck et al., 2023b).”\n23) Line 375 How sensitive is the estimated trend in the terrestrial CO2 sink to the value of\nthe trend assigned to the ocean CO2 sink?\nThe estimated trend in the terrestrial CO sink is directly influenced by the trend assigned to\n2\nthe ocean CO sink when estimating the land sink using the global carbon budget residual.\n0.20 Pg C yr-1 decade-1. Based on the other terms of the global carbon budget, it is possible\nto derive another estimate of the land CO sink and its associated trend using the residual of\n2\nthe other terms. This methodology for obtaining an estimate of the land CO sink (S ) with\nthe carbon budget residual has been added in the supplementary information document\n“Equations of the Global Carbon Budget analysis:\nWithin the Global Carbon Budget analysis, the carbon sinks (atmosphere = G , ocean =\nATM\nS , and land = S ) and emissions (fossil fuel = E , and land use change = E ) are\nOCEAN LAND FOS LUC\nestimated,\n(G + S + S ) = E + E (1)\nATM OCEAN LAND FOS LUC\nBy using the estimates of G , E , and E from the Global Carbon Budget analysis\nATM FOS LUC\npublished in 2023, with our estimate of S from the hybrid approach, an estimate of the\nOCEAN\nS term can be obtained:\nLAND\nS = E + E – (G + S ) (2)”\nLAND FOS LUC ATM OCEAN\nIn the context of the global carbon budget, uncertainties related to land-use change (E )\nLUC\nand the ocean sink for CO (S ) have the greatest influence on the S estimate.\nTherefore, our estimate target the uncertainty in the S related to S as mentioned in\nLAND OCEAN\nthe manuscript by comparing our estimate derived from our hybrid approach with that which\nwould be obtained with fCO -products: “Our estimated trend in the land CO sink is between\n2 2\nthe 0.34 ± 0.20 Pg C yr-1 decade-1 trend estimated by the Dynamic Global Vegetation Model\nensemble used in the global carbon budget analysis (Friedlingstein et al., 2023) and the\ntrend of -0.01 ± 0.13 Pg C yr-1 decade-1 that would be obtained with the ocean CO sink\n2\nestimate from the fCO2-products alone.”\n24) Line 385 The authors miss the opportunity to advocate for more and more extensive\nfCO2 observations, particularly at the high latitudes. Strongly urge that a final paragraph be\nadded to Section 4 that addresses research needs unresolved by the current study.\nWe have modified and strengthened our second to last paragraph of section 4: “Moreover, it\nconfirms the importance of high-density fCO observations, which are notably lacking in the\n2\nfCO2-products could be further improved by scrutinising the extrapolation of observations in\nthe 2010s, which are evolving over the different versions released, in order to understand\ndifferences among fCO -products, and then help improve them. ”\n2\n10\n25) Lines 449-450 “Model simulations are too short to fully represent the input of river fluxes\nand their subsequent outgassing of CO2 in the open ocean” This is a potentially serious\nmodel deficiency as riverine inputs will significantly alter coastal CO2 uptake and will have a\nmeasurable impact on the global ocean CO2 sink\nThis treatment of river fluxes is conformed to the definition of the ocean CO sink and the\n2\nprotocol used by the research teams providing ocean CO sink estimates from GOBMs for\n2\nthe global carbon budget analysis. In the global carbon budget analysis, the river fluxes are\nderived from specific published estimates. The river flux influences the mean air-sea CO\n2\nflux, but would not directly influence its variability, which is the focus of the current analysis.\nResponse to Reviewer #2:\n1) Why was the comparison only done with NEMO-PLANK model? There are other\n(higher-resolution) ocean circulation models for physics-biogeochemistry coupling, for\nexample, ECCO-DARWIN (Carroll et al, 2020, JAMES). This model is higher resolution than\nNEMO-PLANK and is data-assimilating. I suggest the authors compare how the global\nestimates from their models compare with ECCO-DARWIN and whether their hybrid\napproach performs better than ECCO-DARWIN.\nIn our study, we compare the ocean CO sink estimate of the hybrid approach with that of\nfCO2-products used in the global carbon budget analysis. In the manuscript, we highlight the\ndifference between the hybrid approach and NEMO-PlankTOM12.1, since the hybrid\napproach constrains the ocean CO sink simulated by NEMO-PlankTOM12.1. It is therefore\n2\nimportant to note the changes made to the NEMO-PlankTOM12.1 model output by the\nhybrid approach.\nBelow, we made a comparison between ECCO-Darwin and the fCO -products, the GOBMs,\nTable R1 | Decadal trends of the ocean CO sink between 1995 and 2018 estimated by different\n2\nmethods.\n2000-2009 2010-2018 2000-2018\n(PgC/yr/decade) (PgC/yr/decade) (PgC/yr/decade)\nECCO-Darwin 0.17 0.96 0.67\nfCO2-products 0.71 ± 0.38 0.67 ± 0.20 0.61 ± 0.18\nGOBMs 0.40 ± 0.10 0.44 ± 0.12 0.37 ± 0.05\nNEMO-PlankTOM12.1 0.27 0.60 0.46\nHybrid approach 0.80 ± 0.21 0.42 ± 019 0.51 ± 0.07\n2\napproach, with a very low value for the period 2000-2009 and a high value for the period\n2010-2018. Overall, over the period 2000-2018, ECCO-Darwin suggests a high decadal\ntrend, within the range of fCO -products, but with unique temporal patterns that is not\n2\ncoherent with results from model ensembles presented here.\nDeVries et al (2023) reported that ECCO-Darwin exhibits significant temporal variability in\nthe ocean CO sink and in the spatial variability of dissolved inorganic carbon (DIC)\n2\naccumulation rates. They mentioned that ECCO-Darwin results are subject to caveats\nconcerning the possibility of assimilation-induced model drifts affecting their results. For\nexample, DeVries et al (2023) mentioned that changes in biogeochemical model parameters\nduring assimilation can result in DIC gradients that are advected by the mean ocean\ncirculation and could affect air-sea CO fluxes and DIC accumulation. It is plausible that\n2\nthese potential model drifts explain the differences in decadal trends mentioned above\nbecause the data coverage evolves through time, but we are not able to make a firm\nassessment here because we cannot explore the robustness of the decadal trends to\nchanging data coverage and/or to different model choices within ECCO-Darwin.\nFurthermore, the ECCO-Darwin model results presented here (from the RECCAP2 project)\nonly cover the period 1995-2018, while our analysis covers three decades during the period\n1990-2022. For all these reasons, we prefer not to include the results of the ECCO-Darwin\nmodel in our analysis and to limit ourselves to comparisons with other types of estimates\nused in the global carbon budget analysis, for which several different estimates are provided\nwhich can be analysed as ensembles.\nHowever, we agree that more estimates from data-assimilation models are needed and can\nprovide additional information to constrain the ocean CO sink. For example, the Global\nprediction systems, which are based on assimilating physical atmospheric and oceanic data\nproducts into the ESMs. These models were used as a new line of evidence in predicting the\n2023 ocean and land CO sinks, and the atmospheric CO growth. For this reason, we have\n2 2\nnow added some references about ongoing research activity on the development of models\nusing data-assimilation: “Results from models using data-assimilation (Brasseur et al., 2009;\nVerdy and Mazloff, 2017), also suggest an underestimated decadal variability in GOBMs\nair-sea CO flux (Carroll et al., 2020; Li et al., 2023).”\n2\n2) As Nature Communications is for a broader audience, the authors should discuss the\nmeaning of fCO2 and the limitation of this metric in relation to ocean carbon storage in the\nintroduction. Specifically, fCO2 is an ocean-surface measurement and does not represent\nthe changes in ocean carbon content at depth, importantly below the permanent pycnocline.\nIn this context, the authors should discuss recent studies done for the ocean carbon storage\nin 3D space, e.g., Gruber et al (2019, Science), Zemskova et al (2022, Nat Comm).\nWith have added this information in the introduction:\n“This has been made possible by the annual release of quality-controlled\nobservations of CO2 fugacity (fCO2) at the sea surface — analogous to the partial\npressure of CO2 — compiled within the Surface Ocean CO2 Atlas (SOCAT)...”\n12\n“It should be noted that these fCO2-products only assess the air-sea CO2 flux, and\nnot where anthropogenic CO2 is ultimately stored in the ocean, which would require\nadditional measurements of carbon in the water column, as well as more\nassumptions (Gruber et al., 2019; Zemskova et al., 2022; Keppler et al., 2023; Müller\net al., 2023).”\nThe discussion related to the ocean carbon content at depth, which provides a constraint\nmostly on the mean ocean CO sink value than its yearly variations, was already present in\n2\nthe discussion section. We have extended it and included the references mentioned by the\nreviewer:\n“Finally for the mean ocean CO sink, the hybrid approach returns a higher mean\nconsistent bias of overestimation of the surface ocean fCO2. Studies based on\nemergent constraint properties have also suggested…\n[...]\nHowever, this means that our hybrid approach is less robust in estimating the mean\nocean CO2 sink than the variability and trend of this sink, because the mean ocean\nCO2 sink also depends on mixing between the surface ocean and the deep ocean\n(Orr et al., 2001), a process that is weakly constrained when only using surface fCO\n2\nobservations, as is the case in the hybrid approach. Further work, in particular the\nuse and/or assimilation of ocean interior carbon data (Carroll et al., 2020; Li et al.,\n2023; Zemskova et al., 2022), would be better suited to constrain the mean ocean\nCO sink. ”\n2\n3) The manuscript contains many acronyms and numbers, which deteriorates the readability\nof the manuscript. For example, are the acronyms A-IV and A-DV really necessary? The\nauthors could just say \"interannual\" and \"decadal\" variability in words. I also suggest that the\nauthors comb through the manuscript to limit the number of numbers that they present to\nthose that are crucial to the central points of the paper. These points need to be\nemphasized, as they are currently lost in the plethora of values presented in text. This\ntrimming would make the paper more suitable for a broader higher-impact journal like Nature\nCommunications rather than a methodology-oriented journal.\nWe have significantly reduced the use of acronyms and numbers to improve the readability\nof our manuscript. As suggested, we have removed the acronyms A-IV and A-DV and used\nthe terms 'interannual' and 'decadal' variability throughout the manuscript.\nResponse to Reviewer #3:\nMy major comments are as follows (in order of my perception of priority):\nA) The authors chose to exclude the coastal ocean from their analysis--everything within 3\ndegrees of the coast and shallower than 1000 m. Given the outsize role of these regions for\nmarine carbon fluxes, I believe the quantitative consequences of that choice should be\nclarified and visualized for readers.\nfrom open-ocean and coastal zones. When doing so, we are using 94% of the SOCAT\nobservations. The remaining 6% of SOCAT observations not being used are located in areas\nof shallow water not simulated by NEMO-PlankTOM12.1. Below, we have tested the\nrobustness of our results to the inclusion, or not, of coastal observations, by building two\ndifferent open-ocean masks.\nA first mask is made out by using the monthly fCO fields from the fCO -products and\n2 2\nselecting 1°x1° grid cells covered by all fCO -products. In fact, some fCO -products only\n2 2\ncover 90% of the ocean surface, because they do not resolve the ocean CO sink in\n2\ncontinental shelves and high latitude regions. We therefore compare our results using\nSOCAT observations in the open-ocean areas common to all fCO -products. Using this\n2\nmask, we use on average 88% of the SOCAT observations.\nA second mask is made out by using the definition of the coastal zone from Laruelle et al.\n(2018). This mask was used in the original submitted version of our manuscript. By using\nthis mask we are using on average only 69% of SOCAT observations.\nResults show (Table R2 and Figure R1) that whether or not coastal observations are\nincluded, this has a limited impact on the conclusions of our analysis concerning the decadal\ntrend for the period 2000-2022. Global and regional variations are similar with and without\ncoastal observations, with the decadal trend of the global ocean CO sink for the period\n2\n2000-2022 constrained with the hybrid approach remaining unchanged (0.42 Pg C / yr /\ndecade). The only exception is that the global ocean CO sink is lower between 2005 and\n2\n2011 in the results of the hybrid approach that used only open ocean observations as\ndefined in Laruelle et al. (2018). However, using this definition of coastal areas limits the\nnumber of SOCAT observations used to 69%, and removes observations that are\nconsidered as open ocean observations by the fCO -products. Therefore, to provide a better\n2\ncomparison with fCO -products, and simplify the explanation of our hybrid approach to\n2\nreaders, we decided to use all available observations.\nHowever, we note that the ocean CO sink in coastal regions, and its evolution, remains\n2\nuncertain and add this information in the discussion:\n“In addition, the ocean CO sink in the northern region is also more influenced by\n2\ncoastal processes than the southern region, which despite their importance remain\nuncertain (Dai et al., 2022; Regnier et al., 2022; Laruelle et al., 2018; Mathis et al.,\n2024; Resplandy, 2024). Consequently, in the northern region, the way in which\ncoastal fCO2 observations are taken into account by the various fCO2-products\ncould induce some of the discrepancies among fCO2-products. This would partly\nexplain the lack of coherence between the GOBMs and the fCO2-products over this\nrecent decade (Friedlingstein et al., 2023; Gruber et al., 2023; Keppler et al., 2023;\nMüller et al., 2023).”\n2\nnot coastal observations.\nInterannual Decadal 1990-1999 2000-2009 2010-2019 2000-2022 PercentageofSOCAT\nvariability variability (PgC/yr/decade) (PgC/yr/decade) (PgC/yr/decade) (PgC/yr/decade) observationsbeingused\n(PgC/yr) (PgC/yr) (annualaverage)\nHybridapproachwithall 0.22±0.02 0.13±0.02 -0.19±0.17 0.80±0.21 0.44±0.15 0.42±0.06 94%\navailableobservations\nHybridapproachwithall\nobservationsinthearea\n0.23 0.13 -0.17 0.80 0.46 0.42 88%\ncoveredbyall\nfCO-products\ncoastalzones,asdefined 0.24 0.12 -0.14 0.47 0.45 0.42 69%\nbyLaruelleetal(2018)\nFigure R1 | The global and regional hybrid approach results when including or not coastal\nobservations. In red, the hybrid approach results with all available observations. In blue, when using the\nopen-ocean area covered by all fCO -products. In green, when using the open-ocean area as defined by\nB. In my opinion, this paper would be more impactful and robust if the authors included the\nresults of sensitivity analyses used to inform their approach. Specifically:\ni. Several possible mechanisms are proposed to drive discrepancies between models and\ndata-products. Varying subsampling and model parameters would help identify the most\nlikely of these mechanisms.\nOur main focus and results are about the discrepancy between the GOBMs and\nfCO -products for the decadal trend of the global ocean CO sink between 2000 and 2022.\n2 2\nBased on our analysis, we suggest an overestimation of the decadal trend in the 2010s by\nthe fCO -products. The proposed mechanism to explain this overestimation is a lack of fCO\n2 2\n15\nmeasurements, and some methodological issues in the interpolation and extrapolation of the\navailable data. In order to support this suggestion, we report the evolution of the\nfCO -product estimates within the last three global carbon budgets (see text below and the\n2\nnew figure 6). This analysis is more relevant than performing a subsampling experiment with\nour hybrid approach because it directly demonstrates that the trend of the 2010s from\ndata-products has been sensitive to revisions in data and methods, while that from GOBMs\nhas varied very little: “In addition to the extrapolation problems in the under-sampled\nnorthern regions mentioned above, on a global scale, the estimate of the positive trend in the\n2010s from the fCO2-product ensemble has been revised downwards in successive\npublications of the global carbon budget analysis between 2021 and 2023, while their trends\nfor the 1990s and 2000s have remained relatively similar (Fig. 6). For each of the global\ncarbon budget analyses published between 2021 and 2023, the fCO2-product ensemble\naverage has always been produced from seven estimates. However, two of the seven\nfCO2-products were introduced, replacing previously submitted products that were not\nupdated, and five were slightly updated. Among these five fCO2-products, on average, the\n2010 trend between the 2021 and 2023 publications decreased by -0.05 Pg C yr-1\ndecade-1. Thus it was mainly the turnover in the last two fCO2-products between 2021 and\n2023 that led to a visible decrease in the ensemble average of -0.25 Pg C yr-1 decade-1 for\nthe 2010 trend. Consequently, the downward revision observed for the data-based estimate\nwas mainly due to a change in two fCO2-product methodologies and, to a lesser extent, to\nthe annual updates of the SOCAT database and fCO2-product methods, suggesting that the\ntrend of the 2010s estimated with the data-products is not robust at this stage.”\n“Our results suggest that the estimate of the ocean CO2 sink trend in the 2010s by the\nfCO2-product ensemble is overestimated and sensitive to the availability and distribution of\nfCO2 observations. Over the last three annual updates of global carbon budgets\n(Friedlingstein et al., 2022a, b, 2023), although the individual members of the fCO2-product\nensemble underwent very few changes, the 2010s trend estimated from the ensemble of\navailable products has decreased by 14% each year. This supports our finding of an\noverestimated trend in the 2010s ocean CO2 sink from the fCO2-products ensemble, which\nis adjusted downwards as new data become available. In addition, the replacement of two\nmembers of the fCO2-product ensemble by a hybrid approach along the same lines as\npresented here (Gloege et al., 2022; Bennington et al., 2022b) and by a revised\nfCO2-product aimed at improving the retrieval of the ocean CO2 sink trend (Zeng et al.,\n2022) have led to this downward revision of the 2010s trend in the latest ensemble.\nConsequently, our analysis demonstrates the importance of regular updates and efforts to\ncollect fCO2 observations as part of SOCAT (Bakker et al., 2016), as well as evaluations of\nfCO2-product methodologies (Hauck et al., 2023a; Gloege et al., 2021).”\n2\ncarbon budget analysis. The last update of the fCO -product estimate was in 2023 (black\n2\nline and grey shade). The previous two estimates were in 2022 (orange) and 2021 (blue).\nThe left y-axis represents the global CO flux anomaly (for each time series the long-term\n2\nmean between 1990 and 2020 was removed). GOBM estimates for the same three global\ncarbon budget updates are also shown in dotted lines. The number of annual observations in\nthe SOCAT database, for each annual version of SOCAT (v2021, v2022 and v2023), is\nshown in the bottom right-hand corner using the same colour code as for the line. Years prior\nto 2005 are not shown.\nWe also report a discrepancy between GOBMs and fCO -products for the decadal trend of\n2\nthe ocean CO sink in the northern region. To highlight this discrepancy we subsampled the\n2\nΔfCO data from fCO -products and GOBMs, and compared them to the SOCAT database.\n2 2\nThis analysis suggested potential methodological issues within four fCO -products. We also\n2\nsuggest that the way in which coastal fCO observations are taken into account by the\n2\ndifferent methodologies could be important. However, no obvious influence of the coastal\nobservations in the North region was detected with our coastal observation sub-sampling\nexperiment carried out for our response to the previous comment (A).\nFinally, we report a discrepancy between GOBMs and fCO -products in decadal variability of\n2\nthe ocean CO sink. The proposed mechanisms are those cited in the existing literature.\nperformance of the methods used by fCO -products to reproduce the decadal variability of\n2\nthe fCO fields of GOBMs (Gloege et al., 2021; Denvil-Sommer et al., 2021; Hauck et al.,\n2\n2023a). Therefore, performing such an analysis with our hybrid approach would not provide\nadditional information on this discrepancy in the decadal variability of the ocean CO sink\n2\nbetween GOBM and fCO -products, which has already been addressed in several published\n2\nanalyses.\nii. Presenting the results of the sensitivity analysis used to select and tune key parameters\nfor optimization and perturbation--within the constraints identified in the literature--would\nincrease confidence in the results.\ncomment C below. We explain our choice to remove the paragraphs devoted to model\noptimisation (which was made in a broader context of model development rather than\nfocused on this specific study) and to develop model validation instead which is directly\nrelevant to this paper.\nFor the selection of the parameters to perturb NEMO-PlankTOM12.1: we have tested\nparameters related to the bacteria (as in the original submission and used in the main\nmanuscript) and the phytoplankton functional types (see Table R3 and Figure R2 below).\nBacteria parameters were selected for our main results because bacteria are ubiquitous in\nthe ocean, and therefore modifications of their parameters affect carbon fluxes in all ocean\nregions in a consistent and uniform manner. Phytoplankton are already more patchy than\nbacteria. Tests done on phytoplankton parameters provide an idea of robustness but would\nnot be the preferred parameters to use because the response is not as uniform as when\nusing the bacteria parameters due to the inhomogeneous distribution of phytoplankton\nbiomass. We have also tested two different model configurations using two different\natmospheric forcings: NCEP and ERA5. Regardless of the parameter and forcing used, the\nkey results of our analysis remain unchanged about the variations of the decadal trends in\nocean CO sink over the period 1990-2022, with the lowest trend observed in 1990s, and a\n2\nhigher trend in the 2000s than in the 2010s. We now detail the results of the sensitivity\nanalysis in our manuscript and supplementary material (supplementary table 3, and\nsupplementary figures 5 and 6):\n“However, both GOBMs and fCO -products ensembles suggest similar trends between\n2\n2000s and 2010s, while the hybrid approach (including the sensitivity analyses) always led\nto a higher trend in the 2000s compared to the 2010s.”\n“Although our hybrid approach always suggests an underestimation of the decadal variability\nby GOBMs, the exact value is sensitive to the specific model configuration.”\n“Finally, we performed sensitivity analyses to test the robustness of our results to the choice\nof perturbed model parameters and model configurations. The perturbed simulations were\nrepeated with parameters of phytoplankton respiration, and with a combination of both\nbacterial half-saturation and phytoplankton respiration. The model configuration was\nchanged by using ERA5 reanalysis as weather forcing data. In total, we thus have applied\nthe hybrid approach to six different set ups, with three choices of perturbation parameters\nand two choices of forcing configurations (Supplementary Figure 5). Regardless of the\nparameter and configuration used, the results of our analysis remain unchanged about the\nvariations of the decadal trends in ocean CO2 sink over the period 1990-2022, with the\nlowest trend observed in the 1990s, and a higher trend in the 2000s than in the 2010s\n(Supplementary Figure 6 and Supplementary Table 3). We show here results of the model\nforced with NCEP, which has a lower RMSE (38.5 μatm) compared to the configuration\nforced with ERA5 (40.0 μatm), and the perturbation of the half-saturation constant of\nbacterial remineralisation which produces perturbations that are more uniform across the\nocean.”\nconfiguration and by perturbing the bacteria. Because when the model is forced with NCEP,\nit has a lower RMSE (38.5 μatm) than when forced with ERA5 (40.0 μatm). As mentioned\nabove, the biomass of bacteria varies less in GOBMs (and in reality) than the biomass of\nphytoplankton. Finally, the range of perturbations of the bacterial parameter that we used\n(5-18.10-6 mol/L) is in agreement with the literature, e.g., Mentges et al. (2019, Scientific\nReport). Note that Mentges et al., (2019) also mentioned this parameter as “being the least\nconstrained by published referenced values”. Therefore, our analysis carried out with\nperturbed bacteria could be more easily reproduced by other modelling teams. It should be\nnoted that the lower 2000-2022 trends (0.27-0.24 PgC/yr/decade) associated with the\nperturbed phytoplankton members could be explained by a large number of yearly values\nconstrained but outside the bounds of the perturbation experiments (uncertain values)\nobtained with these members.\nTable R3 | Temporal variations of ocean CO2 sink from hybrid approach results with different\nperturbed parameters.\nInterannual Decadal 1990-1999 2000-2009 2010-2019 2000-2022\nvariability variability (PgC/yr/dec (PgC/yr/dec (PgC/yr/dec (PgC/yr/dec\n(PgC/yr) (PgC/yr) ade) ade) ade) ade)\nHybrid approach\n0.22 ± 0.02 0.13 ± 0.02 -0.19 ± 0.17 0.80 ± 0.21 0.44 ± 0.15 0.42 ± 0.06\nNCEP\nbact\nHybrid approach\n0.28 0.16 0.17 0.97 0.48 0.27\nNCEP\nphyto\nHybrid approach\nNCEP 0.18 0.09 -0.37 0.73 0.37 0.41\nphyto + bact\nHybrid approach\n0.23 0.25 -0.93 1.27 0.35 0.36\nERA\nbact\nHybrid approach\nERA 0.19 0.18 -0.60 0.87 0.23 0.24\nphyto\nHybrid approach\nERA 0.17 0.14 -0.43 1.14 0.23 0.36\nphyto + bact\nthe grey shading, represents the results and error estimate from the main manuscript (using\nNEMO-PlankTOM12 forced with NCEP, and perturbed bacteria). The other coloured lines represent hybrid\napproach results obtained with a different forcing (i.e., ERA), and/or when perturbing phytoplankton and/or\nbacteria. The white dots represent uncertain results as defined in the manuscript.\nC. I think the model manipulations described here constitute a sufficiently different approach\nand parameter set that validation should be presented. e.g., an evaluation of how major\nvariables like temperature and salinity, as well as fCO2, differ between this and prior\nPlankTOM configurations and climatological hydrography.\nAs mentioned earlier, in this revised version of the manuscript we are using the\nNEMO-PlankTOM12.1 model that has been submitted to the Global Carbon Budget 2023.\nAs part of the submission to the Global Carbon Budget 2023, the NEMO-PlankTOM12.1\nmodel and the other GOBMs undergo a standard protocol for model validation, which has\nbeen greatly expanded in the Global Carbon Budget 2023 and includes, in addition to a\ncomparison with fCO observations from SOCAT, comparisons with observed :\n2\n● mean ocean CO sink in the 1990s\n2\n● Atlantic Meridional Overturning Circulation\n● Southern Ocean sea surface salinity\n● Southern Ocean stratification index\n● surface ocean Revelle factor\nThe physical model setup from the NEMO-PlankTOM12.1 model used here is unchanged\nsince. It was validated within the 2023 publication of the Global Carbon Budget. The model\nversion provided here differs from the 2023 submission to the Global Carbon Budget 2023\nthrough its forcing only (NCEP instead of ERA5). We have reconducted the model validation\nas performed by the Global Carbon Budget analysis (Figure R3).\nThis version of NEMO-PlankTOM12.1, forced with NCEP, has a better representation of\nfCO than the version forced with ERA5. We also show that this model is within the range of\nnow fully explained and described within the section 5.2: “Description and validation of the\nNEMO-PlankTOM12.1 model”\nD. Various results are compared, but statistical tests are not presented. Including these\nwould help support interpretations of whether or not values are similar or meaningfully\ndifferent.\nThe text makes multiple comparisons of the hybrid approach results with the GOBMs and\nfCO -products results, as for example with the Table 1. Such comparisons take into account\n2\nthe uncertainty around the hybrid approach results, as well as the ±1σ (standard deviation)\nrange associated with the results for GOBMs and fCO -products. As a conservative\n2\napproach, if the hybrid approach value, with its uncertainty, is within the range of the GOBMs\nor fCO -products, we consider that the hybrid approach is not significantly different from the\n2\nvalues provided by the GOBMs or fCO -products. Moreover, we have added results from\n2\nsensitivity analysis and mentioned them whether they strengthen or not our results.\nOtherwise, when statistical comparisons or tests were performed (for example with Figure\n4), we reported the associated p-value in the text.\nE. I think the nudging/data assimilation approach could be better referenced to illustrate how\nit draws on or differs from prior literature.\nThis comment refers to this specific comment:\n“242. My perception is that statistical approaches to 'nudging' earth, ocean, and\nclimate model results to incorporate observational constraints has a long history,\nwhich should probably be referenced here. Examples and context include:\nhttp://dx.doi.org/10.5065/D60Z716B,\nhttps://doi.org/10.1175/1520-0493(2000)128>3664:ASNTFD<2.0.CO;2,\nhttps://www.myroms.org/forum/viewtopic.php?t=320,\nhttps://doi.org/10.5194/gmd-16-1857-2023.\nThe author's will note that these sources point out that nudging to observations is\nsometimes considered a suboptimal data assimilation strategy. I recognize that the\nmodel optimization used here is more complicated, and may not be well summarized\nas \"nudging\".\nOther works have used more data assimilation techniques to incorporate marine\ncarbonate system variables in concert with physical parameters:\nhttps://doi.org/10.1016/j.jmarsys.2016.02.011,\nhttps://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JC012650,\nhttps://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2009GB003531,\nhttps://www.jstor.org/stable/24860162, https://www.jstor.org/stable/24861004.\nGiven this context, I think the authors need to not just reference related concepts and\nworks of this nature, but more clearly delineate how this work is novel in approach\nand focus. I acknowledge that the length-limitations and format of this journal may\nhave obscured significant advances in the approach used by the authors--if so and I\n21\nmissed these, I take that as further evidence that these need to be clearly and\nconcisely delineated in the text.”\nWe realised that the principle of our hybrid approach is not well summarised by the word\n“nudging”. We therefore removed this expression from the text and rephrased our\ndescription (see below).\nWe agree that more estimates from model assimilation techniques have been published and\nprovide more information about the ocean CO sink. For example, the Global Carbon Budget\n2\nanalysis in 2023 used four Earth System Models (ESMs) prediction systems, which are\nbased on assimilating physical atmospheric and oceanic data products into the ESMs. They\nwere used as a new line of evidence in predicting, for 2023, the ocean and land CO sink,\n2\nand the atmospheric CO growth. In addition, such assimilation methodology can provide\n2\ninformation not only on the air-sea CO flux, but also on the sequestration and long-term\n2\nstorage of the CO . Although, there are also some limitations related to the way the\n2\nassimilation of data is performed, as mentioned in DeVries et al (2023) for the ECCO-Darwin\nmodel. Nonetheless, we have now added some references about this ongoing model\nassimilation effort. We also point out to another recent study that also used an hybrid\napproach between GOBMs and fCO observations, and mentioned the similarities in our\n2\nresults:\n“Results from models using data-assimilation (Brasseur et al., 2009; Verdy and Mazloff,\n2017), also suggest an underestimated decadal variability in GOBMs air-sea CO2 flux\n(Carroll et al., 2020; Li et al., 2023).”\n“A similar hybrid approach was recently published (Gloege et al., 2022), but with a machine\nlearning algorithm used to derive the factors influencing the fCO2 variability. Here, the\nmechanism as represented in the NEMO-PlankTOM12.1 model remained unchanged and\nthus also constrained the results.”\n“In addition, the replacement of two members of the fCO2-product ensemble by a hybrid\napproach along the same lines as presented here (Gloege et al., 2022; Bennington et al.,\n2022b) and by a revised fCO2-product aimed at improving the retrieval of the ocean CO2\nsink trend (Zeng et al., 2022) have led to this downward revision of the 2010s trend in the\nlatest ensemble.”\nF. The authors appear to use different sources and calculations for water and air fCO2,\ninstead of the concurrently and self-consistently calculated values from SOCAT. I\nrecommend the authors evaluate whether it is possible to consolidate their data sources in\nthis case.\nThis comment refers to this specific comment:\n“348. SOCAT reports both the wet and dry atmospheric fCO2 measured concurrently\nwith the seawater fCO2 selected by the authors, eliminating the need for this\ncalculation and the potential errors in introduces because of discrepancies in choices\nof equation and virial coefficients used by different groups, differences in the spatial\ngridding, etc. I strongly recommend the authors revisit their data source and this\ncalculation considering the alternative SOCAT availability--if this is not possible,\n22\nplease explain why (e.g., are there no comparable interpolated products for the\nSOCAT air fCO2?).”\nThe SOCAT gridded product does not have xCO2, nor fCO2air values. When measurements\nare submitted to SOCAT, it accepts fCO2air measurements made in parallel to fCO2water\nmeasurements, and archives these fCO2air measurements for future use. However, SOCAT\ndoes not publish these measurements.\nHowever, SOCAT publishes Sea Surface Temperature (SST) and Sea Surface Salinity\n(SSS) measurements. So, we have done the subsampled side of figure 4 with the SST and\nSSS from SOCAT, and the obtained results are identical (see Figure R3). However, for\nconsistency with the non-subsampled side of figure 4, we have decided to keep our\nmethodology.\nFig. R3 | Decadal trends in ΔfCO in the North (>30°N, excluding the Arctic Ocean).\n2\nfCO -product ensemble (blue boxplots) and of the GOBM ensemble (grey boxplots). The\n2\nfCO -products and GOBMs were sub-sampled at SOCAT locations. The trends are\n2\ncalculated from median annual values for the nord region. On the left, the SST values are\nfrom OISST1.2 and the SSS values are from EN4 (as in the manuscript), while on the right\nSSS and SST values are from SOCAT.\nDetailed comments by line number (because this is a long list, I don't expect responses to all\nof these):\n1) Line 53. It would be helpful to further distinguish what is meant by fCO2-products here in\ncontrast to observational data only. Are the authors only referring to gridded, interpolated,\nand climatological datasets that use different methods but all fundamentally rely on SOCAT?\nAre there additional substantive differences readers should be aware of?\nmentioned in the text and in Table 1), and all rely on SOCAT. We have rewritten the\nsentence:\n“Estimates from these observation-based products (fCO -products), that all used the\n2\nsink…”\n2) 61. What are the specific fCO2 products being considered as part of the ensemble? And\nthe specific GOBMs included in the ensemble? I think each product considered should be\nreferenced and included in a table associated with this manuscript. If both ensembles are\nidentical to those used in the Global Carbon Budget 2022 publication, the I think that this\npoint should further be clearly stated in the inline text with an associated reference number.\nThe ensembles are exactly the same as the ones used in the Global Carbon Budget analysis\n2023. We have added this information in the text,\n“Despite recent progress, the two types of estimates within the global carbon budget\nanalysis diverge over the last two decades (2000-2022), with the fCO -product\n2\nensemble suggesting a decadal rate of growth of the ocean CO sink almost twice as\n2\nhigh as that simulated by the GOBMs ensemble (Table 1).”\n3) 72. a. The result in the 2010s trend column--in which the result is significantly lower than\nthe model ensemble, preferred model, or fCO2-product ensemble and thus seems to\ncontradict either observational or model evidence--strikes me as a red flag that requires\nclearer explanation in the text. I'm not sure section 3.4 is sufficient to convince me. Any\nproposed resolution would be more convincing with a model sensitivity study demonstrating\na plausible mechanistic basis.\nb. The larger interannual variability in the presented approach than the other models and\ndata products is also concerning, though less so given the reported uncertainties.\nc. I suggest a concise table title, rather than expanding this into a caption.\nIn the revised manuscript, values within Table 1 have been modified. The value of the ocean\nCO sink trend, between 2000 and 2022, obtained from the hybrid approach now lies\n2\nbetween the GOBM and fCO -product estimates, while still slightly closer to the GOBMs. We\n2\nexplain these changes (between the submitted and revised manuscript) by modifications in\nthe fCO -products results over the last three updates of the global carbon budget analysis\n2\n(see section 3.3 and figure 6), and also by a modification of the NEMO-PlankTOM12.1\nmodel. The description of this result is done in sections 3.3 and 3.4, and discussed in section\n4.\nWe have performed sensitivity analysis about the interannual variability in section 3.1 to\nexplain the observed discrepancy.\nThe table title was changed.\n4) 87. Given that this approach is so central to the analysis, I think a brief overview here (a\ncouple more sentences), prior to the trailing methods section, would be appropriate.\nthe hybrid approach (see response to comment #6)\n5) 90. The trends in the preferred model are 20-60% lower than the ensemble in each\ndecade. Is it fair to call these similar trends? Or perhaps it makes sense to restate this as a\nsimilar multidecadal trend, since the decadal trends diverge. I also suggest double checking\nthe overall 2000-2019 NEMO trend--based on Fig. S1 it appears that either the overall trend\nor detrended decadal component trend should be a bit smaller than that of the ensemble\n(but I may be mistaken).\nWe have modified this sentence as requested:\n“The global annual values of the ocean CO sink estimate of the current model\n2\nversion approximately matches the model ensemble average used within the latest\nglobal carbon budget analysis, and with similar variability and multidecadal trend\n(Table 1, Supplementary Fig. 2)”\nWe have also double checked the values.\n6) 96. A caveat to this comparison to is that these prior studies performed sensitivity\nanalyses across a wide range of parameters and included observational constraints (of\nparameters other than the objective variable) on the the tuned parameters. If this study does\neither, it was not obvious to me. So I think it is appropriate to include some comment on how\nthose prior studies guided the authors approach, and any limitations.\nWe have now added results from sensitivity analysis in the manuscript (see our response to\ncomment B). We have also rewrite the sentence to be more specific about which part of the\nprevious analysis was used in the current study,\n“Note that the hybrid approach used here, which estimates an optimised value of a\ntarget variable based on several model simulations and a cost function, has been\nused in previous studies to constrain global ocean primary production and air-sea\nfluxes of N O and CCl .”\n2 4\n7) 103. Similar to my comment about line 85, I think the authors should introduce a concise\nsummary of the perturbation and tuning approaches here, beyond stating that they occurred.\nThis context seems important for helping readers understand why state estimates\nsometimes arise outside the bounds established by the perturbations.\nDone:\n“This was obtained by perturbing model parameters. Perturbed simulations provided\na range of possible values for the ocean CO sink around the estimate from the\n2\nstandard model simulation. Then, for each year, the optimal CO sink was found\n2\nwithin this range of possibilities by optimising the calculated Mean Square Error\n(MSE) between the simulated fCO2 and the SOCAT observations. The hybrid\napproach also provides a quantitative estimate of uncertainty (see methods).”\n8) 113. This is a matter of personal preference, but now that the authors have defined this\nvariability in terms of amplitude, could the acronyms be shortened to IV and DV? Or, see\ncomment for line 133.\n9) 113. \"Very close\" = indistinguishable? Considering these as a one sided test of difference\n(with table of tests results in supplement?) might clarify the meaning of statements about\nsimilarity ('very close', 'close to', 'marginally overlapping', 'in line with', 'within the range of',\netc.). For consistency of interpretation by readers, the authors might consider choose a\nsingle phrase to indicate not significantly different comparisons.\nAs mentioned and explained in our response to comment D, we have decided to replace this\nstatement by “within the range”.\n10) 120. I think it would be useful to label the perturbed simulation lines with a descriptive\nlabel of how they were perturbed (e.g., microbial K1/2 = ...). Perhaps this is semantic, but in\nthe final sentence of the caption I think \"constrained but outside the bounds of the\nperturbation experiments\" might be more descriptive than 'uncertain constrained'.\nWe have added the range of the perturbed parameter in the figure caption, and modified the\nfinal sentence as requested:\n“Empty red dots are years with an ocean CO sink value constrained but outside the\n2\nbounds of the perturbation experiments (i.e., uncertain values, see methods). The\nperturbed simulations are produced by varying the half-saturation constant of\nbacterial remineralisation (from 5.10-6 mol L-1 to 18.10-6 mol L-1).”\n11) 131. Wouldn't a moving-average as described here often result in lower variability even\nwithout changes in the data distribution? I suggest that smoothing the data in this way is no\nlonger a comparison of interannual variability in the same sense as it was previously\nintroduced. To be very clear, the authors might decide to replace A-DV and A-IV, with a\nnotation that clearly displays the timescale considered. e.g., V_{10y}, V_{3y}, V_{1y}.\nYes, the moving-average lowers the variability. Similarly, applying the hybrid approach over a\nthree year period also lowers the variability. Therefore, comparing the variability from the\nhybrid approach applied to 3 years, with the model and observation-based variability after a\n3-year moving average smooth is reasonable. We decide to not employ the suggested\nnotation to not surcharge the manuscript with acronyms.\n12) 154. a. A mismatch may help diagnose mechanisms, but it is not itself a cause of\nanything. Nor does the presence of a mismatch--absent a clearly aligned expectation\ngrounded in a particular\nmechanism--conclusively indicate that the southern mid-latitude\nCO2 sink is underestimated in models as opposed to overestimated in sparse observations.\nI think this statement needs to be revised--perhaps the authors meant instead that the\noverall regional and global disagreements between data sources reflect these latitudes in\nparticular rather than elsewhere.\nWe have rewritten this sentence:\n26\n“...shows that the mismatch between the model and the fCO observations in the\n2\n40-60°S band could be the main cause of the underestimation of the ocean CO sink\n2\nin the South (Fig. 2). Note that it is the Southern region that has the most influence\non the global ocean CO sink.”\nthe southern high latitudes.\n13) 159. Are the uncertainty bounds smaller than the symbols on most of these bands?\nDoes it make sense that they are this small given the large interannual and interdecadal\nvariability previously discussed?\nYes, most of them are smaller than the symbol. It makes sense as this analysis is based on\nclimatology and not on a specific year.\n14) 194. See comment for line 72a.\nSee our response to comment #3.\n15) 219. If relying on significance for argument, I suggest including test result(s) (null, p,\ntest-statistic).\nThis section has been modified, and see our response to comment D.\n16) 208. I think that this statement implicitly assumes a constant gas transfer velocity. I'm not\nfamiliar enough with global windspeed or momentum flux trends to know whether that is\nlikely but there are observed and modeled changes in these parameters over the last few\ndecades in at least some regions within this latitude band (e.g., Eastern boundary current\nregions). The assumptions here probably need some reinforcing citations.\nThe surface ocean ΔfCO is the driving force behind the air-sea CO flux. The CO flux is\n2 2 2\ngenerally estimated using a bulk formula in which ΔfCO is multiplied by a gas transfer\n2\nvelocity, the solubility of CO in seawater and the ice cover. The gas transfer velocity is\n2\nfrequently parameterized as a function of wind speed.\nIf the air-sea flux of CO were driven primarily by variations in gas transfer velocity, all fCO\n2 2\nproducts would have similar decadal trends for the air-sea flux of CO . Because all\n2\nfCO -products use similar atmospheric forcing datasets for wind speed. For example, see\nstandardised approach for calculating fluxes from fCO products, and still observed large\n2\ndifferences in estimated CO fluxes between different fCO products that are induced by the\n2 2\nvarious fCO mapping techniques.\n2 2\nWe modified the sentence:\n“The strong decadal trend in the ocean CO flux estimated by the fCO -products is\n2 2\nprimarily driven by diverging trends between the CO2 fugacity at the surface of the\nocean compared to that in the atmosphere (ΔfCO ) (Fay et al., 2021).”\n2\n27\n17) 219. a. Please include trends and relevant statistical tests to reinforce this conclusion.\nb. To me, this distinction is a bit arbitrary (or at least reflect a choice of emphasis). Visually,\nany given decade seems unlikely to diverge greatly from the overall 1990-2020 trend. i.e., if\ncomparing each decade to null hypothesis of no difference from the secular trend, are any of\nthese decadal variations actually significant?\nSee our response to comment D.\n18) 238. a. I think this point about spatial bias in the interpolated products is valuable and it\nwould be good to start a new paragraph here to visually highlight this idea.\nb. This point (bifurcation of the interpolation products into high and low clusters) would be\nstrengthened by a supplemental figure.\nWe split this paragraph in two. We did not add a supplementary figure as we believe the text\nis clear enough.\n20) 258. Please see comment regarding line 156a. Without testing one of the mechanisms\nproposed, I think that this statement is not fully supported as written.\nThis result has changed in the revised manuscript\n21) 378. This is unclear to me as written. I suggest including an equation and calculated\nfluxes (and their data sources) in the supplement.\nDone:\n“EquationsoftheGlobalCarbon Budgetanalysis:\nWithin theGlobal Carbon Budgetanalysis,the carbonsinks(atmosphere =G ,\nATM\nocean=S ,andland= S )andemissions(fossil fuel =E ,andlanduse\nOCEAN LAND FOS\nchange=E ) are estimated,\nLUC\n(G +S + S )= E + E (1)\nATM OCEAN LAND FOS LUC\nByusingtheestimates ofG ,E ,andE fromthe GlobalCarbon Budget\nATM FOS LUC\nanalysispublishedin 2023,with ourestimateofS fromthe hybridapproach, an\nOCEAN\nestimateoftheS termcanbeobtained:\nLAND\nS =E + E –(G +S ) (2)\nLAND FOS LUC ATM OCEAN\nNotethat thedifferencebetween theterms inequation 1isequal tothe Budget\nImbalance(B ) ofthe globalcarbon budget:\nIM\nBIM= E + E –(G + S + S ) (3)”\nFOS LUC ATM OCEAN LAND\n28\n23) 417. Randomly selected with what distribution? I ask because the authors previously\nchose to use tests that are suitable for non-normal data. If the data is not well approximated\nby a normal distribution, that should be noted here and an appropriate distribution chosen for\nthis Monte Carlo approach.\nWe have use a uniform distribution, and specify it in the text:\n“...that had been randomly selected (from a uniform distribution) within…”\n24) 444. Some additional details about the treatment of major inorganic carbon components\nwould be appropriate here, even if covered in the reference.\nWe have added:\n“Simulated dissolved inorganic carbon and alkalinity are influenced by air-sea\nexchanges of CO , calcification (production and dissolution), primary production, and\n2\nremineralisation of organic matter (grazing by zooplankton and remineralisation by\nbacteria). The alkalinity is also influenced by denitrification.”\n25) 398. a. If re-optimized as described, this is a new model variant that really requires\nindependent validation from prior published works. Did this process knock any parameters\nback out of reasonable ranges? Are basic properties like T & S well described in at least a\nclimatological sense?\nI acknowledge model validations are a lot of work, but as a reviewer it is difficult for me to have\nconfidence in the results presented here without validation of parameters beyond fCO2 of this new\nproduct presented in the supplement.\nb. This optimization scheme is empirical but it isn't clear to me that the optimized parameters\nare realistically constrained. For example, if I am interpreting the temperature sensitivities\ncorrectly (e.g., based on parameter description by Le Quere 2016), that's equivalent to a Q10 values that are\nlower than typically measured or modeled on either the organismal or ecosystem level (usually 2-3,\nthe values here are equivalent to 1-2). Another example is the use of a constant e-folding scale for\nlight, even though this varies widely in the ocean. Things like this make me question whether\nappropriate constants are being used and if the model optimization is able to distinguish between\ncompensatory drivers.\nc. Also, Table S1 does not include the 47 parameters that the text says were optimized. I think that if\nother variables were unchanged (I think only changed parameters have been tabulated?), then the\ntext should be updated and the table clear on what is included and why.\nPlease, see our response to comment C about the validation of the PlankTOM12.1 model\nbeing used within our submitted manuscript. We also would like to mention that a paper is in\npreparation within our research team that will precisely compare the different versions of\nPlankTOM12 that our group developed over the past few years. This evaluation is heavy and\nout of the scope of this manuscript.\nOur objective for this manuscript was to use the PlankTOM12 version used within the latest\nGlobal Carbon Budget analysis, and that the model performances are within the range of the\nother GOBMs being used within the global carbon budget. This is now clearly explained in\n29\nsection 5.2. Therefore, in this section we have removed the optimization paragraphs and\nelaborate on the model validation instead.\n26) 403. Please elaborate on why these parameters were chosen, at least in a supplemental\ntext. That makes these choices more transparent, and helps other modeling teams\nunderstand which differences are important.\n410. I would prefer representative results to back up this claim were included.\n422. But Table S1 appears to show changed physical parameters as well...please resolve\nthis (appearance of) inconsistency.\nIn the revised manuscript, and in this section, we have removed the optimization paragraphs\nand elaborate on the model validation instead, for the reasons previously mentioned.\n29) 425. Is this different from the A-IV discussed earlier? If so please consider removing the\nnew acronym, as it is only used a couple times in the methods.\nYes, it was different. However, we have modified this section and now used the amplitude of\nthe interannual and decadal variability as described earlier.\n30) 483. What is the order of magnitude of observations required to do this robustly?\nThere is no rule for knowing exactly how many observations per year are needed for the\nhybrid approach to be robust. However, the number of observations before 1980 was\nextremely low (below the 25th percentile = 499 observations), which is why in the revised\nmanuscript we decided to use the hybrid approach from 1980 onwards.\n31) 460. a. This is a choice that I think should be clearly communicated in the main text, and\nas a mask or outline in the maps in Fig. 5 and Fig. S3. It is unclear from this whether the\nmodel outputs and data products are similarly masked out in the same coastal band, so\nplease clarify that these analyses are comparing like to like.\nb. Areas of <1000 m depth and within 3 degrees of shore are disproportionately important for\nmarine carbon cycling. Ignoring them in order to characterize more pelagic regions may\nmake it easier to identify trends and present more constrained appearing values, but at the\ncost of considering real, evnironmentaly relevant variability and significant contributions to\nthe global carbon flux. Reading this, I have diminished confidence in claims about the\nchanged magnitude of the climatological flux (and inferred land flux), trends (which may be\ncompensated or exacerbated by coastal processes), and the relative agreement of variability\nestimates. I think that this is too important a caveat to be buried in the Methods and not\ndiscussed in the Discussion.\nc. How do the results change if coastal areas are considered as well, or separately? I would\nlike to see even a quick discussion of this and followup plots in the supplement.\nuncertainty of ocean CO fluxes associated with coastal areas to the discussion, and how\n2\nthis might explain some of the discrepancies observed in our study.\n32) 494. Is this the most relevant variable to change?\nCould the authors present sensitivity analysis supporting the selection of this parameter\ninstead of/without others?\nWhat are the actual values of the perturbations used? This doesn't seem to be reported\nanywhere, but needs to be.\nIs it reasonable to assume a constant value (as shown in Table S1) or are biogeographic\nvariations in this (and other constants) expected?\nWhat are the observational constraints on this parameter and how were the magnitudes of\nchanges selected within those constraints?\nI think this choice needs significantly more explanation. There many fit results (in time and space)\nplotted outside the range of these perturbed models, which leads me to ask if the perturbations\nreally are providing useful bounds or not--and if they are, whether the hybrid approach would\nbenefit from additional constraints or flexibility in parameters. I understand that this is a poorly\nconstrained parameter, but there does seem to be a finite range that is consistent with measured\nocean DOC values, which could help justify bounds on the perturbations.\nPlease see our response to comment B.\n33) 498. Please provide a representative graphical example of this approach in the\nsupplement (or a couple, a good fit and a bad fit).\nDone\n34) 505. Related to comment on line 467, I'd like to see some more quantitative assessment\nof the statement here.\nLine 505, we mentioned years that have constrained values outside the range of\nperturbation (as understood by the reviewer, see comment #10): “Years with a constrained\nocean CO sink not within the range of the ocean CO sink from the perturbed simulations\n2 2\nare kept but considered as uncertain values”\nWe believe it is a fair approach, and we do not understand how quantitative assessment\nfrom this sentence could be provided.\n35) 514. There is no Fig. S4 in the supplement provided for review, but this would be\nwelcome and would address my request for line 470 as well.\nDone\n36) 528. What quadratic function? Is there a mathematical justification for a different\ncalculation? i.e., is this an arbitrarily larger confidence bound?\nquadratic function. The justification is that the quadratic function provides a larger confidence\ninterval, which is a reasonable approach considering the larger uncertainty related to a\nregional analysis.\n37) Fig. S1. The letters above the last two panels are repeated.\nCorrected\n38) Fig. S3. I suggest making these panels larger.\nDone\n39) In the data link, the hybrid results have a mislabeled column 1 in the tab for results by\nlatitude.\nCorrected\n40) The code looks plausible (I did not not attempt to run it), and additional annotation would\nbe appreciated to help users looking to translate this to an open source programming\nlanguage for testing.\nWe have followed the “data and code availability” policy of the journal.\n32", "reviewer_decision_letter": "Dear Editor,\nThank you for inviting a revision of our work for further consideration in your journal. We\nsincerely thank the reviewers for their invaluable comments on our manuscript. Our revision\naddresses all the points raised by the reviewers and reduces the use of acronyms and\ntechnical language to enhance accessibility.\nIn particular, following your recommendation and the reviewer’s comments, we have made\nthe following important improvements:\n● We conducted and detailed a sensitivity analysis of the hybrid approach, which\nshows that results on decadal trends are robust to the choices of parameters and\nmodel set up.\n● We have included a validation of our model based on the model evaluation process\nused in the Global Carbon Budget 2023, which demonstrates that the\nNEMO-PlankTOM12.1 model is appropriate for estimating the evolution of the\noceanic CO sink, and that it is largely in line with the results of other models.\n2\n● We have improved the description of the hybrid approach, both in the Methods and in\nthe main text.\n● We have used all available fCO observations including those from the coastal zones\n2\n(which were excluded in the original submission).\n● We have clarified the presentation of the results about the North and South regions.\n● We have included a new figure that illustrates the low robustness of the 2010s trend\nin the ocean CO sink estimated with fCO products, giving even more weight to our\n2 2\nown conclusions using our hybrid approach.\nFurthermore, since the submission of our manuscript, the Global Carbon Budget 2023 (GCB\n2023) has been published. We thus updated our analysis to use the latest available\nestimates for the NEMO-PlankTOM12.1, the Global Ocean Biogeochemistry Models\n(GOBMs) and the fCO -products published in GCB 2023, and expanded the period of\n2\nanalysis to year 2022. Our main focus is now on constraining the ocean CO sink trend\n2\nduring 2000-2022 (instead of 2000-2019 in the original submission). Finally, we limited the\nuse of the hybrid approach to start in 1980 as the number of fCO observations prior to that\n2\nwere extremely low and discontinued.\nAll these modifications strengthen our results and improve the readability of our manuscript.\nSee below our point-by-point response to each reviewer, with the reviewer's comments in\nblue, the new text in red and our response in black."}, "author_response": "Dear Editor,\nThank you for inviting a revision of our work for further consideration in your journal. We\nsincerely thank the reviewers for their invaluable comments on our manuscript. Our revision\naddresses all the points raised by the reviewers and reduces the use of acronyms and\ntechnical language to enhance accessibility.\nIn particular, following your recommendation and the reviewer’s comments, we have made\nthe following important improvements:\n● We conducted and detailed a sensitivity analysis of the hybrid approach, which\nshows that results on decadal trends are robust to the choices of parameters and\nmodel set up.\n● We have included a validation of our model based on the model evaluation process\nused in the Global Carbon Budget 2023, which demonstrates that the\nNEMO-PlankTOM12.1 model is appropriate for estimating the evolution of the\noceanic CO sink, and that it is largely in line with the results of other models.\n2\n● We have improved the description of the hybrid approach, both in the Methods and in\nthe main text.\n● We have used all available fCO observations including those from the coastal zones\n2\n(which were excluded in the original submission).\n● We have clarified the presentation of the results about the North and South regions.\n● We have included a new figure that illustrates the low robustness of the 2010s trend\nin the ocean CO sink estimated with fCO products, giving even more weight to our\n2 2\nown conclusions using our hybrid approach.\nFurthermore, since the submission of our manuscript, the Global Carbon Budget 2023 (GCB\n2023) has been published. We thus updated our analysis to use the latest available\nestimates for the NEMO-PlankTOM12.1, the Global Ocean Biogeochemistry Models\n(GOBMs) and the fCO -products published in GCB 2023, and expanded the period of\n2\nanalysis to year 2022. Our main focus is now on constraining the ocean CO sink trend\n2\nduring 2000-2022 (instead of 2000-2019 in the original submission). Finally, we limited the\nuse of the hybrid approach to start in 1980 as the number of fCO observations prior to that\n2\nwere extremely low and discontinued.\nAll these modifications strengthen our results and improve the readability of our manuscript.\nSee below our point-by-point response to each reviewer, with the reviewer's comments in\nblue, the new text in red and our response in black."}]

Decision letter
Dear Editor,
Thank you for inviting a revision of our work for further consideration in your journal. We
sincerely thank the reviewers for their invaluable comments on our manuscript. Our revision
addresses all the points raised by the reviewers and reduces the use of acronyms and
technical language to enhance accessibility.
In particular, following your recommendation and the reviewer’s comments, we have made
the following important improvements:
● We conducted and detailed a sensitivity analysis of the hybrid approach, which
shows that results on decadal trends are robust to the choices of parameters and
model set up.
● We have included a validation of our model based on the model evaluation process
used in the Global Carbon Budget 2023, which demonstrates that the
NEMO-PlankTOM12.1 model is appropriate for estimating the evolution of the
oceanic CO sink, and that it is largely in line with the results of other models.
2
● We have improved the description of the hybrid approach, both in the Methods and in
the main text.
● We have used all available fCO observations including those from the coastal zones
2
(which were excluded in the original submission).
● We have clarified the presentation of the results about the North and South regions.
● We have included a new figure that illustrates the low robustness of the 2010s trend
in the ocean CO sink estimated with fCO products, giving even more weight to our
2 2
own conclusions using our hybrid approach.
Furthermore, since the submission of our manuscript, the Global Carbon Budget 2023 (GCB
2023) has been published. We thus updated our analysis to use the latest available
estimates for the NEMO-PlankTOM12.1, the Global Ocean Biogeochemistry Models
(GOBMs) and the fCO -products published in GCB 2023, and expanded the period of
2
analysis to year 2022. Our main focus is now on constraining the ocean CO sink trend
2
during 2000-2022 (instead of 2000-2019 in the original submission). Finally, we limited the
use of the hybrid approach to start in 1980 as the number of fCO observations prior to that
2
were extremely low and discontinued.
All these modifications strengthen our results and improve the readability of our manuscript.
See below our point-by-point response to each reviewer, with the reviewer's comments in
blue, the new text in red and our response in black.