Validity of Five Satellite-Based Latent Heat Flux Algorithms for Semi-arid Ecosystems

Abstract

Accurate estimation of latent heat flux (LE) is critical in characterizing semiarid ecosystems. Many LE algorithms have been developed during the past few decades. However, the algorithms have not been directly compared, particularly over global semiarid ecosystems. In this paper, we evaluated the performance of five LE models over semiarid ecosystems such as grassland, shrub, and savanna using the Fluxnet dataset of 68 eddy covariance (EC) sites during the period 2000–2009. We also used a modern-era retrospective analysis for research and applications (MERRA) dataset, the Normalized Difference Vegetation Index (NDVI) and Fractional Photosynthetically Active Radiation (FPAR) from the moderate resolution imaging spectroradiometer (MODIS) products; the leaf area index (LAI) from the global land surface satellite (GLASS) products; and the digital elevation model (DEM) from shuttle radar topography mission (SRTM30) dataset to generate LE at region scale during the period 2003–2006. The models were the moderate resolution imaging spectroradiometer LE (MOD16) algorithm, revised remote sensing based Penman–Monteith LE algorithm (RRS), the Priestley–Taylor LE algorithm of the Jet Propulsion Laboratory (PT-JPL), the modified satellite-based Priestley–Taylor LE algorithm (MS-PT), and the semi-empirical Penman LE algorithm (UMD). Direct comparison with ground measured LE showed the PT-JPL and MS-PT algorithms had relative high performance over semiarid ecosystems with the coefficient of determination (R²) ranging from 0.6 to 0.8 and root mean squared error (RMSE) of approximately 20 W/m². Empirical parameters in the structure algorithms of MOD16 and RRS, and calibrated coefficients of the UMD algorithm may be the cause of the reduced performance of these LE algorithms with R² ranging from 0.5 to 0.7 and RMSE ranging from 20 to 35 W/m² for MOD16, RRS and UMD. Sensitivity analysis showed that radiation and vegetation terms were the dominating variables affecting LE Fluxes in global semiarid ecosystem.

1. Introduction

Latent heat flux (LE) is an important component of water and energy exchange in terrestrial ecosystems, and is especially critical for semiarid ecosystems [1,2,3]. However, current research showed that changes in land use and global climate had a direct impact on the hydrological cycle in water-limited regions [4,5,6,7,8]. Peichl et al. [4] found that the choice of grassland management regime was the dominating factor on grassland ecosystem carbon, water, and energy exchanges in the maritime grassland. Brunsell et al. [5] assessed the impact of land cover change in North American mesic grassland on carbon and water cycling. Sun et al. [6,7] found that precipitation and LAI were the dominating factors on LE fluxes in semiarid ecosystem. A number of studies [9,10,11,12,13,14,15] showed semiarid ecosystems such as grassland, shrub, and savanna were extremely susceptible to the consequences of climate change. Considering the vital role of semiarid ecosystems LE in energy, water, and carbon cycle exchanges, accurately estimating LE over semiarid ecosystems is crucial to understanding and modeling ecosystem production, water balance, atmospheric circulation, droughts, and irrigation.

Many models to estimate LE on regional or global scales exist. These LE models fall into three categories: (1) Statistical and empirical models [16,17,18,19]. Wang et al. [17] developed an empirical model based on a regression method. (2) Physical mechanistic processes based models [20,21,22,23,24,25,26,27,28,29]. Norman et al. [20] built a two layer model with directional brightness temperature and its angle of view, vegetation index, net radiation, air temperature and wind speed. MODIS LE product algorithm was also a typical processes based model [25]. (3) Data assimilation models [30,31,32]. Pipunic et al. [30] showed that the extended Kalman assimilation filter had the potential to improve LE predictions. However, previous studies have shown that model input errors, spatial scale mismatch among different data sources, and limitations of the algorithms themselves all affect LE accuracy [33,34,35,36]. Ground measured variable uncertainties also had a great impact on LE estimation [37,38,39]; similarly, the uncertainty in satellite-derived LAI retrieval and NDVI translate into errors in the estimation of constraint factors, which often led to reduced performances of these LE algorithms. The unclosed energy problem of the EC method [40,41] and spatial scale mismatch between flux tower sites and satellite derived NDVI and LAI also influenced the accuracy of LE estimation. Finally, empirical parameters in the LE algorithm structure may lead to large uncertainties. Besides the uncertainties, there was the diversity in model structures, inherent assumptions, driving data and parameter values for different models. Therefore, it is difficult to select an appropriate LE algorithm for a certain biome.

Researchers have focused on comparisons between a limited numbers of models. Zhu et al. [42] compared the Katerji and Perrier (KP), Todorovic (TD) and the Priestley–Taylor (PT) models for an alpine grassland site during the growing seasons, and found the KP model performed well after a simple calibration. Nichols et al. [43] compared six traditional LE models over at the North Tower on the Bosque site. Jiménez et al. [44] compared 12 global monthly mean land surface LE products and found larger LE differences for rain forest land surfaces. Comparison of actual LE algorithms mainly focused on small regions near few EC sites [19,44,45], and large scale studies of different LE algorithms over global semiarid ecosystems were rare. Previous studies indicated that the model performances varied with different land surfaces such as irrigated grass, crops, natural prairies, and forests [33,46,47,48,49], and different conclusions about the performances have been obtained for different land surfaces [49]. Thus, comparisons incorporating in a wide range of vegetated surface (especially grassland and shrub land) under various climatic conditions were still required.

