Update on the Temperature Corrections of Global Air-Sea CO 2 Flux Estimates

The oceans are a major carbon sink. Sea surface temperature (SST) is a crucial variable in the calculation of the air-sea carbon dioxide (CO 2 ) flux from surface observations. Any bias in the SST or any upper ocean vertical temperature gradient (e.g., the cool skin effect) potentially generates a bias in the CO 2 flux estimates. A recent study suggested a substantial increase (∼50% or ∼0.9 Pg C yr −1 ) in the global ocean CO 2 uptake due to this temperature effect. Here, we use a gold standard buoy SST data set as the reference to assess the accuracy of insitu SST used for flux calculation. A physical model is then used to estimate the cool skin effect, which varies with latitude. The bias-corrected SST (assessed by buoy SST) coupled with the physics-based cool skin correction increases the average ocean CO 2 uptake by ∼35% (0.6 Pg C yr −1 ) from 1982 to 2020, which is substantially smaller than the previous correction. After these temperature considerations, we estimate an average net ocean CO 2 uptake of 2.2 ± 0.4 Pg C yr −1 from 1994 to 2007 based on an ensemble of surface observation-based flux estimates, in line with the independent interior ocean carbon storage estimate corrected for the river induced natural outgassing flux (2.1 ± 0.4 Pg C yr −1 ).

