Accurate estimation of satellite-derived ocean latent heat flux (LHF) at high spatial resolution remains a major challenge. Here, we estimate monthly ocean LHF at 4 km spatial resolution over 5 years using bulk algorithm COARE 3.0, driven by satellite data and meteorological variables from reanalysis. We validated the estimated ocean LHF by multiyear observations and by comparison with seven ocean LHF products. Validation results from monthly observations at 96 widely distributed buoy sites from three buoy site arrays (TAO, PIRATA, and RAMA) indicated a bias of less than 7 W/m2 with R2 of more than 0.80 () and with a King–Gupta efficiency (KGE) of over 0.84. Our estimated ocean LHF also performs well in simulating annual variability and predicting between-site variability, as indicated by a bias of lower than 6 W/m2 and an R2 of more than 0.84 (). Overall, the average KGE for estimated ocean LHF increased by 18%–23% compared to other LHF products, indicating robust LHF estimation performance. Importantly, our estimated annual ocean LHF has similar global spatial distribution compared to other LHF products, although there are general differences in LHF values due to the difference in the models and the spatial resolution.

1. Introduction

Ocean latent heat flux (LHF) plays an important role in exchanges of energy and water at the interface of the atmosphere and ocean. Knowledge of these fluxes is essential for understanding the water and energy budget and the linkage between ocean heat fluxes and large-scale climate change [13]. Since the 1980s, moored buoys and ships have been widely used to measure ocean LHF [47]. However, sparse observations hamper the accurate characterization of spatiotemporal global ocean LHF patterns over large spatial scales.

Remote sensing can provide spatially and temporally continuous information on ocean state variables over large scales and has been considered to be the most viable method for estimating global distributed LHF [810]. Currently, most existing satellite-based global ocean LHF products are derived from bulk method with meteorological quantities, such as the Goddard Satellite-Based Surface Turbulent Fluxes (GSSTF), the Hamburg Ocean-Atmosphere Parameters and Fluxes from Satellite Data (HOAPS), and the Japanese Ocean Flux Data Sets with Use of Remote Sensing Observations (J-OFURO). To evaluate these turbulent flux products, Bourras [11] compared five satellite-based products of ocean LHF with the observed data from 75 moored buoys, but validation results demonstrated that the J-OFURO, HOAPS-2, Bourras–Eymard–Liu (BEL) dataset, and GSSTF-2 satellite products had moderate systematic errors; in addition, the Jones Fluxes product had a large systematic deviation with respect to Tropical Atmosphere Ocean (TAO) data. Moreover, the existing satellite-based products, such as HOAPS-3, have high accuracy but a rather coarse spatial resolution of 1° [12, 13]. Later studies [1416] showed the demand of accurate ocean surface LHF with higher spatial resolution (e.g., 0.25°) for oceanography research. Other existing global LHF products (including reanalysis and objectively analyzed products), such as the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-I), provided us with global spatial coverage and high temporal resolution, as well as a large amount of error caused by algorithms and parameterization schemes of the boundary layer [17, 18].

During the last 40 years, there has been a lot of effort to develop models for estimating global ocean LHF. In general, the methods used to derive ocean LHF can be divided into four categories.: (1) eddy covariance methods [1921], physically based models, estimated surface fluxes using direct measurements of temporally continuous information with high-frequency instruments equipped by cruises or specially designed buoys [22]; (2) inertial dissipation methods [2325], based on the turbulent kinetic energy theory and budget equations, estimated surface turbulent fluxes by calculating the dissipation variables determined from vertical wind speed gradients and vertical temperature gradients, but they can only be applied to experiment sites and airborne data [24, 26, 27]; (3) data assimilation methods [2830] are designed based on the assumption of linearity or near-linearity, and most of them have been applied to numerical weather prediction (NWP) models; although they provide reasonable simulations of global coverage of surface fluxes, there are also notable problems caused by their prediction schemes [31]; and (4) bulk aerodynamic methods [27, 3234] generally related turbulent fluxes to satellite-derived meteorological variables by defining turbulent exchange coefficients. Currently, significant effort has been made to develop bulk aerodynamic methods that address the problem of absence of direct measurements. Although these methods are widely used to estimate ocean LHF, a limitation with both the eddy correlation methods and inertial dissipation methods is the essential requirement for high-frequency equipment. Meanwhile, bulk models have been successfully used to estimate ocean LHF by building relationships among meteorological variables. The Coupled Ocean-Atmosphere Response Experiment (COARE) version 3.0 [35] has been considered to be the most reliable bulk model [36, 37].

The bulk turbulent algorithm, COARE-3.0, can be successfully used to estimate ocean LHF with high accuracy. Currently, it has been applied to various existing satellite-based ocean LHF products, such as the South China Sea (SCS), HOAPS, J-OFURO, and objectively analyzed air-sea fluxes (OAFlux). Existing ocean LHF products based on the COARE 3.0 model have accurate estimates and high temporal resolution (e.g., daily) but rather coarse spatial resolution (e.g., >1°). Moreover, some COARE-based LHF estimations (e.g., SCS) are limited to regional coverage. For example, commonly used products such as SCS provide an accurate estimate of LHF, but its research area is limited to the South China Sea region. Therefore, as articulated by numerous groups within the global climate community, there is still a need for high-resolution, accurate ocean surface LHF estimation [38].

