Validation and Application of the Modified Satellite-Based Priestley-Taylor Algorithm for Mapping Terrestrial Evapotranspiration

Satellite-based vegetation indices (VIs) and Apparent Thermal Inertia (ATI) derived from temperature change provide valuable information for estimating evapotranspiration (LE) and detecting the onset and severity of drought. The modified satellite-based Priestley-Taylor (MS-PT) algorithm that we developed earlier, coupling both VI and ATI, is validated based on observed data from 40 flux towers distributed across the world on all continents. The validation results illustrate that the daily LE can be estimated with the Root Mean Square Error (RMSE) varying from 10.7 W/m 2 to 87.6 W/m 2 , OPEN ACCESS Remote Sens. 2014, 6 881 and with the square of correlation coefficient (R 2 ) from 0.41 to 0.89 (p < 0.01). Compared with the Priestley-Taylor-based LE (PT-JPL) algorithm, the MS-PT algorithm improves the LE estimates at most flux tower sites. Importantly, the MS-PT algorithm is also satisfactory in reproducing the inter-annual variability at flux tower sites with at least five years of data. The R 2 between measured and predicted annual LE anomalies is 0.42 (p = 0.02). The MS-PT algorithm is then applied to detect the variations of long-term terrestrial LE over Three-North Shelter Forest Region of China and to monitor global land surface drought. The MS-PT algorithm described here demonstrates the ability to map regional terrestrial LE and identify global soil moisture stress, without requiring precipitation information.


