Harmonization of global surface ocean pCO2 mapped products and their flux calculations; an improved estimate of the ocean carbon sink

Harmonization of global surface ocean pCO2 mapped products and their flux calculations; an improved estimate of the ocean carbon sink Amanda R. Fay, Luke Gregor, Peter Landschützer, Galen A. McKinley, Nicolas Gruber, Marion Gehlen, Yosuke Iida, Goulven G. Laruelle, Christian Rödenbeck, Jiye Zeng 5 1 Columbia University and Lamont Doherty Earth Observatory, Palisades NY, USA 2 Institute of Biogeochemistry and Pollutant Dynamics, ETH Zurich, Zürich, Switzerland 3 Max Planck Institute for Meteorology, 20146 Hamburg, Germany 4 Laboratoire des Sciences du Climat et de l'Environnement, Institut Pierre Simon Laplace, Gif‐Sur‐Yvette, France 5 Atmosphere and Ocean Department, Japan Meteorological Agency, 1-3-4 Otemachi, Chiyoda-Ku, Tokyo 100-8122, Japan 10 6 Dept. Geoscience, Environment & Society (DGES), Université Libre de Bruxelles, Bruxelles, CP16002, Belgium 7 Biogeochemical Signals, Max Planck Institute for Biogeochemistry, P.O. Box 600164, Hans-Knöll-Str. 10, 07745 Jena, Germany 8 National Institute for Environmental Studies (NIES), 16-2 Onogawa, Tsukuba, Ibaraki, 305-8506, Japan

global mean uptake of CO2 from the atmosphere of 1.92 +/-0.35 PgC yr -1 . This work aims to support the community effort to perform model-data intercomparisons which will help to identify missing fluxes as we strive to close the global carbon budget.