In this study, we used a bulk model (COARE 3.0) to estimate global ocean LHF with 4 km spatial resolution, driven by the Moderate-Resolution Imaging Spectroradiometer (MODIS) sea surface temperature (SST) product, Advanced Microwave Scanning Radiometer–EOS (AMSR-E) wind speed data, and reanalysis meteorological data. We had two major objectives. First, we estimated and mapped global ocean LHF using the COARE 3.0 model driven by the MODIS SST product and AMSR-E wind speed data. Second, we validated our estimated ocean LHF using buoy measurements from the moored buoy array of TAO, the Research Moored Array for African-Asian-Australian Monsoon Analysis and Prediction (RAMA), and the Prediction and Research Moored Array in the Tropical Atlantic (PIRATA) and then compared our estimated ocean LHF with the results of eight ocean LHF products (three reanalysis products, three satellite products, and two combined products).

2. Data and Methods

2.1. Data
2.1.1. Satellite and Reanalysis Inputs to the LHF Model

To estimate monthly ocean LHF, the forcing data (Table 1) for COARE 3.0 model include (1) satellite-based data (i.e., the MODIS SST product and AMSR-E wind speed data) and (2) reanalysis variables (Ta and q dataset). The monthly SST (SST4) product [40, 41] with 4 km spatial resolution was obtained from MODIS on board the Terra satellite provided by the National Oceanic and Atmospheric Administration (NOAA). The monthly wind speed data of the AMSR-E provided by the Japan Aerospace Exploration Agency (JAXA) [39] with a spatial resolution of 0.25° is available from May 2002 to October 2011. Input datasets for the COARE 3.0 model also included ERA-Interim [14] 2 m air temperature (Ta) data and 2 m air specific humidity (q) data with a spatial resolution of 0.125°, both of which are commonly used as accuracy datasets [50].

Considering the high correlation between SST and ocean LHF [51, 52], we decided to maintain the high spatial resolution of SST4 dataset to increase spatial heterogeneity. To generate the monthly global ocean LHF products from 2003 to 2007, we used the bilinear interpolation method to interpolate other monthly input data (U, Ta, q) to 4 km spatial resolution over 2003–2007.

2.1.2. Buoy Observations

The monthly ocean state variable observations from different moored buoy arrays were used as the reference data to evaluate the performance of the ocean LHF estimation (Figure 1). The data used in this study mainly include monthly ocean surface LHF. Ninety-seven moored buoys selected here can be divided into three groups due to various environmental conditions: 67 buoys were collected from the TAO/TRITON array (Tropical Atmosphere Ocean/Triangle Trans-Ocean Buoy Network), 18 buoys were collected from the PIRATA array, and 12 buoys were collected from the RAMA array. It is worth noting that the ocean LHF observations of the moored buoy sites used in this article were not direct observations but were computed by buoy-measurement meteorological variables using the COARE 3.0 bulk model. More details about the measurements and data postprocessing are given in https://www.pmel.noaa.gov/gtmba/.

The TAO/TRITON array [46, 47] located in the equatorial Pacific was supported primarily by NOAA of the United States (USA) and Japan. The PIRATA array [48] uses Autonomous Temperature Line Acquisition System (Atlas) moored buoys supported by France, Brazil, and the USA (NOAA), and the RAMA moored buoys array [49] located in the tropical Indian Ocean consists of 30 operational moored buoys, although only 12 buoys were available until 2007.

2.1.3. Comparison with Other Global Ocean LHF Products

To evaluate the performance of the COARE 3.0 model, seven global ocean LHF products were chosen as the comparison datasets (Table 1). These products include reanalysis datasets (MERRA-2, ERA-Interim, National Centers for Environmental Prediction (NCEP-2), and Japanese 25-year Reanalysis (JRA-25)), satellite-derived products (GSSTF version 3 and J-OUFRO version 2), and objectively analyzed air-sea fluxes product (OAFlux).

For comparison with the estimated product, we used the bilinear interpolation method to interpolate all LHF products from 2003 to 2007 to a spatial resolution of 4 km.

(1) MERRA-2. The assimilation system of Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) [42] has a key component that is the GEOS-5 (Goddard Earth Observing System) atmospheric model [53], which uses the finite-volume dynamical core of Lin [54]. To correct the background errors, the incremental analysis update (IAU) method [55] was applied to the analysis for adjusting the background state based on observations. It has global coverage from 1980 with a spatial resolution of 0.5° latitude × 0.625° longitude.

