Satellite Retrieval of Surface Evapotranspiration with Nonparametric Approach : Accuracy Assessment over a Semiarid Region

Surface evapotranspiration (ET) is one of the key surface processes. Reliable estimation of regional ET solely from satellite data remains a challenge.This study applies recently proposednonparametric (NP) approach to retrieve surface ET, in terms of latent heat flux (LE), over a semiarid region. The involved input parameters are surface net radiation, land surface temperature, near-surface air temperature, and soil heat flux, all of which are retrievals or products of the Moderate-Resolution Imaging Spectroradiometer (MODIS). Field observations are used as ground references, which were obtained from six eddy covariance (EC) sites with different land covers including desert, Gobi, village, orchard, vegetable field, and wetland. Our results show that the accuracy of LE retrievals varies with EC sites with a determination of coefficient from 0.02 to 0.76, a bias from −221.56W/m to 143.77W/m, a relative error from 8.82% to 48.35%, and a root mean square error from 67.97W/m to 239.55W/m. The error mainly resulted from the uncertainties from MODIS products or the retrieval of net radiation and soil heat flux in nonvegetated region. It highlights the importance of accurate retrieval of the input parameters from satellite data, which are the ongoing tasks of remote sensing community.


Introduction
Evapotranspiration (ET) includes evaporation from various land surfaces and transpiration from vegetation [1].ET is interchangeable to the associated latent energy (LE) [2].It is a key land surface process in regulating regional hydrological and climatic characteristics.As the only way back to the atmosphere, global land ET returns about 60% of annual land precipitation to the atmosphere [3].Accurate estimation of ET is important for regional water resources management, especially in arid regions.
ET is intrinsically difficult to measure and predict especially at a large spatial scale.Various approaches are proposed to estimate ET in succession after the 18th century, including the Penman evaporation equation [4], the Penman-Monteith (P-M) combination equation [5], and the Priestley and Taylor approach [6].Traditional measurement or estimation of LE is mostly applicable at a point-scale [7][8][9], including eddy covariance (EC) techniques [10].However, representativeness of the point-scale measurement for a large area is generally problematic, and the dense coverage of point measurements is not feasible [11].
Alternatively, remote sensing technology can efficiently solve the representation limitation of the point measurement.However, it cannot observe LE directly, whereas it provides retrievals of geophysical parameters to estimate LE at a regional scale [12].Several retrieval algorithms appeared in the last two decades [13,14], including the triangle approach [15][16][17], the simplified surface energy balance index (S-SEBI) [18], the surface energy balance system (SEBS) [19], the three-temperature model [20][21][22], the MOD16 algorithm [23,24], and other P-M remote sensing methods [25,26].They are widely applied to estimate regional or global ET from remotely sensed data [25][26][27][28][29][30][31].The existing algorithms have 2 Advances in Meteorology a relative error from 10% to 40% [32,33] in their validation sites.For example, Gillies et al. [34] used aircraft scanner data and the tringle approach to retrieve LE in an area covered with grasslands, steppe-shrub, and tall-grass prairie.The LE retrievals had a root mean square error (RMSE) value of 22∼55 W/m 2 and a relative error (RE) of 10% ∼30%, in reference to field measures by EC techniques and Bowenradio approaches.Verstraeten et al. [28] used the advanced very high resolution radiometer (AVHRR) and the S-SEBI approach to retrieve LE in forestland.This achieved a RMSE (RE) of 35 W/m 2 (24%), compared to EC measures.Su [19] used the Moderate-Resolution Imaging Spectroradiometer (MODIS) data and the SEBS approach to estimate LE in wheat, corn, and rainforest areas.The LE retrievals had a RE of 25%, compared to EC measures.Xiong and Qiu [22] used the three-temperature approach and Landsat Thematic Mapper (TM) data to retrieve instantaneous LE in grassland and hills.Relative to Bowen ratio systems, LE had a RE of 4.65% ∼100% or 0.02∼0.20 mm h −1 during the satellite overpassing time.Mu et al. [23] evaluated the MOD16 products and reported an average bias (RE) of 0.31 mm day −1 (24.1%) for daily ET at FLUXNET-EC sites.Zhang et al. [25] applied normalized difference vegetation index-(NDVI-) based ET algorithm to assess global terrestrial LE using AVHRR GIMMS NDVI data.The daily results had a favorable accuracy with RMSE of about 10∼40 W/m 2 at 34 flux tower sites.Similarly, Leuning et al. [26] introduced a remotely sensed leaf area index-(LAI-) based P-M algorithm to calculate regional daily average evaporation using MODIS LAI products.At 15 flux sites globally, the systematic RMSE in daytime mean evaporation was in the range of 0.09∼ 0.50 mm day −1 , whereas the unsystematic component was in the range of 0.28∼0.71mm day −1 .
A nonparametric approach (NP) has been recently proposed for estimating surface evapotranspiration [35].It uses net radiation, surface air temperature, land surface temperature, and soil heat flux as the inputs, without the need of parameterizing surface resistance.All the necessary inputs are measurable, offering a novel but simple approach for practical use.The approach has been validated at 24 EC sites, yet it was not tested with remote sensing application.This paper applies the NP approach to estimate regional LE covering different surfaces, evaluates the accuracy of LE retrievals from MODIS data only, and identifies the error sources which are useful for improving retrieval accuracy.