Global network of EC data across semiarid ecosystems provided a plausible base to validate performances of different LE algorithms. Moreover, satellite-based observations of land surface properties provided a useful technology for monitoring LE at large spatial scales [50]. Therefore, it had great advantages to compare different LE algorithms based on EC data and satellite data. In this study, we compared five satellite based LE algorithms, including MODIS algorithm (MOD16), revised remote sensing-based Penman–Monteith LE algorithm (RRS), Priestley–Taylor LE algorithm of Jet Propulsion Laboratory (PT-JPL), Modified satellite-based Priestley–Taylor algorithm (MS-PT) and Semi-empirical algorithm of the University of Maryland (UMD) across a variety of semiarid ecosystems types. The objectives of this study were to: (1) evaluate LE algorithms performance over grassland, shrub land, and savanna ecosystems based on 68 flux tower sites provided by EC towers network; and (2) assess the performances over time and space of three major semiarid biomes such as grassland, shrub land, and savanna.

2. Data

2.1. Data at EC Sites

Field data from a total of 68 EC towers of the FLUXNET, Aisaflux, Ameriflux and Africa [51] were used to validate the five satellite based LE algorithms. The EC towers provided a continuous half hour to yearly measurement of heat flux with a range of hundreds meters. And it has been considered to be the best way to directly measure LE [52,53]. These towers were equipped for half-hourly or hourly measurements at 2 m height of incident solar radiation (Rs, W/m²), relative humidity (RH), air temperature (Ta, K), diurnal air temperature range (ΔT, K), wind speed (Ws, m/s), vapor pressure (e, Pa), sensible heat flux (H, W/m²), surface net radiation (Rn, W/m²), and ground heat flux (G, W/m²), were used as inputs for the algorithms. Ground measured LE provided the validation. When available, data sets were gap-filled by site principal investigators, daily data were aggregated from half-hourly or hourly data without additional quality control. The tower sites covered three major global land surface semiarid biomes: grassland (GRA; 42 sites), savanna (SAW; 9 sites) and shrub land (SHR; 17 sites) (Figure 1). The in situ LE measurements from the towers were aggregated into daily measured data, filtered for a series of quality control constraints, Gap-filling techniques were employed to ensure high quality data sets. The data were selected to cover the period 2000–2009 with at least one year of reliable data (Table 1).

Table 1. Characteristics of validation data at the EC sites. Sites name, Latitude (Lat), longitude (Lon) and available years are shown here.

Site Name	Country	Biome	Lat	Lon	Available Years
AT-Neu	Austria	GRA	47.12	11.32	2002–2015
CA-Let	Canada	GRA	49.71	−112.94	1998–2015
CH-Oe1	Switzerland	GRA	47.29	7.73	2002–2008
CN-Du2	China	GRA	42.05	116.28	2005–2015
CZ-BK2	China	GRA	49.49	18.54	2004–2012
DE-Gri	Germany	GRA	50.95	13.51	2004–2013
DE-Meh	Germany	GRA	51.28	10.66	2003–2006
DK-Lva	Denmark	GRA	55.68	12.08	2004–2008
DS	China	GRA	44.09	113.57
ES-VDA	Spain	GRA	42.15	1.45	2002–2008
FI-Sii	Finland	GRA	61.83	24.19	2004–2012
FK	China	GRA	44.28	87.92
FR-Lq1	France	GRA	45.64	2.74	2002–2015
FR-Lq2	France	GRA	45.64	2.74	2002–2015
HU-Bug	Hungary	GRA	46.69	19.60	2002–2008
HU-Mat	Hungary	GRA	47.85	19.73	2004–2008
IE-Dri	Ireland	GRA	51.99	−8.75	2002–2013
IT-Amp	Italy	GRA	41.90	13.61	2002–2008
IT-Be2	Italy	GRA	46.00	13.03	2006–2012
IT-MBo	Italy	GRA	46.01	11.05	2002–2015
IT-Mal	Italy	GRA	46.11	11.70	2002–2007
KBU	Mongolia	GRA	47.21	108.74
NL-Ca1	Netherlands	GRA	51.97	4.93	2000–2015
NM	China	GRA	42.93	120.70
PT-Mi2	Portugal	GRA	38.48	−8.02	2004–2008
RU-Ha1	Russia	GRA	54.73	90.00	2002–2004
RU-Ha2	Russia	GRA	54.77	89.96	2002–2003
RU-Ha3	Russia	GRA	54.70	89.08	2004–2004
UK-EBu	United Kingdom	GRA	55.87	−3.21	2004–2012
UK-Tad	United Kingdom	GRA	51.21	−2.83	2000–2001
US-ARb	United States	GRA	35.55	−98.04	2005–2006
US-ARc	United States	GRA	35.55	−98.04	2005–2006
US-Aud	United States	GRA	31.59	−110.51	2002–2015
US-Dk1	United States	GRA	35.97	−79.09	2001–2015
US-FPe	United States	GRA	48.31	−105.10	1999–2015
US-Fwf	United States	GRA	35.45	−111.77	2005–2015
US-Goo	United States	GRA	34.25	−89.87	2002–2007
US-Var	United States	GRA	38.41	−120.95	2000–2015
US-Wkg	United States	GRA	31.74	−109.94	2003–2015
Xfs	China	GRA	44.13	116.33	2004–2006
Xi2	China	GRA	43.55	116.67	2006–2006
YZ	China	GRA	35.95	104.13	2008–2009
ZY	China	GRA	39.09	100.30
AU-How	Australia	SAW	−12.50	131.15	2001–2015
BW-Ghm	Botswana	SAW	−21.20	21.75	2003–2012
BW-Ma1	Botswana	SAW	−19.92	23.56	1999–2001
ES-LMa	Spain	SAW	39.94	−5.77	2004–2012
PR		SAW	38.09	−122.96
US-FR2	United States	SAW	29.95	−98.00	2004–2015
US-SRM	United States	SAW	31.82	−110.87	2004–2015
US-Ton	United States	SAW	38.43	−120.97	2001–2015
ZA-Kru	South Africa	SAW	−25.02	31.50	2000–2013
CA-NS6	Canada	SHR	55.92	−98.96	2000–2006
CA-NS7	Canada	SHR	56.64	−99.95	1999–2006
CA-SF3	Canada	SHR	54.09	−106.01	2002–2006
CG-Tch	Kinshasa	SHR	−4.29	11.66	2006–2009
CN-Ku2	China	SHR	40.38	108.55	2005–2015
ES-LJu	Spain	SHR	36.93	−2.75	2004–2013
IT-Pia	Italy	SHR	42.58	10.08	2002–2006
NI	Niger	SHR	13.48	2.18
RU-Cok	Russia	SHR	70.62	147.88	2003–2013
SP	China	SHR	37.53	105.80
US-Bn3	United States	SHR	63.92	−145.74	2003–2015
US-Los	United States	SHR	46.08	−89.98	2000–2015
US-SO2	United States	SHR	33.37	−116.62	1993–2015
US-SO3	United States	SHR	33.38	−116.62	1993–2015
US-SO4	United States	SHR	33.38	−116.64	2004–2015
US-Wi6	United States	SHR	46.62	−91.30	2002–2003