(2) ERA-I. ERA-Interim utilized an advanced atmospheric model with a spectral resolution of T255 on 60 vertical levels from the surface up to 0.1 hPa, and the data assimilation system uses 4D-Var analysis with a 12-hour analysis window to produce fields at 00, 06, 12, and 18 Z each day. ERA-I has a high spatial resolution at 0.125° and a high temporal resolution of hours.

(3) NCEP-2. NCEP-2 [29] is a reanalysis of NCEP, which is considered to be a revised version of the NCEP-NCAR (NCEP-1), with an approximate horizontal resolution of 1.9°. Assimilation was completed using spectral statistical interpolation (SSI) that is a 3D-Var assimilation system.

(4) JRA-25. JRA-25 [56] is a second generation reanalysis product produced by the Japan Meteorological Agency (JMA) and the Central Research Institute of Electric Power Industry (CRIEPI) covering the period from 1979 to 2004. The JRA-25 reanalysis applied a three-dimensional variational (3D-Var) data assimilation system and global spectral model with a spatial resolution of 1.125°.

(5) J-OFURO. J-OFURO is a third generation product produced by the School of Marine Science and Technology at Tokai University [44]. J-OFURO applied a wider variety of satellite measurements and has updated algorithms for retrieving variables. J-OFURO 3 has also adopted the COARE 3.0 method [35] for LHF estimation and has improved spatial resolution to 0.25°.

(6) GSSTF-3. GSSTF-3 dataset, produced by GES DISC, was released in October 2011 [43]. GSSTF-3 has further adopted an improved model [57] to retrieve air humidity directly related to the corrected Temperature Brightness (Tb), instead of using the two-step approach used in the previous product, which thus improved the accuracy of the ocean LHF estimation. GSSTF-3 used here has the advantage of a high horizontal resolution at 0.25°, and it is available from July 1987 to December 2009.

(7) OAFlux. The OAFlux version 3 [45] produced by WHOI improves the estimates of ocean heat fluxes with COARE 3.0 [35]. It uses an adopted objective analysis method to retrieve surface variables from an optimal blending of satellite-based estimations and three atmospheric reanalyses. The monthly OAFlux has a coarse spatial resolution at 1° from an optimal blending of satellite-based estimations and three atmospheric reanalyses.

2.2. Methods
2.2.1. COARE 3.0 Model

In the COARE 3.0 model [35], Monin–Obukhov [58] similarity theory (MOST) was applied to estimate ocean turbulent fluxes of Hl, sensible heat Hs, and stress τ, with standard bulk expressions as follows:where ρa is the density of air; cp is the specific heat of air at constant pressure; Le is the latent heat of vaporization; and CH, CE, and CD are the transfer coefficients for sensible heat, latent heat, and stress, respectively. The input data require mean quantities, including SST, near-surface potential temperature (θ), average value of wind speed (U) at reference height, and near-surface specific humidity (q).

As described in MOST, the transfer coefficients are partitioned into individual profile components, , , , , and the parameters (, , and ) are the bulk transfer coefficients for temperature, humidity, and wind speed, respectively. They have a dependence on surface stability ζ:where κ is the von Kármán constant at 0.4, the subscript n of cxn refers to neutral stability (where ζ = 0), z is the reference height of variable measurement, and zox (zoq, zot, zo) are the roughness lengths that characterize the transfer properties for variable x. ψx (ζ) is the MOST profile function that explains the stability dependence of the profile. The bulk stability parameter ζb has replaced the stability parameter ζ in the COARE 3.0 model, which reduces the iterations from 20 to 3 by improving the initial stability based on a Richardson number, Rib [59].

Considering the fact that some progress has been made in wave conditions related to bulk flux, an improvement of the COARE model version 3 is that it calculated velocity roughness length under specified wave conditions by implementing different parameterization schemes. Three parameterization schemes (YT96, TY01, and Oo02) have been applied in the COARE 3.0 model. The parameterization schemes of Yelland and Taylor [60] were also used in a previous version and have modified the Charnock parameter from a constant at 0.011 to a parameter increasing with the increment of wind speed. With the improvement in the Charnock parameter (α), the available range of roughness lengths for momentum has been extended to U below 20 m/s. For both TY01 [61] and Oo02 [62], the velocity roughness length was retrieved from wave properties to estimate ocean heat fluxes for diverse regions. Parameterization schemes have been described in more detail elsewhere [2, 35].

2.2.2. Evaluation Metrics

We used metrics of the squared correlation coefficient (R2), mean bias, and root-mean-square error (RMSE) to evaluate the performance of the COARE 3.0 model and different ocean LHF products. The matching degree between two datasets ({xi} {yi}) can be judged by the evaluation metrics of R2, bias, and RMSE, and they are calculated given as follows:

The King–Gupta efficiency (KGE) [63] is also considered to be a comprehensive evaluation metric, as follows:where R is a correlation coefficient between an observations dataset and an estimates dataset, SDe and SDo represent the standard deviation of estimates datasets and observations datasets, respectively, and Me and Mo are the mean values of the estimates dataset and observations dataset, respectively. KGE reaches a maximum of 1 without any simulation errors.

3. Results