Methodology
2.1.The Nonparametric Evapotranspiration Approach.Surface net radiation ( n ) is the net amount of radiation entering and leaving the Earth's surface.A part of  n is transformed into surface soil heat flux ( s ), and another part controls LE and sensible heat flux ( s ).In NP approach, a homogeneous terrestrial ground surface layer is assumed for a macrostate system, and Hamiltonian (potential energy plus kinetic energy) is the total energy of this system. n serves as the potential energy, whereas  s ,  s , and LE serve as kinetic energy.The land surface temperature (LST) serves as a generalized coordinate in this system.The approach calculates the partial differential equations of Hamiltonian with LST ( s ).The final forms are [35] where  is land surface emissivity (LSE),  a is near-surface air temperature (AT), Δ is the slope of saturated vapor pressure at temperature  a ,  is the psychometric constant, and  is the Stefan-Boltzmann's constant (5.67 × 10 −8 Wm −2 K −4 ). can be estimated by the near-surface pressure ().

Retrieval Algorithms for
where  0 is the solar constant at the top of the atmosphere (about 1367 W/m 2 ),  is the surface albedo,  is the solar zenith angle (SZA),  0 is the water vapor pressure,  a is the air emissivity, and  31 and  32 denote the emissivity in bands 31 and 32 of MODIS, respectively. sd is the downwell shortwave radiation. v is the latent heat of vaporization (2.5 × 10 6 J/kg),  v is the gas constant for water vapor (461 J kg −1 K −1 ), and  d is the dew point temperature at screen level.For the surface albedo, it is derived by the following equation [37]: where  1 ,  2 ,  3 ,  4 ,  5 , and  7 are the nadir BRDF-adjusted albedos in bands 1, 2, 3, 4, 5, and 7 of MODIS, respectively.At long time scales,  s is commonly assumed to be negligible.But at subdaily scale, the  s varies with the time of day, and the values of  s are not always negligible.It can be parameterized with the following equation [38]: Obviously,  s is regarded as a function of normalized difference vegetation index (NDVI) and  n .

Correction of EC Measures Used as Reference for Accuracy
Assessment.Although EC is the most accurate technique to measure the turbulent fluxes of sensible and latent heat, the energy balance can not be closed with EC data at the Earth surface with a closure of the energy balance of approximately 80% [39,40].In addition, relative to the fluctuations of LE accuracy,  s is measured in more reliable accuracy [41][42][43].So it is necessary to correct LE directly measured by EC for validation.
A preferred method can be derived from energy balance [44,45].On the large-scale homogeneous surface and steadystate condition, the corrected LE (ERLE) is where  EC is  s measured by EC.The ERLE is the reference of validation for remote sensing retrieved LE (RSLE).

Metrics for Accuracy Assessment.
The linear regression approach is used to describe the accordance between RSLE and ERLE.The determination of coefficient ( 2 ), slope, and intercept of the linear fit between RSLE and ERLE are subsequently obtained.The accordance is more satisfactory when the regression line is nearer to 1 : 1 line and  2 is higher.Their definitions are described as follows [46]: where   means retrieve values,   means reference values,  means the average of   ,  = 1, . . ., , and  is the number of pair data for comparison.
In addition, bias, relative error (RE), and RMSE are used to quantify errors of RSLE.Bias quantifies the average absolute difference between retrieve values and reference values.RE is the absolute value of the bias divided by the magnitude of the reference values.RMSE is the standard deviation of the retrieve values around reference values.Basically, the RMSE represents a combination of standard deviation and bias.
Their definitions are described as follows [47]:  Gobi).The locations of sites imply the similar climatic conditions.The underlying surfaces are homogeneous at all sites excepted for orchard and village sites.Based on the field visit, the fruit trees grow with bean seeding at orchard site.At village site, the underlying surface is composed of the bare soil, house, road, and trees.All instruments were intercompared over the Gobi between 14 and 24 May, 2012 [43].The intercompared and well agreed instruments, accompanied by the uniform data processing steps and standards, ensured data consistency, which guaranteed the reliability of validation [50].All selected sites were covered by different land cover types.Thus, for convenience, the site names were replaced by the types of underlying surfaces here.