Table 1. Characteristics of validation data at the EC sites. Sites name, Latitude (Lat), longitude (Lon) and available years are shown here.

Site Name	Country	Biome	Lat	Lon	Available Years
AT-Neu	Austria	GRA	47.12	11.32	2002–2015
CA-Let	Canada	GRA	49.71	−112.94	1998–2015
CH-Oe1	Switzerland	GRA	47.29	7.73	2002–2008
CN-Du2	China	GRA	42.05	116.28	2005–2015
CZ-BK2	China	GRA	49.49	18.54	2004–2012
DE-Gri	Germany	GRA	50.95	13.51	2004–2013
DE-Meh	Germany	GRA	51.28	10.66	2003–2006
DK-Lva	Denmark	GRA	55.68	12.08	2004–2008
DS	China	GRA	44.09	113.57
ES-VDA	Spain	GRA	42.15	1.45	2002–2008
FI-Sii	Finland	GRA	61.83	24.19	2004–2012
FK	China	GRA	44.28	87.92
FR-Lq1	France	GRA	45.64	2.74	2002–2015
FR-Lq2	France	GRA	45.64	2.74	2002–2015
HU-Bug	Hungary	GRA	46.69	19.60	2002–2008
HU-Mat	Hungary	GRA	47.85	19.73	2004–2008
IE-Dri	Ireland	GRA	51.99	−8.75	2002–2013
IT-Amp	Italy	GRA	41.90	13.61	2002–2008
IT-Be2	Italy	GRA	46.00	13.03	2006–2012
IT-MBo	Italy	GRA	46.01	11.05	2002–2015
IT-Mal	Italy	GRA	46.11	11.70	2002–2007
KBU	Mongolia	GRA	47.21	108.74
NL-Ca1	Netherlands	GRA	51.97	4.93	2000–2015
NM	China	GRA	42.93	120.70
PT-Mi2	Portugal	GRA	38.48	−8.02	2004–2008
RU-Ha1	Russia	GRA	54.73	90.00	2002–2004
RU-Ha2	Russia	GRA	54.77	89.96	2002–2003
RU-Ha3	Russia	GRA	54.70	89.08	2004–2004
UK-EBu	United Kingdom	GRA	55.87	−3.21	2004–2012
UK-Tad	United Kingdom	GRA	51.21	−2.83	2000–2001
US-ARb	United States	GRA	35.55	−98.04	2005–2006
US-ARc	United States	GRA	35.55	−98.04	2005–2006
US-Aud	United States	GRA	31.59	−110.51	2002–2015
US-Dk1	United States	GRA	35.97	−79.09	2001–2015
US-FPe	United States	GRA	48.31	−105.10	1999–2015
US-Fwf	United States	GRA	35.45	−111.77	2005–2015
US-Goo	United States	GRA	34.25	−89.87	2002–2007
US-Var	United States	GRA	38.41	−120.95	2000–2015
US-Wkg	United States	GRA	31.74	−109.94	2003–2015
Xfs	China	GRA	44.13	116.33	2004–2006
Xi2	China	GRA	43.55	116.67	2006–2006
YZ	China	GRA	35.95	104.13	2008–2009
ZY	China	GRA	39.09	100.30
AU-How	Australia	SAW	−12.50	131.15	2001–2015
BW-Ghm	Botswana	SAW	−21.20	21.75	2003–2012
BW-Ma1	Botswana	SAW	−19.92	23.56	1999–2001
ES-LMa	Spain	SAW	39.94	−5.77	2004–2012
PR		SAW	38.09	−122.96
US-FR2	United States	SAW	29.95	−98.00	2004–2015
US-SRM	United States	SAW	31.82	−110.87	2004–2015
US-Ton	United States	SAW	38.43	−120.97	2001–2015
ZA-Kru	South Africa	SAW	−25.02	31.50	2000–2013
CA-NS6	Canada	SHR	55.92	−98.96	2000–2006
CA-NS7	Canada	SHR	56.64	−99.95	1999–2006
CA-SF3	Canada	SHR	54.09	−106.01	2002–2006
CG-Tch	Kinshasa	SHR	−4.29	11.66	2006–2009
CN-Ku2	China	SHR	40.38	108.55	2005–2015
ES-LJu	Spain	SHR	36.93	−2.75	2004–2013
IT-Pia	Italy	SHR	42.58	10.08	2002–2006
NI	Niger	SHR	13.48	2.18
RU-Cok	Russia	SHR	70.62	147.88	2003–2013
SP	China	SHR	37.53	105.80
US-Bn3	United States	SHR	63.92	−145.74	2003–2015
US-Los	United States	SHR	46.08	−89.98	2000–2015
US-SO2	United States	SHR	33.37	−116.62	1993–2015
US-SO3	United States	SHR	33.38	−116.62	1993–2015
US-SO4	United States	SHR	33.38	−116.64	2004–2015
US-Wi6	United States	SHR	46.62	−91.30	2002–2003

