Evaluation of Different Heat Flux Products Over the Tropical Indian Ocean

Net heat flux (Qnet) and its components from four reanalysis (NCEP‐2, CFSR, ERA5, and MERRA) and two blended products (OAFlux & TropFlux) are compared with in situ observation (two Research Moored Array for African‐Asian‐Australian Monsoon Analysis and Prediction buoys and one Woods Hole Oceanographic Institution buoy) over the north Indian Ocean to quantify their uncertainties in daily, seasonal, and annual scales. These comparisons provide the present status of Qnet error in most state of the art reanalysis/blended products. The root‐mean‐square error (RMSE) remains similar to the RMSE a decade earlier, despite more observation and improved models and reanalysis methods. However, there is a clear separation of flux quality from the older generation of reanalysis (NCEP‐2) to the newer production of reanalysis (MERRA, CFSR, and ERA5). While individually ERA5 provides the best estimate, the ensemble mean (i.e., average of ERA5, CFSR, MERRA, TropFlux, and OAFlux) is very close to ERA5 both in terms of correlations and RMSE and provides the most reliable estimate by virtue of removal of some of the uncertainties in estimation of flux by each of the flux products. A significant reduction of RMSE in Qnet estimates from 100 W/m2 (in NCEP‐2) to 45 W/m2 (in the ensemble mean) is considerable progress. It is noteworthy that all the recent flux products estimate the increasing trend of Qnet in the north Bay of Bengal and subseasonal fluctuations with significant fidelity. Also, in the south of equator location vigorous subseasonal fluctuations in boreal winter are well captured. We believe that this is significant progress in the estimation of Qnet over the Indian Ocean.