However, two temperature issues might generate bias in the CO 2 flux estimates by using the SOCAT SST. The first issue is the ship's intake depth (∼5 m instead of micrometers) and the other is the location of the SST sensor (within the warm hull of the ship instead of in the unperturbed seawater).
First, the SOCAT SST represents the bulk seawater temperature, which might not be equal to the temperature at the MBL because many processes can generate vertical temperature gradients in the upper ocean. There is a temperature gradient (red line in Figure 1) in the thermal boundary layer (TBL and gray shaded area) relating to air-sea heat exchange. Infrared radiometer measurements indicate that the skin temperature at ∼10 μm depth (T Skin ) is on average ∼0.17 K (Donlon et al., 2002) lower than the subskin temperature (T Subskin , at ∼0.1-1 m depth) because the ocean surface generally loses heat through longwave radiation and latent and sensible heat fluxes (the so-called cool skin effect; e.g., Donlon et al., 2002Donlon et al., , 2007Minnett et al., 2011;Robertson & Watson, 1992;Zhang et al., 2020). Another process that might create an upper ocean temperature gradient is the diurnal warm layer effect. Water close to the surface (e.g., at 0.5 m depth) is sometimes warmer than deeper water (e.g., at 5 m depth) due to daytime solar insolation, especially under conditions of clear sky and low wind speed (Gentemann & Minnett, 2008;Prytherch et al., 2013;Ward et al., 2004). The warming leads to stabilization of the surface layer and thus helps maintain a layered upper ocean structure. The diurnal warm layer effect is not as ubiquitous as the cool skin effect (Fairall et al., 1996), and the warm layer is complex to characterize. In the absence of the warm layer effect, the bulk seawater temperature (T Bulk ) is approximately equal to T Subskin , and T Thermal (temperature at the base of the TBL) because the water below the TBL is well-mixed by turbulence.
The second issue is the potential warm bias in the SOCAT SST. The SST community has identified a warm bias in shipboard SST measurements in the ICOADS (International Comprehensive Ocean-Atmosphere Data Set; Huang et al., 2021;Kennedy et al., 2011Kennedy et al., , 2019Reynolds & Chelton, 2010). This might be because ship SST measurements are affected by engine room warming because the SST sensor is often located in the engine room or somewhere in the ship interior (Kennedy et al., 2019). The SSTs in SOCAT were almost exclusively measured by shipboard systems (98%), meaning that a warm bias also likely exists in the SOCAT SST data set. It is worth noting that the percentage of the SST data measured by research vessels in SOCAT is likely higher than in the ICOADS shipboard SST data set. The SST measured by research ships (typically external to the ship's hull) is expected to have a higher accuracy than the SST measured by commercial ships (often in the ship's interior/within the engine room), so the warm bias in SOCAT SST may well be different with the warm bias in ICOADS ship SST.
Satellite observation of SST represents a consistent estimate of subskin temperature and avoids the diurnal warm layer effect and any potential warm bias issue. Satellite SST thus has been proposed as an alternative to calculate the bulk air-sea CO 2 flux (Goddijn-Murphy et al., 2015;Shutler et al., 2019;Watson et al., 2020;Woolf et al., 2016). Results based on a satellite SST data set suggest a ∼25% increase (i.e., warm bias correction; the cool skin correction results in another ∼25% increase) in ocean CO 2 uptake compared to the flux estimate based on the SOCAT SST (Watson et al., 2020). However, the satellite SST is not measured concurrently with the fCO 2w . Colocating the 1° × 1°, monthly gridded satellite SSTs with individual fCO 2w in SOCAT might introduce extra uncertainties. In addition, various issues in satellite SSTs (e.g., cloud masking, impact of aerosol, diurnal variability, uncertainty estimation, and validation) have not been fully resolved, especially at high latitudes and in coastal and highly dynamic regions (O'Carroll et al., 2019). A comparison of eight global gap-free satellite/ blended SST products showed that their global mean ranged from 20.02°C to 20.17°C for the period 2003-2018 . Therefore, the current accuracy of satellite SST means that it probably does not allow an optimal estimate of the global air-sea CO 2 flux.
SST observations from drifting buoys are unaffected by engine room warming, and are expected to provide the best quality reference temperature to assess bias in the ship SST, and satellite SST retrievals (Huang et al., 2021;Kennedy et al., 2011Kennedy et al., , 2019Kent et al., 2017;Merchant et al., 2019;Reynolds & Chelton, 2010). This work utilizes drifting buoy SST as the reference temperature to determine the accuracy of the SOCAT SST and to correct for any bias in the SOCAT SST data set.
Subskin temperature with a cool skin correction represents the skin temperature, which can be used to calculate air-sea CO 2 flux. Watson et al. (2020) reported a ∼25% increase in ocean CO 2 uptake by considering a constant cool skin effect (−0.17 K, Donlon et al., 2002) from 1982 to 2020. In this study, the cool skin effect estimated by a physical model (Fairall et al., 1996) and by an empirical model (Donlon et al., 2002) are compared at a global scale. The updated temperature corrections are then used to estimate their impact on the global air-sea CO 2 10.1029/2022GB007360 3 of 13 flux. The revised global air-sea CO 2 flux based on an ensemble of CO 2 flux products (Fay et al., 2021) is then compared with the ocean carbon inventory (Gruber et al., 2019).

Global Air-Sea CO 2 Flux Estimates
The bulk air-sea CO 2 flux equation is: where F (mmol m −2 day −1 ) is the air-sea CO 2 flux and K 660 (cm h −1 ) is the gas transfer velocity (e.g., Wanninkhof, 2014) normalized to a Sc (Schmidt number) of 660. The Sc is defined as the ratio of the kinematic viscosity of water (m 2 s −1 ) and the molecular diffusivity of CO 2 (m 2 s −1 ). The CO 2 solubility (mol L −1 atm −1 ) at the base of the MBL and at the air-sea interface is represented by α w and α i , respectively ( Figure 1). Sc and α are calculated from seawater temperature and salinity Weiss, 1974). Sc is equal to 660 for CO 2 at 20°C and 35 psu seawater. The CO 2 fugacity (μatm) at the base of the MBL and just above the air-sea interface is represented by fCO 2w and fCO 2a , respectively.
To calculate the global air-sea CO 2 flux, fCO 2w measured at the equilibrator temperature is first corrected to the in situ bulk temperature (SOCAT SST). Seawater at ∼5 m depth (ranging from 1 to 10 m depth depending on the ship or sampling platform) is sampled from the ship's underway water intake and pumped through an equilibrator. The equilibrated CO 2 mole fraction in the air of the headspace (χCO 2w ) is measured in a gas analyzer. χCO 2w is then converted to equilibrator fugacity (fCO 2w_equ ) (Text S1 in Supporting Information S1). fCO 2w_equ is further corrected by the chemical temperature normalization (Takahashi et al., 1993) to obtain fCO 2w in the bulk seawater: where T w_bulk is the seawater temperature measured concurrently with fCO 2w at the ship's water intake at typically 5 m depth. Seawater fCO 2w measurements are then interpolated to obtain a global gap-free fCO 2w product (at 1° × 1°, monthly resolution, e.g., Landschützer et al., 2013). A global gap-free SST data set is generally one of the independent input variables for the fCO 2w interpolation process. Other variables in Equation 1 are calculated using a global gap-free SST product and related data sets (e.g., mole fraction of atmospheric CO 2 for the calculation of fCO 2a ). Finally, globally mapped fCO 2w , fCO 2a , Sc, α w , α i , and gas transfer velocity (K 660 , estimated using a global gap-free wind speed data set) are used for the CO 2 flux calculation via Equation 1. Table 1 summarizes the SST types that should be used to calculate variables in Equation 1. Sc should be calculated from the temperature utilized to derive K 660 (e.g., T Bulk for the K 660 derived from the dual-tracer method; e.g., Ho et al., 2006;Nightingale et al., 2000). The air-sea interface temperature (T Interface ) should be used for the calculation of fCO 2a and α i , while the temperature at the base of the MBL (T Mass ) should be employed to calculate fCO 2w (via Equation 2) and α w . However, Woolf et al. (2016) suggested that T Thermal might be a better temperature for calculating fCO 2w and α w . The seawater carbonate system creates a unique situation for air-sea CO 2 exchange, which does not exist for other gases. Seawater temperature changes cause chemical repartitioning of the carbonate species (CO 2 , carbonic acid, bicarbonate, and carbonate; Zeebe & Wolf-Gladrow, 2001). We find that the timescale of this repartitioning equilibration (e-folding time >10 s for typical seawater; Johnson, 1982;Zeebe & Wolf-Gladrow, 2001) is much longer than the timescale (∼1 s) of water mixing below the MBL but within the TBL, where viscous dissipation dominates the water mixing (Jähne, 2009;Jähne et al., 1987;Woolf et al., 2016). The explanation of the timescales is detailed in Text S2 in Supporting Information S1. Although Figure 1. A schematic of the upper ocean (0-10 m depth) using an example where temperature is influenced by a positive (ocean heat loss) sensible heat flux and carbon dioxide (CO 2 ) is being taken up by the ocean. The gray shaded area represents the thermal boundary layer (TBL), and the red line represents the temperature gradient in the TBL. The mass (in this case, CO 2 ) boundary layer (MBL) is embedded within the TBL. The blue line corresponds to the CO 2 concentration gradient within the MBL. The TBL is characteristically 10 times thicker than the MBL because heat is transferred about an order of magnitude quicker than CO 2 (Jähne, 2009). Sea surface temperature is a general term for all temperatures mentioned in the figure. T Interface : the temperature at the air-sea interface; T Skin : the skin temperature at ∼10 μm depth measured by an infrared radiometer; T Mass : the temperature at the base of the MBL (20-200 μm depth); T Thermal : the temperature at the base of the TBL (0.1-2 mm depth); T Subskin : the temperature of seawater below the TBL at a depth of ∼0.1-1 m such as measured by drifting buoys; T Bulk : the temperature at 1-10 m depth as measured at the typical depth of a ship's seawater intake. T Interface , T Mass , and T Thermal are conceptual (black text), whereas T Skin , T Subskin , and T Bulk are from actual measurements (practical, blue text).