Figure 1. The 68 eddy covariance (EC) flux tower sites used in this study. GRA: grassland; SAW: savanna; SHR: shrubland. The base map is Ecoregion map provided by The Nature Conservancy [54].

Figure 1. The 68 eddy covariance (EC) flux tower sites used in this study. GRA: grassland; SAW: savanna; SHR: shrubland. The base map is Ecoregion map provided by The Nature Conservancy [54].

Although the EC technique was regarded as a good measurement for heat flux, the sum of H and LE, as measured by the EC method, was generally less than the available energy [55]. Dolman et al. [53] found that the energy closure ratio was about 0.8 averaged from more than 50 sites in Europe and North America. Previous studies [56,57,58] showed that energy budget correction substantially increased LE estimates and improved the energy closure ratio. Twine et al. [41] found EC sites appeared to underestimate LE and H fluxes systematically by 10%~30% based on the Southern Great Plains 1997 Hydrology Experiment and the closure issue became more significant upon consideration of the long-term water balance. Twine et al. [41] provided a correcting method which was used original LE flux divided by the energy closure ratio. And the ratio was calculated by original LE, sensible heat flux, net radiation and ground heat flux. Therefore, we used the method of Twine et al. [41] to correct the LE from the EC towers due to the unclosed energy problem [41,53].

$$LE=\frac{(R_n-G)}{(L{E}_{ori}+H_{ori})}\times L{E}_{ori}$$

(1)

where LE is the corrected latent heat flux, H_ori and LE_ori are the uncorrected sensible heat flux and latent heat flux, respectively.

NDVI and FPAR for the EC sites were determined from the MODIS products [59]. Leaf area index (LAI) was determined from the GLASS product [60]. The 8-day GLASS LAI (GLASS01A01), 8-day MODIS FPAR (MOD15A2) and 16-day MODIS NDVI (MOD13A2) data at 1 × 1 km spatial resolution were used for models verification at the EC sites. We used an average of surrounding pixels around the tower to get the NDVI and LAI values. The MODIS data were downloaded by the ftp link [61]. The GLASS data can be found at web site [62].

2.2. Data at Regional Scale

We employed five LE algorithms including one statistical and empirical model (UMD) and four physical mechanistic processes based models (MOD16, RRS, MS-PT and PT-JPL) to estimate regional LE based on satellite forcing data. We used the MODIS products with 0.05 degree spatial resolution including monthly NDVI (MOD13C2) [63], monthly FPAR (MOD15C2) and annual land cover type (MCD12C1). We also used monthly LAI product (GLASS01B01) with 0.05 degree spatial resolution from GLASS datasets. The digital elevation model (DEM) product was obtained from SRTM30 dataset. The DEM data can be found at web site [64].

Daily radiation and meteorological data including Rs, RH, Ta, ΔT, Ws and Rn were collected from MERRA datasets [65] for 2000–2009. Previous studies [66,67] showed that MERRA provided comparable results for the global and energy cycle research. MERRA data were collected from Goddard Earth Science Data and Information Services Center web site [68]. Considering the temporal and spatial resolution (½ degree latitude × ⅔ degree longitude) of the MERRA data, we averaged the daily MERRA data to monthly data sets, and we used spatial linear interpolation to match the spatial resolutions of the MODIS data. After that all the variables required in LE algorithms had consistent temporal and spatial resolutions.

3. Methods

3.1. LE Algorithms

3.1.1. MODIS Algorithm (MOD16)

Mu et al. (2011) [25] improved the MOD16 algorithm from Cleugh et al. [69] using the Penman–Monteith model by: (1) simplifying the calculation of the vegetation cover fraction with FPAR; (2) calculating LE as the sum of daytime and nighttime components; (3) improving calculations of aerodynamic, boundary-layer, and canopy resistance; (4) estimating the soil heat flux using available energy and simplified NDVI; (5) dividing the canopy into wet and dry components; and (6) separating moist from saturated wet soil surfaces. Their validation using 46 AmeriFlux showed that MOD16 enhanced the LE estimation [11]. The MOD16 algorithm has been successfully extended to generate a MODIS global terrestrial LE product driven by MODIS land cover, albedo, LAI/FPAR, and a GMAO daily meteorological reanalysis data set.

3.1.2. Revised Remote Sensing-Based Penman–Monteith LE Algorithm (RRS)

Considering advances in the physical structure of the PM equation, Mu et al. [34] developed the MOD16 algorithm with constraint parameters of air temperature and vapor pressure deficit (VPD) for different biomes. To reduce the effects of misclassification, Yuan et al. [70] applied consistent model parameters across different vegetation types to develop the RRS by revising the algorithm parameters, modifying the air temperature constraint for vegetation conductance, and improving calculation of the vegetation cover fraction using EVI (enhanced vegetation index). Chen et al. [35] found that the RRS algorithm showed improved performance compared to the MOD16 algorithm, validating on 23 EC flux tower sites in China.

3.1.3. Priestley–Taylor Algorithm (PT-JPL)

Priestley and Taylor [24] developed a simple LE algorithm by reducing the atmospheric control term in the PM equation and adding an empirical factor to avoid the complexity of parameterizations of both aerodynamic and surface resistance. Based on this algorithm, Fisher et al. [26] proposed a PT based LE algorithm (the Priestley–Taylor LE algorithm of Jet Propulsion Laboratory, Caltech (PT-JPL)) with atmospheric (RH and VPD) and eco-physiological (FPAR and LAI) to downscale potential Evaportranspiration (ET) to actual ET (ET, equivalent to LE). The total ET is calculated as the sum of ETc (canopy transpiration), ETs (soil evaporation) and ETi (interception evaporation). Each component is calculated using the Priestley–Taylor equation and the corresponding eco-physiological condition.