3.1. Validation of the Estimated Global Ocean LHF with Buoy Observations

We used the COARE 3.0 model driven by the MODIS SST product, AMSR-E U data, ECWMF ERA-Interim Ta, and MERRA-2 q data to estimate global monthly ocean LHF at 4 km spatial resolution from 2003 to 2007.

The estimated LHF results using the COARE 3.0 model were directly compared with buoy observations at different buoy sites over 5 years. To reduce the uncertainties from the comparison between the estimation and observations, the average value of the observations for more than one buoy site within one pixel was considered as a reference.

Figure 2 shows scatter plots between the monthly observed ocean LHF and eight ocean LHF estimates for three buoy site arrays. For 67 buoy sites among the TAO array, our estimated LHF using the COARE 3.0 model performs the best at estimating monthly ocean LHF, with the highest KGE (0.83) and R2 (0.80, ), the lowest RMSE (16.0 W/m2), and lower bias (6.7 W/m2). Similarly, it has a good performance for the PIRATA buoy site array for ocean LHF, as indicated by a KGE = 0.89 and bias <6 W/m2. For the RAMA array, our estimated LHF performs unsatisfactorily compared to other buoy site arrays, with an R2 = 0.57 and RMSE = 19.4 W/m2; this result may have been caused by insufficient efficient observations.

Almost all estimated ocean LHF values from all products tend to present poor results for the RAMA buoy site array, with R2 values ranging from 0.25 to 0.57 and bias ranging from 19.4 to 66.8 W/m2. For the PIRATA buoy site array, however, almost all ocean LHF estimations showed good agreement with observations, as indicated by 0.03–0.2 higher R2 values. It is clear that the ocean LHF estimations for PIRATA were superior to that of TAO, as indicated by R2 being 0.03–0.21 higher. Among all ocean LHF products, the J-OFURO product is also produced based on the COARE 3.0 model. Thus, secondary to our estimated LHF, J-OFURO also performs unsatisfactorily for RAMA buoy site array, but it still showed better performance than that of other ocean LHF products and exhibited a higher KGE ranging from 0.50 to 0.83. It is worthwhile to note that the estimated ocean LHF from the ERA-I product outperforms other reanalyses (MERRA-2, JRA-25, and NCEP-2) at most buoy site arrays, as indicated by KGE being 0.13–0.24 higher, R2 being 0.02–0.18 () higher, RMSE being 6.58–26.04 lower, and bias being 4.57–25.05 lower. Almost all products overestimate ocean LHF for all buoy site arrays, especially for the JRA-25 reanalysis, which has biases that are more than 44 W/m2 higher than all buoy site arrays. Overall, our estimated LHF yields the best ocean LHF with the highest R2 and KGE and lower RMSE and bias in comparison to other ocean LHF products.

Figure 3 shows the comparison of the estimated annual ocean LHF and observed ocean LHF at all 96 buoy sites, and the results illustrated that our estimated LHF derived from the COARE 3.0 model has a good ability to accurately estimate ocean seasonal and annual LHF. Overall, the KGE values for our estimated seasonal and annual LHF versus observations were approximately 0.80 and 0.84, the RMSEs were 16.2 W/m2 and 9.9 W/m2, and the bias values were 6.6 W/m2 and 5.0 W/m2, respectively. Similarly, we can also note the ability of our estimated ocean LHF to estimate the among-site variability, where the R2 of site-averaged estimates compared to the observed LHF is approximately 0.87, which is slightly higher than that of J-OFURO (0.84). Overall, our estimated LHF driven by the COARE 3.0 model has a high accuracy according to the validation of temporal and spatial variation in ocean LHF estimates.

Compared to estimated LHF values derived from other products, our estimated ocean LHF is also satisfactory for reproducing interannual variability at the site scale with at least five years of data. Moreover, the bias between our estimated LHF anomaly and the buoy site observations anomaly is significantly lower than those of other LHF products, with values of 0.4 W/m2, and has the highest R2 of 0.61. For all buoy site observations, although the medians of the KGE for eight ocean LHF products were close (Figure 4), JRA-25, GSSTF-3, and NCEP-2 showed large differences among sites, implying inconsistent model performance for buoy site observations. In contrast, our estimated monthly LHF outperforms other ocean LHF products at all buoy sites, as indicated by an average KGE of 0.84, followed by J-OFURO, ERA-I, and OAFlux with KGEs of 0.82, 0.74, and 0.73, respectively.

3.2. Mapping of the Estimated Global Ocean LHF
3.2.1. Seasonal LHF

Figure 5 shows the seasonal patterns of ocean LHF distribution averaged over years (2003–2007), where positive values indicate heat loss from ocean to atmosphere and negative values represent heat loss from atmosphere to ocean. The global average LHF for MAM (March, April, and May), JJA (June, July, and August), SON (September, October, and November), and DJF (December, January, and February) is 220.2 W/m2, 206.5 W/m2, 236.3 W/m2, and 252.8 W/m2, respectively, indicating that maximum heat loss occurs in DJF and that minimum heat loss occurs in JJA. Generally, in both SON and DJF, ocean LHF showed maximum heat loss. In addition, the peak value of DJF ocean LHF appeared in the Northern Hemisphere. However, when JJA comes, the peak value decreased and moved to the Southern Hemisphere (Figure 5). As a result, the peak value of the Northern Hemisphere has disappeared.