Introduction
Evapotranspiration (LE) is a major component of the earth's climate system and global water cycle, and it represents a crucial link between global water, energy and carbon exchanges [1][2][3][4].Although the current Eddy Covariance (ECOR) or Bowen Ratio (BR) systems at flux towers have provided point measurements of terrestrial LE, LE is inherently difficult to measure and predict especially at large spatial scales because sufficient ground observations will never be available [3][4][5].In contrast, remotely sensed data can be used as proxies for retrieving important controlling variables.Therefore, satellite-based estimate of temporal and spatial variations of LE is crucial for improving hydrological and agricultural management [6][7][8][9][10].
An alternative approach, the Priestley-Taylor (PT) algorithm, can be accurate where aerodynamic and surface resistance is not available and a coefficient multiplier, alpha, sets the equation equal to potential LE (PE) [17,21,[33][34][35][36].The general form of the Priestley-Taylor algorithm is: where LE is evapotranspiration in W/m 2 , ∆ is the slope of the saturated vapor pressure curve (kPa/°C), and γ is the psychrometric constant (kPa/°C ).R n and G represent the surface net radiation and the soil heat flux in Wm −2 .ɑ is the PT coefficient.LE cannot exceed R n − G without significant advection and convection.ɑ has a limited range between 0 and (Δ + γ)/Δ and in the standard application of the Priestley-Taylor method, ɑ equals to 1.26 over water or wet surfaces.ɑ affects the partition of the sensible and latent heat flux because the variation of air temperature can lead to the changes of both (Δ + γ)/Δ and sensible heat flux.Generally, (Δ + γ)/Δ can vary from −2.6% to −0.7% when air temperature increases from 10 °C to 40 °C [10,37].
Based on the Priestley-Taylor algorithm, a large number of revised approaches to estimate terrestrial LE have been built.Three types of schemes have been developed to parameterize Priestley-Taylor coefficient (ɑ) from remotely sensed data: (1) two-step interpolation scheme from the dry and wet edges in the Land Surface Temperature (LST)-Normalized Difference Vegetation Index (NDVI) triangular space [10,[37][38][39][40][41][42][43]; (2) eco-physiological constraints derived from vegetation indices (VIs) or vegetation fraction cover (f c ) [2,[11][12][13]17,21,34,[44][45][46]; and (3) the parameterization of key variables characterizing soil moisture using meteorological or remotely sensed data [17,18,21,35].However, deriving surface parameters from the vertex of LST-NDVI scatter-plots requires a continuum of soil moisture and vegetation status to provide a range of surface conditions and this method can not be applied under bare soil or full vegetation cover conditions.In addition, many parameterization schemes still need many variables, such as relative humidity (RH) or precipitation, to improve the complexity of the Priestley-Taylor algorithms, and errors derived from many input variables introduce the uncertainty of LE estimations.
To overcome the difficulty of the satellite-based estimation of relative humidity (RH) and vapor pressure deficit (VPD) for traditional Priestley-Taylor algorithm, Yao et al. [21] developed a modified satellite-based Priestley-Taylor (MS-PT) algorithm using the Apparent Thermal Inertia (ATI) derived from the temperature change and NDVI derived from remote sensing products.This model was validated over 16 eddy covariance flux towers in China with an average R 2 of 0.86 and 9% bias, and applied to analyze the variations of terrestrial LE in China [21].However, this algorithm has not been analyzed and validated its efficacy in LE estimation at other regions.Moreover, this algorithm has only used to detect the variations of LE in China, between 2001 and 2010, due to the short-term sparse regional measurement datasets.
In this study, we present an overview and validation of the modeling algorithm, and describe two current international projects involving hydrological impacts of Three-North Shelter Forest region of China and global surface drought monitoring.It has three major objectives.First, we validate the modified satellite-based Priestley-Taylor algorithm (MS-PT), based on ground-observed data from 40 flux towers distributed across the world on all continents.Second, we evaluate both MS-PT algorithm and PT-JPL algorithm based on flux towers data.Finally, we calculate decadal variation in LE of the Three-North Shelter Forest region of China during 1984-2010 and detect the trends of global land surface drought from 1984-2007.

MS-PT Algorithm
The modified satellite-based Priestley-Taylor algorithm (MS-PT) [21] was specifically designed to minimize the need for ancillary meteorological data while maintaining a physically realistic representation of evapotranspiration process.It only needs four variables, as follows: surface net radiation (R n ), air temperature (T a ), diurnal air temperature range (DT), and NDVI.
MS-PT algorithm estimate LE by calculating the sum of the unsaturated soil evaporation (LE s ), the canopy transpiration (LE c ), the saturated wet soil surface evaporation (LE ws ), and the canopy interception evaporation (LE ic ) (Figure 1).The total LE can be expressed as: Unsaturated soil evaporation can be calculated using an index of soil water deficit (f sm ) and f sm can be acquired from an exponential algorithm of the Apparent Thermal Inertia (ATI), namely, where f wet is the wet surface fraction and can be derived from f sm , DT max is the maximum diurnal air temperature range (DT max = 40 °C ), R ns is the surface net radiation to the soil and can be calculated using both R n and vegetation cover fraction . G is also derived from R n and f c (a g (1 − f c )R n , a g = 0.18).In MS-PT algorithm, we calculated f c using NDVI: where NDVI min and NDVI max were the minimum and maximum NDVI during the study period, set as constants of 0.05 and 0.95 [7,21], respectively.Canopy transpiration can be estimated using the modified Linear Two-Source Model (N95) and can be described as: where f T is the plant temperature constraint and can be calculated using T a and an optimum T opt set as 25 °C.R nv is the surface net radiation to the vegetation and can be calculated using both R n and f c (R nv = R n f c ).Both saturated wet soil surface evaporation and vegetation interception evaporation can be calculated from the following two equations, respectively.

PT-JPL Algorithm
Based on the Priestley-Taylor algorithm, Fisher et al. [17] put forward PT-JPL algorithm by downscaling potential LE from Priestley-Taylor to actual LE.PT-JPL algorithm includes many eco-physiological constraint functions by introducing atmospheric moisture (VPD and RH) and vegetation indices.The PT-JPL algorithm can be expressed as: ) ) ( ) where f c is the green canopy fraction (f APAR /f IPAR ), f M is a plant moisture constraint (f APAR /f APARmax ), f sm is a soil moisture constraint (RH VPD ) and f wet is the relative surface wetness (RH 4 ), f APAR is the absorbed photosynthetically active radiation (FPAR), f IPAR is the intercepted PAR.

Drought Index and Potential LE Calculations
To apply the MS-PT algorithm for global surface drought monitoring, we select the Evaporative Drought Index (EDI), incorporating actual LE and potential LE (PE), to characterize the results of soil moisture response to surface dryness [47].The EDI can be described as: Yao et al. [47] described the physical implications of EDI and monitored surface drought over the conterminous United States using Moderate-Resolution Imaging Spectroradiometer (MODIS) and National Centers for Environmental Prediction (NCEP) Reanalysis-2 Data.Additionally, we also adopt the Hargreaves method to estimate PE.Although the Hargreaves method mainly performs effectively for well-cropped grass, the results of EDI calculation from this method can be accepted on a global scale in this study [48].Based on the Hargreaves model, the PE can be easily estimated from NCEP-2 data.The Hargreaves model is expressed as follows: Here, PE is the potential evapotranspiration (mm/day) and R a is the extraterrestrial solar incident radiation (MJ/m 2 per day).T max is monthly maximum air temperature (°C ), and T min is monthly minimum air temperature (°C ).

Trend Analysis
Linear trend analysis is used to analyze the regional long-term LE trend of Three-North Shelter Forest Region of China and to explore the variation of global land surface EDI, respectively.A simple linear regression equation is used to calculate the annual values and to obtain the long-term LE and EDI trends.bt y y t   0 (17) Here, y t represents the annual value of LE or EDI, t is the year and coefficient b represents the trend of long-term annual LE or EDI.Moreover, the Student's t-test distribution with n − 2 degrees of freedom is adopted to calculate the confidence levels of the derived tendencies [49].We also apply the linear trend analysis pixel-wise to calculate the trend for each pixel and calculate regional average time series, and, then, apply the linear trend analysis to quantify the regional trends.

Eddy Covariance Flux Towers
To evaluate the performance of both MS-PT algorithm and PT-JPL algorithm, we use the ground-observed data from 40 flux towers distributed the world: ten AmeriFlux sites, five AsiaFlux sites, one Atmospheric Radiation Measurement (ARM) site, one Chinese Ecosystem Research Network (CERN) site, one Asian Automatic Weather Station Network Project (ANN) site supported by the GEWEX (Global Energy and Water cycle EXperiment) Asian Monsoon Experiment (GAME AAN) and twenty-three other flux towers sites (Table 1 and Figure 2).These data sets include the longest continuous worldwide multi-site measurements of LE, R n , and corresponding meteorological observations, and each tower provides at least one year of reliable data.The land cover types of the flux towers sites include grasslands, crop, shrub, savanna, wetlands, evergreen forest, deciduous forest, and mixed forests (Table 1).LE collected from all flux towers is measured by the Eddy Covariance (ECOR) method.Although the ECOR method has been widely used in global measurement experiment, this method does not conserve energy [19,20].In this study, we have selected the method developed by Twine et al. [50] to correct the LE from all flux towers.Table 1.A description of site conditions.Land cover types, Latitude (Lat), longitude (Lon), Elevation (Elev, meter), time-period, and network names are shown here.For global land surface drought monitoring, we use global monthly surface downward and upward shortwave and long-wave radiation products at a spatial resolution of 1° × 1° from 1984 through 2007 that are derived from GEWEX Surface Radiation Budget (SRB) products.The monthly air mean temperature, and the maximum and minimum air temperature data are extracted from NCEP-2 data, which are acquired from the NCEP/NCAR Reanalysis Project (CDAS).These datasets have a spatial resolution of 1.875° longitude by approximately 1.9° latitude, and are interpolated into 1° × 1° using bilinear interpolation.The monthly GIMMIS NDVI products at a spatial resolution of 8 km are also interpolated into 1 degree.We also use the monthly Palmer Drought Severity Index (PDSI) products derived from the NCAR CGD's Climate Analysis Section dataset with a 2.5-degree spatial resolution, for the period from 1984 through 2007.To highlight the obvious features of both PDSI and EDI, we have interpolated the PDSI products into 1 degree from a 2.5-degree spatial resolution.

Validation and Comparison
To evaluate the ability of the MS-PT method to predict the spatial variation in LE, we have validated both MS-PT algorithm and PT-JPL algorithm based on the collected ground-measured data from all flux towers sites.Table 2 shows the Root Mean Square Error (RMSE), the bias, and the square of correlation coefficients (R 2 ) of the comparison between the ground-measured and estimated daily LE from 40 flux towers.One can observe that the RMSE of the estimated daily LE using MS-PT algorithm varies from 10.7 W/m 2 to 87.6 W/m 2 , the bias varies from −23.7 W/m 2 to 48.6 W/m 2 , and R 2 varies from 0.41 to 0.89 (p < 0.01).Similarly, the RMSE of the estimated daily LE using PT-JPL algorithm varies from 11.3 W/m 2 to 89.1 W/m 2 , the bias varies from −21.3 W/m 2 to 56.2 W/m 2 , and R 2 varies from 0.40 to 0.88 (p < 0.01).Overall, as compared with PT-JPL algorithm, MS-PT algorithm improves the LE estimates at most flux towers sites.Figure 3 shows an example of the eight-day time series of modeled LE, using both MS-PT algorithm and PT-JPL algorithm driven by tower daily meteorological measurements versus the corresponding tower LE measurements.The results illustrate MS-PT algorithm performs better than PT-JPL algorithm.The MS-PT algorithm provides the favorable agreement with the tower observations and captures observed LE seasonality and associated differences among the major land cover types.Figure 4 illustrates scatter plots of a comparison between annual estimated and ground-measured LE using MS-PT algorithm driven by the ground observations.We notice that the bias of the estimated LE at all 40 flux towers sites is 2.5 W/m 2 .The RMSE is 28.4 W/m 2 and the R 2 is 0.68 (p < 0.01).The accuracy of the LE simulation can be used for estimating the regional or global land surface LE.As the flux towers sites have different land cover types and different climate regimes, the comparison of the site-averaged LE demonstrates the ability of this method to predict the spatial variation in LE.
To fairly assess how well the model predicts long-term variations in LE, we have validated the estimated annual LE anomalies with the observed annual LE anomalies.We only use flux towers sites where five years of data are available.As shown in Figure 5, the bias of estimated annual LE anomalies deviating from ground-based observations is −2.3 W/m 2 , the RMSE is 11.2 W/m 2 and R 2 is 0.42 (p = 0.02).This illustrates that the annual variation of LE is slightly smaller than that expected and perhaps the missing NDVI caused by bad weather can explain this small bias [19,20].In general, MS-PT algorithm works well, indicating that this method may be a good tool for detecting the long-term variation of LE.Rigorous validation of the surface latent heat flux derived from remote sensing data is a challenging scientific problem as soil surface evaporation and plant transpiration involve complicated physical processes The differences in the performance of our MS-PT algorithm among different land cover types or locations are partly caused by errors of ground-observations over flux towers, landscape heterogeneity over flux sites, and the limitations of the MS-PT algorithm.One source that may influence the performance of LE modeling is the flux measurement itself.These flux measurements use the ECOR method to obtain surface radiation fluxes and ECOR suffers energy imbalance.Although we have selected the method proposed by Twine et al. [50] to correct LE, the uncertainty of observed LE still exists.In addition, the quality of MODIS NDVI products also affects the biases of validations.Heterogeneity within the subset around flux towers may influence the correlation between field measurements and remote sensing simulation due to mismatch in spatial representative areas [52,53].MODIS NDVI products are still average values in a given pixel (large area) and often mix with signals of lower vegetation, which will lead to errors in LE modeling among different ecological categories.Another source that can account for the differences between the observed LE from flux towers and the modeled LE is the limitations of the MS-PT algorithm, which neglects the differences of parameters in different biome types.Thus, there may be small biases between the observed and the modeled LE.   [54,55].Therefore, if the fraction of bare soil is miscalculated, the evaporation from the soil fraction based on thermal inertial will be off, as will transpiration from the canopy fraction.It is possible these are compensating errors that reduce overall error.Similarly, for the canopy, it is assumed that the vegetated fraction will transpire at a rate determined by R n with a fixed alpha value, moderated by the vegetation cover and optimal temperature.In fact, stomatal conductance is sensitive to a number of environmental factors and many variables can not be acquired only using remotely sensed data.In this study, we use vegetation cover derived from NDVI to represent the variations in vegetation state and canopy response to changes in environmental conditions such as water availability in the vegetation root zone, the plant water potential, FPAR and CO 2 concentration [56,57].Large error and uncertainty may be introduced by these approximations and assumptions.

Sensitivity Analysis
To test the change in LE from the change in key input variables (R n , NDVI, T a and DT), sensitivity analysis of major parameters for the MS-PT algorithm is also conducted (Figure 6).As illustrated for the MS-PT algorithm, the largest change of LE is caused by the variation of R n .Sensitivity of R n is the highest in the MS-PT algorithm because the Priestley-Taylor equation is calculated as the sum of the surface energy balance term.The averaged LE varies at all flux towers up to ±20% for the MS-PT algorithm by changing R n with ±20%.In response to the change in NDVI with ±20%, LE varies by ±10% at all flux towers.Similarly, LE can increase by 6% for T a change of 20%.However, DT shows different sensitivity with the MS-PT algorithm and LE relatively increases up to 3% for DT change of −20%.Overall, LE estimation by the MS-PT algorithm shows the obvious sensitivity orders: R n > NDVI > T a > DT.In the MS-PT algorithm proposed in this study, R n acquires more highly dependent and this is generally consistent with the previous literature [58,59].For example, Jang et al. [58] considered the sensitivity orders of estimated LE using satellite-based models with R n (±12% change of LE for ±20% change of the variable) > LAI (approximate ±10%) > VPD and T a (less than ±5%) at different land cover types.Hwang and Choi [59] found the sensitivity orders: R n > LAI > P a for the revised remote sensing-based PM model.Besides R n , the Dependency of NDVI in the MS-PT algorithm is higher in the all flux towers as vegetation amount quantified by vegetation index (NDVI) and LAI affect the vegetation photosynthesis and transpiration [18,56].Similarly, DT, reflecting the ATI, is an important factor controlling soil evaporation after surface soil moisture is deficient [60,61].The variations of DT can lead to slight changes in soil evaporation.Compared with R n , NDVI, and T a , DT plays an insignificant role.

Application I: Mapping Terrestrial Evapotranspiration of the Three-North Shelter Forest Region of China
Three-North Shelter Forest Programme of China (TNSFP) is a large ecological restoration project aiming to restore and protect regional vegetation over northwestern, northern, and northeastern regions of China [62,63].Three-North Shelter Forest Region of China (TNSFR) includes more than 550 Chinese counties and covers an area of 4,069,000 km 2 mostly in arid and semi-arid regions [63].Although great progress has been made to monitor the long-term variations of vegetation cover over TNSFR, better understanding of the degree of the variation of LE caused by increasing vegetation cover on decadal scales is still critical due to the deficiency of needed ground-observations data [62][63][64].To this end, the MS-PT algorithm driven by GMAO-MERRA and GIMMIS-NDVI products has been applied to generate monthly LE with at a 0.1-degree spatial resolution in this region during 1982-2009.Here, we do not directly use the GMAO-MERRA LE products as these products exist large uncertainty [5]. Figure 7 displays the map of annual composites of monthly LE for 1982-2009.The largest annual LE occurs in the boreal forest regions of Northeast China, followed by temperate regions of North China and the smallest LE occurs in the arid and semi-arid regions of Northwest China.The highest LE over Northeast China is closely linked to the higher forest cover and healthy forests can maintain long-term LE because forest roots can acquire moisture from deeper soil layers [5,21,45,46,65].The higher LE over North China is mainly caused by large-scale crops.The lowest LE over Northwest China may be attributed to the fact that abnormal precipitation deficiency has raised the likelihood of increased droughts and reduced terrestrial LE.Despite these general similarities for these annual LEs, slight inter-annual variability in spatial LE distributions does exist, driven mainly by differences in climate, hydrology, and vegetation status. Figure 8a demonstrates the spatial pattern of trend in annual LE over Three-North Shelter Forest Region of China.Overall, the actual LE has increased on average over the entire regions and the positive linear slope is about 0.15 W/m 2 per year (p = 0.03) from 1982 through 2009.Regionally, the LE increased over large areas in both Northeast China and Northwest China (except for the Tarim Basin) while decreasing in both North China and the Tarim Basin.The variations of both air temperature and precipitation can explain these differences.In west arid and semi-arid regions, LE is dominated by the precipitation because most of the region is a precipitation-limited environment [5,66].However, in east temperate regions, temperature plays a dominant role in controlling vegetation growth [5,21,67,68].The rising temperature has prolonged the growing season to improve plant productivity and increased LE [63,69].As shown in Figure 8, the spatial variation of the LE and the precipitation trends show the same general trends in west regions of China, while there are good agreement between the spatial variations of the LE and the air temperature trends in east regions of China.These spatial correspondences support the increasing air temperature in the East China and increasing precipitation in the West China, associated with climate warming during the past 30-year period, can be account for the variations of LE over Three-North Shelter Forest Region of China.

Application II: Monitoring Global Land Surface Drought
Drought is the most severe natural disaster causing global environmental changes and has attracted widespread attention from global scientists [70,71].Recently, many satellite-based drought indices (e.g., Temperature-Vegetation Dryness Index, TVDI, and Evaporative Drought Index, EDI), have been widely used to monitor global or regional surface drought [45][46][47]72,73].Particularly for the LE-based EDI, it has integrated the actual LE and potential LE to isolate the effects due to spatially varying soil moisture availability, and has reflected a good correspondence with other meteorological drought indices, but at a significantly higher spatial resolution [32,[45][46][47]73].
To monitor long-term global surface drought, the MS-PT algorithm driven by GEWEX radiation products, NCEP-2 data, and GIMMIS-NDVI products has been applied to generate global monthly LE at a 1-degree spatial resolution during 1984-2007.GEWEX radiation products are used in this study as satellite-based radiation products are the main focus in global surface drought monitoring.Annual standardized anomalies in both EDI and PDSI from 1984 through 2007 are compared in Figure 9. Generally, the EDI successfully reproduces patterns evident in the PDSI, indicating the value of the satellite-derived NDVI and ATI signal as a surface moisture proxy.For example, MS-PT algorithm-based EDI captures the major global surface drought events occurring in 2000, 2002, 2003, 2005, and 2007, even in the Amazon Basin where there is dense vegetation cover and little exposure of the dry soil surface.Similarly, EDI also captures the wet events occurring in 1998 due to the last large El Niño event, which is consistent with the findings of other studies [21,24].For instance, Jung et al. [24] reported that 1998 marks a transition period in which the global land LE trend decreases.The spatial pattern of trends in annual total EDI and PDSI from 1984 through 2007 has been examined (Figure 10).Substantial spatio-temporal variability appears in the dryness trends of EDI derived from GEWEX datasets that is almost consistent with that of PDSI.However, in the Eastern US, Western Australia, and the western regions of South America, the drought trend in EDI differs from the drought variability in PDSI.The missing and contaminated NDVI caused by cloud may lead to retrieval error in the EDI.Detailed speculations on the possible causes of this difference are challenging because drought is affected by variability in local atmospheric conditions (vapor pressure deficit, wind speed, air temperature, and relative humidity), moisture availability (precipitation), radiative forcing (cloud cover and sun angle), and vegetation amount.
Incorporating the optical remote sensing with a higher spatial resolution, we can produce daily time series required to calculate EDI at regional or field scales.For example, Yao et al. [47] generated daily EDI at 4-km resolution for April-September of 2003-2005 across the continental United States and found good performance of EDI in assessing drought at continental scales.These will facilitate drought assessment at field, regional and global scales, which will be valuable for drought monitoring and distribution of drought induced crop yield loss compensation.