3.1.4. Modified Satellite-Based Priestley–Taylor Algorithm (MS-PT)

Yao et al. [28] developed the MS-PT, modifying the PT-JPL algorithm by parameterizing of vegetation transpiration, introducing soil moisture constraints calculated from the apparent thermal inertia (ATI), derived from diurnal temperature range (DT) and the revised linear two-source model (N95) [20,22,28]. Their model needed few inputs variables including Rn, Ta, DT, and NDVI to calculate saturated wet soil evaporation, unsaturated wet soil evaporation, vegetation transpiration, and interception of evaporation from vegetation. MS-PT performance was verified using 16 EC flux tower sites across China with good validation.

3.1.5. Semi-Empirical Algorithm of the University of Maryland (UMD)

Wang et al. [29] used the Penman equation [71] to develop a semi-empirical LE algorithm which was used to analyze long term variation of global LE because of a lack of long term satellite and ground measured data [29,32]. The input variables include the incident solar radiation, air temperature, VPD, RH, vegetation indices, and wind speed. Unlike other LE methods, wind speed was used to linearly parameterize aerodynamics conductance [72,73]. Validating the UMD model for 16 day average daily LE at 64 globally distributed flux tower sites, they found an average R² of 0.94 and average RMSE 17 W/m².

3.2. Data Analysis

To compare the five LE models at site scale and large spatial scale, there were several preprocessing steps required. For site scale comparison, we calculated the daily LE by collecting the meteorological variables from EC sites and NDVI and LAI from satellite data (MOD13A2, MOD15A2 and GLASS01A01). We used an average of surrounding pixels to get the vegetation index to match the EC site footprint. The five models were also compared based on MERRA meteorological data. Considering the coarse resolution of MERRA data, we used the values simulated in the pixels which the towers are placed. We calculated the coefficient of determination (R²), RMSE, Nash–Sutcliffe efficiency coefficient (NSE), and average bias for estimated LE to evaluate model performances [17,74]. Six EC sites were randomly selected to show the seasonal trend of LE by five LE models for each biome. Based on the daily LE results, we filtered the time series by averaging the daily results to 8-days results.

The LE models have different algorithm structures and input variables. Therefore, when comparing the LE models, it was crucial to examine the sensitivity of input variables for different models. Therefore we conducted a sensitivity analyses to assess the sensitivity of input variables for the five LE models [17,74]:

$$SC=\frac{\partial LE}{\partial I(x)}$$

(2)

where SC is the sensitivity coefficient, I(X) is the input variable (e.g., Rn, RH, Ta, and NDVI). The sensitivity coefficient provided a significance of correlation between relative error in input variables and relative error in LE. Hence, we firstly drew scatter plots of relative error between input variables and LE. Then we calculated sensitivity coefficient for each input variable to get the main factors of each LE model.

For large spatial scale comparison, we used MERRA daily data, monthly (FPAR)/LAI with 0.05 degree spatial resolution, monthly MODIS NDVI product with 0.05 degree spatial resolution; monthly global land surface satellite (GLASS) LAI product with 0.05 degree spatial resolution, annual land cover type (MCD12C1) product with 0.05 degree spatial resolution and DEM product. Before comparison, the spatial and temporal resolution of input variables should be identical. Thus, MERRA data were spatially interpolated to attain the spatial resolution matching MODIS data. We averaged the daily MERRA data to monthly data. After the data preprocessing, we calculated the five models LE over grassland, shrubland and savanna biomes based on four years data. And we also computed the bias which was calculated by LE values based on one model minus others.

4. Results

4.1. Validation

At the site scale, the model performance differed among the five models during the period 2000–2009. Generally, the LE models explained 60%~75% of LE variability over all measurements (Figure 2). The MS-PT and PT-JPL models showed good overall performance for all three biomes with relatively high R² values and low bias and RMSE values across grassland, shrub, and savanna ecosystems (Figure 3), which may be due to the simple characterization of LE variation limiting the uncertainties. Other model performances varied over the biomes. The median LE simulations of MOD16, RRS, MS-PT, PT-JPL and UMD varied across different ecosystems, which implied the inconsistent model performance over the study sites (Figure 3). According to NSE, the performance of MS-PT and PT-JPL were better (Figure 3).

We analyzed the daily LE results cover the period 2000–2009. The comparative results for GRA indicated that PT-JPL had good performance with the highest R² (0.76), lowest RMSE (19.16 W/m²) and highest NSE_avg (0.6). The next best model was MS-PT with R² = 0.73, RMSE = 20.29 W/m², and NSE = 0.59 followed by MOD16, RRS and UMD. As expected, MOD16 and RRS had almost similar R² (0.68 and 0.62) and RMSE (25.61 and 24.34 W/m²). However, both models had low NSE (0.40 and 0.35). UMD showed reduced performance with R² = 0.58, RMSE = 27.02 W/m² and NSE = 0.39. RRS and MS-PT had negative bias (−8.07 and −5.53 W/m², respectively) compared with MOD16, PT-JPL and UMD (8.14, 19.16 and 5.79 W/m²).

Figure 2. Comparison of statistical indicators including latent heat flux (LE) values, Nash–Sutcliffe efficiency coefficient (NSE), the coefficient of determination (R²) and root mean squared error (RMSE) by daily LE observations and simulations of five LE algorithms at 68 EC sites for GRA (grassland), SAW (savanna), SHR (shrubland) during the period 2000–2009.