For DJF in the Northern Hemisphere, the maximum heat loss occurs in the western boundary current regions (i.e., the Kuroshio and Gulf Stream and its extension). In the western boundary current regions of the Northern Hemisphere, the seasonal maximum heat loss in DJF exceeds 350 W/m2, while it falls to below 150 W/m2 in JJA. Similarly, the large seasonal change of LHF in the Southern Hemisphere occurs in the boundary current regions, especially the western boundary current regions. One notices that the LHF difference between JJA and DJF in the Southern Hemisphere is much less than that in the Northern Hemisphere. This outcome partly reflects the differences in the distribution of land and sea in the Southern and Northern Hemispheres, allowing the air masses to move southward and reduce the temperature difference and air-sea humidity difference [13].

3.2.2. Annual LHF

We mapped annual global LHF averaged over the period of 2003 to 2007 from our estimated LHF, OAFlux, J-OFURO, GSSTF-3, MERRA-2, ERA-I, NCEP-2, and JRA-25 (Figure 6). They have a similar global distribution of ocean LHF, though general differences exist in the spatial LHF distributions due to the difference between models and their spatial resolution. All the LHF estimates yield higher annual LHF in western Pacific area, and Eastern Pacific Cold Tongue region has lower ocean LHF owing to the decreased SST and moisture limitations. The most striking differences among ocean LHF estimations are the heat loss occurring in low-latitude regions, especially in the Northern Hemisphere. In contrast, the reanalysis (MERRA-2, ERA-I, NCEP-2, and JRA-25) yields higher average annual ocean LHF than that of other ocean LHF products in low-latitude regions, as indicated by the bias being more than 30 W/m2.

3.3. Comparison with Other Global LHF Products
3.3.1. Spatial Differences among Ocean LHF Products

Figure 7 is a scatter plot between annual ocean LHF products and our estimated LHF for all pixels averaged from 2003–2007, and it shows that seven ocean LHF products are highly correlated to our estimated ocean LHF, as indicated by an R2 of more than 0.89 (). It is clear that the LHF results derived from ERA-I and J-OFURO products are similar with our estimated LHF, with an R2 of more than 0.92, a KGE higher than 0.81, and bias lower than 14 W/m2. The OAFlux product also provides a similar estimation to our estimated LHF, as indicated by the highest KGE of 0.89 and the lowest bias of less than 6 W/m2. Among these ocean LHF products, the JRA-25 and NCEP-2 have the greatest differences from our estimated LHF, with KGEs less than 0.65 and biases more than 24 W/m2. We note that the strip distribution of pixels from NCEP-2 and JRA-25 is mainly caused by the coarse spatial resolution.

Figure 8 shows the multiyear (2003–2007) average spatial distribution of the differences between our estimated LHF and other LHF products. The result illustrated that the spatial distribution of our LHF estimates is similar to those of most ocean LHF products, such as ERA-I, MERRA-2, OAFlux, and J-OFURO, especially the satellite-based product J-OFURO, which is also estimated using the COARE 3.0 model. In tropical regions, our estimated LHF is slightly higher than OAFlux and similar to ERA-I and J-OFURO; however, in mid-latitude and high-latitude regions, our estimated ocean LHF is slightly lower than OAFlux, ERA-I, and J-OFURO product.

The most striking differences between our estimated LHF and the satellite-based GSSTF-3 product are the positive differences (more than 30 W/m2) that occur in the continental boundary region and high-latitude areas, which may be attributed to large uncertainties from forcing data. Both NCEP-2 and JRA-25 reanalyses show a significant increase in LHF in most areas, with average differences of 24.2 W/m2 and 26.5 W/m2, which are probably due to the combination effect of bulk variables, such as the low humidity difference, the weak wind speed, and the low SST. Overall, our LHF estimates have a similar spatial distribution to those of other products, but there are substantial differences in some areas due to the differences from different models and forcing data.

3.3.2. Seasonal and Zonal Averages

Figure 9 compares the estimated monthly global average LHF with other ocean LHF products. They showed strong seasonality and similar temporal variations among all ocean LHF estimates. Most ocean LHF estimates decrease from January to July and increase afterwards, with maximum values occurring in DJF and minimum values occurring in JJA. Meanwhile, there are still differences in the monthly LHF among all products; the multiyear monthly average LHF of GSSTF-3, JRA-25, and NCEP-2 is higher than those of other LHF estimates. Among them, the overestimation of LHF in GSSTF product is caused by the lack of effective value in high-latitude area, which is the location of low ocean LHF, resulting in the overestimation of average LHF. Besides, it is found that the temporal pattern of ERA-I is similar to reanalysis, like JRA-25 and MERRA. Meanwhile, the temporal pattern of our estimated LHF is similar to satellite-based products (J-OFURO and GSSTF-3) and OAFlux product, which indicated that the LHF difference between JJA and DJF in reanalysis is much less than that in satellite-based products.