of 13
there is a temperature gradient in the TBL due to the cool skin effect, the carbonate species are not expected to have time to thermally adjust, which suggests that T Thermal is the optimal temperature for calculating fCO 2w and α w .
T Thermal , T Mass , and T Interface are conceptual temperatures, which can be approximated by practical temperatures (Figure 1). Satellite SST, which represents the subskin temperature, is a good approximation for T Thermal Watson et al., 2020;Woolf et al., 2016). A satellite T Subskin product can be used to calculate α w and Sc, and to map fCO 2w for the global ocean. T Subskin with a cool skin correction can then be utilized to calculate global fCO 2a , and α i . In situ T Subskin should ideally be used to correct fCO 2w from the equilibrator temperature to the subskin seawater temperature. However, the in situ temperature measured concurrently with the fCO 2w in SOCAT is T Bulk , and in situ T Subskin measurements are unavailable to exactly match the SOCAT space and time stamp. Using in situ T Bulk (i.e., SOCAT SST) to correct fCO 2w is reasonable in the absence of a warm layer effect, but it is important to account for the potential warm bias in the SOCAT SST. Table 1 also summarizes the influence of SST and the corresponding importance for the variables used to make air-sea CO 2 flux estimates (after Woolf et al., 2016). The Sc and fCO 2a variations due to the bias in the SST product have a small influence on the global air-sea CO 2 flux. However, any bias in the SST data used for the calculation of α w , α i , and especially fCO 2w can result in a considerable bias in the flux. The temperature influence on the fCO 2w mapping should be significantly dampened by the interpolation process. The most significant influence on the CO 2 flux due to temperature bias comes from individual fCO 2w (∼160% K −1 , Table 1). An average bias of 0.1 K could result in a bias in fCO 2w of ∼1.6 μatm, which corresponds to ∼16% of the net air-sea CO 2 flux for the last decade .
The skin temperature should be used for the calculation of α i and fCO 2a . The T Skin can be obtained from T Subskin with a cool skin correction. If T Subskin is used rather than T Skin for the calculation of α i , and fCO 2a , the ocean CO 2 uptake is in theory underestimated by ∼19% for the last decade, with a mean cool skin effect of 0.17 K (Donlon et al., 2002).

Bias Assessment
The in situ bulk SST in SOCAT is generally used to correct individual fCO 2w observations from the equilibrator temperature to the seawater temperature (e.g., studies in Table S1 in Supporting Information S1). However, a warm bias might exist in the SOCAT SST due to heating in the engine room. Watson et al. (2020) colocated the DOISST v2.0 (NOAA Daily Optimum Interpolation SST data set; Reynolds et al., 2007; representing the subskin temperature) with individual in situ SST measurements in SOCAT. They found that the SOCAT SST is on average 0.13 ± 0.78 K higher than the colocated DOISST v2.0. However, Huang et al. (2021) pointed out that there might be a cold bias in the DOISST v2.0 and DOISST v2.1 products (the difference between DOISST v2.0 and v2.1 can be seen in Text S4 in Supporting Information S1).