Meteorological and Surface Flux Data.
Acquired from the HiWATER, the EC and AWS data span from June 25 to September 15 in 2012, the overlapping time span of EC data at all sites (Table 2).All AWS and EC data were produced, archived, and made available to the scientific community by the Cold and Arid Regions Science Data Center at Lanzhou [43,51].They were used for validation and error source analysis of RSLE.All LE validation references were obtained by  s (directly measured by EC),  n , and  s (measured by AWS) (see (5)).
The parameters obtained by AWS were averaged every 10 minutes, whereas the temporal resolution of EC was 30 where  lu ( ld ) is the surface upwelling (downwell) longwave radiation.Atmosphere Archive and Distribution System (LAADS), the National Aeronautics and Space Administration (NASA).In our study, both MYD and MCD (MYD07, MYD11, MYD13, and MCD43) were selected to retrieve LE.MYD07 provided the profile of air temperature and moisture to obtain the near-surface atmospheric temperature and dew temperature [53].MYD11 provided the land surface temperature and emissivity [54].MYD13 provided the 16-day NDVI [55,56].MCD43 provided 8-day nadir BRDF-adjusted albedos [57,58], including the albedo in band 1 to band 7. ASTER radiometer has 5 thermal infrared (TIR) bands to provide TIR spectral emissivity variations at 90 m spatial resolution [59].LSE product is produced by the Temperature and Emissivity Separation (TES) algorithm [60].In our study, ASTER product was also provided by the Cold and Arid Regions Science Data Center at Lanzhou [61,62].ASTER product was used to estimate LSE.LSE can be represented by the ASTER narrowband emissivities using the following linear equation [63]:

Remote Sensing
where  10 ∼  14 are the five ASTER narrowband emissivities.It was regarded as the real value of LSE because of the high accuracy and spatial resolution.

Data Processing.
Aiming to ascertain the applicability of LE retrieve algorithm, the validation and error source analysis of RSLE were made at the different sites (Figure 2).Firstly, the LE obtained by EC was corrected by  n ,  s (obtained by AWS), and  s (obtained by EC) in the way of energy residual correction (see (5)).Secondly,  n and  s were retrieved by MODIS product in Bisht's and Moran's algorithm (( 2) and ( 4)), respectively.Then, RSLE can be estimated in NP approach (see (1)).The errors of RSLE were derived from the discrepancies between RSLE and ERLE.Thirdly, to qualify the error contribution due to the input error of one parameter, that parameter (derived from MODIS) and the other actual parameters (measured by AWS/ASTER) were brought in NP approach to get LE estimation (see (1)).Similarly, all actual input parameters (measured by AWS/ASTER) were bought in NP approach to get LE estimation, too.The difference among these two estimations was regarded as the error contribution of the parameter at the retrieval moment [64].Fourthly, according to the error analysis, we searched for the probable ways to improve the retrieve accuracy.3 showed energy fluxes and environmental parameters measured by surface observations at the six sites.In the order of decreasing surface moistures, these sites were listed as wetland, vegetable, orchard, village, Gobi, and desert sites.The decreasing order matched pretty well with the vegetation abundances at these sites.Generally, there were higher LE in vegetated region with higher  n /LSE and lower  s /LST.In detail,  n was higher at vegetated than nonvegetated sites.Mean amount of  n was 625∼735 W/m 2 for vegetated sites, and it decreased to Mean LSE was in the range of 0.978∼0.981and 0.932∼0.975for vegetated and nonvegetated sites.On the contrary, LST was considerably lower at vegetated sites with typical values of 300∼305 K, compared to 315∼320 K at nonvegetated sites.The near-surface pressure was higher at vegetated sites than at nonvegetated sites (up to 89.11 kpa for wetland, and down to 83.28 kpa for desert).The difference of AT was little, and there was about 299 K of AT at all sites.These environmental parameters were the background of LE retrieve.

