Evaluation of Daytime Evaporative Fraction from MODIS TOA Radiances Using FLUXNET Observations

In recent decades, the land surface temperature/vegetation index (LST/NDVI) feature space has been widely used to estimate actual evapotranspiration (ETa) or evaporative fraction (EF, defined as the ratio of latent heat flux to surface available energy). Traditionally, it is essential to pre-process satellite top of atmosphere (TOA) radiances to obtain LST before estimating EF. However, pre-processing TOA radiances is a cumbersome task including corrections for atmospheric, adjacency and directional effects. Based on the contextual relationship between LST and NDVI, some studies proposed the direct use of TOA radiances instead of satellite retrieved LST products to estimate EF, and found that use of TOA radiances is applicable in some regional studies. The purpose of the present study is to test the robustness of the TOA radiances based EF estimation scheme over different climatic and surface conditions. Flux measurements from 16 FLUXNET (a global network of eddy covariance towers) sites were used to validate the Moderate Resolution Imaging Spectro radiometer (MODIS) TOA radiances estimated daytime EF. It is found that the EF estimates perform well across a wide variety of climate and biome types—Grasslands, crops, cropland/natural vegetation mosaic, closed shrublands, mixed forest, deciduous broadleaf forest, and savannas. The overall mean bias error (BIAS), mean absolute difference (MAD), root mean square difference (RMSD) and correlation coefficient (R) values for all the sites are 0.018, 0.147, 0.178 and 0.590, respectively, which are comparable with published results in the literature. We conclude that the direct use of measured TOA radiances instead of LST to estimate daytime EF can avoid complex atmospheric corrections associated with the satellite derived products, and would facilitate the relevant applications where minimum pre-processing is important. OPEN ACCESS Remote Sens. 2014, 6 5960


Introduction
Mapping of the land surface heat fluxes and modeling the mass and energy interactions between land and atmosphere are significant for better understanding the mechanism of climate change, and plays a crucial role in hydrological, agricultural and meteorological studies [1,2].The complicated physical mechanisms such as turbulent transport, the feedback in the soil-plant-atmosphere continuum and the heterogeneity of land surface all combine to make estimation of energy balance components a challenge [3,4].Nevertheless, much work has been done in surface turbulent fluxes and partitioning among energy balance components [5][6][7].As one fundamental parameter of surface heat fluxes, the evaporative fraction (EF, defined as the ratio of latent heat flux to available energy) represents the surface control on latent heat and sensible heat fluxes portioning [8,9].Considering the spatial and temporal variability of the EF characteristic, satellite remote sensing is recognized as a promising technique to provide reasonable EF estimates over large areas and continents [10].The approaches proposed in the literature use either empirical or physically based schemes [11][12][13].
One popular method among them is using the relationship between land surface temperature (LST) and normalized difference vegetation index (NDVI) [14].If the satellite derived LST and NDVI over heterogeneous areas are plotted, the shape of the pixel envelope resembles a physically meaningful triangular or trapezoidal feature space (Figure 1).The physical properties encapsulated in the LST/NDVI space are simple: The LST has low sensitivity over vegetated areas, but increased sensitivity over bare soil regions.The wet edge presents high EF because it has high thermal inertia and strong evaporative cooling.In contrast, the lowest EF occurs at dry edge due to its weakest evaporative cooling [15,16].The LST/NDVI method is unique in interpreting the spatial variations of LST and NDVI to infer EF without largely depending on information from ground-based observations.Due to its simplicity and relatively high accuracy, this type of approach has already been widely accepted and used [17][18][19].The main limitation of the triangle method is that a large number of pixels over a flat area with a wide range of soil wetness and fractional vegetation cover are required to make sure that the dry and wet limits exist in the triangular space, and at the same time relatively uniform atmospheric forcing [20][21][22].For a detailed description and discussion of the LST/NDVI triangle feature space method, the reader is referred to the reviews given by Carlson et al. [23] and Petropoulos et al. [24].
The notable advantage of this method is that it needs only satellite data to retrieve EF with minimum requirements of model inputs and ancillary data [25].However, when satellite data are used to represent traditionally ground-based measurement, such as LST, correcting satellite data based on radiative transfer theory are necessary.Because the information received by the satellite sensors is TOA (Top of Atmosphere) radiances that are affected by viewing angles and atmospheric variables [26].This actually poses great challenges for the remote sensing community to develop different correction procedures to eliminate the atmospheric attenuation effects [27].Although much effort has been devoted to establish sophisticated and streamlined data correction procedures, these procedures are troublesome and need independent observation data to satisfy the mathematical and physical constraints.These limitations practically increase the operational difficulty.Under this background, several studies focus on the EF estimation from the TOA radiances rather than satellite products, and found that use of TOA radiances is adequate to estimate EF [28][29][30].The feasibility of estimating EF using TOA radiances stems from the contextual relationship between LST and NDVI.For applications that depend on contextual information from remote sensing, absolute radiometric calibration to remote sensing data is thought to be unnecessary [31,32].A recent study by Peng et al. [33] investigated the general possibility of estimating Normalized Difference Temperature Index (NDTI, a key parameter for EF estimation) from TOA radiances through a physical understanding of Planck radiation law and radiative transfer equation, and performed a detailed sensitivity analysis of NDTI on surface and atmosphere variability.They concluded that the TOA radiances based estimates have similar level of accuracy as obtained using atmospherically corrected data products.However, as indicated by Peng et al. [28,33], more validation work of TOA radiances retrieved EF against ground-based measurements over different climatic and surface conditions still needs to be carried out.The main objective of this study is to evaluate the applicability and robustness of the TOA radiances based daytime EF estimation scheme through comparison with measurements from 16 FLUXNET (a global network of eddy covariance towers) sites, which have already been widely used for validation of estimates from remote sensing [34,35].These sites represent a wide range of climates and biome types-grasslands, crops, cropland/natural vegetation mosaic, closed shrublands, mixed forest, deciduous broadleaf forest, and savannas.

