Advances in the Estimation of Global Surface Net Heat Flux Based on Satellite Observation: J-OFURO3 V1.1

The reliability of surface net heat flux data obtained from the latest satellite-based estimation [the third-generation Japanese Ocean Flux Data Sets with Use of Remote Sensing Observations (J-OFURO3, V1.1)] was investigated. Three metrics were utilized: (1) the global long-term (30 years) mean for 1988–2017, (2) the local accuracy evaluation based on comparison with observations recorded at buoys located at 11 global oceanic points with varying climatological characteristics, and (3) the physical consistency with the freshwater balance related to the global water cycle. The globally averaged value of the surface net heat flux of J-OFURO3 was −22.2 W m−2, which is largely imbalanced to heat the ocean surface. This imbalance was due to the turbulent heat flux being smaller than the net downward surface radiation. On the other hand, compared with the local buoy observations, the average difference was −5.8 W m−2, indicating good agreement. These results indicate a paradox of the global surface net heat flux. In relation to the global water cycle, the balance between surface latent heat flux (ocean evaporation) and precipitation was estimated to be almost 0 when river runoff from the land was taken into consideration. The reliability of the estimation of the latent heat flux was reconciled by two different methods. Systematic ocean-heating biases by surface sensible heat flux (SHF) and long wave radiation were identified. The bias in the SHF was globally persistent and especially large in the mid- and high latitudes. The correction of the bias has an impact on improving the global mean net heat flux by +5.5 W m−2. Furthermore, since J-OFURO3 SHF has low data coverage in high-latitudes areas containing sea ice, its impact on global net heat flux was assessed using the latest atmospheric reanalysis product. When including the sea ice region, the globally averaged value of SHF was approximately 1.4 times larger. In addition to the bias correction mentioned above, when assuming that the global ocean average of J3 SHF is 1.4 times larger, the net heat flux value changes to the improved value (−11.3 W m−2), which is approximately half the original value (−22.2 W m−2).


INTRODUCTION
Surface net heat flux, defined as the total heat exchange between the atmosphere and oceans, affects both atmospheric and oceanic processes. In addition, the surface net heat flux determines the actual state of atmospheric-ocean interaction and the climate system. Therefore, the surface net heat flux is an essential climate variable (ECV) and an essential observational variable (EOV). Accurate and observational estimations are required globally (Cronin et al., 2019). Global estimates based on observations are necessary to understand long-term climate change and related responses, in addition to validating climate model results. Consequently, estimating surface net heat flux using satellite observations and improvements are of vital importance.
Several efforts have been made to estimate based on satellite observations (e.g., Pinker et al., 2014) and the data products are available, but how reliable is the satellite estimation of surface net heat flux? This question is not self-evident. This is because the satellite-based surface net heat flux estimation is obtained by combining the output of the turbulent flux estimation and the radiation flux estimation, which are being promoted as separate research projects. Although each product has been previously evaluated in several studies (e.g., Andersson et al., 2011;Rutan et al., 2015;Bentamy et al., 2017), there are few studies related to net heat flux. Therefore, it must be evaluated as a surface net heat flux.
A recent study evaluated the estimation of the surface net heat flux resulting from ocean reanalysis as well as atmospheric reanalysis and satellite-based estimations (Valdivieso et al., 2017). The results showed that the ocean reanalysis gave close to 0 for the global mean surface net heat flux, while the atmospheric reanalysis and satellite-based estimates indicated that the global mean had a bias of heating the ocean. In addition, a comparison with observations from local buoys indicates that the satellitebased estimation was in good agreement, but ocean reanalysis estimates had a bias of cooling the ocean. Yu (2019) reported that modifying the bulk equation of the turbulent heat fluxes improved the global heat balance in satellite-based estimations (OAFlux-HR). However, it was indicated that a large physical inconsistency regarding freshwater balance occurs when using the modified equation.
These two studies highlight the "paradox" of surface net heat flux estimations. This occurs because of the poor agreement between the global heat balance and the local accuracy in addition to similar inconsistencies between the heat and freshwater balances.
In this study, the latest satellite-derived surface net heat flux dataset, J-OFURO3 (Tomita et al., 2019) is evaluated using the following three metrics: (1) the long-term (30 years) mean, (2) local accuracy, and (3) physical consistency. In addition, the advancement in the satellite data will be estimated by comparing the current data with the previous generation dataset, J-OFURO2 . The number of buoys used for comparison has also increased from those in past studies because of the inclusion of buoys in mid-and high-latitude areas. Through these efforts, the state of the latest satellite-based surface net heat flux estimations is better understood. Finally, suggestions for future improvements are provided.