Although all eight products showed similar latitudinal variability of ocean LHF, there are still substantial differences in eight ocean LHF estimates. The latitudinal distribution pattern of latent heat is bimodal with the maximum in the subtropical region and the minimum at the pole, suggesting a decrease in LHF from subtropical to high latitudes (Figure 10). The highest annual ocean LHF occurs in the subtropical trade wind zone due to high winds, leading to incremental changes in LHF. Comparing all LHF products, we found that the peak values of LHF products varied from 130 W/m2 to 175 W/m2, while the products of NCEP-2 and JRA-25 simulated the highest LHF values, especially in low-latitude areas. Our estimated LHF is slightly lower than J-OFURO, and the pattern of our estimated LHF is similar to the OAFlux product, with a difference of less than 10 W/m2.

4. Discussion

4.1. Performance of the COARE 3.0 Model in Estimating Global Ocean LHF

Validation for 96 globally distributed buoy sites (TAO, PIRATA, and RAMA) for the period of 2003 to 2007 illustrated that the bulk aerodynamic model COARE 3.0 is reliable for estimating ocean latent heat LHF and is robust for low-latitude regions. Comparing buoy observations reveals that biases do exist among different LHF estimates, whereas Figures 2 and 3 both demonstrate that the ocean surface LHF estimates derived from the COARE 3.0 model (OAFlux, J-OFURO, and our estimated LHF) have no significant LHF bias and yield LHF values close to buoy site observations relative to other LHF estimates (such as GSSTF-3, JRA-25, and NCEP-2). For example, compared to the average statistical metrics of five other ocean LHF products, the COARE 3.0 model shows better performance for all buoy site observations, with KGE being 11%–30% higher and RMSE being 30%–35% lower. A number of studies have shown that the COARE bulk method has improved the performance of ocean LHF estimation [64, 65]. Brunke et al. [36] evaluated and ranked 12 bulk aerodynamic methods and demonstrated that COARE 3.0 has the least problematic algorithms for estimating ocean latent heat flux. A similar conclusion has been proposed by Iwasaki et al. [37], who evaluated four bulk methods (COARE 3.0, Chou, Kondo, and Zeng) based on observed LHF during 15 cruises obtained using direct eddy correlation and inertial dissipation flux methods. The result illustrated that the bias of COARE 3.0nowc (COARE 3.0 method without warm layer and cool skin models) is lower than 10 W/m2, while those of other methods are higher than 10 W/m2 in all wind speed regions, and thus they consider COARE 3.0 to be the best bulk algorithm for calculating ocean LHF. Error from bulk methods can lead to 30–50% of total error of ocean LHF; thus, using an applied robust COARE 3.0 algorithm to estimate ocean LHF contributed to the high accuracy of our estimated LHF in this study.

We also noticed that most products overestimate ocean LHF according to validations at the buoy site scale, but our estimated LHF in this study has reduced biases by 4–25 W/m2. In summary, these comparison results (Figures 3, 6, 9, and 10) provide confidence regarding COARE 3.0-derived estimates for mapping of ocean LHF, and thus we concluded that the COARE 3.0 model is a better choice for the calculation of LHF.

However, even if the model is less problematic on the whole, it is highly possible that the COARE 3.0 model does not perform well for some regions. For instance, although the OAFlux product has lower RMSE (18.9 W/m2, 38.0 W/m2) and bias (5.9 W/m2, 27.7 W/m2) than ERA-I at some buoy arrays (PIRATA and RAMA), ERA-I has a better performance with higher KGEs. As illustrated by Fairall et al. [35], the bulk COARE 3.0 model is applicable to wind speeds below 20 m/s. Andreas et al. [66] documented the same conclusion; that is, the bulk method could successfully estimate ocean LHF in moderate wind (U < 20 m/s). Furthermore, Andreas et al. [66] also reported that the method is not suitable in high wind speeds due to the nonlinear relationship between high wind speed and ocean turbulent flux. It is necessary to improve the parameterization schemes for unstable conditions included in COARE 3.0.

4.2. Global Ocean LHF Estimation

Another goal of this study is to assess the quality of our estimated ocean surface LHF, which is derived from a bulk aerodynamic method utilizing both reanalyzed meteorological variables and satellite-derived parameters. Due to the lack of spatial continuous observations for ocean latent heat flux, it is difficult to apply the rigorous validation of ocean LHF. Therefore, we demonstrated the reliability of our estimated ocean LHF values outside of the tropics by comparing them with those of other LHF products (Figures 7, 9, and 10).