Conclusions
We have presented a modified satellite-based Priestley-Taylor (MS-PT) algorithm to estimate terrestrial LE using remote sensing products and meteorological data.This algorithm is physically based, requiring no subjective parameter calibration as employed by many other traditional LE methods.We have also validated it using ground observations collected from 40 flux towers distributed across the world on all continents.The validations show that the RMSE of the estimated daily LE using MS-PT algorithm varies from 10.7 W/m 2 to 87.6 W/m 2 , the bias varies from −23.7 W/m 2 to 48.6 W/m 2 , and R 2 varies from 0.41 to 0.89 (p < 0.01).Similarly, the RMSE of the estimated daily LE using PT-JPL algorithm varies from 11.3 W/m 2 to 89.1 W/m 2 , the bias varies from −21.3 W/m 2 to 56.2 W/m 2 , and R 2 varies from 0.40 to 0.88 (p < 0.01).The average daily LE can be estimated reasonably in terms of the Root Mean Square Error and correlation coefficients.As compared with the PT-JPL algorithm, the MS-PT algorithm improves the LE estimates at most flux towers sites.
To evaluate how well the model predicts long-term variations in LE, we have validated the estimated annual LE anomalies with the observed annual LE anomalies.We only use flux tower sites where five years of data are available.The results illustrate that the bias of estimated annual LE anomalies deviating from ground-based observations is −2.3 W/m 2 , the RMSE is 11.2 W/m 2 and R 2 is 0.42 (p = 0.02).The MS-PT algorithm is also satisfactory in reproducing the inter-annual variability at sites with at least five years of data, indicating this method may be a good tool for analyzing the long-term variation of LE.The sensitivity analysis of major parameters for the MS-PT algorithm demonstrates that the sensitivity orders of estimated LE with R n (±20% change of LE for ±20% change of the variable) > NDVI (approximate ± 10%) > T a (±6%) and DT (±3%) for different land cover types.
This algorithm has been applied in mapping terrestrial LE of Three-North Shelter Forest Region of China and applied in monitoring global land surface drought.The decadal variation in LE of Three-North Shelter Forest region of China, during 1982-2009, illustrates that increasing air temperature in East China and increasing precipitation in West China, associated with climate warming during the past 30-year period, can explain the variations of LE over this region.Additionally, Evaporative Drought Index (EDI) has provided useful surface moisture proxy information without requiring precipitation data for global land surface drought monitoring.MS-PT algorithm-based EDI captures the major global surface drought events occurring in 2000, 2002, 2003, 2005, and 2007, even in the Amazon Basin where there is dense vegetation cover and little exposure of the dry soil surface.Similarly, EDI also captures the wet events occurring in 1998 due to the last large El Niño event, which is consistent with the findings of other studies.MS-PT algorithm-based EDI derived from optical remote sensing with a higher spatial resolution will facilitate drought assessment at field, regional and global scales, which will be valuable for drought monitoring and distribution of drought induced crop yield loss compensation.