Introduction
Surface ocean partial pressure of CO2 (pCO2) observations play a key role in constraining the global ocean carbon sink. This is because variations in surface ocean pCO2 is the driving force governing the exchange of CO2 across the air-sea interface, which is commonly described through a bulk formula (Garbe et al. 2014;Wanninkhof 2014): wind speed data products. Utilizing the MPI-SOMFFN pCO2 product (Landschützer et al. 2020a) and a quadratic parameterization (Wanninkhof 1992) they find flux estimates that diverge by 12% depending on the choice of wind speed 70 products. Additionally, they find regional discrepancies to be much more pronounced than global differences, specifically highlighting the equatorial Pacific, Southern Ocean, and North Atlantic as regions most impacted by the choice of wind product. Roobaert et al. (2018) stress that to minimize the uncertainties associated with the wind speed product chosen, the global coefficient of gas transfer must be individually calculated for each (Wanninkhof 1992(Wanninkhof , 2014. In this work, we assess the impact of wind speed product choice and scaling on six pCO2 products' calculated air-sea flux estimates. By applying a 75 consistent flux calculation methodology to each pCO2 product, we minimize the methodological divergence of fluxes within the ensemble. SeaFlux provides a more consistent approach specifically targeting the most commonly used pCO2 data products to deliver an end product for consistent intercomparisons within assessment studies such as the Global Carbon Budget (Friedlingstein 80 et al. 2020;Hauck et al. 2020). By first addressing differences in spatial coverage between the observation-based products we are able to better present a true global pCO2 estimate for each product. This SeaFlux package also provides a means to normalize the gas transfer velocity to a consistent 14 C inventory. By calculating fluxes using multiple scaled gas transfer velocities for different wind products, we present a methodologically consistent database of air-sea CO2 fluxes. The SeaFlux package is an ensemble data product along with documented code allowing the community to reproduce consistent flux 85 calculations from various data-based pCO2 reconstructions.

Methods
The SeaFlux method is based on six observation-based pCO2 products and spans years 1988-2018 (Table 1). These six include three neural network derived products (MPI-SOMFFN, CMEMS-FFNN, NIES-FNN), a mixed layer scheme product 90 (JENA-MLS), a multiple linear regression (JMA-MLR), and a machine learning ensemble (CSIR-ML6). These select products are included as they have been regularly updated to extend their time period and incorporate additional data that comes with each annual release of the SOCAT database.
All of these methods provide full three-dimensional fields (latitude, longitude, time) of the sea surface partial pressure of 95 CO2 (pCO2) and the air-sea CO2 flux. In their original form each product may utilize different choices for the inputs to Equation 1 (Table A1). In this work recompute the fluxes using the following inputs to the bulk parameterization approach Equation 1: kw is the gas transfer velocity (further discussed in Sect. 2.3), sol is the solubility of CO2 in seawater, in units mol m -3 uatm -1 , calculated using the formulation by Weiss (1974), EN4 salinity (Good et al. 2013), Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) sea surface temperature (Good et al. 2020), and European Centre for Medium-100 Range Weather Forecasts (ECMWF) ERA5 sea level pressure (Hersbach et al. 2020); ice is the sea ice fraction from OSTIA  (Good et al. 2020); pCO2 is the partial pressure of oceanic CO2 in atm for each observation-based product after filling as discussed in Sect. 2.1, and pCO2 atm is the dry air mixing ratio of atmospheric CO2 (xCO2) from the ESRL surface marine boundary layer CO2 product available at https://www.esrl.noaa.gov/gmd/ccgg/mbl/data.php (Dlugokencky et al. 2017) multiplied by ERA5 sea level pressure (Hersbach et al. 2020) at monthly resolution, and applying the water vapor correction 105 according to Dickson et al. (2007).
Flux is defined positive upward, i.e., CO2 release from the ocean into the atmosphere is positive, and uptake by the ocean is negative. In the following sections we discuss the three steps that have the greatest impact on the inconsistencies between unadjusted flux calculations in the six pCO2 products and the approach that we utilize for the SeaFlux ensemble product. 110

Step 1: Area filling
Machine learning methods aim to maximize the utility of the existing in situ observations by extrapolation using various proxy variables for processes influencing changes in ocean pCO2. Extrapolation with these independently observed variables is possible due to the nonlinear relationship between pCO2 in the surface ocean and the proxies that drive these changes.
However, not all of the proxy variables have complete global ocean coverage for all months, so the resulting pCO2 products 115 are limited by the extent of the proxy variables ( Figure 1). Additionally, in coastal regions there is the potential that different relationships of pCO2 are expected than in the open ocean, thus limiting the extrapolations. In contrast, the mixed layer scheme (utilized by the JENA-MLS product) does not suffer from such missing areas but does not distinguish between coastal and open ocean. While the area extent of the available air-sea flux estimates varies between products, there are consistent patterns; nearly all products cover the open ocean, whereas larger differences exist in the coverage of coastal 120 regions, shelf seas, marginal seas and the Arctic Ocean.
To account for differing area coverage, past studies (Friedlingstein et al. 2019(Friedlingstein et al. , 2020Hauck et al. 2020) have adjusted simply by scaling based on the percent of the total ocean area covered by each observation-based product. This does not account for the fact that some areas have CO2 flux densities that are higher or lower than the global average (Table 1,3). 125 Thus, the magnitude of the adjustment by area-scaling is likely an underestimate (McKinley et al. 2020). One specific example is the northern high latitudes where coverage by the six products varies substantially. Similarly, three products provide estimates in marginal seas such as the Mediterranean while the other three products have no reported pCO2 values here. Shutler et al (2016) report that subtle differences in regional definitions can cause differences of >10% in the calculated net fluxes. year changes in pCO2 in the missing areas (given that the extended MPI-ULB-SOMFFN product is a monthly climatology centered on year 2006) and is obtained as follows.
To extend the open and coastal merged monthly climatology (MPI-ULB-SOMFFN) to 1988-2018, we calculate a global 140 scaling factor based on the product-based ensemble mean pCO2 for regions which are covered consistently by all six pCO2 products. We first mask all pCO2 products to a common sea mask before taking an ensemble mean (pCO2 ens ). Next, we divide this ensemble mean by the MPI-ULB-SOMFFN climatology (pCO2 clim ) at monthly 1° by 1° resolution (Equation 2).
The monthly scaling factor (sfpCO2) is calculated by taking the mean over the spatial dimensions.

145
The scaling factor calculation can be represented as where 2 is the one-dimensional scaling factor (time dimension), 2 is the ensemble mean of all pCO2 products at 150 three-dimension, monthly 1° by 1° resolution, 2 is the MPI-ULB-SOMFFN climatology, also at three-dimension but limited to just one climatological year. The x and y indicate that we take the area-weighted average over longitude (x) and latitude (y) resulting in the monthly scaling value. If a product mean is exactly equal to the climatology mean, the scaling factor is 1. Value ranges from 0.91 to 1.06 over the 31-year time period.

155
The one-dimensional scaling factor is then multiplied by the MPI-ULB-SOMFFN climatology for each spatial point resulting in a three-dimensional scaled filling map. These values are then used to fill in missing grid cells in each observation-based product.
Globally, the adjustments are all less than 20% of the total flux, with the mean adjustment for the six products at 9%. In the Northern Hemisphere however, the filling process can drive adjustments of up to 35% (Table 3). As expected, the 160 observationally-based products with more complete spatial coverage tend to have smaller flux adjustments, however the impact on the final CO2 flux depends on the ∆pCO2 and wind speed of the areas being filled (Figures 2-3, Table 1,3). The only product that does not change during this adjustment process is the JENA-MLS mixed layer scheme-based product (Rödenbeck et al. 2013) which is produced with full spatial coverage and therefore needs no spatial filling.

165
Our approach is not without its own assumptions and limitations. We rely on a single estimate of the missing pCO2 in coastal ocean regions, marginal seas, and the full Arctic Ocean, given that this is the only publicly available product currently existing. Nevertheless, the fact that common missing areas along coastal regions and marginal seas are reconstructed using specific coastal observations provides a step forward from the linear-scaling approach currently used by the Global Carbon Budget (Friedlingstein et al. 2019(Friedlingstein et al. , 2020. Further confidence is provided by previous research showing that climatological 170 relevant signals, i.e. mean state and seasonality, are well reconstructed by the MPI-SOMFFN method (Gloege et al. 2021).
Furthermore, our scaled filling methodology assumes that pCO2 in the missing ocean regions is increasing at the same rate as the common area of open-ocean pCO2 used to calculate the scaling factor. Research from coastal ocean regions and shelf seas reveal that, in spite of a large spatial heterogeneity, this is a reasonable first order approximation (Laruelle et al. 2018). 175 While our approach has a constant scaling factor for the missing ocean areas regardless of latitude we acknowledge that this could be improved with increased understanding.

Step 2: Wind product selection
Historical wind speed observations (including measurements from satellites and moored buoys) are aggregated and 180 extrapolated through modeling and data assimilation systems to create global wind reanalyses. These reanalyses are required to compute =air-sea gas exchange. Air-sea flux is commonly parameterized as a function of the gradient of CO2 between the ocean and the atmosphere with wind speed modulating the rate of the gas exchange (Equation 1). Each of these wind reanalyses has strengths and weaknesses, specifically on regional and seasonal scales ( Table A2.

Step 3: Calculation of gas exchange coefficient 190
We employ the quadratic windspeed dependence of the gas transfer velocity (Wanninkhof 1992) and calculate the piston velocity (kw) for each of the wind reanalysis products as where the units of kw are in cm h -1 , Sc is the dimensionless Schmidt number, and 〈 2 〉 denotes the second moment of average 10-m height winds (m s -1 ). We choose the quadratic dependence of the gas transfer velocity as it is widely accepted and used in the literature (Wanninkhof, 1992 Ocean confirm that even in this high wind environment, a quadratic parameterization fits the observations best (Butterworth & Miller 2016). Future updates of the SeaFlux product will include options for other parameterizations.
We calculate the square of the wind speed at the native resolution of each wind product and then average it to 1° by 1° monthly resolution (see Table A2). The order of this calculation is important as information is lost when resampling data to lower resolutions because of the concavity of the quadratic function. For example, if the second moment were calculated 205 from time-averaged wind speeds, it would result in an underestimate of the gas transfer velocity (Sarmiento and Gruber 2006;Sweeney et al. 2007). The resulting second moment is equivalent to <U 2 > = Umean 2 + Ustd 2 where Umean and Ustd are the temporal mean and standard deviation calculated from the native temporal resolution of U.

210
In addition to the choice of wind parameterization, large differences in flux can result due to the scaling coefficient of gas transfer (a) that is applied when calculating the global mean piston velocity. This constant originates from the gas exchange process studies (Krakauer et al. 2006;Sweeney et al. 2007;Müller et al. 2008;Naegler 2009) which utilize observations of radiocarbon data from the GEOSECS and WOCE/JGOFS expeditions (Key et al. 2004). The 14 C released from nuclear bomb testing (hence bomb-14 C) in the mid twentieth century has since been taken up by the ocean. The number of bomb-14 C atoms 215 in the ocean, relative to the pre-bomb 14 C, can thus be used as a constraint on the long-term rate of exchange of carbon between the atmosphere and the ocean. A probability distribution of wind speed is used to optimize the coefficient of gas transfer based on these observed natural and bomb 14 C invasion rates. This coefficient must be individually calculated and is not consistent for each wind product. Further, the gas transfer velocity used by the different pCO2 mapping products are not scaled to the same bomb-14 C estimate (Table A1). The range of the different bomb-14 C estimates is within the range of the 220 uncertainty from the associated studies (Naegler, 2009), but the choice would introduce inconsistency that is easily addressed here.
We scale the gas transfer velocity to a bomb-14 C flux estimate of 16.5 cm hr −1 as recommended by Naegler (2009). The 225 coefficient (a) is calculated for each wind product via a cost function which optimizes the coefficient of gas transfer where parameters are as defined in Equation 3 40% range (Wanninkhof 2014). By determining the optimal a coefficient for each of the reanalysis winds, uncertainty in the global fluxes can be decreased. Our scaled coefficients (Table 2) correspond well with the estimate of Wanninkhof (2014) who uses the CCMP wind product to estimate a as 0.251. Differences in the coefficient will also result from the time period considered and definition of global area and ice fraction applied in the calculation. 235 This scaling of the gas exchange coefficient (a) for each wind product is an essential, and an inconsistently applied step (Table A1), that has large implications for air-sea flux estimates (Figure 4). Without individual scaling, and instead utilizing a set value for the gas transfer coefficient (a) regardless of wind product, our results show that calculated global fluxes could be as high as 9% different depending on which pCO2 and wind reanalysis product considered (Roobaert et al. 2018). 240

Further parameters for flux calculation
The remaining parameters of Equation 1 are the solubility of CO2 in seawater (sol), the atmospheric partial pressure of CO2 (pCO2 atm ), and the area weighting to account for sea ice cover. While the choices of products used for these parameters can also result in differences in flux estimates, the impacts are much smaller as compared with the parameters discussed above. 245 Atmospheric pCO2 is calculated as the product of surface xCO2 and sea level pressure corrected for the contribution of water vapor pressure. The choice of the sea level pressure product, or absence of the water vapor correction can have small, but not insignificant, impact on the calculated fluxes. Additionally, some products utilize the output of an atmospheric CO2 inversion product (e.g. CarboScope, Rödenbeck et al. 2013;CAMS CO2 inversion, Chevallier, 2013) which can introduce differences 250 in the flux estimate outside of the sources related to a product's surface ocean pCO2 mapping method. Importantly, we do not advocate that our estimate of pCO2 atm is an improvement over other estimates thereof; rather we provide an estimate of pCO2 atm that has few assumptions and leads to a methodologically consistent estimate of ∆pCO2. We maintain the same philosophy in our estimates of solubility of CO2 in seawater and sea-ice area weighting and therefore we do not elaborate on them here. 255

SeaFlux air-sea CO2 flux calculation
Following Equation 1, CO2 flux is calculated individually for each of the six observation-based products with each available wind product (CCMPv2, ERA5, JRA55) as discussed in Sect. 2.2 (Table 4). Since we account for spatial coverage 260 differences via our filling method (Sect. 2.1), taking a global mean flux for each of the data products is now straight forward. one example (other products in Figure A2). The three wind products show very consistent fluxes throughout the time series, however the importance of appropriate scaling of the gas exchange coefficient (a) is evident by the significant differences between global mean fluxes calculated with unscaled and scaled a value (Figure 4). It is clear that the impact of applying the 265 appropriate gas exchange coefficient through proper scaling has a larger impact on the resulting flux time series than solely the choice of wind product.

SeaFlux ensemble flux
By calculating each product's flux using these consistent methods, we permit for a more accurate comparison of fluxes and increase confidence in the SeaFlux product ensemble mean flux estimate of -1.92 +/-0.35 PgC yr -1 (Table 4). Here, the stated 270 uncertainty represents 2 as calculated from the 18 realizations of flux included in the SeaFlux ensemble (six pCO2 products and three wind products). This result is further strengthened by the use of multiple wind products which we consider to be independent estimates for the purpose of the uncertainty calculation.
These flux values will be different from those produced by the observation-based pCO2 product's original creator, both spatially and on the mean ( Figure 5, Table A1, A3). However, by calculating fluxes using this standardized approach we 275 have higher confidence in the uncertainties and in the ensemble mean of global fluxes.

Issues not addressed by SeaFlux
While the SeaFlux data set allows us to standardize much of the calculation of air-sea carbon flux, the community is still working towards consensus on other issues that impact this estimate. One source of uncertainty has been raised by Watson et al. (2020) who contend that a correction should be applied to pCO 2 observations to account for the vertical temperature 280 gradient between the ship water intake depth and the surface skin layer where gas exchange actually takes place. A further correction should be applied when calculating fluxes to account for the "cool skin" effect caused by evaporation (Woolf et al. 2016;Watson et al. 2020). Applying these corrections results in an increasing CO2 sink by up to 0.9 PgC yr −1 (Watson et al. 2020). Here, we do not take such corrections into account for two reasons. Firstly, the skin temperature correction to pCO2 needs to be applied directly to the measurements and not the final interpolated pCO2 from the data products. Hence, it 285 is up to the developers of the SOCAT dataset and the developers of the pCO2 mapping products to decide on the inclusion of this correction. It would then be up to the developers of the data products to update their mapped products. Secondly, the cool skin correction would be equally applied to all methods and would not contribute to the inconsistencies that we are trying to address here. As the ocean carbon community moves towards consensus on such issues, the SeaFlux product will be updated to include revised protocols. Fnet itself, which is our focus.

Conclusions
We introduce a standardized approach for flux calculations from observationally-based pCO2 products. The SeaFlux approach for flux calculations from available surface ocean pCO2 estimates enhances consistency and comparability for this ensemble of products. Specifically, we address the two largest sources of divergence, namely the differences in spatial 310 coverage between the products, and the scaling of the gas transfer velocity for available wind speed products based on global 14 C-based constraints. The area adjustment is the largest contributor to the methodological discrepancies, resulting in an increase in CO2 uptake of 0-20% relative to the original, possibly incomplete coverage (depending on pCO2 product). The global scaling of the gas transfer velocity can change the CO2 flux on average by 6% relative to non-standardized flux calculations. The impact of applying the appropriate gas exchange coefficient through proper scaling has a larger impact on 315 the resulting flux time series than solely the choice of wind product. By accounting for these sources of differences, the global mean calculated air-sea carbon flux calculated from the six available products is adjusted by up to 24%. The ensemble mean air-sea carbon flux is estimated to be -1.92 +/-0.35 PgC yr -1 with the uncertainty representing 2σ as calculated from the 18 realizations. applications, such as coverage over longer time periods, higher spatial or temporal resolution, or runs incorporating further auxiliary data sets or pCO2 data (e.g., SOCCOM float data, Bushinsky et al. 2019).
Along with the ensemble of CO2 flux fields, we also provide a public-use coding package allowing users to apply the presented standardized flux calculations to own data-based pCO2 reconstructions. 330

Product
Area     Table A3: Mean fluxes, PgC yr -1 , 1988-2018 for each observational pCO2 product. Mean flux calculated from unfilled (filled) coverage pCO2 map and unscaled (scaled) gas exchange coefficient; calculated for 3 wind products (CCMP2, ERA5, JRA55) with the average shown here. Percent change is calculated as the difference between the unfilled/unscaled and filled/scaled as a fraction of the filled/scaled; 690 does not indicate an error in the product's flux but is a representation of the impact the filling and scaling can have on the end flux estimate. The mean flux as reported in the original pCO2 product is included for comparison ( Figure 5)