The COARE 3.0-based estimate of annual average global (90S°-90N°) LHF is 69.8 W/m2 from 2003 through 2007; this value is comparable to other LHF estimations. Chou et al. [67] compared four annual products and reported that global (60S°-60N°) average LHF varied from 88.5 W/m2 derived from the HOAPS product to 108.2 W/m2 from the GSSTF-3 product during 2002 to 2003, while the COARE 3.0-derived annual ocean LHF (60S°-60N°) was 92.8 W/m2 from 2003 to 2007. L Yu and Weller [13] inferred that the global (60S°-60N°) average estimates of ocean LHF from OAFlux product ranged from 86 W/m2 to 95 W/m2 during the period of 1982 to 2004. Despite general differences among different products, Figure 6 shows the common characteristics of all LHF estimations; that is, the annual average maximum ocean LHF occurs in the western ocean due to increased SST, and lower ocean LHF occurs in the Eastern Pacific Cold Tongue region owing to the lower sea surface temperature and lower air-sea humidity difference. We then concluded that our estimation successfully captures the spatial and temporal information of ocean LHF (Figures 9 and 10) with an average difference of less than 10 W/m2 compared to most products.

Compared with other LHF products, our estimated LHF has the better performance for the tropics, according to observation validations. In addition to the method, we also attributed the better performance of our estimated LHF to high-resolution SST input data. The most important strength of SST4 data is that high spatial resolution can increase the spatial heterogeneity [40]. To evaluate the effect of SST4 data on the accuracy of ocean LHF estimation, we recalculated OAFlux and J-OFURO products based on the same COARE 3.0 model driven by high-resolution SST4 data to replace the SST data used in LHF products. As shown in Figure 11, the new products (OAFlux_new and JOFURO_new) perform well in estimating monthly ocean LHF, with an R2 increase of approximately 0.04, a bias reduced by 2–5 W/m2, and KGE values increased by approximately 0.04. Compared to previous products for all buoy sites, JOFURO_new performs well with higher KGE at 0.84, lower average bias at 8.0 W/m2, and higher R2 at 0.78. Therefore, SST4 data improves the performance of observations validation compared to the other coarser LHF products at the buoy site scale.

4.3. Uncertainties in Ocean LHF Estimation

Compared with buoy site observations, validation results illustrated that the uncertainty in monthly LHF estimates based on the COARE 3.0 model varies from 18% to 35%. We have attributed the bias of ocean LHF estimation to the uncertainties from buoy site observations, errors in forcing bulk data (e.g., reanalysis variables and satellite-derived data), algorithm limitations, and spatial scale mismatches among different data sources.

First, the buoy site observations were used as reference data to determine the accuracy of ocean LHF. However, the errors of observations in buoy sites cannot be ignored. The uncertainty in the buoy estimate of LHF is approximately 10 W/m2 [68], which would have an impact on the accuracy of the estimated LHF. Another noticeable problem is the method used in calculating ocean surface LHF for buoy site observations. The ocean surface turbulent flux is calculated using the COARE 3.0 bulk model rather than the eddy correlation (EC) method. In fact, no EC systems were installed at the buoy sites; consequently, most buoy site observations were applied to the COARE 3.0 bulk model to estimate ocean surface LHF using observed meteorological variables. This may be a reason for the high correlation of observed data to the LHF estimates based on the COARE 3.0 model (in our estimated LHF, OAFlux, and J-OFURO), by eliminating model errors between LHF products.

Second, biases of forcing data for the COARE 3.0 model are another factor contributing to the uncertainties for the global ocean LHF estimation. For all forcing inputs, the correlation between ocean LHF and the humidity difference in the tropics can reach 0.87, and for the temperature difference the correlation can reach 0.74 [69]. Therefore, the accuracy of bulk variable data impacts the estimation of ocean LHF. Many studies have illustrated that there are large errors in bulk variables; for instance, ERA-I data tends to underestimate air temperature (Ta) when compared to observations from Kwajalein Experiment (KWAJEX) [70]. In addition, in our study, it was found that continuous variations of wind speed occur in high-latitude areas for individual months, resulting in anomalous LHF in the same area. This indicates that biases in wind speed data can introduce substantial uncertainties into LHF estimates. These results suggest that it is necessary to minimize variable-caused biases in the retrieval to improve the accuracy of the input data and thus improve the estimation of ocean LHF, for example, through the increment analysis update (IAU) technique used in MERRA-2 [70].

Third, spatial scale mismatch among the different datasets may also lead to uncertainty of ocean LHF estimates. The horizontal resolution of reanalysis and wind speed data was greater than 12.5 km, which was much greater than the SST data with the spatial resolution at 4 km. In addition, the pixel average for LHF products is larger than 0.25° (e.g., MERRA-2), and pixel size in some products reaches approximately 2° (e.g., in NCEP-2), whereas buoy site observations can only represent several hundred meters. The spatial scale mismatch may lead to uncertainties in the validations.