Remote Sensing Data
Moderate Resolution Imaging Spectro radiometer (MODIS) is the primary satellite sensor in the NASA Earth Observing System (EOS) for land, ocean and atmosphere research [36].Specifically, the Terra-MODIS Collection 5 data products used in this study mainly include MOD021KM, MOD03, MOD09GA and MOD35_L2.The MOD021KM (band 31) and MOD03 datasets are used to provide geolocated and calibrated thermal band TOA radiance.In addition, the MOD09GA product contains surface reflectance that is used to calculate NDVI in this study.It should be noted that the MODIS data used in this work were transformed from Hierarchical Data Format-Earth Observation System (HDF-EOS) swath format to a Universal Transverse Mercator (UTM) projected GeoTIFF image and resampleed for 1 km pixel size.A prerequisite for the applicability of the LST/NDVI method are homogeneous atmospheric conditions and flat topography.A square domain was therefore defined around each flux tower site and was used as the study area for the subsequent processing.The size of each domain is about 60,000 km 2 , representing a wide range of fractional vegetation cover and soil wetness.On the basis of MOD35_L2 cloud mask product, the clear sky condition was identified when larger than 85% of the study domain being clear.The number of clear sky days (clear sky at MODIS overpass time) for each flux site is shown in Table 1.

FLUXNET Observations
The publicly available FLUXNET observations were used to validate TOA radiances estimated EF. Through a suit of instruments, the flux towers can measure half-hourly averaged fluxes (net radiation, soil heat flux, latent heat flux and sensible heat flux) and meteorological data (e.g., air temperature, precipitation).The relevant FLUXNET methodologies and summaries could be found in the research of Aubinet et al. [37] and Baldocchi et al. [38].In this study, 16 FLUXNET sites were selected mainly according to the following criteria: (1) the land cover for the 1 km 2 area centered on the flux tower is homogeneous; (2) the land cover for a 10,000 km 2 area encompassing the flux tower is heterogeneous with a range of fractional vegetation cover, and the terrain of the area is also flat.More information about these sites is given in Table 1 [39][40][41][42][43][44][45][46][47][48][49][50][51][52].
These sites are located across Europe, North America and Africa and cover a broad range of land cover types (Table 1).According to the IGBP (International Geosphere-Biosphere Program) classification scheme, the sites are classified into seven main groups: grasslands, crops, cropland/natural vegetation mosaic, closed shrublands, mixed forest, deciduous broadleaf forest, and savannas.Further details about these sites are provided by the corresponding publications and references therein.
The eddy covariance technique is known to have problems with energy balance closure [53][54][55].Therefore, the eddy covariance measured heat fluxes need to be corrected for energy balance closure .One commonly used approach is the Bowen ratio method, which repartitions the mismatch in the energy budget closure in accordance to the Bowen ratio observed.Thus, the basic assumption is that the error in the budget closure can be distributed proportional to the ratio of the turbulent fluxes [54].Since this method is relatively simple and accurate, and has been successfully used by many applications [56][57][58], we also applied the method in the present study to correct the eddy covariance measured latent and sensible heat fluxes.