Cool Skin Effect Estimate
The cool skin effect is ubiquitous in the ocean (Donlon et al., 2002) and should be considered when estimating air-sea CO 2 fluxes. Watson et al. (2020) used a constant value (−0.17 K) to account for the impact of the cool skin effect on air-sea CO 2 fluxes. However, the cool skin effect is affected by many environmental processes. Donlon et al. (2002) proposed a wind speed-dependent cool skin effect based on skin and bulk temperature measurements (Donlon02, hereafter). A physical model for the cool skin effect proposed by Saunders (1967) and developed by Fairall et al. (1996) considers wind speed, longwave radiation, heat flux, and solar radiation (Fairall96, hereafter). Fairall96 has been included in the COARE 3.5 model (Edson et al., 2013) and recent studies (Alappattu et al., 2017;Embury et al., 2012;Zhang et al., 2020) suggest that Fairall96 better accounts for the cool skin effect than the parameterization dependent upon a single variable (wind speed).
We employ the ERA5 wind speed data (Hersbach et al., 2020) (Kalnay et al., 1996), ERA5 monthly averaged reanalysis data sets (Hersbach et al., 2020) for wind speed, 2 m above mean sea level (AMSL) air temperature, relative humidity (calculated from 2 m AMSL air temperature and dew point temperature using the August-Roche-Magnus approximation), downward shortwave radiation, downward longwave radiation, and boundary layer height.

Global Air-Sea CO 2 Flux Estimates With the Temperature Correction
We use two different methods to account for the bias in the SOCAT SST for the global air-sea CO 2 flux estimates. For the first method, we use the buoy SST as the reference temperature to assess the bias in SOCAT SST (bias_buoy, hereafter). We correct the 1° × 1°, monthly fCO 2w in SOCAT v2021 via Equation 2 (i.e., fCO 2w_corrected = fCO 2w e −0.0423 * ΔSST ) by the temperature difference (ΔSST) between SOCAT SST and buoy SST. The ΔSST varies with latitude (with a 10° latitude running mean, see the orange line in Figure 2b) but we do not consider the variation of ΔSST over time. The number of matched data points between SOCAT SST and buoy SST is small in most years, so ΔSST is averaged from 1982 to 2020. In addition, only fCO 2w data within 70°S to 70°N are corrected because of the small number of measurements in the polar oceans. For the second method, the colocated DOISST v2.1 replaces SOCAT SST in Equation 2 to reanalyze fCO 2w (bias_OI, hereafter; Watson et al., 2020). The reanalyzed fCO 2w is used for the flux calculation (see Goddijn-Murphy et al., 2015;Holding et al., 2019 for the reanalysis process).
We employ the MPI-SOMFFN neural network technique (Landschützer et al., 2013) to interpolate the fCO 2w_corrected and the reanalyzed fCO 2w to the global ocean from 1982 to 2020, using a set of input variables. We use the same data sets as Landschützer et al. (2014) for the neural network inputs, except for the SST product. The CCI SST (Merchant et al., 2019) represents the subskin temperature and is independent of in situ SST measurements, so we utilize the 1° × 1°, monthly CCI SST v2.1 for the neural network training process. The CCI SST v2.1 is also used to calculate Sc and α w , while the CCI SST v2.1 with a cool skin correction is employed to calculate α i and fCO 2a .
We use two models (Fairall96 and Donlon02) to estimate the cool skin effect. Both Fairall96 and Donlon02 cool skin effect estimates are applied to the CCI SST v2.1 to calculate α i and fCO 2a , respectively. The quadratic wind speed-dependent formulation (K 660 = a U 10 2 ; Ho et al., 2006;Wanninkhof, 2014) is used to calculate gas transfer velocity. The 1° × 1°, monthly ERA5 wind speed data from 1982 to 2020 is utilized to scale the transfer coefficient a to match to a global mean K 660 of 18.2 cm hr −1 (equal to 16.5 cm hr −1 for K) from the 14 C inventory method (Naegler, 2009). It is worth noting that the cool skin effect and the warm layer effect do not impact the global mean K 660 calculated from the 14 C inventory because the air-sea 14 C concentration difference (Δ 14 C) is very large (Naegler, 2009;Sweeney et al., 2007), and the upper ocean temperature gradients only result in a minor change in Δ 14 C. In the end, we substitute all the variables above into Equation 1 to calculate the global air-sea 10.1029/2022GB007360 6 of 13 CO 2 flux. This study typically adopts 1 standard deviation (i.e., 1 sigma) as a representation of uncertainty unless specified otherwise.