Finally, the structure of the COARE 3.0 model will also reduce the accuracy of ocean LHF estimates due to the different parameterization schemes. The parameterization of important parameters (i.e., CH, CD, and roughness length) is the key to reducing algorithm uncertainty [71, 72]. Specifically, variations in the exchange coefficients have substantial differences between moderate wind and high wind [69, 73]. Additionally, the parameterization schemes of exchange coefficients and roughness length are also influenced by wave state. It has been a popular topic to develop more applicable parameterization schemes for exchange coefficients; these schemes could improve the accuracy of ocean LHF estimates [71, 72]. Zhang et al. [69] demonstrated that the neutral drag coefficients derived from the Charnock parameter become larger [63, 74] at high wind speeds (U > 30 m/s); therefore, they have combined a new parameterization scheme (DREMAKIN) of roughness length with the COARE 3.0 model to suit various wind speed conditions. Pan et al. [72] developed a regression model of roughness length (PS07) that is based on the highly sensitive relationship between wave age and dimensionless roughness length. In comparison with other parameterization schemes consistent with the COARE 3.0 model (TY01, YT96, and Oo02), they found that the PS07 model had the highest accuracy in the North Sea Platform experiment and the Zhanjiang Seaport experiment. Evaluation of LHF product performance at the site scale in this study was only performed in the tropics due to the lack of available and time-continuous buoy observations over a wide area outside the tropics. To analyze the performance of estimated LHF globally, it is urgent to develop more and longer observations distributed in wide regions.

5. Conclusions

We applied the bulk aerodynamic COARE 3.0 model driven by many bulk variables (i.e., MODIS SST with 4 km spatial resolution, AMSR-E wind speed, specific humidity, and air temperature obtained from ERA-Interim) to estimate monthly ocean LHF from 2003 to 2007. To evaluate the performance of the COARE 3.0 model, we validated our estimated global ocean LHF using buoy data collected from 96 buoy sites from 3 buoy site arrays (TAO, PIRATA, and RAMA). Additionally, we also compared our estimated global ocean LHF values with seven ocean LHF products: 4 reanalysis products (MERRA-2, ERA-I, JRA-25, and NCEP-2), 2 satellite-derived products (GSSTF-3 and J-OFURO), and a combined product (OAFlux).

Compared to the observations from buoy sites, the results for ocean LHF estimation driven by reanalysis variables and satellite-derived data showed that our LHF estimates yield better performance than the other seven ocean LHF products mentioned above. At the buoy site scale, the validation results show that monthly ocean LHF estimates based on the COARE 3.0 method performed better, with the highest R2 values, lower bias and RMSE values, and the highest KGE values for all buoy site arrays. Compared to their monthly products, all evaluated annual ocean LHF products are closer to all 96 buoy sites, and the R2 of the results from NCEP-2 and JRA-25 increased by 20% and 15%, respectively. Moreover, the mean annual COARE 3.0-based global (60S°-60N°) ocean LHF was 92.8 W/m2 from 2003 to 2007, which is in good agreement with other studies.

These eight products showed strong and similar seasonality and have similar temporal variation of LHF. Specifically, high heat loss occurs in DJF, while lower ocean LHF occurs in JJA. The spatial distribution of the average annual ocean LHF values shows high values in the western boundary current regions and the subtropics, while low average annual LHF values are observed in the high-latitude areas. Our estimated ocean LHF performs well with high accuracy; we attributed this to the advantaged model COARE 3.0 and the high spatial resolution input data which would improve the performance of ocean LHF estimation. Future work could focus on generating long-term LHF estimations using SST4 data and reanalysis wind speed data.

Data Availability

The data used to support the findings of this study are included within the article and listed in the reference list. Buoy site observations (TAO, PIRATA, and RAMA) were obtained from Pacific Marine Environmental Laboratory (PMEL) of the National Oceanic and Atmospheric Administration (https://www.pmel.noaa.gov/gtmba/). MERRA-2 and GSSTF-3 were obtained from the Earthdata website (https://earthdata.nasa.gov/). ERA-I was obtained from public dataset of ECMWF website (https://www.ecmwf.int/). AMSR-E, JRA-25, J-OFURO, NCEP-2, and OAFlux were obtained from the Asia-Pacific Data-Research Center (APDRC) (http://apdrc.soest.hawaii.edu/index.php). SST4 was obtained from NASA’s Ocean Color Web, which is supported by the Ocean Biology Processing Group (OBPG) at NASA’s Goddard Space Flight Center (https://oceancolor.gsfc.nasa.gov/).3

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors thank Prof. Shaomin Liu, Dr. Tongren Xu, Dr. Zhongli Zhu, and Dr. Linna Chai from the Faculty of Geographical Science, Beijing Normal University, China, and Prof. Xin Wang and Dr. Rongwang Zhang from the State Key Laboratory of Tropical Oceanography, South China Sea Institute of Oceanology, Chinese Academy of Sciences, China, for their suggestions to improve this manuscript. The authors also thank the International Science Editing group (http://www.international.scienceediting.com) for editing this manuscript. This work was partially supported by the National Key Research Development Program of China (2016YFA0600102) and the National Natural Science Fund of China (41671331).