Satellite-Derived Air-Sea Heat Flux Datasets
J-OFURO (https://j-ofuro.isee.nagoya-u.ac.jp) is a research project on estimating surface heat, momentum, and freshwater fluxes based on satellite remote-sensing observations. The project also provides the global dataset for the research community. Although the first dataset (Kubota et al., 2002) did not cover the entire global region, the second-generation dataset, J-OFURO2  provided global surface net heat flux data for 1988-2008, with their own turbulent heat flux estimation and surface radiations obtained from ISCCP (Rossow and Schiffer, 1991). From here, we refer to the J-OFURO2 dataset as J2.
The third-generation dataset: J-OFURO3 (Tomita et al., 2019) was first released as V1.0 for 1988-2013. The J-OFURO3 is characterized by the use of multi-satellite, multi-microwave sensors, and the state-of-the-arts estimation algorithm (e.g., Tomita et al., 2018). During the data period of 1988-2013, data from the satellite microwave radiometer sensors: SSMI/SSMIS series, TMI, AMSR-E, and AMSR2 were used to estimate atmospheric specific humidity, which is essential for estimating latent heat flux. In addition to above microwave radiometers, data from microwave scatterometers: ERS AMI series, QuikSCAT, ASCAT, and OSCAT series were used to estimate ocean surface winds. In addition to these estimates of surface specific humidity and surface winds, using ensembles obtained from 12 types of global SST products including satellite observation data and surface air temperature data obtained from atmospheric reanalysis data, the turbulent heat fluxes were calculated. J-OFURO3 provides a global net heat flux using its own surface turbulent heat flux estimates and surface radiation flux estimates utilizing ISCCP FD and CERES SYN1D (Doelling et al., 2013(Doelling et al., , 2016Loeb et al., 2018). Please note that the J3 upward long wave radiation flux has been recalculated using J-OFURO3 ′ s ensemble sea surface temperature data for consistency with other fluxes in J-OFURO3. This procedure is also the same as that of J2, despite the sea surface temperature data being different. Furthermore, J-OFURO3 calculates the evaporation from the ocean surface based on the J-OFURO3 surface latent heat flux and provides data on the freshwater flux in combination with the data for precipitation obtained from GPCP (Adler et al., 2003, version 2.3). These data were also used to confirm the consistency of the surface heat flux with the hydrological cycle. More details on J-OFURO3 V1.0 can be found in Tomita et al. (2019) and the official data documentation (Tomita, 2017).
The latest version, J-OFURO3 V1.1 with some updates including source data version updates, minor algorithm changes, and extended data periods covering 30 years for 1988-2017 have been released. For surface radiation data in V1.0, we found a temporal discontinuity in 2000. This temporal discontinuity is caused by changing input source data (i.e., from ISCCP to CERES products). Therefore, in the V1.1 we adjusted the radiations of ISCCP to CERES, assuming CERES is well calibrated. From here, we refer to the J-OFURO3 V1.1 dataset as J3.
The surface net heat flux (NHF) is calculated as the sum of the following components: net shortwave radiation (SWR), net long wave radiation (LWR), surface latent heat flux (LHF), and sensible heat flux (SHF), that is, NHF = SWR + LWR + LHF + SHF. In this study, all heat fluxes assumed to be positive when they are directed upward, away from the ocean surface to the atmosphere.
In this study, evaluation of J3 for 1988-2017 was the main focus, but to confirm the progress from J2, a comparison for 1988-2008 was also conducted. For J2, the monthly data of the 1-degree grid was used, and for J3, the monthly data of the 0.25-degree grid was used.
Furthermore, we have used another satellite-based product for comparison, namely: HOAPS-4.0 (Andersson et al., 2017). HOAPS is characterized by the unique development of both precipitation and evaporation (LHF) using SSMI and SSMIS series observations. The EUMETSAT Satellite Application Facility on Climate Monitoring (CM SAF) provides monthly global data from July 1987 to December 2014, with a spatial resolution of 0.5.

Calculation of Global Long-Term Average
Because the satellite-derived air-sea net heat flux datasets used in this study are gridded data, the globally averaged value is indicated as the area-weighted average value obtained from the data of each original grid size. The "global" means the region 0-360E, 90S-90N. However, it should be noted that the J2 and J3 data do not include data over land and sea ice areas. The global long-term average is calculated by the arithmetic mean over time after calculating the global area-weighted average.

In situ Observation Data
Buoy data were used to obtain the sea truth of the surface net heat flux. To obtain the surface net heat flux, the buoy measurements must provide a dataset of all components to estimate surface heat fluxes. Although there are few such buoys having sensors for radiation measurement, there are 11 in the global oceans that capture varying climatological characteristics (Figure 1, Table 1). These buoys are part of the following observation networks: ocean climate stations and the global tropical moored buoy array in NOAA/PMEL, the Ocean Reference Stations in the WHOI. The KEO, PAPA, and NTAS buoys are in the North Pacific region, and three TAO buoys  are in the tropical Pacific Ocean. There are two RAMA buoys (McPhaden et al., 2009) in the Indian Ocean. The PIRATA (Servain et al., 1998) and WHOTS buoys are in the Atlantic Ocean. STRATUS (Weller, 2015) is the only buoy in the Southern Hemisphere. The Southern Ocean Flux Station (SOFS, Schulz et al., 2012) and Agulhas Return Current (ARC) buoys are located in the Southern Hemisphere, but because they do not provide sufficient observational data, they were excluded from the main comparison of this study.
The four surface heat flux components (SWR, LWR, LHF, and SHF) were calculated from the hourly observation data of each buoy. Subsequently, the daily average values were derived after the flux calculations. Furthermore, the monthly averaged value was calculated from the daily averaged value of the flux data. The NHF was calculated from the monthly average value, and if any components were missing, all the components were set as the missing values.
The flux calculation was performed according to the method described by Tomita et al. (2010). The net SWR was calculated from the observed downward SWR according to Equation (1): where α is the surface albedo, and the climatological monthly mean values on each grid obtained from the ISCCP have been used in this study. The net LWR was calculated from the observed downward LWR and the calculated upward LWR value from the sea surface temperature (SST) according to Equations (2) and (3), where ε is the emissivity at the ocean surface, set as 0.984 following Konda et al. (1994), and σ is the Stefan-Boltzmann constant (5.679 × 10 −8 W m −2 K −4 ). For the LHF and SHF, the bulk flux calculation algorithm, COARE 3.0  was used. The input parameters required in the flux calculation using the algorithm are air temperature, humidity, winds, SST, and sea level pressure. For all parameters, the observed values at each buoy were used. The algorithm also requires the observation height of each parameter.  Supplementary Table 1 was used. For the SST, a skin/warm layer correction was not conducted.

The observation height information for each buoy listed in
The accuracy of the data on surface net heat flux obtained from in situ buoys is estimated to be 8 W/m 2 on average (Colbo and Weller, 2009;Cronin et al., 2019), while the values are slightly higher on a daily scale. For each flux component, the long-term averaged accuracy for SWR, LWR, LHF, and SHF are estimated 5.0, 3.9, 1.5, and 5.0 W/m 2 , respectively. It can also be confirmed that the mid-latitude buoys (KEO) almost exhibit the same range .

Comparison
The buoy data are point values while the satellite data are gridded. Therefore, we compared the values on gridded satellite data that include the locations of the buoys with the values calculated from the buoy measurements. All comparisons were conducted monthly. The statistics: bias, RMS, and correlation coefficient, r, were calculated for each flux component using the following equations: where s and b are the satellite gridded value and buoy point value, respectively, and n is the number of monthly data at each buoy (see Table 1). It should be noted that the RMS is defined as a form in which the bias is removed from the difference between s and b (Taylor, 2001). All statistics values are available as Supplementary Data.

RESULTS
Global Long-Term Mean Figure 1 shows the distribution of the long-term (1988-2017) mean for the global NHF obtained from J3. The figure represents the climatological features of the distribution for the NHF. In the tropical zone, a net heat flux exists from the atmosphere to the ocean, and in mid-and high-latitudes, there is a net heat flux from the ocean to the atmosphere. Examining regional features, there is a larger heat flux from the atmosphere to the ocean in the eastern tropical Pacific and at the equator. These areas contain upwelling ocean currents. In addition, there is large net heat flux from the ocean to the atmosphere at the western boundary current region for both hemispheres. Moreover, there is a strong flux contrast corresponding to the ocean fronts in these areas.
The global long-term average value calculated from J3 is −22.2 W m −2 . This indicates that the net heat flow is to the ocean surface. Although the characteristics of the qualitative distribution are not significantly different from common knowledge, this value is more than one order of magnitude larger from the viewpoint of global surface heat balance (NHF → 0).
The results are similar when compared with the global average value obtained from the previous generation dataset (J2). J2 data    Comparison With Buoys Figure 4 shows the bias between the J3 estimates and buoy observations. At all buoy points except RAMA (15N, 90E), the bias is negative, indicating a significant ocean heating bias in J3. In particular, large negative biases were found in buoys in the western tropical Pacific (0N, 165E), central tropical Pacific (0N, 170W), and central tropical Indian Ocean (0N, 80.5E). The overall averaged bias is −5.8 W m −2 , and the averaged bias at mid-high latitudes excluding buoys in the tropical zone is −3.1 W m −2 , which shows good agreement.
There are cases in which the positive and negative biases of each component cancel each other out ( Figure 4B). For example, at KEO, SW, and LH show positive biases, while LW and SH show negative biases. The sum of absolute biases of the components is 27.6 W m −2 , which is significantly larger than the absolute NHF bias of 1.5 W m −2 . Similar canceling out of biases was seen in Stratus, RAMA (15N, 90E), and TAO (0N, 140W).
Unexpectedly, Figure 4B indicates that the SH bias is always negative for these data, while biases of the other components show both negative and positive biases depending on the buoy locations. A comparison using more comprehensive global buoy data that can calculate turbulent heat flux confirms these negative SH biases, especially over the open ocean area (Tomita et al., 2019). The influence of SH biases on the global long-term mean of NHF is discussed in the "Discussion" section. In addition, a pattern of characteristic bias was also observed in the LW. Relatively large negative biases in LW were found in the KEO, STRATUS, and WHOTS buoys located in the subtropics and mid-latitudes. The negative LW biases were relatively small at buoys in the tropics. The KEO and STRATUS buoy networks correspond to areas that have significant cloudiness consisting of low-level clouds.
The same comparisons were made using both J2 and J3 data for the period up to 2008. An improvement from the previous generation data was confirmed. Figure 5 shows the bias, RMS, and the correlation for the NHF for J2 and J3. From the data in this figure, significant improvements in statistics from J2 were confirmed. At most buoy points, the RMS and correlation coefficients are improved. For the bias, on average, the absolute bias of J3 is slightly higher than that of J2. However, J2 has a small bias due to canceling out of the large positive and large negative biases. Such large positive and negative biases were improved in J3, while a slight negative bias was observed.

Consistency With Freshwater Flux
The LHF, which is one of the surface heat flux components, is proportional to the evaporation rate from the ocean surface and is a part of the freshwater flux as a counterpart for precipitation. Therefore, by checking the consistency between the latent heat flux (evaporation) and the surface freshwater balance one can evaluate the surface heat balance from a different perspective.
In general, the global ocean freshwater flux defined as (evaporation minus precipitation) is positive because a significant amount of water evaporates from the ocean. The evaporated water is transported to land by atmospheric advection, mainly in the form of water vapor. On a long-term average, the changes in atmosphere disappear and the net positive freshwater flux over the ocean is balanced by the runoff of river water from the continent into the ocean (i.e., freshwater flux-runoff = 0).
The global long-term mean value of ocean evaporation in J3 was 3.4 mm/day, while the precipitation over the ocean obtained from GPCP V2.3 was 3.0 mm/day. Therefore, the global longterm mean of freshwater flux was calculated as 0.4 mm per day. The result showed a good balance after considering the runoff from land. Various estimates have been obtained by studies on river runoff. These values range from approximately 0.27 to 0.34 mm/day (Schlosser and Houser, 2007). More recent studies estimate river runoff as 0.29 (Ghiggi et al., 2019) and 0.31 (Wilkinson et al., 2014) mm/day. According to these previous studies, if we assume a value of 0.3 mm/day of river runoff, the freshwater balance estimated from J3, GPCP, and the river runoff is 0.1 mm/day, which is a reasonable result. An improvement was confirmed compared to the estimation using J-OFURO2 (Iwasaki et al., 2014).
Although there are various global precipitation datasets (Kidd and Huffman, 2011), GPCP is used as the standard dataset in numerous studies (e.g., Andersson et al., 2011;Tapiador et al., 2017;Yu, 2019;Gutenstein et al., 2021). However, most studies suggest that much of the uncertainty in water balance lies in precipitation products as well as evaporation. To confirm the differences in the results that depend on the satellite products, we reconfirmed the results using another satellite precipitation/evaporation product, HOAPS-4.0. Consequently, it was confirmed that the long-term mean precipitation of HOAPS-4.0 for 1988-2014 (2.9 mm/day) was slightly smaller than that of GPCP V2.3 for the same period (3.0 mm/day). For evaporation, the long-term mean HOAPS-4.0 value for 1988-2014 was 3.4 mm/day. Therefore, the global ocean freshwater flux value (0.5 mm/day) was slightly higher than that estimated by J3 (0.4 mm/day, for 1988-2014), while being sufficiently comparable.

DISCUSSION
The global long-term mean value of the NHF from the J3 data was consistent with the previous generation data. The tendency of ocean heating is similar to other estimates such as satellite and ocean reanalysis. However, its value was large in magnitude, −22.2 W m −2 , indicating a significant negative imbalance. In contrast, by comparison with local 11 buoys, the average bias was found to be −5.8 W m −2 and the negative largest value of −15.7 W m −2 was found at buoy in the western tropical Pacific. Therefore, the relationship between the global mean value of J3 and the local bias was inconsistent, and the paradox of surface heat flux was confirmed. In previous studies (Pinker et al., 2017;Valdivieso et al., 2017), several tropical and subtropical buoys data (Stratus) were used for the evaluation of surface net heat flux. In contrast, for this study, the comparison was performed using more buoy data which included mid-high latitude buoys (KEO and PAPA). These data contained longer time series, and the paradox was confirmed. The cause is discussed in the following text.
In contrast to the excessively large negative imbalance for the NHF in the global long-term mean, the globally averaged ocean surface freshwater flux estimated by surface latent heat flux (evaporation) in J3, GPCP, and runoff was almost 0. This result was consistent with a previous estimate (Trenberth et al., 2007). Because the sum of these independently estimated components was close to 0, J3 the LHF is considered to be very reliable. Therefore, the cause of the imbalance might be other than the LHF. The comparison of the J3 LHF with more comprehensive global buoys ( Figure 6A) revealed that the total bias was fairly small (<1 W m −2 ) while there are some regional biases. This fact also strongly suggests that the cause of the excessively large global long-term mean imbalance of NHF is likely to be other than LHF.
In contrast to the LHF, the J3 SHF had a persistent bias. As shown in the Results (Section "Comparison with Buoys"), the SHF showed negative biases in the comparisons with all of the 11 buoys. Negative biases were also confirmed by comparison with more comprehensive buoy observations ( Figure 6B). There are negative biases in almost all open ocean areas except for the coastal area, and a larger negative bias occurs especially in midhigh latitudes (Tomita et al., 2019). Figure 7 shows the bias of the SHF as a function of latitude. In order to investigate the effect of this SHF bias characteristic on the global averaged value, this bias was corrected by using a fitting curve and the global averaged value was recalculated. The global long-term mean value of the SHF without bias correction was +8.1 W m −2 , while the global long-term average value after bias correction was 13.3 W m −2 . This bias correction has an impact of improving the global mean of the NHF by +5.5 W m −2 , but a large imbalance of −16.7 W m −2 still remains. However, the number of buoy observations on which this bias correction was based does not completely cover global oceans (as seen in Figure 6B). It is necessary to consider a more robust correction method in the future. As shown in Figure 7B, the cause of this SHF bias is in the air temperature. The J3 uses atmospheric reanalysis data instead of satellite retrieval for air temperature estimation, and it is desirable to refer to better air temperature estimates or develop a satellite-based retrieval method in the future.
Furthermore, the data coverage of J3 SHF over high-latitudes is small. In the presence of sea ice, J3 cannot calculate the turbulent heat flux; therefore, the estimation of turbulent heat flux over regions with sea ice is overlooked. In a simple test performed using ERA5 (Hersbach et al., 2020), which has complete global coverage, the global ocean average value is approximately 30% smaller than the original ERA5 value when the sea ice area is excluded for simulating the coverage of J3. When including the sea ice region, the global ocean average SHF value is approximately 1.4 times larger. This is a reasonable result considering the large air-sea temperature difference (i.e., large SHF) with sea ice at high latitudes. The same test for LHF does not give the same result. In the case of LHF, the value corresponding to the sea ice region does not have a large influence on the global ocean averages, and by including the sea ice region, the global averaged value decreases slightly. In addition to the bias correction mentioned above, when assuming that the global ocean average of J3 SHF is approximately 1.4 times larger, the NHF value changes to the improved value (−11.3 W/m 2 ), which is approximately half the original value (−22.2 W/m 2 ). This indicates the limits of microwave satellite-based flux products such as J3 and the importance of considering the value of SHF over the sea ice region in the global ocean heat balance.
The LWR also had a notable bias characteristic. The LWR bias was relatively small in the tropics, while it was comparatively large in the subtropics and mid-latitudes. For example, the largest biases were found in the KEO and Stratus data. These were −5.7 and −4.7 W m −2 , respectively. A detailed comparison was performed to investigate the cause in detail. Figure 8 shows the bias in the upward and downward components of the long wave radiation described by Equation (2). As shown in Equation (3), the upward component of the LWR is not completely independent of the downward component, but the bias shows small negative values (<1 W m −2 ). The major factor of the LWR bias is the downward component. Validation of the downward LWR assessments was performed using more comprehensive buoy observations (Rutan et al., 2015;Kato et al., 2018). Our results are consistent with their results, indicating that the bias is <5 W m −2 ; however, the spatial characteristics of the error have never been investigated thoroughly.
The buoy locations of Stratus and KEO are known as oceanic areas frequently covered by low-level clouds (e.g., Klein and Hartmann, 1993). In contrast, for the high latitude area of the North Pacific, PAPA, which is also characterized by low-level clouds, the downward LWR shows a good agreement. A more detailed investigation of the bias and the relationship with clouds, air temperature, and sea surface temperature will be needed better understand this phenomenon.
In the above, we discussed the possibility of large biases outside the 11 buoys (which showed relatively good agreement on the global heat balance). As another possibility, we discuss the effect of the difference due to the bulk formula and the associated calculation method on the global heat balance. In general, the selection of a bulk formula has a major influence on the estimation of the global turbulent heat and momentum fluxes. Based on a comparative study (Brunke et al., 2003;Iwasaki et al., 2010), COARE 3.0 is used in J3 and other satellite products. Brodeau et al. (2017) estimated that changes in the bulk formula will affect the global heat balance by 10%, and the use of different bulk formulae may significantly change the global heat balance. However, it is necessary to pay attention to consistency with other physics by changing the bulk formula. Yu (2019) confirmed that although the global heat balance was improved by changing the bulk formula, the change in LH caused a freshwater imbalance.   (Tomita and Hihara, 2017). Note that the data in near coastal region (the distance from coastline < 200 km) were removed in this comparison.
Similar results are expected for J3. Improvements in the bulk formula for turbulent heat fluxes are needed, while maintaining consistency with other physics.
In this study, we focused on the long-term mean surface net heat flux. The daily satellite-derived data set is very useful for analyzing the flux variation over time-scales varying from several days to inter-annual or decades. This type of analysis was not in the scope of this research. Weller (2015) showed that precise and long-term buoy observation revealed long-term flux trends. It will be useful for the verification of satellite data in the future. However, the number of buoys will still be small to understand the overall characteristics of these fluctuations.
In general, the uncertainty of precipitation and evaporation from the ocean is a major challenge in understanding of the water cycle. Improving satellite-based products should address this challenge. In this study, the state-of-the-art satellite-based products, J-OFURO3 (Tomita et al., 2019) and HOAPS-4.0 (Andersson et al., 2017), were confirmed to be consistent with each other in the estimation of freshwater flux, which confirms the improvement of satellite products. Estimating and providing uncertainty information is also important. An approach that combines satellite and ocean observations with the estimation of atmospheric energy transport derived from atmospheric reanalysis data is also a powerful tool to better estimate global surface fluxes (e.g., Liu et al., 2015Liu et al., , 2017Carton et al., 2018). In the future, it will be necessary to combine multiple approaches, while improving satellite products by their comparison with such approaches.

AUTHOR CONTRIBUTIONS
HT designed the study, led the research, interpreted the data, wrote the manuscript, and managed the project. HT analyzed global satellite-based surface heat flux datasets. KK, TH, and HT collected in situ buoy observation data and calculated surface heat fluxes. KK and HT compared satellite estimates with buoy observations. TH contributed some useful scripts for this study. All authors contributed to improving the manuscript.