Introduction
Accuracy in estimation of the air-sea fluxes over the ocean surface is essential as they represent a medium of communication between the ocean and atmosphere, leading to air-sea interaction processes on diurnal to interdecadal variability. Surface net heat flux (Qnet) plays a crucial role in controlling feedback between the ocean and atmosphere . Accurate turbulent flux estimates are essential to assess the energy budget (Dong & Kelly, 2004) as it modulates the mixed layer temperature (Niiler & Kraus, 1977), sea surface temperature (Yu et al., 2006, heat content (Garternicht & Schott, 1997;Godfrey, 1995), and ocean circulation (Schott et al., 2009). Typically, changes in the Ocean Heat Content of the upper ocean layers can be quantified through the estimation of imbalance of surface flux components. Thus, accurate surface net heat flux through the heat budget in the upper ocean can give insight into the roles of the atmosphere and ocean in the evolution of sea surface temperature (SST), heat content (Dong & Kelly, 2004), and ocean-atmosphere feedback with changing climate. Therefore, the correct representation of surface fluxes is essential for improving climate models and their forecast skills.
Qnet is a combination of radiative fluxes (shortwave and longwave radiation) and turbulent heat fluxes (latent and sensible flux). So any uncertainties in the estimation of net heat flux are contributed by the uncertainties in each of the components. In situ flux measurements acquired from buoys and ships have long been used as a reference base to quantify the accuracy of surface heat flux products from reanalysis and satellite-based products (e.g., Brunke et al., 2011;Cronin et al., 2006;Josey, 2001;Pinker et al., 2009;Smith et al., 2001;Yu et al., 2004a;Yu et al., 2006). Previous studies have reported significant uncertainty in the estimation of Qnet products over the global oceans (e.g., Large & Yeager, 2009;Schott et al., 2009;Yu et al., 2006). These uncertainties are mainly due to the errors in the input parameters (e.g., Large & Yeager, 2009;Rahaman & Ravichandran, 2013) and the bulk algorithms used for estimation of flux component (Brunke et al., 2003;Kalnay et al., 1996;Uppala et al., 2005). Several studies addressed the uncertainties in the estimation of Qnet in various flux products. For example, a recent paper by Valdivieso et al. (2017) reviewed the surface heat fluxes from an ocean reanalysis perspective over the globe using 16 monthly flux products. Bentamy et al. (2017) accessed various turbulent heat flux derived from satellite-based and atmospheric reanalysis products in spatial and temporal scales and found that all of them exhibit similar space and time patterns with significant differences in magnitude in some specific regions. More recently, Yu (2019) addresses the dominant source of uncertainties of surface flux products and the reliability of using these products for the budget closures.
However, this kind of study of systematic evaluation is lacking over the Indian Ocean. It has been reported in many studies that Indian Ocean SST is crucial for ISMR prediction and is mostly governed by Qnet on seasonal and intraseasonal time scales (e.g., Parampil et al., 2016;Sengupta & Ravichandran, 2001). Hence, in this study, we focus on quantifying the uncertainty (error) in estimating the Qnet over the Indian Ocean (IO) by the latest hybrid flux products and latest reanalysis flux products to ascertain the progress made on the unacceptably large errors known to exist in earlier products . This is important as oceanatmosphere interactions are critical in maintaining the annual cycle of SST and precipitation (Webster et al., 1998), interannual variability (Saji et al., 1999), and the monsoon intraseasonal oscillations (Sengupta & Ravichandran, 2001). Further, uncertainties of the order of 10 Wm −2 in Qnet estimate may result in uncertainty of 0.5°C in SST over three months, implying a strong need for reduction of the uncertainties in flux to less than 10 Wm −2 (Roberts, 2011). Besides, seasonally reversing differences in Qnet over the northern and southern tropical IO drive a shallow meridional circulation (Schott & Mccreary, 2001) in the tropical IO. While models do simulate the shallow meridional circulation and its variability, it is poorly constrained due to the absence of adequate current observations. For confidence in the simulated shallow meridional circulation by the models, their ability to reproduce Qnet correctly is critical. Hence, a comparison of model-simulated Qnet with "observation" and reliable estimate of Qnet in flux products over the IO region is essential. Over the Indian Ocean basin, since the study by  who showed that the uncertainties in net heat flux could be as significant as 60-100 W/m 2 over the region, no comprehensive evaluation studies have been carried out over the Indian Ocean afterward. In the backdrop of availability of a new hybrid flux product, namely, the TropFlux and latest reanalysis products like CSFR and ERA5 together with some new in situ observations from Research Moored Array for African-Asian-Australian Monsoon Analysis and Prediction (RAMA) moorings and the Woods Hole Oceanographic Institution (WHOI) buoy (see section 2), we ask: is the uncertainty in Qnet estimation reduced in new individual flux products? What is the remaining interproduct uncertainty? This study is motivated to answer these questions.
In this study, we have used six widely used flux products in their common time frame (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009). Four flux products are based on reanalysis, and the rest two are blended/hybrid products. The reanalysis consists of NCEP-2 (Kanamitsu et al., 2002), Modern-Era Retrospective analysis for Research Applications (MERRA; Rienecker et al., 2011), NCEP CFSR (Saha et al., 2010), and ECMWF ERA5 (Olauson, 2018) together with two hybrid products are OAFlux  and TropFlux (Kumar et al., 2012). The MERRA, CFSR, and ERA5 are more recent reanalysis data sets employing improved model and analysis systems, and NCEP-2 is a relatively older reanalysis. Inclusion of NCEP-2 allows us to quantify the level of improvements achieved in the recent reanalysis. Both the hybrid flux products (OAFlux and TropFlux) use International Satellite Cloud Climatology Project (Zhang, 2004) for radiation component (shortwave and longwave radiations) and surface meteorological fields from a combination of buoy, ship, reanalysis, and satellite data for estimating turbulent flux component (latent and sensible). Although previous studies have shown the superiority of hybrid flux products as compared to the reanalysis flux products over the globe (Yu et al., 2006Farmer, 2012;Liu et al., 2015) for the estimation of net heat flux, still the uncertainties are quite significant. In recent times, when all the flux products have considerably improved due to ((1)) enhanced sampling from moored buoys and another in situ platform, (2) significant increase in the coverage of satellite data from multisatellite observations, (3) improved retrieval algorithms for satellite humidity and air temperature with much-reduced uncertainties, and (4) better assimilation methods combined with availability of large number of quality data for the reanalysis model. Therefore, it is necessary to reassess the uncertainties in the latest flux data and assess their consistencies using the most recent in situ data. The paper is organized as follows. Section 2 describes the data sets used in the study. The annual and seasonal mean features of the different flux products are discussed in section 3. The contribution of heat flux and ocean processes toward SST variability is elaborated in section 4, and a summary of the study is given in section 5.

Data and Methods
The present study utilizes six flux data from three different sources. Two are hybrid flux products (OAFlux and TropFlux), four are from reanalysis (NCEP-2, MERRA, CFSR, and ERA5), and remaining two flux data are from in situ moorings (RAMA and WHOI buoys) which are described in detail below.

Blended Flux Products
The OAFlux products ( Yu et al., 2004a( Yu et al., , 2004b use the COARE bulk algorithm 3.0  to give the daily estimates of surface latent heat and sensible heat fluxes over the global ice-free ocean. The input parameter required for the bulk algorithm is calculated from the combination of International Satellite Cloud Climatology Project (ISCCP) for radiative downward fluxes (shortwave and infrared fluxes) and surface meteorological state variables (10-m winds, 2-m air and sea temperature, 2-m air relative humidity) are objectively synthesized using buoy, ship, satellite, and reanalysis data. The ISCCP flux data (ISCCP-FD) are calculated at the surface using a complete radiative transfer model from the GISS GCM with improved observations of the physical properties of the surface, atmosphere, and clouds based on the ISCCP data sets. The ISCCP-FD data are available on 3-hr interval over the globe covering the period July 1983 through December 2009. The data were daily averaged and linearly interpolated onto a 1°grid to combine with OAFlux latent and sensible heat flux components. The resulting net heat flux product covers the period from July 1983 to December 2009. We have used a monthly product of net shortwave radiation flux (NSWR), net longwave radiation flux (NLWR), latent heat net flux (LHF), sensible heat net flux (SHF), and Qnet for our study.
The TropFlux product (Kumar et al., 2012) has been developed under collaboration between the Laboratoire d'Océanographie: Expérimentation et Approches Numériques from Institut Pierre Simon Laplace, France and National Institute of Oceanography/CSIR, India and distributed through Indian National Center for Ocean Information Services, India. The TropFlux heat fluxes are also estimated using the COARE 3.0 algorithm. TropFlux uses bias and amplitude corrected ERA-I (10-m winds, 2-m air and sea temperature, 2-m air relative humidity, and downward radiative fluxes) and ISCCP (shortwave radiation) fluxes. All bias corrections are derived based on comparisons with the Global Tropical Moored Buoy Array data (Kumar et al., 2012). At present TropFlux provides data of 1°× 1°spatial resolution for the entire 30°N-30°S region. Daily and monthly data of NSWR, NLWR, LHF, SHF, and Qnet variables at the surface are used during the period of 2000 to 2009.

Reanalysis Flux Products
NCEP-2 (NCEP-NCAR reanalysis; Kanamitsu et al., 2002) is an update to NCEP-1 (Kalnay et al., 1996) that rectified errors and updated parameterizations of physical processes. It is generated from a frozen forecastanalysis platform that consists of a T62 model (equivalent to a horizontal resolution of about 210 km) with 28 vertical levels and 6-hr intervals. It is using a state-of-the-art analysis/forecast system to perform data assimilation using past data from 1979 through the previous year, currently available from January 1979 to August 2019. Components of radiative heat flux on monthly scale for DSWR (downward shortwave radiation flux), USWR (upward shortwave radiation flux), DLWR (downward longwave radiation flux), ULWR (upward longwave radiation flux), and turbulent heat flux (LHF and SHF) are directly obtained from the APDRC site (http://apdrc.soest.hawaii.edu/datadoc/ncep_mon.php). Qnet from NCEP-2 is calculated as where NSWR = DSWR − USWR and NLWR = ULWR − DLWR.
The MERRA data set (Rienecker et al., 2011) was released in 2009. It is based on a version of the Goddard Earth Observing System Data Assimilation System version 5 atmospheric data assimilation system that was frozen in 2008. MERRA data spanned the period 1979 through February 2016 and were produced on a 0.5°× 0.66°grid with 72 layers. MERRA was used to drive stand-alone reanalysis of the land surface (MERRA-Land) and atmospheric aerosols (MERRAero). The shortwave radiation scheme of Chou and Suarez (1999) and the long-wave radiation scheme of Chou et al. (2001) are utilized in MERRA products. Roberts et al. (2011) and Brunke et al. (2011) analyze surface turbulent fluxes over the ocean from MERRA and other data products. Daily and monthly data of NSWR, NLWR, LHF, and SHF are downloaded (http://apdrc.soest.hawaii.edu/datadoc/merra.php), and the Qnet is calculated using equation (1).
The NCEP CFSR (Saha et al., 2010) is the coupled assimilation reanalysis which uses the NCEP coupled forecast system (CFSv2) model. The atmospheric model has a spectral resolution of T382 (38 km) and 64 hybrid vertical levels. The ocean model consists of the Geophysical Fluid Dynamics Laboratory modular ocean model (version 4p0d; Griffies et al., 2004), which is a finite-difference model at a resolution of 0.5°C with 40 levels in the vertical. The atmosphere and ocean models are coupled with no flux adjustment. CFSR uses simplified Arakawa-Schubert convection with momentum mixing. CFSR implements orographic gravity wave drag based on the approach of Kim and Arakawa (1995) and subgrid-scale mountain blocking by Lott and Miller (1997). It uses rapid radiative transfer model shortwave radiation with maximum random cloud overlap (Clough et al., 2005;Iacono et al., 2000). In CFSR implementation, both shortwave and longwave radiations are invoked at 1-hr interval. It is also coupled to a four-layer Noah land surface model (Ek et al., 2003) and a two-layer sea ice model (Wu et al., 2005). The NCEP CFSR uses the gridded statistical interpolation data assimilation system for the atmosphere. Flow dependence for the background error variances is included as well as first-order time interpolation to the observation (Rancic et al., 2008). Variational quality control of observations (Anderson & Järvinen, 1999) is also included. An ocean analysis for SST is also performed using optimum interpolation. A full range of observations is used as in the other reanalyses, which are quality-controlled and bias-corrected, including satellite radiances. Observations on ocean temperature and salinity are also used. Further description of the CFSR is available in the study by Saha et al. (2010).
Data from recently released fifth generation of ECMWF reanalysis, ERA5 (Hersbach & Dee, 2016), are also used in the study. ERA5 is a climate reanalysis data set, covering the period 1950 to present. Data processing for ERA5 is carried out by ECMWF, using ECMWF's Earth System model IFS, cycle 41r2. A significant improvement is observed in the performance of ERA5 over ERA-Interim. ERA5 is highly accurate, representing the magnitude and variability of near-surface air temperature and wind regimes. The higher spatial and temporal resolution provided by ERA5 reduces significantly the cold coastal biases identified in ERA-Interim and increases the accuracy representing the wind direction and wind speed in the escarpment (Tetzner et al., 2019).

In Situ Mooring Buoys
The RAMA was designed to study the Indian Ocean's role in the monsoons (McPhaden et al., 2009 Daily and monthly products of NSWR, NLWR, LHF, SHF, and Qnet are downloaded from their website http://uop.whoi.edu/projects/Bengal/QCData.html. We have included this buoy data to check the consistency of quality of the flux products and to make the statistical comparison more robust.

Taylor Diagram
The Taylor plot provides the statistics of correlation, a root-mean-square error (RMSE), and the ratio of variances together (Taylor, 2001). The distance from the origin is the standard deviation of the field, normalized by the standard deviation of the observational climatology. The distance from the reference point to the plotted point gives the RMSE. The correlation between the model and the climatology is the cosine of the polar angle. Thus, the data which have the most significant correlation coefficient, smaller RMSE, and comparable variance will be near to the observation and is considered the best among all others.

Validation With In Situ Buoy Data
The validation of average daily Qnet and its all components with independent buoys at three different locations using Taylor plot brings out the comparative assessment of flux products with respect to in situ values (Figures 1a-1c). Two buoys are located in the Bay of Bengal basin and one over the equatorial Indian Ocean ( Figure 4d). The Taylor Figure 1a). The radiative components (SWR = solid square; LWR = solid delta) of ERA5 have a higher correlation (SWR = 0.80; LWR = 0.96) as compared to other flux products, except for CFSR-LWR which performs equally good with CC values of 0.96. The variance of the SWR component is overestimated (underestimated) in NCEP-2 (CFSR) by almost quarter (half) of the total observed variance and marginally underestimated by less than 10th of the observed variance in case of MERRA, ERA5, and TropFlux. For LWR, all products overestimate the observed variance, and among them, TropFlux has marginally higher overestimation as compared to ERA5, CFSR, and MERRA. In the case of the turbulent flux (LHF = solid star) too, ERA5 performs the best in terms of both correlation and variance. Thus, considering both correlation and variance, ERA5 stands out as the best match with RAMA buoy not only in terms of Qnet but also in all components of it.

Earth and Space Science
MERRA underestimates it by about 25%, ERA5, and TropFlux by 18-20% and NCEP-2 overestimates it by 30% (Figure 1b). For the two major components (SWR and LHF), ERA5 has the largest correlations but underestimates the variance of SWR by almost 35%, whereas in the case of LHF, the variance is very close to that of the buoy.
However, for a robust statistical comparison of the flux products we have included all the collocated available data from the RAMA buoys during 2008 to 2015 irrespective of long gaps in between the data continuity. For buoy B1 (15°N, 90°E) and B2 (8°S, 67°E) the number of valid points is 1,879 and 915, respectively. For B3, we have considered one-year data which is 365 points. Tables 1-3 show this comparison statistics generated from five flux products (ERA5, CFSR, NCEP-2, MERRA, and TropFlux) together with the ensemble mean (EM) of the four data sets (ERA5, CFSR, TropFlux, and MERRA) with respect to three buoy locations. The values are almost similar to the results of Taylor plots with slight differences. Here too ERA5 (NCEP-2) shows the least (highest) RMSE and highest (least) CC in Qnet among all products at all the buoy locations.
All the new products (MERRA, CFSR, ERA5, and TropFlux) have improved significantly as compared to older reanalysis products (NCEP-2), which is very clear in the scatterplot of correlation and RMSE at these buoy location as shown in Figures 2a-2c. The position of scatter Qnet data points at B1 buoy (15°N, 90°E) corresponding to NCEP-2, MERRA, CFSR, TropFlux, ERA5, and Ensemble is diagonally arranged from down-right toward top-left corner of the plot indicating the improvement of product quality in terms of increased correlation and decreased RMSE, respectively ( Figure 2a). Qnet at other two buoy location also has the similar diagonal spread except that ERA5 takes over Ensemble in terms of correlation but RMSEwise Ensemble still fares well and sometimes different product switch position at different buoy locations (Figures 2b and 2c). SWR and LHF data points too are arranged in similar way at all buoy locations except LWR and SHF which show improvements in terms of correlation but not in term of RMSE. Further from Tables 1-3 the large underestimation bias of NCEP-2, often going up to more than 100 W/m 2 during some part of the year, is reduced significantly with all the four products being clustered around the buoy values. Interestingly the ensemble mean has improved further and outperformed ERA5 in most of the cases in terms of RMSE. It demonstrates that while individually ERA5 provides the best estimate of all the flux Note. Root-mean-square-error (RMSE) of the data products with respect to buoy are in brackets. components, the ensemble mean being close to ERA5 both in correlations and RMSE for all flux components provides the most reliable estimate by virtue of removal of some of the uncertainty in estimation of flux by each of the flux products as also evident by the lowest RMSE among the products in almost all the variables. Notwithstanding significant error remaining in the estimation of Qnet by flux products when compared to buoy locations, a considerable reduction of RMSE of Qnet estimates from~100 W/m 2 in NCEP-2 to~45 W/m 2 in the ensemble mean is significant progress since last few decades. These improvements can also be seen in the statistical values shown in Tables 1-3.
To further substantiate the Taylor plots (Figures 1a-1c) and Tables 1-3, we have also plotted the daily time series of Qnet from all the products and compared with RAMA and WHOI buoys at Bay of Bengal and southern western Indian Ocean location (Figures 3a-3c). The running mean of 15 days is used to smooth daily data by removing the high-frequency oscillations, which is also ideal for identifying the persistent biases in the Qnet values of the various products. The ensemble mean of all products except NCEP-2 is also plotted  for the consolidated view. We have not included NCEP-2 for the ensemble mean generation since NCEP-2 performs as an outlier in comparison to observations and other individual products. Inclusion of NCEP-2 degrades the ensemble mean. All the individual product could able to capture the observed daily variation with good accuracy. Among all, ERA5 shows the closest match as compared to others.
Two essential characteristics of the Qnet in the north Bay of Bengal (BoB; Figures 3a and 3c) are the steady increase from about −50 W/m 2 in June to about +100 W/m 2 followed by large amplitude (trough to peak often spanning 150 W/m 2 ) subseasonal oscillations of roughly one-month period during the summer monsoon season. In the south equatorial buoy location, there is no such increasing trend from winter to premonsoon months, and the subseasonal oscillations are vigorous in winter (December-April), but weak during the summer monsoon season. The causes behind this vigorous subseasonal oscillation during winter are discussed in the study by Saji et al. (2006) and Jayakumar and Gnanaseelan (2012). It is noteworthy that all the recent flux products estimate the increasing trend in the north BoB and subseasonal fluctuations with considerable fidelity. Also, in the south of the equator location, the lack of trend during the premonsoon months and vigorous subseasonal fluctuations in boreal winter are well captured. We believe that this is also major progress in the estimation of Qnet over the Indian Ocean.

Mean Qnet Distribution
As mentioned earlier, the reversing north-south dipole pattern of the Qnet is critical in driving a shallow meridional circulation in the tropical IO that plays a vital role in maintaining the annual cycle of the SST distribution. Therefore, in Figure 4, we examine the yearly mean Qnet distribution over the tropical Indian Ocean from all the products. Ideally, we would need in situ observations on 0.25 × 0.25 grid boxes to verify the accuracy of the spatial pattern from the products. In the absence of such observations, we shall use the ERA5 pattern as a reference and compare the patterns of other products with respect to it. Annually there is a net heat gain over most of the study domain (40°E-120°E, 30°S-30°N) in all the data sets except NCEP-2, which has very high heat loss particularly over the Southern Oceans (Figure 4a). We note that ERA5 has the lowest basin averaged annual heat gain (9.21 Wm −2 ) and is likely to be closer to the truth with CFSR, MERRA, TropFlux, and OAFlux overestimating this gain by on an average 15 Wm −2 (gain of these products ranging between 23 and 26 Wm −2 ) with respect to ERA5. The gain may be partly because of the

Earth and Space Science
POKHREL ET AL.
limited domain considered here and partly true responsible for basin-wide increasing trend of SST in the region (Roxy et al., 2014). The similarity among the spatial patterns of the annual mean Qnet in five products except for NCEP-2 is rather striking. To quantify the similarity, correlation/RMSE between the ERA5 Qnet pattern and those of CFSR, MERRA, TropFlux, and OAFlux are calculated and indicated in Figure 4. With high pattern correlation (0.94) and low RMSE (14.64 Wm −2 ), the spatial pattern of annual mean Qnet in CFSR and ERA5 are quite close. The other three products (MERRA, OAFlux, and TropFlux) are also close to ERA5 with correlations ranging from 0.83 to 0.89 and RMSE ranging from 15.37 to 19.46 Wm −2 . The intermonthly standard deviation is more or less identical across the products, with minimal appreciable changes indicating the temporal variability to be better. Figure 5 shows the basin-averaged mean annual cycle of Qnet, shortwave radiation (SWR), longwave radiation (LWR), and latent heat flux (LHF). The yearly cycle was computed from the monthly data during the common time frame (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009). We note that the asymmetry between summer and winter Qnet (Figure 5a) is minimum for ERA5 consistent with minimum annual mean heat gain of 9.21 Wm −2 . The systematic negative bias of NCEP-2 Qnet of the order of 30-50 Wm −2 throughout the year causing underestimation of heat gain in winter and overestimation of heat loss during summer. We further note that this bias in Qnet in NCEP-2 arises due to systematic underestimation of the SWR and systematic overestimation of LHF throughout the year (Figures 5b and 5d). Improvement of model and analysis system in CFSR has overcome the systematic negative bias in NCEP-2 to slight positive bias throughout the year as compared to ERA5

10.1029/2019EA000988
Earth and Space Science ( Figure 5a). This improvement has come mainly from improved LHF estimation by CFSR (Figure 5d), while the positive bias seems to have come from about 15 Wm −2 of SWR throughout the year (Figure 5b). It is interesting to note that the OAFlux has almost identical positive bias in Qnet like CFSR (Figure 5a), but unlike CFSR, it arises due to underestimation of LHF (Figure 5d) rather than SWF as in CFSR. On the other hand, the Qnet in TropFlux is also close to CFSR (Figure 5a) but is a result of a more complex combination of underestimation of SWR during summer (Figure 5b) and systematic negative bias of~15 Wm −2 in LHF throughout the year (Figure 5d). Our analysis indicates that while the annual mean Qnet in the four products (ERA5, CFSR, TropFlux, and OAFlux) are close to each other, there is a need for further reduction of the uncertainty of approximately 20 Wm −2 still remaining among the products.

Latitudinal Variation of Annual and Seasonal Qnet
What is the latitudinal dependence of the uncertainty in the estimation of Qnet? In pursuit of an answer to this question, we plotted the latitudinal variation of annual zonal average (longitude averaged from 40°E to 120°E) Qnet from all the products and found that NCEP-2 remains an outlier and the other five products clustered together. Therefore, in Figure 6a, we present the latitudinal profile of zonally averaged annual mean Qnet from NCEP-2 (red) along with the latitudinal profile of ensemble mean of zonally averaged annual mean Qnet (black). The spread among the five products is shown by the shaded envelope with the latitudinal mean spread indicating the uncertainty in the zonal mean Qnet among the products. With Qnet from ERA5 agreeing closely at buoy locations and the other four products clustering around it with uncertainty of~10 Wm −2 provided confidence that all these products are close to observations. The ensemble mean clearly shows the heat gain zone of tropics is from 15°S to 20°N, and thereafter it reduces at higher latitudes ( Figure 6a). The significant heat loss zone over southern latitudes in NCEP-2 data is visible in zonal mean values too, where NCEP-2 zonal average Qnet is consistently below its compatriots (CFSR, ERA5, MERRA, TropFlux, and OAFlux) from 30°S to 20°N and the highest difference is around 20°S where NCEP-2 values are almost less than 60 Wm −2 from other flux products.
The zonally averaged net heat flux in four different seasons shows the nature of latitude-wise pattern among the flux products (Figures 6b-6e). The close clustering of five products around ERA5 is evident in all seasons, and so is the NCEP-2 being an outlier in all four seasons. The reversal of the heat loss (gain) in the north (south) IO from December-January-February (DJF) (Figure 6b) to June-July-August (JJA) (Figure 6d) is noteworthy. The zonal average of Qnet in general shows that the regions in northern (southern) part of the equator gain (loss) most of the heat during spring (March-April-May (MAM); Figure 6c) and summer (JJA; Figure 6d) and loss (gain) heat during autumn (September-October-November; Figure 6e) and winter (DJF; Figure 6b), this is as per the solar position during the ensuing season. The uncertainty among the products is the least (10.80 Wm −2 ) during MAM, while during the other three seasons (DJF, JJA, September-October-November), it is around~15 Wm −2 . This kind of uncertainty among various flux products was also shown by  almost a decade earlier, indicating that the improvements among flux products have not reduced the uncertainty to a great extent, and still, there is an urgent need for the development.

Annual Mean Qnet Component Distribution
It is insightful to examine the role of each component on Qnet to pinpoint the region-wise relative importance of each of them on surface heat budget (Scott & Alexander, 1999;Tomita & Kubota, 2004). Figures 7a-7x show the annual mean distribution of net shortwave radiation (SWR), net longwave radiation (LWR), latent heat (LHF), and sensible heat (SHF) from all the six flux products over the north Indian Ocean. The positive values of SWR represent heat gain by the ocean surface; however, in the case of LWR, LHF, and SHF positive (negative) values represents heat loss (gain) by the ocean surface. The cloud-free zone over the north-western Indian ocean has the highest amounts of SWR (close to 220-250 Wm −2 ) in all the flux products (Figures 7a, 7e, 7i, 7m, 7q, and 7u), which is only limited to north-western area of the Arabian Sea in NCEP-2 ( Figure 7a) but extensively extended to the western coast of Africa and equatorial Indian Ocean in case of rest all products. The Bay of Bengal and eastern equatorial Indian Ocean region have the least SWR values due to persistent clouds throughout the year. In general SWR distribution is almost similar in all products except older generation reanalysis (i.e., NCEP-2) where a meagre value of SWR is evident at equatorial Indian Ocean regions, and that has improved significantly in its newer

10.1029/2019EA000988
Earth and Space Science

10.1029/2019EA000988
Earth and Space Science generation of reanalysis (CFSR). This may be one of the probable reasons for highly negative net heat flux in NCEP-2 as compared to other flux products.
Net longwave flux at the sea surface is the difference between the longwave radiated upward from the sea surface, and the longwave radiated downward from the atmosphere. It is contributed by the three components, the sea surface, clouds, and the atmosphere. Maximum heat loss at ocean surface due to LWR (close to 70-90 Wm −2 ) is concentrated over the northern tip of Arabian Sea, in all data (Figures 7b, 7f, 7j, 7n, 7r, and 7v), whose spatial extent is minimum in NCEP-2 data ( Figure 7b) and maximum in OAFlux data (Figure 7j). This region is the cloud-free zone, hence show large LWR in all products. MERRA (Figure 7f) data also show a high loss of LWR over the southern India Ocean, which is not very prominent in other data. Furthermore, this loss exists throughout the year, as evident in the monthly climatology of LWR, which is almost invariant (see Figure 5c). Equatorial Indian Ocean basin (15°S to 10°N) has the least loss of LWR in all the data, which is extended till far south (30°S) in NCEP-2 data. Overall, NCEP-2, followed by TropFlux, has the minimum loss of LWR over most of the basin, which is due to the less gain of SWR in the same data. TropFlux (Figure 7n), CFSR (Figure 7r), and ERA5 (Figure 7v) have almost identical LWR distribution and may be considered close to the observation.
Latent heat flux distribution clearly shows that the highest evaporative zone is situated at the south of the equator (between 10°S and 20°S) in all the data sets (Figures 7c, 7g, 7k, 7o, 7s, and 7w), which are also trade wind regions with very high wind speed. This also confirms the previous studies (e.g., Pokhrel et al., 2012). The very high overestimation of LHF in NCEP-2 (180-200 Wm −2 ) is considerably reduced in the present-day reanalysis (MERRA, CFSR, and ERA5) and the blended products (OAFlux and TropFlux). The low (high) values of SWR (LHF) over the southern equatorial region in NCEP-2 product lead to highly negative net heat flux, as shown in Figure 4a. The sensible heat flux depends on the difference in temperature between the sea surface and the surrounding air. All the data products show almost similar distribution of SHF over the Indian Ocean basin (Figures 7d, 7h, 7l, 7p, 7t, and 7x) with ocean gaining (losing) sensible heat from northern (southern) Indian Ocean basin.

Spatial Pattern of Seasonal Mean Qnet and Their Components
To see whether the spatial distribution pattern of Qnet and its component is persistent throughout the year, we examine the spatial distribution of seasonal mean Qnet and its component over our study domain. For brevity, only the winter season (DJF) and summer season (JJA) of individual flux products are shown. Figure 8 shows the winter season distribution of Qnet and its major components from all the flux products. In general, the DJF season is characterized by the dipole kind of net heat flux distribution, where region south (north) of 10°N of the Indian Ocean basin gains (losses) heat flux (Figures 8a, 8e, 8i, 8m, 8q, and 8u). The heat gain region is much more prominent as compared to the heat loss region of the Arabian Sea and Bay of Bengal combined. The gain in net heat flux is least (highest) in NCEP-2 (CFSR), and rest have almost similar magnitude. The major components (SWR, LWR, and LHF) confirm the share of each, resulting in Qnet distribution. Spatial distribution pattern of Qnet follows the spatial pattern of SWR, which is the most dominant component among all. The SWR is more (less) over the southern (northern) Indian Ocean, this is due to the distribution clouds during the winter season and most likely is associated with the seasonal migration of Sun, which in general has very less spatial extent over the southern Indian Ocean during winter season (Figures 8b,8f,8j,8n,8r,and 8v). This is also shown by Bony et al. (2000) using INSAT data, where the percentage of clear sky pixels is higher than 60% over the southern Indian Ocean. In case of LWR, the heat loss on the order of 40-60 Wm −2 is predominant over most of the Indian Ocean, except the northern stretches of both Arabian Sea and Bay of Bengal where the heat loss is much higher (~90-110 Wm −2 ) in all the products (Figures 8c, 8g, 8k, 8o, 8s, and 8w). NCEP-2 has the least heat loss due to LWR (Figure 8c). However, this loss is higher (60 Wm −2 ) in CFSR (Figure 8s), indicating the improvement in the reanalysis model from the previous generation. The higher or lower heat loss of LWR is due to higher or lower heat gain from SWR. In case of LHF too, NCEP-2 has the highest evaporative loss (~180-200 Wm −2 ) at the southern and northern ocean again (Figure 8d), which is almost less by the margin of 40-50 Wm −2 in rest of the flux products (Figures 8h, 8l, 8p, 8t, and 8x). Thus, comparatively, the fewer values of Qnet in NCEP-2 data result from less gain in terms of SWR and more loss in terms of LHF. Thus, CFSR shows remarkable improvement in all flux components as well.
In the JJA season, Qnet distribution is roughly opposite to that of the DJF season with southern (northern) zones of the Indian Ocean domain having heat loss (gain) with few exceptions (Figures 9a, 9e, 9i, 9m, 9q, and 9u). Similar to the DJF season, NCEP-2 (Figure 9a) is the outlier as compared to other products in terms of maximum heat loss. Surprisingly the latest reanalysis ERA5 (Figure 9u) shows the 10.1029/2019EA000988

Earth and Space Science
heat loss over the north Arabian and Bay of Bengal region, which is almost absent in other flux products. To understand this, we have compared different Qnet during JJA with RAMA observation in the "B1" location for the year (2009,2010,2013,2014,2015,2017,2018) and WHOI buoy for year 2015 (Figures S1 and S2 in the supporting information). The same is compared with ERA5 for respective locations. It is observed that although the correlation is high (around  Earth and Space Science   Table 4) as compared with the annual time series for any year (Tables 1-3). Further, both reanalysis products (NCEP-2 and MERRA) also have heat loss over the southern Arabian Sea and southern Bay of Bengal region with different magnitude (Figures 9e and 9i). However, this is mostly absent in both the blended flux products as well the newest reanalysis (OAFlux and TropFlux (Figures 9i and 9m) and CFSR and ERA5 (Figures 9q  and 9u)). Instead of heat loss, these two blended products and the newest reanalysis product (ERA5) show heat gain. This may be associated with the lack of representative assimilated data in the reanalysis product, which may be rectified in blended and latest reanalysis products by including the satellite data (Kumar et al., 2012;. The significant improvement in the SWR distribution is visible in CFSR ( Figure 9r) as compared to its predecessor NCEP-2 ( Figure 9b) wherein the heat gain from SWR is drastically increased by almost 100 Wm −2 north of 10°S implying a significant improvement in the physical processes affecting the SWR. In general, the heat gain due to SWR is higher (~180-240 Wm −2 ) over the northwestern part of the Indian Ocean basin as compared to the southern and eastern parts (~120-140 Wm −2 ) in all the flux products (Figures 9b, 9f, 9j, 9n, 9r, and 9v). The low values of SWR (~120-140 Wm −2 ) over the northern Bay of Bengal and eastern equatorial the Indian Ocean is due to the presence of deep clouds during the JJA season over this region, which obstructs the incoming shortwave radiation (Pokhrel & Sikka, 2013). Hatzianastassiou et al. (2005) have shown a strong anticorrelation between incoming shortwave radiation at the surface and cloud amount. Hence, the lower SWR values are associated with the regions where more convections and cloud cover present such as ITCZ, which is reflected in Figures 8 and 9. Further, the area and amplitude of Bay of Bengal heat gain vary in all products, which itself indicates a large number of uncertainties in the calculation of heat flux in various data sets. In LWR distribution NCEP-2 shows basin-wide low values (~30-50 Wm −2 ), unlike other products that show higher values over the southern Indian Ocean. These higher values over the southern Indian Ocean are mainly due to the absence of clouds (Figures 9c, 9g, 9k, 9o, 9s, and 9w). NCEP-2 has the least heat gain over the southern basin ( Figure 9b); correspondingly, it has the least LWR loss over the southern basin (Figure 9c). The LHF loss is the most important contributor toward the Qnet, which is predominantly higher in the NCEP-2 (Figure 9d) product as compared to the other products (Figures 7h, 7l, 7p, 7t, and 7x). This implies that the major contributor toward the highly negative Qnet over the Southern Oceans in the NCEP-2 data is due to higher latent heat loss, as this loss is subdued in other products resulting in very less negative Qnet over the Southern Oceans.
The annual and seasonal mean of Qnet and its major components in the ensemble mean are shown in (Figures 10a-10t). The ensemble mean Qnet distribution is the most reliable spatial distribution since its spatial pattern matches with satellite-based observations (see Paramil et al., 2016, Figure 2). The heat loss region of the southern Indian Ocean is now stabilized in the ensemble annual Qnet distribution (Figure 10a) as it is neither very high as that of CFSR (Figure 4e) nor very less as that of OAFlux (Figure 4c). Similarly, the season-wise distribution of Qnet has realistically captured the north-south dipole from the southern minimum in Qnet in MAM and JJA (Figures 10e and 10i) to northern minimum in September-October-November and DJF (Figures 10m and 10q; Parampil et al., 2016). Qnet Components are now much better represented annually and seasonally due to cancellation uncertainties in individual components. This establishes that the ensemble mean being the most stable product of net heat flux in the present scenario when the individual flux uncertainties are still high.

Lead-Lag Correlation of Qnet and Its Components With SST
As evident from Figure 3, the Intraseasonal Oscillation (ISO) features in Qnet is prominently visible in all the buoy locations. Many previous studies showed that SST ISO is mainly governed by Qnet over BoB (Sengupta & Ravichandran, 2001). Now to test the robustness of the products, it will be worthwhile to check whether the subseasonal relation of Qnet and its components with SST is reproduced accurately? The lead-lag correlation of Qnet, SWR, LHF, and wind speed with SST over the Bay of Bengal RAMA buoy location (15°N, 90°E) is plotted using buoy, and ensemble mean data using daily June to September data (Figures 11a and 11b). The correlation at 99% significance level is marked as black dash lines. Qnet leads SST by approximately six to eight days; this nature of the relationship is precisely captured by Ensemble, but the lead time is less (approximately three to four days). Similarly, two most significant drivers of Qnet, namely, SWR leads SST and LHF lags SST in both Buoy as well as Ensemble data, only the lead/lag days varies by one to two days and moreover the lead/lag has a stronger relationship in buoy data which is weaker in case of Ensemble mean data. Wind speed too lags SST by 10-12 days. This result corroborates a larger number of previous researchers who have studied the relation of SST and Net heat Flux and its components (e.g., Wu et al., 2008;Zhang et al., 1995;Sengupta and Ravichandan, 2001). Thus, the conformity that the Ensemble product can capture the physical relationship in the intraseasonal time scale between SST and other Flux components establishes its robustness.

Conclusion
Surface flux products are known to have significant uncertainties that arise from both the uncertainties in the input variables (u, v, U, qa, Ta, and Ts) and the uncertainties in estimation of transfer coefficients (cd, ce, and ch) in the bulk flux algorithms (Brunke et al., 2003;Isemer et al., 1989;Josey et al., 1999;Valdivieso et al., 2017). In this study, we try to show the uncertainties of widely used flux products over the Indian Ocean region. A comparative assessment of surface net heat flux (Qnet) and its components from four reanalyses (NCEP-2, MERRA, CFSR, and ERA5) and two blended products (OAFlux and TropFlux) is done over the Indian Ocean for the 10 years from 2000 to 2009. The daily time series of Qnet and its components (viz., SWR, LWR, and LHF) from each product is validated using observation from three independent buoys. The analysis revealed that both new generation reanalysis ERA5 and CFSR compares well with buoy observation as compared to the older generation reanalysis NCEP-2. The blended product TropFlux is also close to observation. ERA5, in general, outperform all other products for most of the cases with highest correlation (0.89-15N, 90°E; 0.74-8°S, 67°E; and 0.86-18°N, 89°E), least RMSE, and most explained variance with respect to the in situ values. The time series analysis clearly shows that five flux products (ERA5, CFSR, MERRA, OAFlux, and TropFlux) can capture the observed variation with reasonable accuracy. However, NCEP-2 is unable to obtain the observed variation, particularly during winter and summer. These five fluxes are very close to each other in terms of day-to-day variation; therefore, the ensemble of these five products are created, and the ensemble mean shows the best fit with the observed variation by virtue of removal of some of the uncertainty in the estimation of flux by each of the flux products. A significant reduction of RMSE of Qnet estimates from~100 W/m 2 in NCEP-2 to~45 W/m 2 in the ensemble mean, which is a major progress in recent times. Time series comparisons clearly show the ability of all the products to capture the summer ISO feature at BoB buoy location and wintertime ISO at the southern Indian Ocean buoy location.
The Qnet spatial distribution is strikingly close to each other except NCEP-2, which shows a significant negative bias over the southern Indian Ocean. This is confirmed by the in situ observations from the buoy located at the thermocline dome region (67°E, 8°S). A comparison with buoy observations showed that ERA5 performs best among all other products. It is hence assuming ERA5 as a reference to the pattern correlation between ERA5 and each of the product, which ranges from 0.83 to 0.94. The ensemble mean spatial distribution on the annual and seasonal scale shows the best realization since it gives the closest pattern seen in the previously reported observations (Parampil et al., 2016). The annual cycle of Qnet and its components (SWR, LWR, and LHF) brings out the improvement in newer generation CFSR with respect to older generation NCEP-2 by excellently improving LHF.
The zonal mean of ensemble product shows that the heat gain zone of tropics is from 15°S to 20°N, after that it reduces at higher latitudes, and in the case of NCEP-2 heat loss zone is much bigger than the heat loss zone. The reversal of the heat loss (gain) in the north (south) IO from DJF to JJA is noteworthy. The average zonal uncertainty of 14.86 (10.80) Wm −2 among product is highest (lowest) at JJA (MAM) season. These uncertainties among the product are not reduced considerably, as reported by Yu et al. (2006), thus indicating a need for improvement in observations, assimilation schemes, and model physics to improve the reanalysis products. Except for the uncertainty in the ensemble net heat flux product, the intraseasonal lead-lag relation of SST with Qnet, SWR, LHF, and wind speed captures the known physical relation, thus justifying its robustness.
This study will be useful to both modelers and observation to check the present status of the available observational flux products to improve the flux quality and to have proper flux to force the model.