Figure 2. Comparison of statistical indicators including latent heat flux (LE) values, Nash–Sutcliffe efficiency coefficient (NSE), the coefficient of determination (R²) and root mean squared error (RMSE) by daily LE observations and simulations of five LE algorithms at 68 EC sites for GRA (grassland), SAW (savanna), SHR (shrubland) during the period 2000–2009.

Figure 3. LE observed at the EC sites and predicted by the MOD16, RRS, PT, MS-PT and UMD algorithms during the period 2000–2009.

Figure 3. LE observed at the EC sites and predicted by the MOD16, RRS, PT, MS-PT and UMD algorithms during the period 2000–2009.

For SAW, the PT-JPL and MS-PT had relative lower RMSE (24.35 and 22.96 W/m²) and higher R² (0.72 and 0.71) compared to RRS and UMD. MOD16 has the highest RMSE (28.98 W/m²). However, all models showed low NSE except RRS, and overestimate LE ranging from 7 to 12 W/m². RRS had negative bias (−12.64 W/m²) which may be due to the model parameterization of surface resistance [34].

For SHR, MS-PT have highest R² (0.61) and lowest RMSE (20.98 W/m²). UMD was second with R² = 0.6 and RMSE = 21.28 W/m². Similar to the SAW biome, all models showed low NSE and almost overestimate LE ranging from 3 to 13 W/m².

The five MERRA driven models explained 50%~60% of LE variability over all measurements (Figure 4). The R² of all five models were close with a range from 0.45 to 0.59. Compared with in situ data, the MERRA driven LE models showed low performance over all biomes sites, which was likely a consequence of the uncertainty introduced by MERRA data.

Figure 4. Observed and estimated LE at EC sites driven by MERRA meteorological data during the period 2000–2009.

Figure 4. Observed and estimated LE at EC sites driven by MERRA meteorological data during the period 2000–2009.

Seasonal variations in LE by five models could be found at several sites. In general, the measured and predicted seasonal curves were in good agreement (Figure 5). All models showed high LE in summer and low LE in winter. This was due to the seasonal variation of Rn and vegetation index terms. Compared with other models, the estimated LE by PT-JPL is more close to the ground observed LE for all sites.

We calculated the sensitivity coefficient of the input variables for each model [17]. We drew a scatter plot of relative error between input variable and LE (Figure 6). The calculation of sensitivity coefficients demonstrated that Rn, Rs, LAI and NDVI had bigger impact on LE than the other input variables. Consequently, models energy (Rn or Rs) and vegetation (NDVI or LAI) terms were demonstrated high correlations with LE, accounting for 50%~85% of LE variation for all five models.

Figure 5. Measured and estimated time series for the 8 day LE average during the period 2000–2007. The GRA sites include US-FPe and US-Goo. The SAW sites include US-Ton and AU-How. The SHR sites include CA-NS6 and CA-NS7. Locations include FPe*—(Fort Peck); Goo*—Goodwin Creek; Ton*—Tonzi Ranch; How*—Howard Springs; NS6*—UCI 1989; NS7*—UCI 1998. US*—United States; CA*—Canada; AU*—Australia.

Figure 5. Measured and estimated time series for the 8 day LE average during the period 2000–2007. The GRA sites include US-FPe and US-Goo. The SAW sites include US-Ton and AU-How. The SHR sites include CA-NS6 and CA-NS7. Locations include FPe*—(Fort Peck); Goo*—Goodwin Creek; Ton*—Tonzi Ranch; How*—Howard Springs; NS6*—UCI 1989; NS7*—UCI 1998. US*—United States; CA*—Canada; AU*—Australia.

Figure 6. Sensitivity analysis for estimated LE by the (a) MOD16; (b) MS-PT; (c) PT-JPL; (d) RRS and (e) UMD algorithms with corresponding input variables. The black line is 1:1 line.

Figure 6. Sensitivity analysis for estimated LE by the (a) MOD16; (b) MS-PT; (c) PT-JPL; (d) RRS and (e) UMD algorithms with corresponding input variables. The black line is 1:1 line.

4.2. Spatial Distribution of LE

We analyzed daily EC sites data over the period 2000~2009 and found all algorithms had high correlation coefficients (0.60~0.96) at sites located in Europe and North America, and lower correlations at sites located in central Asia (Figure 7). PT-JPL and MS-PT had low RMSE (32~10 W/m²) compared to the other algorithms. LE bias, estimated by PT-JPL, was slightly positive in North America and central Asia, and slightly negative in Europe. LE bias for RRS was consistently negative at most sites. For MOD16, the bias tended to be positive for North America and Europe sites and negative for central Asia. PT-JPL and MS-PT showed high performance over most sites with high R², low RMSE and low bias.

The LE algorithms were also applied to estimate global semiarid ecosystems for 2003–2006 at spatial resolution of 0.05° using GMAO-MERRA meteorological data and MODIS products. Thus we produced annual average global terrestrial LE for GRA, SHR and SAW, respectively (Figure 8). The estimated LE by UMD model was generally higher than the other models in all regions. For GRA, average annual LE from MOD16 was 18.37 W/m², higher than MS-PT (13.62 W/m²) and RRS (18.30 W/m²), and lower than PT-JPL (27.95 W/m²) and UMD (41.91 W/m²). Despite differences in spatial LE distributions among different models, all of the models predicted high annual LE over temperate grassland in Argentina pampas, South of central North America prairie, and tropical grassland in Kenya, whereas arid and desert regions in Eurasian grassland and central North America prairie have low annual LE due to moisture limitations. For SAW, all models showed higher annual average LE than GRA. The UMD model yielded the highest annual average LE for SAW from (68.07 W/m²), followed by PT-JPL (52.57 W/m²), MOD16 (46.00 W/m²), RRS (41.71 W/m²), and MS-PT (39.20 W/m²). All of the models yielded high annual average LE over the Cerrado in southern Brazil and low annual average LE over tropical savannah in the Serengeti of northern Tanzania and southern Kenya. For SHR, MS-PT average annual LE was 6.24 W/m², which was the lowest among the five methods. UMD showed the highest average annual LE (26.57 W/m²) followed by PT-JPL. The MOD16 and RRS showed similar average annual LE (10.46 and 10.29 W/m²).

