Evaluation of Regional Surface Energy Budget Over Ocean Derived From Satellites

The energy balance equation of an atmospheric column indicates that two approaches are possible to compute regional net surface energy flux. The first approach is to use the sum of surface energy flux components Fnet,c and the second approach is to use net top-of-atmosphere (TOA) irradiance and horizontal energy transport by the atmosphere Fnet,t. When regional net energy flux is averaged over the global ocean, Fnet,c and Fnet,t are, respectively, 16 and 2 Wm–2, both larger than the ocean heating rate derived from ocean temperature measurements. The difference is larger than the estimated uncertainty of Fnet,t of 11 Wm–2. Larger regional differences between Fnet,c and Fnet,t exist over tropical ocean. The seasonal variability of energy flux components averaged between 45°N and 45°S ocean reveals that the surface provides net energy to the atmosphere from May to July. These two examples demonstrates that the energy balance can be used to assess the quality of energy flux data products.


INTRODUCTION
Estimating the surface energy budget is one of the key components of understanding energy flow within the Earth system. Surface energy fluxes determine the energy input to the ocean and energy transfer through the ocean-atmosphere boundary. In addition, surface fluxes often drive processes occurring near the surface. For example, energy fluxes at the surface affect cloud processes occurring in the boundary layer (e.g., Betts, 1985;Betts and Ridgway, 1988;Albrecht et al., 1990;Bretherton and Wyant, 1997;Wood, 2012). In addition, surface radiative fluxes play a key role in determining sea ice melts (Hudson et al., 2013). Both low-level clouds (Soden et al., 2008;Loeb et al., 2018) and sea ice play a critical role determining climate feedback and Earth's energy budget (Hartmann and Ceppi, 2014).
Components of the surface energy flux are radiative flux (irradiance), turbulent flux, and flux associated with mass transfer. Surface irradiance is composed of shortwave and longwave irradiances. Downward surface shortwave irradiance is the solar irradiance transmitted through the atmosphere. Part of the irradiance that reaches the surface is reflected by the surface, which can be reflected back to the surface by the atmosphere, and the rest is absorbed by the surface. Broadband ocean surface albedo is about 0.05 (e.g., Kato et al., 2002) but the albedo depends on solar zenith angle and surface wind speed (Cox and Munk, 1955) and chlorophyll concentration (Jin et al., 2004). Downward surface longwave irradiance is the emitted irradiance by the atmosphere. Similar to the shortwave irradiance, a part of the downward longwave irradiance is reflected and the rest is absorbed by the surface. The surface also emits longwave irradiance, with a magnitude depending on the temperature and emissivity, with the latter being a function of wind speed over ocean (Sidran, 1981;Masuda et al., 1988).
Turbulent heat fluxes consist of sensible and latent heat fluxes. These are enthalpy fluxes and depend on near surface and surface properties (e.g., Cronin et al., 2019). In addition to these enthalpy fluxes, the enthalpy is transferred through the atmosphereocean boundary when water is transported by precipitation and evaporation (Mayer et al., 2017;Trenberth and Fasullo, 2018;Kato et al., 2021). From a regional ocean energy budget perspective, surface fluxes are needed to separate the energy input from that by horizontal energy transport by ocean, provided that the regional ocean temperature tendency is known (Trenberth and Solomon, 1994). Despite the importance of surface energy fluxes, there is a large uncertainty associated with surface fluxes derived from satellite observations. The uncertainty in the regional monthly mean net surface irradiance over ocean estimated by Kato et al. (2020) is 13 Wm −2 . Similarly, a typical error in a long-term mean net energy flux over ocean is 10 Wm −2 (Cronin et al., 2019). In addition, temperatures of raindrops and water vapor are needed to estimate the enthalpy transfer associated with water mass transfer (Mayer et al., 2017;Kato et al., 2021). As a consequence, when satellite derived surface flux data products are integrated to assess the surface energy balance, the residual of annual global surface energy balance is about 10-15 Wm −2 (Kato et al., 2011;Stephens et al., 2012;Loeb et al., 2014;L'Ecuyer et al., 2015;Wild et al., 2015). Therefore, annual net surface energy fluxes averaged over global ocean is one order of magnitude larger (Meyssignac et al., 2019) than the mean ocean heating rate of 0.8 Wm −2 (Johnson et al., 2016;Loeb et al., 2018;von Schuckmann et al., 2020).
Computed surface fluxes have been evaluated using observations. However, the spatial and temporal coverage of surface observations are limited. Comparisons of similar data products are useful to identify outliers and unphysical assumptions. Even though the regional energy budget bias in satellite-based data products is significant, the bias can be explained by the sum of uncertainties of all components . A less explored approach for evaluating energy budget flux products is to use horizontal energy transport by the atmosphere and to compare the resulting net surface energy flux with the sum of surface energy flux components. In this study, we integrate surface energy data products and demonstrate the use of energy divergence in the atmosphere in evaluating surface energy flux data products. In addition, we examine regional water mass balance to understand whether the regional water mass balance residual can explain the energy balance residual. Furthermore, we analyze the seasonal variability of surface energy fluxes to test whether the variability can be used in the evaluation. If the physical processes are robust, the seasonal variability might be used to evaluate the quality of energy data products. Section "Regional Water Mass and Atmosphere Energy Budget Equations" introduces regional energy and water mass budget equations that are used in this study. Data products used in this study are described in section "Data Products." Regional net surface energy flux and mass balance, as well as seasonal variability of surface energy fluxes are discussed in section "Results."

REGIONAL WATER MASS AND ATMOSPHERE ENERGY BUDGET EQUATIONS
One form of energy balance equations of an atmospheric column is the balance between energy tendency in the column and fluxes at boundaries, i.e., top-of-atmosphere (TOA), lateral boundaries, and at surface boundary (Trenberth, 1997). The energy flux at the top boundary is the net TOA irradiance R TOA . R TOA is the sum of absorbed shortwave irradiance, which is the insolation minus reflected shortwave irradiance and emitted longwave irradiance multiplied by −1 (because the net is defined downward positive in this study). Fluxes at the surface include net surface irradiance R sfc , sensible heat flux F SH and latent heat flux F LH , and enthalpy flux associated with precipitation F fallout and evaporation F v . The atmosphere can transport enthalpy, potential energy and kinetic energy through lateral boundaries. The energy budget equation of an atmospheric column is then expressed as (Mayer et al., 2017;Trenberth and Fasullo, 2018;Kato et al., 2021), where U a is the dry air horizontal velocity, c p is the specific heat capacity, T is the temperature, is the geopotential, s is the surface geopotential, k is the kinetic energy, p is the pressure, p sfc is the surface pressure, t is time, and g is the gravitation acceleration at sea level. All energy terms include contributions from dry air, water vapor and hydrometeors and they are added weighted by their respective mixing ratios r, k = k a + r v k v + r l k l + r i k i + r r k r + r s k s (2) c p = c p a + r v c p v + r l c l + r i c i + r r c r + r s c s where c p a is the specific heat capacity of dry air at constant pressure, c p v is the specific heat capacity of water vapor at constant pressure, and c is the specific heat capacity of hydrometeors. Subscripts v, l, i, r, s are, respectively, water vapor, liquid and ice cloud particles, rain, and snow. The latent heat l is the sum of the enthalpy change when water vapor is condensed or hydrometeors are evaporated, l v is the enthalpy of vaporization, l f is enthalpy of fusion, r is the mixing ratio. The latent heat is determined by the enthalpy release when water vapor is condensed or snow and ice crystals are melted. Therefore, l l in this case is the specific heat capacity multiplied by the temperature difference between cloud particles or raindrops and the reference temperature. The F h term on the right side of Eq. 1 represents the divergence of moist static and kinetic energy associated with hydrometeors moving at different velocities from the dry air velocity (Kato et al., 2021), Because F h is not provided by data products used in this study, we assume that velocities of all hydrometeors are equal to the dry air velocity (i.e., F h = 0). The enthalpy fluxes associated with water mass transfer F m are the sum of, and where c w is the specific heat capacity of either ice or liquid water, T w is 2 m wet bulb temperature, T skin is the ocean skin temperature, T 0 is 0 • C (Mayer et al., 2017;Trenberth and Fasullo, 2018),Ṗ is the precipitation rate andĖ is the evaporation rate. The energy balance Eq. 1 suggests two possible ways of computing net surface energy flux. One way is to use all the surface flux components, where F net,c is positive downward (hereinafter the component approach). The second method for deriving the net surface flux suggested by Eq. 1 is to use TOA net irradiance, atmospheric divergence, and tendency terms (hereinafter the transport approach), This approach is used by, for example, Trenberth (1997) Trenberth and Fasullo (2018), and Liu et al. (2020). Similar to the energy balance Eq. 1, the water mass balance equation for an atmospheric column is (Trenberth, 1991), where and the precipitation rateṖ =Ṗ r +Ṗ s is the sum of rain and snow rates. Eq. 11 states that the tendency and convergence of water mass is balanced with the difference of precipitation and evaporation rates at the surface. Because of the assumption of hydrometeors traveling with the dry air velocity, F r = 0. The energy balance Eq. 1 and water mass balance Eq. 11 are related through the diabatic heating term in the total energy equation of the atmospheric column (e.g., Kato et al., 2021).

DATA PRODUCTS
The TOA and surface irradiance data product used in this study is the Edition 4.1 Clouds and the Earth's Radiant Energy System (CERES) Energy Balance and Filled (EBAF). Global net TOA irradiances averaged over 10 years from July 2005 to June 2015 is adjusted to 0.71 Wm −2 based on ocean heating rates, ice warming and melts, and atmospheric and lithospheric warming (Johnson et al., 2016;Loeb et al., 2018). A detailed description of the method to produce TOA and surface irradiances included in the product and their uncertainty are given, respectively, in Loeb et al. (2018) and Kato et al. (2018). Horizontal transport of moist static energy and water vapor, as well as surface turbulent fluxes are derived from the ERA-Interim reanalysis data product (Dee et al., 2011). Although mass is not conserved in the original data product (Mayer et al., 2017;Trenberth and Fasullo, 2018;Liu et al., 2020), it is conserved in the ERA-Interim data product used in this study. Mass correction procedures and the method to compute horizontal transport are discussed in Trenberth and Fasullo (2018). Two other sets of surface turbulence fluxes used in this study are Version 3 SeaFlux data product  and the Version 3 of the Objectively Analyzed Air-sea Fluxes (OAFlux) (Yu et al., 2008). Precipitation rate data product used in this study is Version 2.3 Global Precipitation Climatology Project (GPCP) data product (Adler et al., 2012). EBAF irradiances, SeaFlux turbulent fluxes, and GPCP precipitation rates are derived using primarily satellite observations. The time period used in this study is from January 2001 to December 2016.

RESULTS
In this section, we investigate the consistency of regional surface energy budget derived by two different approaches (Eqs. 9, 10). Seasonal variabilities of surface net energy fluxes averaged over 45 • N to 45 • S ocean are compared to check the consistency of the temporal variability.

Regional Surface Energy Budget and
Water Mass Balance Figure 1 shows four surface energy flux components that appear on the right side of Eq. 1. The net irradiance R sfc is primarily a function of latitude. A positive net irradiance over the tropical ocean is largely compensated by latent heat flux F LH . Sensible heat flux F SH is much smaller than F LH . In addition, the meridional gradient of F SH is much smaller than the meridional gradient of R sfc and F LH . The spatial pattern of the enthalpy flux associated with water mass transfer F m resembles the spatial pattern of precipitation rate. When these components are averaged over the global open water, R sfc , F LH , F SH , and F m are, respectively, 126, −97, −12, and −1 Wm −2 . The sum of these four components is 16 Wm −2 , which is one order of magnitude larger than the annual global mean ocean heating rates. Note that the annual global mean ocean heating rate given in Johnson et al. (2016) of 0.68 Wm −2 is averaged over the entire global area. Therefore, the ocean heating rate averaged over the global ocean is 0.93 Wm −2 , which is estimated as the product of 0.68 Wm −2 and the ratio of the global area to the ocean area of 1.37 (=0.510 / 0.372). We ignore enthalpy transported to the global ocean by river runoff, which is about 10% of water evaporated from the ocean . The regional energy budget computed by Eq. 9 F net,c , expressed as the sum of R sfc , F LH , F SH , and F m , is shown in the left panel of Figure 2. A positive value indicates that net energy is transferred to the ocean. The right panel of Figure 2 also shows the regional energy budget computed by Eq. 10 F net,t . The net surface flux F net,t averaged over the global open water is 2.1 Wm −2 , whereas its F net,c counterpart is 16.3 Wm −2 . Mayer et al. (2021) who computed net surface energy flux over ocean with ECMWF ERA5 (Hersbach et al., 2018) by the transport approach report F net,t of 1.6 Wm −2 . Josey et al. (2013) who computed net surface energy flux over ocean with OAFlux turbulent fluxes (Yu and Weller, 2007) and International Satellite Cloud Climatology Project (ISCCP) irradiances (Zhang et al., 2004) by the component approach report F net,c of 29 Wm −2 . The difference between F net,c and F net,t is the atmospheric energy balance residual ε E , The difference between F net,c and F net,t exists over tropics (bottom left plot of Figure 2). Spatial distribution of F net,t is not necessarily correct, but a casual comparison for F net,c and F net,t showing large ε E over the tropics leads us to conclude that the tropical region is largely responsible for the larger energy budget residual from the component approach. The standard deviation over the tropics (bottom right plot of Figure 2) is nearly equal to ε E over the tropics while it is larger than ε E over mid-latitude. This indicates that seasonal variability of ε E is small over the tropics, but larger seasonal variability exists over mid-latitude. Note that the spatial pattern of ε E shown in the bottom left plot of Figure 2 differs from the spatial pattern of Figure 3 of Kato et al. (2016) because different energy balance equations of an atmospheric column, hence different energy flux data products, are used. Because F LH is equal toĖl v , Eqs. 9, and 11 suggest that a bias in the regional water mass balance can affect regional surface energy balance. In order to assess regional water mass balance, in Figure 3 we plot the divergence of water vapor as a form of the latent heat plus the tendency term on the top left panel and l v0 (Ṗ −Ė) on the top right panel where l v0 is the enthalpy of vaporization at 0 • C. In the tropics, water vapor converges toward the inter tropical convergence zone (ITCZ) where larger precipitation rates are present. For the atmosphere to balance water mass by Eq. 11 regionally, the sum of values shown in left and right plots needs to be 0. The bottom left plot of Figure 3 shows the water mass balance residual expressed in Wm −2 . Regions with negative values are where convergence is too large, precipitation rate is too small, or evaporation rate is too large. Although the size is unknown, the F h term defined by Eq. 5 that are neglected in this study can contribute to both energy and water mass balance residuals (Kato et al., 2021). While the standard deviation is relatively large over the tropics (bottom right plot of Figure 3), the spatial pattern of water mass residual differs from the spatial pattern of net surface energy flux (top left panel of Figure 2). This indicates that the bias is not limited to a common component, i.e., F LH and l vĖ , which affect both balances, but also in other energy flux components likely play a role.
To understand the effect of water mass balance residual to surface net energy flux, however, we assume that all water mass residual is caused by one component of water mass fluxes. If the residual is caused by the bias in F LH , the bias corrected F net,c is, where the water mass balance residual ε M is If the residual is caused by latent heat divergence then the bias corrected F net,t is, 0 U a c p T + + k + l dp−l v0 ε M . (16) The water mass balanced surface net energy flux is F net,c − l v0 ε M and F net,t − l v0 ε M , which are shown, respectively, on the left and right panel of Figure 4. The surface net energy flux shown in Figure 4 is the flux if regional water mass were balanced. We can further combine surface net energy fluxed derived from two approaches and water mass correction to form the water mass corrected atmospheric energy balance residual, Figure 5 shows regional ε E − l v0 ε M . The bias can be associated with any of the energy flux components that appear in Eqs 9, and 10. The bias in F LH can also contribute to ε E − l v0 ε M , provided that biases in water mass flux components are present in a way that the sum is not altering the regional water mass balance. Regions with a positive value need a larger upward surface energy flux or a smaller moist energy and kinetic energy divergence.   Figure 6 shows climatological seasonal variability of energy flux components averaged between 45 • N and 45 • S over ocean. The net TOA irradiance (blue line in the top plot) seasonal variability is similar to the seasonal variability of the climatological global mean shown in Fasullo and Trenberth (2008). The seasonal variability is largely caused by insolation driven by the Earth-tosun distance, which in turn affects absorbed shortwave irradiance. The energy divergence in the atmosphere averaged between 45 • N and 45 • S over ocean is relatively constant. During May through July, the atmosphere transports more energy than net TOA irradiance, indicating that additional net energy is provided from the surface during these months. The bottom plot of Figure 6 shows the surface flux components separately. The net surface irradiance is indicated by the blue line and three estimates of turbulent flux from SeaFlux, OAFlux and ERA-Interim are indicated, respectively, by red, magenta, and green lines. In Figure 6, the net surface irradiance is positive downward and the turbulent fluxes are positive upward. From May to July, ERA-Interim turbulent flux is larger than the net irradiance, while SeaFlux and OAFlux turbulent fluxes are smaller than the net irradiance throughout the year. The larger ERA-Interim turbulent flux than the net irradiance during May through July is interpreted as the surface providing net energy to the atmosphere, which is consistent with the top plot of Figure 6. This result suggests that SeaFlux and OAFlux turbulent fluxes are too small or both the net surface irradiance and divergence are too large during these 3 months. More importantly, annual mean turbulent fluxes averaged over 45 • N to 45 • S over ocean derived from three different products differ by nearly 20 Wm −2 , which is approximately 14% of the annual mean value. The 20 Wm −2 difference of turbulent fluxes is nearly equivalent to the difference between the net surface energy fluxes derived from two methods. Although the net surface energy fluxes derived by the two methods differ by more than 10 Wm −2 , Figure 6 indicates that differences of turbulent fluxes derived from ERA-Interim, SeaFlux, and OAFlux are nearly constant throughout the year. As a consequence, once respective annual mean values are subtracted, seasonal variabilities derived from two data products are similar (Figure 7, bottom plot). After the corresponding annual mean is subtracted, the difference in climatological monthly mean variability is less than ±2 Wm −2 , which is approximately 10% of the amplitude of the seasonal cycle.

DISCUSSION AND SUMMARY
Energy flux data products were integrated by two different methods to assess the regional surface energy balance. The first method is to use all surface flux components and the resulting net surface flux is denoted by F net,c . The second method is to use the TOA net irradiance, horizontal energy transport by the atmosphere, and tendency for which the resulting net surface flux is denoted by F net,t . While the uncertainty in regional F net,c is larger than F net,t , the advantage of the component approach is that it provides all components. The transport approach only provides the net surface energy flux. However, F net,c averaged over the global ocean is 16 Wm −2 while F net,t averaged over the global ocean is 2 Wm −2 . The comparisons of net surface energy fluxes derived from both approaches shed lights into the uncertainties of their components. For example, the net surface energy flux averaged over global ocean is nearly equal to the ocean heating rate provided that the enthalpy transported by river runoff is negligible. Because annual mean ocean heating rate is 0.93 Wm −2 , the size of combined biases from flux components used in estimating the net surface energy flux averaged over global ocean is nearly equal to F net,c . Our results also indicate that the residual of energy balance (16-0.93 Wm −2 ) is larger than the residual of mass balance (4 Wm −2 ) when they are averaged over global ocean. It is not obvious why the energy balance residual is larger than the mass balance residual. One reason might be due to the components achieving energy and mass balances. The surface energy balance is achieved largely among R sfc , F SH , and F LH while the mass balance is achieved among F LH , precipitation rate, and horizontal transport of water.
The water mass residual indirectly affects energy budget residual. Kato et al. (2021) show that regional water mass is not conserved when multiple data products are integrated. This is also illustrated in Figure 3. Horizontal transport of water vapor derived from ERA-Interim is shown on the top left plot of Figure 3. The annual mean water vapor transport from ocean to land is 4 × 10 4 km 3 yr −1 in the units of volume, which is equivalent to 3.2 PW or approximately 10 Wm −2 . This is nearly the same size of the energy flux from ocean to land (Fasullo and Trenberth, 2008), indicating that the ocean to land energy transport is dominated by latent heat divergence. Our value of water mass divergence averaged over the global ocean in the units of energy flux of 11 Wm −2 is consistent with the value given in Trenberth and Fasullo (2013) and Rodell et al. (2015). At a regional scale, water vapor divergence plus tendency in the atmosphere should balance withĖ −Ṗ (Eq. 3). This implies thaṫ P −Ė averaged over the global ocean should be approximately −10 Wm −2 , provided that the tendency term is small. When GPCP precipitation rate and SeaFlux latent heat flux are used, P −Ė averaged over global ocean is −7 Wm −2 . This suggests that biases inṖ andĖ are nearly the same size and partially cancels whenĖ is subtracted fromṖ. This may be the reason for the smaller water mass balance residual than the energy balance residual averaged over the global ocean. Note that the annual mean latent heat flux form OAFlux averaged over global ocean is approximately 5.5 Wm −2 smaller than the SeaFlux counterpart value so thatṖ −Ė averaged over global ocean is −1.5 Wm −2 , which gives a larger ε M than ε M computed with SeaFlux.
TheĖ −Ṗ value also depends on precipitation data product. Generally, precipitation rate estimated over ocean is uncertain due to lack of ground-based direct observations compared to precipitation rate estimated over land (Sun et al., 2018). A study by Behrange and Song (2020) suggests that GPCP precipitation rate over global ocean is underestimated by 9%. Precipitation rates over the ITCZ regions estimated by GPCP are generally smaller than those estimated from a similar multiple observation merged global precipitation data product of the CPC Merged Analysis of Precipitation (CMAP; Xie and Arkin, 1997). In addition, GPCP precipitation over the tropics is smaller than TRMM 3B42 precipitation (Masunaga et al., 2019). In particular, heavy rain rates larger than 30 mm day −1 over ocean occur less frequent in GPCP precipitation than TRMM 3B42v7 precipitation (Masunaga et al., 2019). The spatial pattern of water mass residual (lower left of Figure 3) suggests that the residual is not limited over the ITCZ, indicating that the water mass balance residual is not entirely caused by our choice of data products.
The surface net energy flux F net,c averaged over the global ocean of 16 Wm −2 is smaller than the maximum standard deviation of regional net fluxes of approximately 30 Wm −2 derived from 12 data products shown in Yu (2019). Therefore, regional F net,c shown in Figure 2 could change depending on the data product. Nevertheless, given the regional water vapor balance residual shown in Figure 3, the water mass balance residual contributes a large portion of the energy balance residual for some regions.
We showed that sums of latent heat and sensible heat fluxes from SeaFlux and OAFlux averaged between 45 • N and 45 • S over ocean during May through July are likely to be too small because the atmosphere transports energy more than energy input from the top and bottom boundary. If the moist static divergence is too large, the net surface irradiance is also too large. The difference between annual mean F net,c and F net,t averaged over 45 • N to 45 • S is 18 Wm −2 . The 20 Wm −2 difference of turbulent fluxes found between three products is nearly equivalent to the  difference of F net,c and F net,t . In addition, if we assume that the 2.5 Wm −2 uncertainty in regional outgoing TOA shortwave and longwave irradiances  and the 10 Wm −2 uncertainty in horizontal transport (Trenberth and Fasullo, 2018) are independent, then the uncertainty in regional F net,t is 11 Wm −2 [=(2.5 2 +2.5 2 +10 2 ) 1/2 ]. Therefore, the 18 Wm −2 is likely to be larger than the uncertainty in TOA net irradiance and transport combined.
The root-mean-square (RMS) difference of monthly mean downward irradiances compared with observations at buoy sites that are mostly located in the tropics is 11 and 5 Wm −2 , respectively, for shortwave and longwave irradiances (Kato et al., 2018). The estimated uncertainty in the regional monthly mean net surface irradiance over ocean is 13 Wm −2 (Kato et al., 2020). The uncertainty in regional turbulent flux in F net,c is difficult to estimate. A casual comparison of latent heat flux shown in If there is a random error contribution, the RMS difference in monthly means is expected to be smaller. In fact, climatological mean differences of regional F SH and F LH from SeaFLux and OAFlux are within 10 Wm −2 for most regions (Figure 8), although the differences of F SH and F LH over tropical regions have the same sign. All these results suggest that estimated uncertainties in monthly mean R sfc , F SH , and F LH are of the order of 10 Wm −2 .
As mentioned in Cronin et al. (2019), the uncertainty in satellite derived gridded energy products is larger than the in situ observation uncertainty of ∼10 Wm −2 in a long-term averaged value (e.g., monthly mean) largely due to limitations in observing surface air temperature and humidity. Near surface temperature and humidity affect both downward longwave irradiance and sensible and latent heat fluxes. In addition, ocean skin temperature retrieved from satellite observations is the temperature of the thermal skin layer of about a 0.1 mm thick below ocean surface (Wong and Minnett, 2018). Furthermore, the ocean thermal skin layer temperature is generally lower than ocean temperature below the surface at a depth of ∼2 m because of evaporative cooling at the air-sea interface (Smith et al., 1996). Therefore, satellite derived skin temperature can be different from ocean surface temperature (Cronin et al., 2019).
In addition to these input variables, there is the uncertainty associated with the turbulent flux parametrizations and the use of bulk state variables (Cronin et al., 2019;Yu, 2019). Gustiness of wind speed can alter turbulent fluxes and needs to be taken into account. Wind speed dependent ocean surface roughness length and transfer coefficients need to be improved, and their dependence on the sea state might need to be incorporated into parameterizations. Furthermore, deriving latent heat flux from satellites is challenging because of their diurnal sampling and calibration stability . Also, a larger footprint size and broad passive microwave channel weighting functions are not ideal for deriving near surface humidity with a required spatial resolution .
Advances in remote sensing techniques to derive surface air temperature and humidity will help reduce errors in surface irradiance. In addition, better observations of ocean surface temperature, wind speed at a high temporal resolution, and sea state at a global scale are needed to improve turbulent fluxes. Some promising observations to improve surface fluxes are proposed in Cronin et al. (2019) and are expected to provide indispensable observations toward reducing uncertainties in estimating regional surface energy fluxes. Our results highlight the weakness of the component approach to estimate regional net surface energy flux and identify regions where larger issues remain. Our analysis indicates that large surface energy flux biases exist in the tropics, which is consistent with earlier studies (Loeb et al., 2014;Kato et al., 2016). Improving future observations need to target the regions with large residuals, which could lead to improved understanding and reduction of these biases. In addition, toward achieving the goal of improving surface energy flux, data product integration and physical processes demonstrated in this study can be used to assess the quality of individual data products. Reducing surface energy budget uncertainty needs coordinate efforts among surface radiation and turbulent flux data providers.

AUTHOR CONTRIBUTIONS
SK wrote the first draft of the manuscript. DP contributed to the editing. FR generated the figures. SK, FR, F-LC, DP, and WS contributed to the analysis. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the NASA CERES and the NASA's Energy and Water cycle Study (NEWS) projects.