Methodology
According to the definition of EF, the instantaneous EF (dimensionless) calculated from instantaneous flux tower measurements can be written as: where R n is the surface net radiation (W• m −2 ), G is the ground heat flux (W• m −2 ), LE is latent heat flux (W• m −2 ) and H the sensible heat flux (W• m −2 ) at time t.Meanwhile, the daytime EF is determined using the following equation [10]: where the time difference t 2 − t 1 refers to the time from 8:00 LT to 17:00 LT in the present study.The daytime rather than daily was selected as study period, because the eddy covariance technique is more reliable during daytime [59].The daytime EF can also be estimated by averaging instantaneous EF from different time periods.The two approaches would have similar results if EF is stable during daytime [10,60].In practice, the latter approach is more sensitive to the errors in LE and H measurements when they have low absolute values during early morning and late afternoon.
On the basis of the LST/NDVI feature space obtained from plotting remotely sensed LST against NDVI, the parameterization of instantaneous EF using TOA radiances is given as follows: (3 where Δ is the slope of saturated vapor pressure at the air temperature (kPa• K −1 ) and γ is the psychrometric constant (kPa• K −1 ) [61].ϕ combines the effects of Budyko-Thornthwaite-Mather wetness parameter and Priestley Taylor coefficient, which accounts for aerodynamic and canopy resistances [62].In this study, ϕ is calculated using TOA radiances rather than LST products: (4) where ϕ max is the maximum ϕ without surface water stress and often set to 1.26 [63].L s is the observed TOA radiance for a given pixel whose NDVI value is NDVI i , L max and L min are the corresponding highest and smallest TOA radiance which have the same NDVI i value.In order to obtain the value of ϕ for each pixel, a three-step linear interpolation scheme based on the LST/NDVI triangle (Figure 1) is used in the following manner [14]: (1) determine the dry and wet edges in the triangular space; (2) the global minimum and maximum ϕ are respectively set to ϕ min = 0 for the driest bare soil pixel and ϕ max = 1.26 for the densely vegetated pixel with largest NDVI and lowest TOA radiance, then ϕ i min is linearly interpolated for each NDVI interval (NDVI i ) between ϕ min and ϕ max , and ϕ i max for each NDVI i is obtained from the lowest TOA radiance pixel with that NDVI interval (ϕ i max is generally set to ϕ i max = ϕ max =1.26); (3) ϕ i value within each NDVI interval is interpolated between the lowest TOA radiance pixel and highest TOA radiance pixel.Consequently, the ϕ value for each pixel can be calculated using Equation (4).

Algorithm Evaluation
A number of quantitative indices, including mean bias error (BIAS), mean absolute difference (MAD), root mean square difference (RMSD), relative error (RE) and correlation coefficient (R) are selected in this study to evaluate the model performance [64].Furthermore, the results were also compared with published studies.

Energy Imbalance of Flux Tower Measurements
The energy closure of the flux tower measurements was investigated for the selected clear sky case days.Figure 2 [55] and Foken [65].However, the causes of the lack of energy balance closure are still under discussion and might be related to systematic bias in instrumentation, neglected energy sinks, landscape heterogeneity and mismatch in source areas.

Can Near Noon Instantaneous EF Represent Daytime EF?
The EF has been found to be stable during daylight hours in many studies [11,66,67], making it possible to extrapolate instantaneous EF values to daytime scale.In this study, the instantaneous EF at MODIS overpass time was used to represent daytime EF value.To examine the feasibility of this assumption, the FLUXNET measurements were used to respectively estimate instantaneous and daytime EF using Equations ( 1) and ( 2). Figure 3 shows the comparison between the instantaneous EF at the time of MODIS overpass and the daytime EF for all the FLUXNET sites in our study.Table 2 presents the BIAS, MAD, RMSD, RE and R values for every FLUXNET site.It can be seen that these statistical indices values range from −0.035 to −0.011 for BIAS, from 0.023 to 0.050 for MAD, from 0.031 to 0.084 for RMSD, from −7.19% to −2.82% for RE and from 0.931 to 0.996 for R. On the whole, a good agreement and negligible bias between the instantaneous and daytime EF is obtained with BIAS = −0.020,MAD = 0.031, RMSD = 0.042 and RE (Relative Error) = −4.47%for all the sites.Besides, it can be observed that R = 0.977 (Table 2).These statistics results indicate the reliability of using EF constant assumption to estimate daytime EF without incurring substantial errors.

Evaluation of Daytime EF from MODIS TOA Radiances
The MODIS TOA radiances estimated daytime EF is evaluated with FLUXNET measured daytime EF.The comparison results are illustrated by Figure 4.In general, the derived EF agrees well with tower-measured EF with data points distributed around the dashed 1:1 line without a significant discrepancy.The grasslands, croplands, cropland/natural vegetation mosaic and mixed forest show good accuracy, whereas the savannas have slightly poorer performance.It may be attributed to its lack of full range of vegetation cover and surface temperature.The implicit assumption of the triangle method is that ETa primarily depends on soil moisture and vegetation cover.This assumption requires a heterogeneous area with a full range of possible soil moisture and vegetation fraction values, and at the same time relatively uniform atmospheric forcing [20].The poor performance of savannas further demonstrate the limitation of triangle method.
Table 3 gives a comprehensive summary of the statistical metrics for each FLUXNET site.It can be observed that the BIAS values range from −0.08 to 0.12.The RMSD values vary from 0.103 to 0.224, and the R values appear to be low to high ranging from −0.280 to 0.846.The overall BIAS, MAD, RMSD and R values for all the sites are 0.018, 0.147, 0.178 and 0.590, respectively.This suggests the feasibility of estimating EF with MODIS TOA radiances alone.If we excluded the results of savannas, the BIAS, MAD and RMSD are further reduced to 0.006, 0.138, and 0.168 with a better R of 0.648.The LST/NDVI feature space has already extensively been used to estimate EF.The performance of this method has also been reported in the literature and it has been shown that EF estimates based on TOA radiances result in comparable results than using products of geophysical surface variables [28,68].Table 4 summarizes the statistics of the differences between satellite products derived EF and observed EF from previous published studies.It can be seen that the statistical results found in this study are comparable to those reported previously.It further suggests that using TOA radiances can provide reasonable estimation accuracy for EF, while requiring less input data and preprocessing like in classical LST/NDVI feature space approaches.In terms of error assessment, most discrepancies between estimated and measured EF are less than 0.2, presented in Figure 4. On the basis of semi-empirical error analysis, Jiang et al. [69] found that the upper bounds of absolute error and relative error in LST/NDVI estimated EF are less than 0.25% and 33.3%, respectively.As precipitation could introduce a large error of EF estimation through enhancing soil moisture and leaf interception, we used plus symbols to indicate the days contaminated by precipitation which occurred before satellite overpass.Besides, Wang and Dickinson [7] found that the LST/NDVI method is most suitable for a growing season in middle latitude areas.Because the key assumption of triangle method is that the ETa is negatively correlated with surface temperature.However, the ETa in high latitudes and cold areas is generally positively correlated with temperature [70,71].In other words, the range of vegetation index should be large enough, and the soil moisture rather than air temperature or available energy is the key control of EF.It means the sample days that are not during growing season will have more uncertainties, which is in accordance with our results.In Figure 4, the points with x markers are the samples days outside the growing season.Furthermore, the use of TOA radiances itself can introduce slight errors, as Peng et al. [33] found that NDTI can be estimated from TOA radiances with an accuracy of 90%.Other sources of uncertainty of our results could be related to the scale discrepancy between the FLUXNET footprint measurements and satellite pixel estimates [35,72,73], the relative error between instantaneous EF and daily EF, as well as the linear parameterization of ϕ (without inclusion of wind speed and surface humidity) within the LST/NDVI feature space.Considering the simplicity and acceptable accuracy level of the LST/NDVI method, we conclude that the use of TOA radiances appears to be adequate for the estimation of EF, and would facilitate direct use of remote sensing data for the situations (e.g., multi-sensor studies, data assimilation) where minimum pre-processing is important.

Conclusions
The evaporative fraction (EF, defined as the ratio of latent heat flux to surface available energy) has been estimated from top of atmosphere (TOA) radiances instead of satellite retrieved land surface temperature (LST) products in some regional studies.In order to comprehensively examine the robustness of the TOA radiances based EF estimation scheme over different climate and surface conditions, flux measurements from 16 FLUXNET (a global network of eddy covariance towers) sites were used to validate the daytime EF estimated from Moderate Resolution Imaging Spectro radiometer (MODIS) TOA radiances.It was found that the EF estimated from TOA radiances perform well across a wide variety of climate conditions and biome types.The accuracy level is also comparable with published results in the literature.Besides, using FLUXNET measurements, the instantaneous EF is also found to be applicable for representing daytime EF, without incurring substantial errors.Overall, the present study, together with the work by Peng et al. [28,33] demonstrated that the direct use of measured TOA radiances instead of satellite retrieved LST products to estimate EF is feasible and applicable.The notable advantage of this approach is that no atmospheric corrections are required.This would facilitate data assimilation and multi-sensor studies due to minimal pre-processing requirements.The problem of using polar-orbiting satellites data such as MODIS is the temporal gaps due to cloud cover.One possible solution is utilizing the geostationary satellites data with high temporal resolution.Future work will examine the full capacities of the TOA radiances approach for geostationary satellites data.

Figure 2 .
Figure 2. Comparison of the observed available energy (R n − G) against the sum of latent heat and sensible fluxes (LE + H) after correction for energy balance closure for all the FLUXNET sites.CRm, Em and 1:1 line represent the averaged closure ratio, mean residual energy and perfect agreement, respectively.

Figure 3 .
Figure 3. Comparisons of MODIS overpass time EF and daytime average EF for all the FLUXNET sites.

Figure 4 .
Figure 4. Comparisons of daytime estimated and observed EF at FLUXNET sites for different biome types: (a) Grasslands; (b) Croplands; (c) Cropland/Natural Vegetation Mosaic; (d) Closed Shrublands; (e) Deciduous Broadleaf Forest; (f) Mixed Forest; (g) Savannas.Dashed line is the 1:1 line.Grey area represents that the discrepancies between estimated and observed EF are less than 0.2.The points with plus markers are the samples contaminated by precipitation.The x marker indicates those samples outside the growing season.

Table 1 .
Details about the FULXNET sites used in this study, Elev in the table refers to elevation.
shows the comparisons of observed daytime average available energy R n − G against measured turbulent fluxes LE + H after correction for energy balance closure for all the FLUXNET sites except CASF2 and CASF3, because these two sites both have no G measurements.The closure ratio CR = (LE + H)/(R n − G) varies from 0.57 to 0.98 with residual energy E = R n − LE − H − G ranging from −136.53 to − 8.83 W/m 2 at all the sites.Overall, the S (slope) and R from the linear least squares regression for observed LE + H and R n − G are 0.70 and 0.79, respectively.It can also be seen that the measured LE + H fluxes are generally less than R n − G for all these sites with the averaged closure ratio CR m of 0.78 and mean residual energy E m of 78.14 W/m 2 .Similar results were reported byWilson et al.

Table 2 .
Statistical results for the comparisons between instantaneous EF at MODIS overpass time and daytime average EF.

Table 3 .
Statistical results for the comparisons between estimated and observed EF for each FLUXNET site.

Table 4 .
Accuracy assessment of the LST/NDVI feature space method used to derive EF in the literature.