Figure 2 .
Figure 2. Location of the 40 flux tower sites used in this study.

Figure 3 .
Figure 3. Eight-day time series comparisons of the modeled LE (daily total) estimates based on two PT algorithms and the ground-measured LE using the data collected from the ten flux towers in their respective land cover classes from the validation tower set.DBF: deciduous broadleaf forest; DNF: Deciduous needleleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; MF: mixed forest; SHR: shrubland.All r values are significant with 99% confidence.

Figure 4 .
Figure 4. Comparison of the predicted and ground-measured annual LE collected at all 40 flux towers sites shown in Table1.

Figure 5 .
Figure 5.Comparison of the annual anomalies of predicted LE and ground-measured LE collected at the flux towers sites where at least five years of data are available.

Figure 6 .
Figure 6.Sensitivity analysis of LE with net radiation, NDVI, DT, and air pressure near surface.

Figure 7 .
Figure 7. Map of annual composites of monthly LE over Three-North Shelter Forest Region of China for 1982-2009.

Figure 8 .
Figure 8. Spatial pattern of linear trends in annual (a) LE based on MS-PT algorithm driven by GMAO-MERRA data and GIMMIS-NDVI products; (b) Precipitation from GMAO-MERRA data; and (c) air temperature from GMAO-MERRA data during 1982-2009.

Figure 9 .
Figure 9. Annual anomalies in global land surface EDI and PDSI for 1984-2007.

Figure 10 .
Figure 10.Spatial pattern of linear trends in annual (a) EDI and (b) PDSI, from 1984 through 2007.

Table 2 .
Statistics of estimated daily LE against the eddy-flux tower observations.The bias and Root Mean Square Error (RMSE) are in units of W/m 2 .All r values are significant with 99% confidence.