Warm Bias in the In Situ SOCAT SST
The temperature assessment using the buoy SST suggests a cold bias in the DOISST v2.1 (0.09 K on average, standard error 4.7 × 10 −4 K) and a small warm bias (0.02 K on average, standard error 4.1 × 10 −3 K) in the SOCAT SST, which indicates that while a warm bias exists in the SOCAT SST, using the colocated DOISST would overestimate this bias in SOCAT SST (Figure 2a). Figure 2b) is small and the standard error for the temperature difference (gray shading) is large in the high latitude oceans. Therefore, we only consider data between 70°S and 70°N. The SOCAT SST minus buoy SST (ΔSST, orange line in Figure 2b) shows apparent variation with latitude. ΔSST is on average positive, but is slightly negative at 35°N and 30°S. In the northern hemisphere, ΔSST is +0.04 K near the equator and increases by +0.1 K to a maximum at 25°N and then decreases to −0.05 K at 35°N. ΔSST also increases from 35°N to a maximum of +0.15 K at 50°N and then decreases further north. The ΔSST pattern in the southern hemisphere roughly mirrors that in the northern hemisphere with a 5° northward shift.

Figure 2b shows the latitudinal variation of the bias in SOCAT SST. The number of grid cells with both SOCAT and buoy data (green bars in
It is worth noting that under-sampling affects these bias assessments for SOCAT SST. If we consider all paired cells with both buoy and SOCAT SST measurements, the warm bias is on average +0.02 K. If we only consider cells with at least 10 buoy SST and 10 SOCAT SST measurements, the warm bias is on average +0.03 K ( Figure  S2a in Supporting Information S1). The latitudinal variation of the bias is very similar no matter how many measurements are within a cell ( Figure S2b in Supporting Information S1).

Figure 2.
Latitudinal variation in sea surface temperature (SST) differences, number of matched grid cells, the gas transfer velocity (K 660 ) and the fraction of the globe's surface area covered by ocean: (a) 1° latitude average temperature difference between DOISST v2.1 and buoy SST (red line) ± 1 standard error (gray shading). The input data are from 1982 to 2020 and have a 1° × 1°, monthly resolution. Blue bars show the number of cells (5° latitude bin) containing both DOISST and buoy SST data (b) 10° latitude running mean of the temperature difference between SOCAT SST (from SOCATv2021) and buoy SST (orange line, i.e., ΔSST in the main text) ± 1 standard error (gray shading). Green bars correspond to the number of cells (5° latitude bin) containing both gridded SOCAT and buoy SST; (c) 1° latitude average K 660 (purple line) calculated with a wind speed-dependent parameterization (Ho et al., 2006) using the ERA5 wind speed data (Hersbach et al., 2020) for the global ocean. The blue-shaded area corresponds to the fraction of ocean area in different latitudes (1° latitude average).
It is important to consider latitudinal variation when correcting for bias in SOCAT SST. For instance, SOCAT SST has a relatively large warm bias (thus a large bias in the fCO 2w ) in the Southern Ocean (south of 35°S, Figure 2b), which coupled with a high K 660 and a large surface ocean area (Figure 2c) results in a substantial bias in Southern Ocean CO 2 flux estimates. This study uses a latitude-varying temperature bias (i.e., the orange line in Figure 2b) to correct the air-sea CO 2 flux between 70°S and 70°N (see Section 2.4). Figure 3 shows the cool skin effect estimated by Donlon02 and Fairall96. The Fairall96 estimate of the cool skin effect is stronger than the Donlon02 estimate for low wind speeds (U 10 < 9 m s −1 ) but weaker for high wind speeds (9 m s −1 < U 10 < 16 m s −1 ) (Figure 3a). The monthly wind speed distribution (green bars in Figure 3a) shows that wind speeds less than 9 m s −1 account for 80% of the wind conditions. Therefore, the cool skin effect estimated by Fairall96 is typically stronger than that estimated by Donlon02. The standard deviation of the Fairall96 cool skin effect is much higher at low wind speeds than at high wind speeds, which reflects that the drivers (longwave radiation, heat flux, and solar radiation) can produce substantial variations in the cool skin effect under relatively calm conditions. The Donlon02 cool skin effect only has a slight latitudinal variation that is not substantially different from a constant (−0.17 K) value (Figure 3b), which was used by a previous study for air-sea CO 2 flux correction (Watson et al., 2020). In contrast, the Fairall96 cool skin estimate shows a clear latitudinal variation with two relatively small cool skin effect regions at around 50°S and 50°N where wind speeds are high. The Fairall96 cool skin effect is stable in the tropical zone and decreases toward both poles to ∼50° and then increases at even higher latitudes.

The Cool Skin Effect
In most ocean regions, the Fairall96 cool skin effect follows variations in wind speed. Intriguingly, the Fairall96 cool skin effect is nearly constant within the tropical and subtropical zones, even though the wind speed is much lower near the equator than in the subtropics. Drivers other than wind speed (i.e., latent and sensible heat fluxes, and longwave radiation) might counteract the low wind speed effect in this area.  (Fairall et al., 1996, solid blue line), the Donlon02 wind speed-dependent empirical model (Donlon et al., 2002, dashed blue line) and a constant value (−0.17 K, gray line; Donlon et al., 2002). The light blue-shaded area in both subplots indicates one standard deviation of the bin averages in Fairall96 cool skin estimates. Global ocean 1° × 1° monthly data sets are used to estimate the cool skin effect (see Section 2.3).

Variation in the CO 2 Flux Correction
In this section, we discuss the impact of the warm bias and cool skin effects on global air-sea CO 2 flux estimates. The corrections are applied over time (between 1982 and2020, Figures 4a and4b) and by latitude (Figures 4c  and 4d).
The bias correction using the buoy SST assessment (bias_buoy) leads to an average increase in ocean CO 2 uptake of 0.19 Pg C yr −1 , while the bias correction utilizing the colocated DOISST (bias_OI) suggests an average increase of 0.43 Pg C yr −1 (Figure 4a). Adopting the cool skin correction from Fairall96 and Donlon02 increases the 1982-2020 average ocean CO 2 uptake by 0.39 Pg C yr −1 and 0.43 Pg C yr −1 , respectively (Figure 4b). A constant cool skin correction of −0.17 K increases the flux by an amount similar to using the Donlon02 correction. Zhang et al. (2020) show that the mean difference between the Fairall96 cool skin effect and the observed cool skin effect (7,239 observations) is 0.04 K. If we take this value as the uncertainty of the Fairall96 cool skin estimate, the corresponding relative uncertainty in the Fairall96 flux correction is ∼20% (i.e., 0.08 Pg C yr −1 ). In total, the flux correction using the bias_buoy and Fairall96 is on average ∼0.3 Pg C yr −1 lower than if the bias_OI and Donlon02 are used from 1982 to 2020. The interannual variation in the net air-sea CO 2 flux with different temperature corrections is shown in Figure S4 in Supporting Information S1. Figures 4a and 4c show the change in the air-sea CO 2 flux (ΔFlux) generated by correcting for the warm bias in SOCAT SST. The temporal and latitudinal variation of the two flux corrections (bias_buoy and bias_OI) follow similar patterns, but the magnitude is different. Using bias_OI creates a ΔFlux that is twofold larger (in absolute terms) than that using bias_buoy. The data in Figure 2a suggest that using bias_OI may overestimate the bias in SOCAT SST, which would result in a ∼0.25 Pg C yr −1 overestimation of the air-sea CO 2 flux correction. Therefore, we favor the bias_buoy correction over the bias_OI correction.
While we use the same latitude-varying temperature difference (i.e., bias_buoy) to correct the bias in SOCAT SST every year, the flux correction shows clear interannual variation (green line in Figure 4a). A possible reason is that the number of measurements in each year of SOCAT is different ( Figure S2 in Supporting Information S1), and their spatial distribution differs between years. The latitude-dependent bias correction, when applied to the Negative ΔFlux values represent increased ocean CO 2 uptake. Green and red lines represent ΔFlux due to the bias correction assessed by drifting buoy SST (bias_buoy) and by colocated DOISST (bias_OI), respectively. Blue and purple lines represent ΔFlux due to the Fairall96 and the Donlon02 cool skin corrections, respectively. ΔFlux in (a and b) is the global annual mean, while ΔFlux in (c and d) is the long-term average  in 1° latitude bins. Results are based on the MPI-SOMFFN fCO 2w mapping method (Landschützer et al., 2013) (See Methods). The interannual variation of the global air-sea CO 2 flux with different temperature corrections can be seen in Figure S4 (Supporting Information S1). Our preferred corrections are bias_buoy for warm bias in SOCAT SST and Fairall96 for the cool skin effect (see Section 4.1).
different year-to-year spatial distribution in the SOCAT data, results in a time-varying annual mean bias correction ( Figure S2 in Supporting Information S1).

Figures 4b and 4d
show the change in air-sea CO 2 flux when accounting for the cool skin effect using Fairall96 and Donlon02 models. Figure 4b indicates an increase over time in both flux corrections (absolute value), which is driven by the increase in fCO 2a (see Equation 1 and Table 1). The impact of the cool skin effect on the air-sea CO 2 flux is through α i * fCO 2a . The ever rising atmospheric CO 2 concentration and thus fCO 2a , result in the growing cool skin flux correction.
The flux correction using Donlon02 exceeds that by Fairall96 by ∼0.05 Pg C yr −1 (in absolute terms). The largest difference in flux between the two cool skin corrections occurs in the Southern Ocean (Figure 4d). The Donlon02 cool skin effect has minimal latitudinal variation, so the flux correction is largest at ∼50°S where the gas transfer velocity is maximum and the ocean area is relatively large (Figure 2c). The Fairall96 cool skin effect has an apparent latitudinal variation and a minimum (absolute) value at ∼50°S (Figure 3). This minimum cool skin effect offsets the maximum wind speed and large ocean area, resulting in a smaller flux correction (in absolute terms) at ∼50°S for Fairall96 than for Donlon02. Recent work (Alappattu et al., 2017;Embury et al., 2012;Zhang et al., 2020) has suggested that the Fairall96 cool skin model is better than Donlon02 at capturing the cool skin effect at a global scale and this, coupled with our estimates, indicates that using the Donlon02 model may lead to an overcorrection of the air-sea CO 2 flux, especially in the Southern Ocean.

Implications for Air-Sea CO 2 Flux Estimates
This study deals with the potential bias in the fCO 2w -based air-sea CO 2 flux estimates due to upper ocean temperature effects. A large amount of uncertainty in this fCO 2w -based flux also comes from the gas transfer velocity . The air-sea CO 2 flux estimated from the ocean carbon inventory (Gruber et al., 2019) does not require the gas transfer velocity, is unaffected by upper ocean temperature effects, and provides an independent estimate of ocean CO 2 uptake. To compare the fCO 2w -based net air-sea CO 2 flux with the anthropogenic air-sea CO 2 flux of the ocean carbon inventory, we need to adjust for river-induced CO 2 outgassing. The riverine carbon flux has been estimated as 0.23 Pg C yr −1 (Lacroix et al., 2020), 0.45 Pg C yr −1 , 0.65 Pg C yr −1 (Regnier et al., 2022) and 0.78 Pg C yr −1 (Resplandy et al., 2018). Here, we adopt the mean of these values (0.53 ± 0.21 Pg C yr −1 ).
The net air-sea CO 2 flux derived from the ocean carbon inventory from 1994 to 2007 is −2.1 ± 0.4 Pg C yr −1 (i.e., −2.6 Pg C yr −1 anthropogenic flux plus 0.53 Pg C yr −1 river carbon flux; see the footnote of Table 2 for the propagated uncertainty) (Gruber et al., 2019), which is shown in Table 2 along with the ensemble mean of eighteen fCO 2w -based fluxes (Fay et al., 2021). Fluxes from six fCO 2w products and three wind speed products (three wind products are used for each fCO 2w product) are utilized to generate the ensemble mean flux, where missing fCO 2w has been filled with a scaled climatology and gas transfer velocity (K 660 ) has been calibrated to a global average of 18.2 cm hr −1 over the ice-free ocean based on 14 C-bomb flux estimates (Fay et al., 2021). All six fCO 2w products (which include the MPI SOMFFN method) have been developed from the SOCAT v2021 data set. So the corrections to the ensemble mean flux for the temperature effects should be similar to the corrections in this study based on the MPI-SOMFFN fCO 2w mapping method (Landschützer et al., 2013). Furthermore, an ensemble Net air-sea CO 2 flux estimates (Pg C yr −1 ) Flux without a temperature correction Flux with warm bias correction Flux with warm bias and cool skin correction bias_buoy bias_OI bias_buoy + Fairall96 bias_OI + Donlon02 Ensemble mean of fCO 2w -based fluxes a −1.7 ± 0.4 −1.8 ± 0.4 −2.0 ± 0.4 −2.2 ± 0.4 −2.4 ± 0.4 Ocean carbon inventory b −2.1 ± 0.4 Note. Here, bias_buoy and bias_OI represent the bias correction (to SOCAT sea surface temperature (SST)) using the assessment from buoy SST and colocated DOISST, respectively. Fairall96 (Fairall et al., 1996) and Donlon02 (Donlon et al., 2002) correspond to the cool skin effect estimated by the physical and the empirical models, respectively. We favor the bias_buoy and Fairall96 corrections (see Section 4.1).
a The ensemble mean of the fluxes from six fCO 2 products and three wind speed products (Fay et al., 2021  of different data interpolation methods and different wind products provides a more robust flux estimate than a single interpolation method based on a single wind product. The flux corrections estimated in this study are applied to the ensemble mean flux. The ensemble mean air-sea CO 2 flux without any bias and cool skin corrections (−1.7 ± 0.4 Pg C yr −1 ) is 0.4 Pg C yr −1 lower than the net flux estimate from the ocean carbon inventory. The ensemble mean CO 2 flux with bias_buoy and Fairall96 cool skin corrections is −2.2 ± 0.4 Pg C yr −1 , similar to the ocean carbon inventory derived net ocean CO 2 uptake. The corrections using the bias_OI and the Donlon02 suggested by a previous study (Watson et al., 2020) push the ensemble mean air-sea CO 2 flux (−2.4 ± 0.4 Pg C yr −1 ) toward the lower limit of the ocean carbon inventory flux estimate (Table 2). However, these comparisons depend on the choice of the riverine carbon flux correction. The riverine flux is still an unresolved issue and the flux estimates span from 0.23 Pg C yr −1 to 0.78 Pg C yr −1 Lacroix et al., 2020;Regnier et al., 2022;Resplandy et al., 2018). Without knowing which of the riverine flux estimates is most accurate, an average is simply taken here. Therefore, an accurate estimate of the river flux is required to increase our confidence for the comparison above.
Another question is whether the warm bias and cool skin flux corrections conflict with our understanding of air-sea CO 2 fluxes. One might argue that the preindustrial ocean and atmosphere would have been in a natural equilibrium (i.e., the global total of steady state natural air-sea CO 2 fluxes would have been zero; see Hauck et al., 2020 for details), but the temperature corrections would create a preindustrial ocean carbon sink. However, the warm bias in SOCAT SST is not a natural phenomenon and should not affect the preindustrial flux estimate. Furthermore, while cool skin is a natural phenomenon, the flux correction due to the cool skin effect includes both natural and anthropogenic contributions. Figure 4b shows that the cool skin flux correction decreased almost linearly by ∼0.1 Pg C yr −1 (from −0.34 to −0.43 Pg C yr −1 ) due to the increase in atmospheric CO 2 (∼70 ppm or μmol mol −1 , from 341 to 414 ppm) from 1982 to 2020 (Dlugokencky & Tans, 2018). Preindustrial atmospheric CO 2 was ∼260-280 ppm (Wigley, 1983), which is ∼70 ppm lower than atmospheric CO 2 in 1982. Thus, the preindustrial natural air-sea CO 2 flux correction due to the cool skin effect could be ∼−0.25 Pg C yr −1 , with the remaining correction (∼-0.2 Pg C yr −1 in 2020) due to the increase in atmospheric CO 2 by anthropogenic emissions.
A flux correction for the cool skin effect is only related to the fCO 2w observation-based flux estimate, which is available from the 1980s onwards (Friedlingstein et al., 2022). There were no fCO 2w measurements in preindustrial times, so the total preindustrial air-sea CO 2 flux (the sum of steady state natural flux and river flux) is based on model studies, theory, and lateral transport constraints (Hauck et al., 2020). Although the cool skin effect might result in an ∼−0.25 Pg C yr −1 flux, we can still assume that the ocean and atmosphere were in a natural equilibrium in preindustrial times. Specifically, the cool skin effect has been implicitly included in the preindustrial natural equilibrium assumption. Therefore, this study improves our understanding by suggesting an increasing anthropogenic contribution to the air-sea CO 2 flux while there is no contradiction between the temperature correction and the preindustrial natural equilibrium assumption.
The cool skin effect and its impact on the air-sea CO 2 flux have been discussed for decades. While the cool skin effect itself has been well observed and modeled, its impact on the air-sea CO 2 flux is mainly based on theoretical arguments. We still lack strong observational evidence to confirm the need to include the cool skin effect on estimates of air-sea CO 2 flux-an important topic we urge the community to demonstrate experimentally. The eddy covariance method (e.g., Dong et al., 2021) provides direct flux measurements that could be used as a reference CO 2 flux to assess the accuracy of the bulk CO 2 flux. Long-term eddy covariance measurements at a place with |ΔfCO 2 | ∼ 0 would be insightful because the relative effect of cool skin on the bulk CO 2 flux is in theory more prominent for regions of low |ΔfCO 2 |. Appropriate laboratory experiments may yield further insight.
In summary, this work updates the temperature corrections to the fCO 2w -based air-sea CO 2 flux estimates. It shows that there is a slight warm bias in SOCAT SST and a latitude-varying cool skin effect, resulting in ∼0.6 Pg C yr −1 additional ocean CO 2 uptake from 1982 to 2020. The corrected air-sea CO 2 flux for an ensemble of six gap-filled air-sea CO 2 flux products agrees well with the ocean carbon inventory derived net flux. The extreme sensitivity of the air-sea CO 2 flux to the accuracy of SST means that we should carefully choose the reference temperature to assess any bias in the SOCAT SST. The importance of the Southern Ocean for atmospheric CO 2 uptake, and the strong winds encountered there mean that large scale assessments need a suitable model for the cool skin correction to the air-sea CO 2 flux.