MODIS Retrievals of
n ,  s , and LE. n ,  s , and LE were all derived from satellite retrieve.Table 4 showed input parameters obtained by satellite retrieve at the six sites.The dew point temperature, pressure, and LSE were similar among all sites, with values of about 278 K, 78 kpa, and 0.96, respectively.Other parameters were different among sites.In detail, AT were 296∼299 K at all sites except for orchard (292 K) and desert (289 K) sites.Low LST appeared at vegetated sites (about 301 K), whereas high LST occurred at nonvegetated sites (300∼309 K).Similarly, albedo was slightly lower at vegetated sites with value of about 0.17, compared to about 0.19 at nonvegetated sites.Thus, generally, the retrieved  n was higher in vegetated region (650∼685 W/m 2 ) than in the nonvegetated region (570∼685 W/m 2 ).Considering the lower NDVI at nonvegetated sites (less than 0.2), the higher retrieved  s (more than 250 W/m 2 ) appeared there.On the basis of the retrieved  n and  s , RSLE was higher at vegetated sites (350∼450 W/m 2 ) than at the nonvegetated sites (120∼ 440 W/m 2 ).The instantaneous retrieve results of  n ,  s , and LE in a part of Zhangye region at 05:55 (UTC) in August 20, 2012, were shown in Figure 3.The distribution of LE,  n , and  s was in good accordance with the oasis-desert ecosystem.The desert (the east and south of the region) had lower  n because of the high albedo here.The desert also had higher  s because of bare surface.In addition, the oasis (the middle of the region) and wetland (the north of the region) had more evaporation than desert because of the irrigation.In view of retrieve values, the LE was up to 300∼400 W/m 2 in the region of oasis, whereas the LE decreased to 150∼250 W/m 2 in the region of desert.In general, the distribution of retrieve results was deemed to be reliable.

Accuracy Assessment of MODIS-Retrieved LE. Figure 4 revealed the relationship between ERLE (donated by 𝑥-axis)
and RSLE (donated by -axis) at the retrieval moment.In general, relative to ERLE in 30 minutes, the RSLE was generally accurate and underestimated with bias, RMSE, RE,  To validate the retrieved LE further, the LE directly observed by EC (ECLE) was also compared with RSLE in Figure 6.Similarly, the RSLE was in relatively good agreement with ECLE with  2 of 0.11∼0.36 at vegetated sites.Nevertheless, at the nonvegetated sites, the accordance disappeared with  2 of 0.05∼0.23.Relative to the ECLE in Figure 6, the ERLE matched better with retrieved surface LE in Figure 5, especially at nonvegetated sites.Thus, at least, the RSLE had a satisfactory accuracy at vegetated sites, and it also probably had a relatively good accuracy at nonvegetated sites.

Error Sources and Their Contributions to MODIS-Retrieved LE.
To reveal the error contributions of input error, the input error of each parameter was showed firstly (Figure 7). n was retrieved with low accuracy at village, orchard, and Gobi sites with bias value of more than 80 W/m 2 .There was a low accuracy of retrieved  s at almost all sites, especially at Gobi and desert sites (bias values of about 170 W/m 2 ).The large difference (about 10 kpa) of surface pressure appeared at Gobi and desert sites.At most sites, the biases of AT were 2∼5 K, except for orchard (−7 K) and desert (−10 K) sites.Similarly, the biases of LST were also 0∼4 K at vegetated sites, whereas they were more than 9 K at nonvegetated sites.The LSE difference between MODIS and ASTER products was less than 0.01 at all sites except for Gobi (−0.036) and desert sites (0.031).
On the basis of input errors, the error contribution can be revealed.Figure 8 showed the error contributions of each factor (shown as the line) and the error of RSLE (shown as the columns) at 6 sites.Except for orchard and village sites, the RSLE were in relatively satisfactory accuracy, and the biases were within −90∼50 W/m 2 at the other 4 sites.Based on the analysis of error sources, it was clear that the major error sources (inducing more than 40 W/m 2 RSLE error) were  n ,  s , LST, and AT at nonvegetated sites, with error contributions of 40∼110 W/m 2 , −120∼10 W/m 2 , >60 W/m 2 , and −100∼10 W/m 2 , respectively.At vegetated sites, input errors were not the dominant error sources of RSLE.
In detail, the large  s error contribution (causing more than 100 W/m 2 RSLE error) appeared at Gobi and desert sites.About 100 W/m 2 RSLE error caused by the  n error occurred at village and Gobi sites.The influence of AT and LST error on LE error was below 40 W/m 2 at most of sites, except for the LST at village and Gobi sites, the AT and LST at desert site, and AT at orchard site.The error contribution of LSE accounted for RSLE error as a small part (leading to less than 10 W/m 2 LE error).Similarly, the input errors of pressure affected RSLE error quite little (mostly less than 1 W/m 2 ).In addition, according to validation, the accuracy of RSLE was unsatisfactory at orchard and village with a bias of −221.56W/m 2 and 143.77W/m 2 , respectively.At these two sites, the main errors sources (inducing more than 100 W/m 2 RSLE error) were AT and  n .