Figure 7. Spatial distribution of the correlation coefficients, RMSE and Bias for estimated daily LE calculated by MOD16, RRS, PT-JPL, MS-PT, and UMD over the period 2000–2009.

Figure 7. Spatial distribution of the correlation coefficients, RMSE and Bias for estimated daily LE calculated by MOD16, RRS, PT-JPL, MS-PT, and UMD over the period 2000–2009.

Figure 8. Estimated annual global terrestrial LE for grassland, savanna and shrubland averaged for 2003–2006 at spatial resolution of 0.05° from algorithms driven by GMAO-GERRA meteorology.

Figure 8. Estimated annual global terrestrial LE for grassland, savanna and shrubland averaged for 2003–2006 at spatial resolution of 0.05° from algorithms driven by GMAO-GERRA meteorology.

Our comparisons of five LE models demonstrated that distinct differences of five LE models occurred for different climate regions and land covers (Figure 9). For grassland, MOD16 predicted higher LE than MS-PT in Asia (about 2~11 W/m², 8%~164%), America (about 2~16 W/m², 20%~73%), and Europe (about 4~18 W/m², 4%~105%), and lower in central Africa (about −1~−12 W/m², −5%~−80%). PT-JPL simulated the higher LE than MS-PT model in Asia (about 6~28 W/m², 44%~135%), America (about 11~35 W/m², 82%~264%) and Europe (about 11~23 W/m², 40%~222%), and lower in central Africa (about −7~−3 W/m², −6%~−80%). MOD16 showed mostly lower LE than PT-JPL over nearly all regions (about −1~−30 W/m², −10%~−96%) except south Brazil, West Eurasian and England grassland. RRS and MOD16 showed a similar spatial pattern for grassland with difference −10 to 10 W/m². RRS showed lower estimated LE than PT-JPL (about −4~−57 W/m², −11%~−98%). UMD estimated higher LE than MOD16 (about 8~45 W/m², 64%~336%), MS-PT (about 7~41 W/m², 10%~305%), and RRS (about 8~45 W/m², 42%~288%) in all regions. MS-PT had lower estimated LE than PT-JPL (about −3~−30 W/m², −60%~−210%) and slightly lower estimated LE than RRS (−3~−12 W/m², −26%~−112%), except for central Africa.

Figure 9. Estimated annual global terrestrial LE biases for three biomes (grassland, savanna, and shrubland) averaged for 2003–2006 at spatial resolution of 0.05° from the algorithms driven by GMAO-GERRA meteorology. Bias was computed as the difference of the models.

Figure 9. Estimated annual global terrestrial LE biases for three biomes (grassland, savanna, and shrubland) averaged for 2003–2006 at spatial resolution of 0.05° from the algorithms driven by GMAO-GERRA meteorology. Bias was computed as the difference of the models.

For shrubland (Figure 9), MOD16 had higher values of LE than MS-PT (about 2~11 W/m², 4%~96%) almost in all regions except the Red Sea, and was larger than RRS at low latitudes (about 1~4 W/m², 4%~50%), but smaller at high latitudes (about −1 W/m², −4%~−7%), which was related to the air temperature constraint for vegetation conductance in the RRS algorithm. MOD16 was generally lower than PT-JPL in most regions (about −2~−47 W/m², −17%~−85%), although higher in central Australia (about 3~5 W/m², 20%~44%) and South Africa (about 2~11 W/m², 12%~102%). MS-PT was also larger than RRS at low latitudes (about 5~22 W/m², 25%~291%), but smaller at high latitudes (about −5~−11 W/m², −33%~−77%). MS-PT showed lower LE than PT-JPL in most regions (−1~−30 W/m², −13%~−63%) except central Australia. PT-JPL was generally larger than RRS (about 1~41 W/m², 10%~240%). Overall, UMD model showed the highest estimated LE.

For savanna (Figure 9), MOD16 showed larger estimated LE than MS-PT (about 4~28 W/m², 3%~41%) and RRS (about 5~18 W/m², 10%~28%) in all regions, and PT-JPL (about 3~9 W/m², 4%~28%) in south Brazil, South Zambia and South Chad. MS-PT was lower than PT-JPL (about −1~−18 W/m², −2%~−39%) except in South Zambia and South Chad, and MS-PT and RRS almost had only a slight difference −10 to 10 W/m². RRS was lower than PT-JPL (about −1~−22 W/m², −6%~−20%). Generally, UMD model had the highest estimated LE.

5. Discussion

5.1. Performance of the LE Algorithms

PT-JPL and MS-PT exhibited good performance over GRA, SAW and SHR sites, which could be mainly attributed to these models not requiring aerodynamic and surface resistances to reduce the uncertainties in forcing data. Net radiation and air temperature were the main driving forces for the PT-JPL and MS-PT models, and they generally had lower observation uncertainty [28,75]. This good overall performance of PT-JPL and MS-PT models was reported previously [26,36,76]. Vinukollu et al. [36] evaluated the PT-JPL algorithm for 12 EC tower sites and found that PT-JPL demonstrated the most consistent performance for most sites. Yao et al. [76] also found that simulated latent heat fluxes by MS-PT algorithm showed improved agreement with 40 EC sites compared to the PT-JPL algorithm.