Discussion
The improvement of input data can possibly benefit LE retrieve in the future. n and  s are not the direct input parameters of algorithms but the indirect parameters retrieved by Bisht's and Moran's algorithm [36,38].The improvement of  n and  s retrieval algorithms is helpful to the improvement of RSLE.In our study, the Bisht's algorithm under clear sky [36] is selected as  n retrieval algorithm.Considering the widespread cloud, Bisht's algorithm under all sky [65] can be chosen to broaden applications.Besides, the accuracy of retrieved  n is especially unsatisfactory in arid nonvegetated region [36].It is suitable to replace the retrieved  n by the observation directly measured by the Clouds and the Earth's Radiant Energy System (CERES) project [66].CERES's accuracy meets the demand of  n input in the region of desert [67,68].For  s , because of the subtle temporal variation of NDVI in the arid region,  s retrieved by Moran's algorithm varies slightly.Thus, some more sensitive algorithm can be chosen for  s , such as the  s retrieval algorithm in SEBS model or in the way of thermal inertia [19,27].
AT and the dew point temperature are derived from the atmospheric profile data of MYD07 in our study.They are not the values at the near-surface but the values of the atmospheric profile nearest to the surface.It is evident that the more accurate atmospheric information contributes to the retrieve of  n and LE.So the other Earth observations with the high accuracy (e.g., the Goddard Earth Observing System Model, Version 5 (GEOS-5)) are optional substitutions [69].For LST and LSE, they are derived from the MYD11 in the split-windows algorithm [54].Wan et al. [70] have reported that just the average of bands 31 and 32 emissivities could lead to an overestimation of LSE, especially in arid and semiarid region.Accordingly, LST is underestimated in these regions.The MYD21 C6 is estimated to be published in 2016, and it can supply the more accurate LST and LSE based on TES algorithm [71], especially in the arid and semiarid region.That means the remarkable improvement of  n and LE retrieve if we replace MYD11 by MYD21.
[48,49]dyArea and Ground Sites.As a typical inland river basin in the northwest of China, the Heihe River Basin is located between 97 ∘ 24  ∼102 ∘ 10  E and 37 ∘ 41  ∼42 ∘ 42  N and covers an area of approximately 130 000 km 2 .The selected 6 ground observation sites are parts of multiscale EC observation matrices belonging to Heihe Watershed Allied Telemetry Experimental Research (HiWATER)[48,49], and they are acquired over the Zhangye region (100 ∘ 25  E, 38 ∘ 51 N, 1519 m) in the middle reaches of Heihe River Basin (Figure1).A total of 6 EC sites measuring LE and  s are used for data analysis and accuracy assessment (Table1), accompanied by 6 automatic weather stations (AWS) which are used to measure the near-surface meteorological parameters.In view of spatial homogeneity and underlying representativeness, observations focus on six different areas with landscapes ranging from moist vegetated surfaces (vegetable, orchard, and wetland) to arid nonvegetated surfaces (village, desert, and

Table 1 :
Descriptive information of EC sites.
N Orchard Vegetable Wetland Desert Gobi Village Figure 1: Geolocation of EC observation sites and illustration of land surface.

Table 2 :
Datasets for analysis.
[52].Remote sensing data were obtained from the MODIS and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) products.The pixels located at the sites were obtained to retrieve LE.All chosen images were acquired under clear sky during 13:00 to 15:00 (local time) from June 25 to September 15, 2012.The temporal and spatial resolutions of these products were listed in Table2.MODIS are onboard NASA's Earth Observation System TERRA and AQUA satellites[52].As mentioned above, various products (MOD, MYD, and MCD) are provided by MODIS.They have been produced, archived, and made available to the scientific community by the Level 1 and s (EC) Figure 2: Schematic representation of the procedure for analysis in our study.

Table 3 :
Average of fluxes and environment parameters derived from AWS/ASTER at the six sites in the EC matrices, including net radiation ( n ), soil heat flux ( s ), land surface temperature ( s ), near-surface air temperature ( a ), land surface emissivity (LSE), near-surface pressure (), latent heat (LE), and sensible heat ( s ).The energy residual corrected (ER) LE is also revealed.

Table 4 :
Average of input parameters derived by MODIS products at the six sites, including land surface temperature ( s ), near-surface air temperature ( a ), land surface emissivity (LSE), near-surface pressure (), dew point temperature ( d ), albedo, and NDVI.The retrieved net radiation ( n ), soil heat flux ( s ), and latent heat flux (LE) are also shown.