Compared with PT-JPL and MS-PT, MOD16 and RRS models show reduced performances over GRA, SAW and SHR sites. These two models are both based on the Penman and Penman–Monteith equations, which require accurate estimation of canopy resistance [25,70]. In semiarid ecosystems, canopy resistance was sensitive to soil moisture [77]. The calculation of canopy resistance of MOD16 and RRS was driven by VPD, near-surface air temperature and soil moisture. However, relative humidity was used as a proxy for soil moisture in MOD16 and RRS models. Hence the complicated calculation of canopy resistance introduced larger uncertainty in LE estimation [36,78].

The UMD algorithm provided a simple yet robust method to build functional relationships between LE and predictor variables [29,79]. Despite the fact that it ignored explicit biophysical mechanisms, we found that UMD produced acceptable results in most grassland sites. Similarly, Wang et al. [29,79] found that the 16 day average daily LE can be reasonably predicted with an average correlation coefficient of 0.94 and average RMSE of 17 W/m² using 64 FLUXNET sites. However, large uncertainties might be introduced due to the representativeness of the limited training dataset [18,29,79]. Chen et al. [35] reported that the parameters of the UMD algorithm may have different combination due to the independent environmental factors of ET.

5.2. Spatial Differences of Five LE Algorithms

Although these LE models were developed by mechanistic processes, LE models showed large inter-model differences from −30 to 30 W/m² based on MERRA data for 2003–2006. MS-PT and PT-JPL performed well overall, which may be due to the simple parameterization of bio-physiological process between vegetation and the atmosphere.

MOD16 and RRS models were both developed based on the PM equation, but showed significantly different performances. This difference might attribute to the calculation of surface resistance. Mu et al. (2011) [25] used look-up table to determine the inputs variables for each biome. The stomatal conductance limited by minimum air temperatures, close VPD and open VPD. These variables were different in these two models. For MOD16, the open VPD is all equal to 650 Pa, close VPD varied within the range of 4400~3600 Pa and open minimum air temperatures is ranging from 12.2~8.61 °C. However, these variables in RRS were set to constant for all biomes with open VPD = 930 Pa, close VPD = 2500 Pa and open minimum air temperatures = 8.31 °C [70]. These differences caused the different calculations of constraint parameters in surface resistance estimation. Hence the difference would be more obvious in extreme environment where the values of temperature and VPD near the breakpoints in piecewise function during the calculation of surface resistance.

The UMD model showed the highest estimated LE in all biomes. The empirical parameters in the functional relationships between the LE and input variables could be a factor in the overestimation. For grassland, the largest difference between UMD and the other models occurred in North America, whereas for savanna, the largest difference in LE estimates between UMD and other models occurred in north Brazil and east Tanzania, and for shrubland in South America, West Australia and Somalia in Africa.

5.3. Limitations of LE Estimation Based on Five LE Algorithms

Our results of sensitivity analysis showed that radiation and vegetation indices were the vital meteorological variables related with LE across semiarid ecosystems (Figure 6). Thus, bias or uncertainty in these variables might introduce substantial uncertainty for the LE algorithms. In this study, MERRA reanalysis data were used as meteorological forcing data for LE estimation. However, algorithms performances based on MERRA data were lower than when based on ground measured inputs. One reason for this poor performance was attributed to bias in the MERRA radiation variables For example, Zhao et al. (2013) [80] found that MERRA surface solar radiation has an average bias error of +20.2 W/m² on monthly and annual scales from American FLUXNET sites, Zib et al. (2012) [81] reported an annual mean bias of 3.9 W/m² at two baseline surface radiation network (BSRN) sites for surface solar radiation, and Wang and Zeng [82] found an overestimation of up to 40 W/m² for surface solar radiation.

NDVI and LAI were often used to describe the vegetation photosynthesis and canopy conductance which was closely related to transpiration [28,83,84,85]. In many previous studies, NDVI or LAI were used as a proxy for vegetation moisture and RH for soil moisture to develop satellite-based LE algorithms [25,26,29,70,79]. However NDVI or LAI might fail to capture the vegetation dynamics of certain biomes and uncertainties in NDVI or LAI estimation translate to errors in the constraint factors or canopy conductance.

Bias might lead to LE estimate errors from scale mismatch between the input data resolution and the field measurement footprint. Kustas et al. [86] analyzed different pixel resolution of remote sensing inputs, and showed that variation in ET flux between corn and soybean fields cannot be effectively distinguished when the input is of the order of 1000 m. Zhang et al. (2010) [87] found that coarse NCEP-NCAR reanalysis (NNR) meteorology data (National Centers for Environmental Prediction and the National Center for Atmospheric Research) could introduce bias to match the local tower footprint in some regions. The empirical parameters in the algorithms can also introduce considerable uncertainty in calculating LE estimates.

6. Conclusions

We evaluated five satellite based LE algorithms: MOD16, RRS, PT-JPL, MS-PT and UMD over semiarid ecosystems based on 68 EC flux tower sites from the FLUXNET project. We compared the LE model performances with in situ and MERRA meteorological forcing data. The sensitivity of input variables for these LE algorithms was analyzed. Results of the sensitivity analysis indicated that overall radiation (Rn or Rs) and NDVI (LAI) were the major factors, and accounted for 50%–85% of the variation in LE estimates.

All the LE models produced acceptable results for most grassland, savanna, and shrubland sites. PT-JPL and MS-PT performed better than the other algorithms with lower bias and higher R², lower RMSE over semiarid ecosystems. UMD algorithm showed the highest LE estimates in all biomes compared with other models. All algorithms have high correlation coefficients ranging from 0.60 to 0.96 in Europe and North America, whereas have low correlation in central Asia. The results of evaluation among these algorithms showed that the five satellite-based LE models produced acceptable results in most grassland, savanna and shrubland sites. The uncertainties from the radiation and vegetation terms have great impact on final LE estimation by these algorithms.