Estimating Daily Global Evapotranspiration Using Penman – Monteith Equation and Remotely Sensed Land Surface Temperature

Daily evapotranspiration (ET) is modeled globally for the period 2000–2013 based on the Penman–Monteith equation with radiation and vapor pressures derived using remotely sensed Land Surface Temperature (LST) from the MODerate resolution Imaging Spectroradiometer (MODIS) on the Aqua and Terra satellites. The ET for a given land area is based on four surface conditions: wet/dry and vegetated/non-vegetated. For each, the ET resistance terms are based on land cover, leaf area index (LAI) and literature values. The vegetated/non-vegetated fractions of the land surface are estimated using land cover, LAI, a simplified version of the Beer–Lambert law for describing light transition through vegetation and newly derived light extension coefficients for each MODIS land cover type. The wet/dry fractions of the land surface are nonlinear functions of LST derived humidity calibrated using in-situ ET measurements. Results are compared to in-situ measurements (average of the root mean squared errors and mean absolute errors for 39 sites are 0.81 mm day−1 and 0.59 mm day−1, respectively) and the MODIS ET product, MOD16, (mean bias during 2001–2013 is −0.2 mm day−1). Although the mean global difference between MOD16 and ET estimates is only 0.2 mm day−1, local temperature derived vapor pressures are the likely contributor to differences, especially in energy and water limited regions. The intended application for the presented model is simulating ET based on long-term climate forecasts (e.g., using only minimum, maximum and mean daily or monthly temperatures).


Introduction
Evapotranspiration represents the summation of coupled processes controlling evaporation from impervious, soil, vegetation and open-water surfaces and transpiration from plants.Although it is a major component of the hydrologic cycle, with roughly 60-65% of all precipitation over land being removed from the terrestrial water cycle through ET [1][2][3][4], it is hard to quantify because of the many different factors influencing the many ET processes [5].To further our understanding of ET and support the development and application of process models, in-situ measurement networks encompassing varied landscape, vegetation and meteorological conditions have been operating since 1990's [6] including the AmeriFlux Eddy Covariance Flux Tower network [7] that are used in the present study.Measurements from these networks are generally available from the global FluxNet or AmeriFlux networks [8].To simulate ET, process model input requirements vary from only temperature [9] to radiation [10] to a suite of meteorological, vegetation and soil parameters [11].Additionally, these models generally estimate potential ET (PET) and must be combined with soil moisture models to limit ET based on moisture conditions (i.e., limited water availability) to determine actual ET [12][13][14].
Remote Sens. 2017, 9, 1138 2 of 24 Although significant research has been conducted on ET, much of the modeling focus has been on agricultural applications.For example, the Food and Agriculture Organization (FAO) of the United Nations has published an extensive document to determine crop water requirements based on an array of varying data availability [15].Although the FAO methods can be applied as needed based on data availability, they do not provide global estimates of ET.For global estimates, gridded (i.e., 0.5-1.0degree resolution) ET estimates from fully coupled water-energy balance hydrologic models operating at 3-h to monthly temporal resolutions are freely available from NASA'S Global Land Data Assimilation System (GLDAS) [16] or the National Center for Atmospheric Research's (NCAR) Community Earth System Model (CESM) [17].There is also the relatively new MODIS (Moderate Resolution Imaging Spectroradiometer) ET product, MOD16, which provides ET estimates based on the Penman-Monteith (PM) method using a combination of remotely sensed vegetation measurements and re-analysis derived meteorological inputs [18,19].
Although it is possible to use ET (e.g., estimated from GLDAS or the MODIS ET product) as an independent time series for hydrologic study or as input to hydrologic models focused on event rainfall-runoff dynamics, there are two main challenges: limited ability to couple soil moisture-ET process models, which can impact model calibration, and the application of future climate conditions (e.g., IPCC AR5) for which the above ET estimates are not available due to lack of land surface vegetation characteristics in global scale climate models.The drawback with simulating ET processes within a hydrologic model is the uncertainty resulting from the number of required input parameters for determining available energy, meteorological conditions (humidity, vapor pressure deficit, wind speed, etc.) and resistance terms controlling water transfer.Models that use simulated soil moisture to alter ET resistance terms introduce additional parameter requirements and uncertainties [12].Additionally, most of the ET products are not available on a daily scale, which is needed for hydrological modeling applications.For example, MODIS ET product is only available at 8-day, monthly, and annual scales [19].However, there are some uncertainties and limitations associated with MODIS ET in different land cover types.For example, Liu, et al. [20] investigated eight sites from ChinaFLUX and showed MODIS ET performs better in grasslands as compared to croplands and evergreen forests with a general over-estimation of 6 mm/8-day and RMSE of 11 mm/8-day, where ET errors were linked to LAI estimation errors.Bhattarai, et al. [21] reported that MODIS ET under-estimated ET by 26% as evaluated at 13 AmeriFlux sites inside the United States; with site RMSE ranging from 8.5 to 9.4 mm/8-day.They linked their ET errors to the heterogeneity associated with each MODIS Land Surface Temperature (LST) pixel and the mismatch in spatial scales between MODIS products and ground based measurements.Velpuri,et al. [22] investigated MODIS ET estimations over 58 AmeriFlux sites and reported up to 25% uncertainty in MODIS ET product, suggesting the need for algorithm parameters refinement in the MOD16 product.Long, et al. [23] reported 10-15 mm/month uncertainty in MODIS ET over 3 regions in south central United States, relating the errors to energy and water balance constraints.
While the PM equation (e.g., used in MODIS ET) and associated forcing and parameter estimation methods typically reference air temperature (i.e., temperature measured under shelter at 2 m above the surface [24]), air temperature measurements are available at relatively few meteorological stations around the world providing limited data to capture relevant spatial patterns [24].However, MODIS LST provides global surface temperatures, with no use of reanalysis data, providing consistent information, especially over places where there are either no or very few meteorological stations available.Thus, a central theme of this study is the applicability of MODIS LST for estimating global ET based on the PM equation.
There has been extensive research on the applicability of the MODIS LST.Benali, et al. [24] investigated the accuracy of MODIS LST over 106 meteorological stations in Portugal for a 10 year period and concluded that errors generally fall within 2 to 3 • C. Wan [25] reported that, in clear sky conditions where LST ranges from −10 • C to 58 • C, the accuracy of daily MODIS LST version 5 (V5) is better than 1 • C as compared to the in-situ measurements for most places and the RMSE is less than Remote Sens. 2017, 9, 1138 3 of 24 0.5 • C for all places except for those with apparently heavy aerosol loadings.However, Wang, et al. [26] reported an under-estimation of 2-3 • C during nighttime for daily MODIS LST V4.Zhu, et al. [27] also reported RMSEs of 7.45 • C and 2.97 • C for daytime and nighttime daily LST, respectively, obtained from the Terra satellite and 9.44 • C and 4.41 • C obtained from Aqua satellite, over the Xiangride River basin in the north Tibetan Plateau during 2009 and 2010.Vancutsem, et al. [28] showed that nighttime MODIS LST provides good nighttime temperatures over different ecosystems in Africa, whereas, daytime estimations strongly vary with seasonality, ecosystems, solar radiation and cloud cover.Lin,et al. [29] reported disagreement of MODIS LST estimations and daily measurements in East Africa during the day (MAE = 6.9 ± 5 • C), while suggesting better agreement during the night (MAE = 1.9 ± 1.7 • C).Overall, the accuracy of MODIS LST varies from day to night, seasonally, with land cover, with cloud cover and aerosol loading [25].
Building on MODIS ET concepts, a globally applicable approach using MODIS LST for estimating daily ET is presented and compared to the MOD16 ET product and ET observations at 39 AmeriFlux sites.The novel aspect of this approach is the use of only remotely sensed measurements (i.e., no use of meteorological data, re-analysis data or in-situ measurements).For example, MOD16 ET uses derived meteorological inputs (humidity, air temperature, radiation) from NASA's Global Modeling and Assimilation Office (GMAO, v.4.0.0) re-analysis dataset [19].Here, the MODIS ET approach [18,19,30] is combined with FAO limited data methods for estimating vapor pressure and radiation [31] reducing meteorological input requirements to only temperature and vegetation characteristics, where both quantities are derived from MODIS products.This approach is taken to provide ET for use in large scale hydrologic models (i.e., as ET input or assessment) derived only from remotely sensed measurements (i.e., for regions without in-situ data) and to enable the prediction of future ET using coarse resolution future climate projections produced by Global Climate Models (GCMs).In this study, daily and monthly MODIS LST [25,32] are used and described in Table 1.However, the approach is applicable to any temperature series (e.g., daily or monthly climate model forecasts).For applications using future climate projections, spatial land cover and seasonal vegetation patterns and dynamics must be known or approximated based on current conditions.

Materials and Methods
Central to this study is the Penman-Monteith (PM) method for estimating ET (mm d −1 ) based on the combination of energy and mass transfer concepts [11,15,30,42]: where λ (MJ kg −1 ) is the latent heat of vaporization; ∆ (kPa C −1 ) is the slope of the saturation vapor pressure-temperature; R n (MJ m 2 day −1 ) is net radiation; G (MJ m 2 day −1 ) is the soil heat flux; 86,400 is a unit conversion factor; ρ a (kg m −3 ) is the density of moist air density; C p (MJ kg −1 • C −1 ) is the specific heat capacity of air; e s and e a (kPa) are saturated and actual vapor pressure, respectively; γ (kPa • C −1 ) is the psychometric constant; r s (s m −1 ) is water vaporization resistance from plants, soil or water surfaces; and r a (s m −1 ) is aerodynamic resistance.Here, λ, ∆, ρ a , γ, e s and e a are a function of temperature and/or pressure, with pressure being a function of elevation.The resistance terms (r s & r a ) are largely based on Vegetation Type and LAI [33].For R n , the FAO energy method [15] is used, which results in R n at a given location being influenced only by latitude, Julian day, albedo (α), cloudiness, temperature (T min , T mean , T max ) and elevation (z).The soil heat flux is a minor component in the energy budget (often less than 5%) at daily and coarser scales [43].Given our focus on daily ET, and the soil heat flux is neglected herein.

Datasets
The datasets in Table 1 are used to determine elevation, temperature, albedo as needed to solve the PM equation.Specific details are discussed below.Temperature forcings used in this study are based on four different MODIS derived LST measurements from NASA's Aqua and Terra satellites at roughly 1:30 a.m., 10:30 a.m., 1:30 p.m. and 10:30 p.m. (solar local times) [44].Although the actual acquisition times vary over time, which likely influences the derived temperatures from the two satellites, the MODIS data products for day/night LST are used throughout the studied period.From the four measurements (i.e., day/night from Aqua/Terra), minimum and maximum values are determined, and the mean temperature is determined by averaging the minimum and maximum values.Note, the mean is based on averaging the two temperatures rather than the four to provide consistency for when measurements are available from only one satellite (i.e., two per day or month).Also note that, the level 3 MODIS LST data are used, which masks out bad data due to cloudiness or other reasons.Therefore, our analysis is only based on reliable data with no fill values for the missing days.Thus there is no ET estimation for days with no LST data.
Mean monthly wind speed at 50 m above the surface is used from NASA's Surface Meteorology and Solar Energy (SSE) data product; release 5 dataset (January 2005) consisting of 10-year monthly average (July 1983-June 1993).Wind speed at 50 m is converted to wind speed at 2 m above the surface (based on a method used by Allen,et al. [15]) and only used for estimating open water evaporation.
The datasets are obtained based on their product resolution listed in Table 1.In order to deal with different spatial resolutions, all datasets not gridded with 0.05 degree by 0.05 degree pixels, were either resampled to 0.05 degrees (e.g., Aerosol Optical Depth 1.0 degrees to 0.05 degrees) or upscaled to 0.05 degrees (e.g., LAI 1 km to 0.05 degrees).For resampling, all 0.05 degree pixels within the coarse resolution pixel were assigned the same value as the coarse resolution pixel.For upscaling, all finer resolution pixels within the 0.05 degree resolution pixel were averaged and the mean value was assigned to the 0.05 degree resolution pixel.For temporal resolution, the model was run daily and used the corresponding 8-day or monthly or annual products for each day.

Radiation
Net radiation (R n ) is the difference between the net shortwave radiation (R ns ) and the net longwave radiation (R nl ), where R ns is the fraction of solar shortwave radiation (R s ) not reflected based on surface albedo (α): R ns = (1 − α)R s .Albedo is estimated from black-sky, direct beam, albedo (BSA) and white-sky, completely diffuse, albedo (WSA) based on the fraction of diffuse skylight, S(θ, AOD), which is a function of solar zenith angle (θ) and Aerosol Optical Depth (AOD) [45]: The solar shortwave radiation is determined by: R s = 0.16 which uses the daily temperature range to account for cloudiness without using any cloudiness data [31,46,47], where R a (MJ m 2 day −1 ) is the solar radiation above the atmosphere (extraterrestrial radiation), which is a function of the latitude, Julian day, and solar time (i.e., solar hour angle) [31].R nl is determined by: There are different methods for estimating the resistance terms (r s & r a ) in the PM equation depending on the ET components being considered.Here, four ET water fluxes are represented for a given land area (e.g., 0.05 by 0.05 degree pixel) based on the fraction of vegetated (F c ) or non-vegetated (1 − F c ) land surface and the fraction of wet (F w ) or dry (1 − F w ) land surface Table 2. Studies have shown each component can control ET in different climate conditions.For example Li,et al. [48] showed that the soil water availability in hot-dry seasons and LAI in rainy seasons control ET.The resistance terms for each ET component are estimated based on methods listed in Mu, et al. [19], where each vegetation type has pre-defined resistance terms.Table 3 lists the Biome properties look-up table used here.The values were first derived by Running, et al. [49], and then revised by Heinsch, et al. [50].The values listed in Table 3 are updated (http://files.ntsg.umt.edu/data/NTSG_Products/MOD16/readme) from values listed in Mu, et al. [19].For land cover classified as open-water, evaporation is estimated based on Vallet-Coulomb, et al. [51] Table 2. Total ET for a given location (e.g., 0.05 • grid cell) is the sum of all of ET components proportionally based on the location's land cover fractions.Note, land cover fractions are re-scaled to include only the largest three categories.Transpiration from Vegetation where u (m s −1 ) is wind speed 2 m above surface [51].
* For canopy, vegetation and soil components, ET(r s , r a ) is determined using Equation (1) where the resistance terms vary based on the ET component following Mu, et al. [19] and Table 3.To estimate the four ET components explained in Table 2 for land pixels, F c and F w need to be determined.Here, the wet fraction (F w ) is based on humidity, and the vegetated fraction (F c ) is based on LAI.Each grid cell is divided into wet (F w ) and dry (1 − F w ) fractions.It is important to note that these fraction terms represent indexes used to determine the resistance terms that should be used in Equation (1).Building on Fisher, et al. [52], Mu, et al. [19] used the following equation to determine F w .
saturated vapor pressure at T min [53-55], and RH is calculated by: RH = e s,Tmin (e s,Tmin + e s,Tmax )/2 × 100 (7) where e s,Tmin and e s,Tmax (kPa) are the saturated vapor pressures at T min and T max , respectively, and determined using: Mu, et al. [19] F w estimation method suggests a step at relative humidity of 70%.For some cases, this can create a large jump in ET (Figure 1), which is especially sensitive to temperature derived VPD and RH.As an alternative, we developed a revised approach (Equation ( 9)), based on the assumption of zero wet fraction for humidity conditions below a minimum threshold (RH min ) and a wet fraction of 1 for 100% humidity [52,56].For humidity values from RH min to 100%, the wet fraction varies from 0 to 1 as a linear or nonlinear function depending on the value of the exponent (β).
Values of RH min and β are determined by minimizing the AmeriFlux sites meteorological data driven ET estimate (ET met ) errors at AmeriFlux sites for the period 2000-2006.Similar to surface wetness, each grid cell is divided into vegetated (Fc) and non-vegetated (1 − Fc) fractions.In previous studies, fractionation has been determined using Normalized or Enhanced Vegetation Indices products (NDVI, EVI) [18] or by using the Fraction of Photosynthetic Active Radiation (FPAR) [19].To reduce the input data requirements (e.g., MODIS FPAR or EVI/NDVI products as in Mu, et al. [18] and Mu, et al. [19]) and capitalize on the strong correlation between LAI and FPAR [57][58][59], relationships between LAI and FPAR for each vegetation type [19,60] are developed to approximate Fc.A regression analysis based on a simplified model from the Beer-Lambert law for describing the light transition through various vegetation types [61][62][63] is used to estimate Fc: where the light extinction coefficient (k) can substantially vary for different land covers [64].For example, k ranges from 0.4 to 0.7 for needleleaf forests and 0.5 to 0.8 for broadleaf forests.Aubin, et al. [65] reported a general value of 0.54 for all vegetation types and varying values of 0.37 to 0.98 for different vegetation types.In the present study, vegetation types are determined from MODIS Land Cover (MCD12C1; see Figure 2) and combined with MODIS FPAR (MCD12C1) and LAI (MCD15A2) to determine light extinction coefficients for each land cover type.Similar to surface wetness, each grid cell is divided into vegetated (F c ) and non-vegetated (1 − F c ) fractions.In previous studies, fractionation has been determined using Normalized or Enhanced Vegetation Indices products (NDVI, EVI) [18] or by using the Fraction of Photosynthetic Active Radiation (FPAR) [19].To reduce the input data requirements (e.g., MODIS FPAR or EVI/NDVI products as in Mu, et al. [18] and Mu, et al. [19]) and capitalize on the strong correlation between LAI and FPAR [57][58][59], relationships between LAI and FPAR for each vegetation type [19,60] are developed to approximate F c .A regression analysis based on a simplified model from the Beer-Lambert law for describing the light transition through various vegetation types [61][62][63] is used to estimate F c : where the light extinction coefficient (k) can substantially vary for different land covers [64].For example, k ranges from 0.4 to 0.7 for needleleaf forests and 0.5 to 0.8 for broadleaf forests.Aubin, et al. [65] reported a general value of 0.54 for all vegetation types and varying values of 0.37 to 0.98 for different vegetation types.In the present study, vegetation types are determined from MODIS Land Cover (MCD12C1; see Figure 2) and combined with MODIS FPAR (MCD12C1) and LAI (MCD15A2) to determine light extinction coefficients for each land cover type.

Sheltering Factor
Although transpiration should decrease under water or temperature induced stress, these processes can be difficult to account for without simulating soil moisture or in-situ measurements.
Here, temperature and VPD are used to reduce the potential stomatal conductance in cases of low humidity or cold temperatures to inhibit photosynthesis via mTmin and mVPD multipliers [18,19].However, under some specific conditions (e.g., highly vegetated locations where daily temperature gradients are roughly 10 °C or larger) simulated transpiration results in ET estimates exceeding Priestley-Taylor potential ET (PT ET) (Figure 3).While these conditions may be real, the use of temperature gradient to estimate VPD results in very low humidity, which is not realistic given the large fraction of vegetation.The high daily temperature gradient leads to high VPD, which both decrease the vegetation resistance to transpiration and increase the mass transfer flux.Such conditions provide great potential for transpiration, which especially takes place with introducing a high vegetation density to the model.Figure 3 shows such conditions for which transpiration estimation goes up almost monotonically with LAI leading to high ET.To overcome this challenge, we build on the Dingman [66] approach related to sheltering.Dingman [66] discussed that some leaves transpire at lower rates because they are "sheltered" from the sun and wind and defined a shelter factor (fs) to account for this process.This factor decreases from 1 (open canopy) to 0.5 (closed canopy) as LAI increases [67,68].Based on the mentioned studies, Jang, et al. [69] suggested the below constrain as a third multiplier (i.e., in addition to VPD and Tmin controls) for canopy conductance: Here, the stomatal conductance is multiplied by the shelter factor (fs) in addition to factors for VPD and Tmin to further reduce conductance.As is shown in Figure 3, this factor decreases transpiration estimations.However, for the specific conditions represented in Figure 3, the reduced ET still exceeds potential PT ET as LAI increases above 4.6.Thus, additional research is needed to better understand the sheltering effect on transpiration.

Sheltering Factor
Although transpiration should decrease under water or temperature induced stress, these processes can be difficult to account for without simulating soil moisture or in-situ measurements.Here, temperature and VPD are used to reduce the potential stomatal conductance in cases of low humidity or cold temperatures to inhibit photosynthesis via mT min and mVPD multipliers [18,19].However, under some specific conditions (e.g., highly vegetated locations where daily temperature gradients are roughly 10 • C or larger) simulated transpiration results in ET estimates exceeding Priestley-Taylor potential ET (PT ET) (Figure 3).While these conditions may be real, the use of temperature gradient to estimate VPD results in very low humidity, which is not realistic given the large fraction of vegetation.The high daily temperature gradient leads to high VPD, which both decrease the vegetation resistance to transpiration and increase the mass transfer flux.Such conditions provide great potential for transpiration, which especially takes place with introducing a high vegetation density to the model.Figure 3 shows such conditions for which transpiration estimation goes up almost monotonically with LAI leading to high ET.To overcome this challenge, we build on the Dingman [66] approach related to sheltering.Dingman [66] discussed that some leaves transpire at lower rates because they are "sheltered" from the sun and wind and defined a shelter factor (f s ) to account for this process.This factor decreases from 1 (open canopy) to 0.5 (closed canopy) as LAI increases [67,68].Based on the mentioned studies, Jang, et al. [69] suggested the below constrain as a third multiplier (i.e., in addition to VPD and T min controls) for canopy conductance: Here, the stomatal conductance is multiplied by the shelter factor (f s ) in addition to factors for VPD and T min to further reduce conductance.As is shown in Figure 3, this factor decreases transpiration estimations.However, for the specific conditions represented in Figure 3, the reduced ET still exceeds potential PT ET as LAI increases above 4.6.Thus, additional research is needed to better understand the sheltering effect on transpiration.

Evaluation Method
The above methods enable ET to be estimated globally at daily temporal resolutions with a spatial resolution of 0.05 degrees.To evaluate simulated ET (ETLST), in-situ measurements from the AmeriFlux network are used (ETobs).In order to be consistent with the evaluation approach presented by Mu, et al. [19,39] Ameriflux tower sites (Figure 4) with consistent vegetation types reported between the flux network and MOD12 land cover type 2, and measurements for more than half a year during 2000-2006 are used [19].Following Mu, et al. [19] methods, these in-situ daily measurements are determined based on the Level 4 data consisting of half hour air temperature, VPD, incoming global solar radiation (shortwave radiation, 0.15-4.0μm waveband including both direct radiation and diffuse radiation), and the latent heat flux (LE) at canopy level.To determine daily averages, only days containing 40 or more half hour measurements with quality control flags of "most reliable" were used with no daily fill values.However, note that the number of days at each site reported by Mu, et al. [19] is different than what is used in this study.This might be because the AmeriFlux network now provides more data at some sites and unlike the above study no fill value is used in the present study.
Note that, there is a mismatch in spatial resolution between tower measurements and simulated ETLST.For example, tower data represent a footprint area, which varies with wind speed/direction and measurement heights, while remote sensing measurements and the resulting ETLST estimates represent spatially averaged values over a sizeable landscape.Wang,et al. [70] assessed the representativeness of tower measurements and suggested that every eddy covariance tower footprints have an irregular spatial pattern similar to spatial distribution of prevailing wind direction and the flux sites with zonal vegetation has a better representativeness than site with non-zonal vegetation.Studies have suggested different up-scaling methods (e.g., [71]) to account for spatial mismatches.However, to provide error measures consistent with those reported for the MODIS ET product, the tower measurements (ETobs) are compared to gridded ETLST estimates directly [18,19].Although the tower footprint areas (e.g., ∼100's m 2 to ∼1's km 2 ) are less than the simulated ETLST pixel area (∼25 km 2 ) and there are often energy balance closure problems associated with tower measurements [33,72,73], the tower-based estimates are applied uniformly over the pixel area.Thus, in this study, the flux tower measurements of ETobs (i.e., LE/λ yielding ETobs in units of mm/d for the associated footprint area) are compared directly to simulated ETLST (i.e., ETLST in mm/day for the pixel containing the tower location), for consistency with Mu, et al. [19].Using the available measurements, three ET errors and two correlation-based terms are used to assess model performance:

Evaluation Method
The above methods enable ET to be estimated globally at daily temporal resolutions with a spatial resolution of 0.05 degrees.To evaluate simulated ET (ET LST ), in-situ measurements from the AmeriFlux network are used (ET obs ).In order to be consistent with the evaluation approach presented by Mu, et al. [19,39] Ameriflux tower sites (Figure 4) with consistent vegetation types reported between the flux network and MOD12 land cover type 2, and measurements for more than half a year during 2000-2006 are used [19].Following Mu, et al. [19] methods, these in-situ daily measurements are determined based on the Level 4 data consisting of half hour air temperature, VPD, incoming global solar radiation (shortwave radiation, 0.15-4.0µm waveband including both direct radiation and diffuse radiation), and the latent heat flux (LE) at canopy level.To determine daily averages, only days containing 40 or more half hour measurements with quality control flags of "most reliable" were used with no daily fill values.However, note that the number of days at each site reported by Mu, et al. [19] is different than what is used in this study.This might be because the AmeriFlux network now provides more data at some sites and unlike the above study no fill value is used in the present study.
Note that, there is a mismatch in spatial resolution between tower measurements and simulated ET LST .For example, tower data represent a footprint area, which varies with wind speed/direction and measurement heights, while remote sensing measurements and the resulting ET LST estimates represent spatially averaged values over a sizeable landscape.Wang,et al. [70] assessed the representativeness of tower measurements and suggested that every eddy covariance tower footprints have an irregular spatial pattern similar to spatial distribution of prevailing wind direction and the flux sites with zonal vegetation has a better representativeness than site with non-zonal vegetation.Studies have suggested different up-scaling methods (e.g., [71]) to account for spatial mismatches.However, to provide error measures consistent with those reported for the MODIS ET product, the tower measurements (ET obs ) are compared to gridded ET LST estimates directly [18,19].Although the tower footprint areas (e.g., ~100's m 2 to ~1's km 2 ) are less than the simulated ET LST pixel area (~25 km 2 ) and there are often energy balance closure problems associated with tower measurements [33,72,73], the tower-based estimates are applied uniformly over the pixel area.Thus, in this study, the flux tower measurements of ET obs (i.e., LE/λ yielding ET obs in units of mm/d for the associated footprint area) are compared directly to simulated ET LST (i.e., ET LST in mm/day for the pixel containing the tower location), for consistency with Mu, et al. [19].Using the available measurements, three ET errors and two correlation-based terms are used to assess model performance: where M i and O i are, respectively, ET estimates and observations for day i, and n is the total number of days between 2000 and 2006 that have both ET observations and estimates at a particular AmeriFlux site.In addition, Correlation (R) and Taylor Skill Score (S) [74] are used: where, σ M and σ O are the standard deviations in ET estimates and observations at a given site, and R 0 is the theoretical maximum correlation; assumed to be equal to one at all sites.Unlike RMSE and MAE which can be impacted by the larger errors, R and S are more affected by trends/patterns.The range for R is −1 to 1 and S is 0 to 1 with higher values indicating better agreement between simulated and measured values for both.These four metrics are also used to provide consistency with previous MOD16 assessments [19].
Remote Sens. 2017, 9, 1138 10 of 24 where Mi and Oi are, respectively, ET estimates and observations for day i, and n is the total number of days between 2000 and 2006 that have both ET observations and estimates at a particular AmeriFlux site.In addition, Correlation (R) and Taylor Skill Score (S) [74] are used: where, σM and σO are the standard deviations in ET estimates and observations at a given site, and R0 is the theoretical maximum correlation; assumed to be equal to one at all sites.Unlike RMSE and MAE which can be impacted by the larger errors, R and S are more affected by trends/patterns.The range for R is −1 to 1 and S is 0 to 1 with higher values indicating better agreement between simulated and measured values for both.These four metrics are also used to provide consistency with previous MOD16 assessments [19].

Surface Fractions
First, results are presented for determining the wet and vegetated fractions.For F w , values of RH min and β in Equation ( 9) were optimized.The ET model was implemented with β values ranging from 1 to 10 stepping by 0.1 and RH min values ranging from 0 to 1 stepping by 0.01.As shown in Figure 5a, the lowest RMSE results from an RH min of 0 and β of 4. Figure 5e shows the average AmeriFlux sites meteorological data driven ET estimates (ET met ) for all sites with different RH min and β values.The constant plain represents the average ET obs value of 1.25 mm day −1 measured for all sites during the period of 2000-2006.The intersection of the average measurement and simulated surface is where the model matches the measurement mean.The resulting equation does not generate a step in the wet fraction and eliminates sudden changes in ET (Figure 1).

Surface Fractions
First, results are presented for determining the wet and vegetated fractions.For Fw, values of RHmin and β in Equation ( 9) were optimized.The ET model was implemented with β values ranging from 1 to 10 stepping by 0.1 and RHmin values ranging from 0 to 1 stepping by 0.01.As shown in Figure 5a, the lowest RMSE results from an RHmin of 0 and β of 4. Figure 5e shows the average AmeriFlux sites meteorological data driven ET estimates (ETmet) for all sites with different RHmin and β values.The constant plain represents the average ETobs value of 1.25 mm day −1 measured for all sites during the period of 2000-2006.The intersection of the average measurement and simulated surface is where the model matches the measurement mean.The resulting equation does not generate a step in the wet fraction and eliminates sudden changes in ET (Figure 1).To determine Fc using Equation ( 5), light extinction coefficients (k) for each LC type were determined.Twelve tiles from the MCD15A2 product were manually selected to capture representative vegetation regions distributed throughout the globe while ensuring all MODIS vegetation types were included for a large number of pixels.For sub-sampling within each tile, FPAR and LAI values (from MCD15A2 product) from one thousand randomly selected points for each LC type were selected.To capture seasonal and annual variability, points were selected from: January, April, July, and October, in the years 2003, 2008 and 2013.Note that, the FPAR/LAI correlation varies little from year to year.The main differences occur within a given year between the different seasons.For example, LAI increases dramatically during the spring growing season in North America.For each LC, k values were determined using all the selected points and Equation (5).To account for point selection uncertainty, the above process was performed three separate times (i.e., different randomly To determine F c using Equation ( 5), light extinction coefficients (k) for each LC type were determined.Twelve tiles from the MCD15A2 product were manually selected to capture representative vegetation regions distributed throughout the globe while ensuring all MODIS vegetation types were included for a large number of pixels.For sub-sampling within each tile, FPAR and LAI values (from MCD15A2 product) from one thousand randomly selected points for each LC type were selected.To capture seasonal and annual variability, points were selected from: January, April, July, and October, in the years 2003, 2008 and 2013.Note that, the FPAR/LAI correlation varies little from year to year.The main differences occur within a given year between the different seasons.For example, LAI increases dramatically during the spring growing season in North America.For each LC, k values were determined using all the selected points and Equation (5).To account for point selection uncertainty, the above process was performed three separate times (i.e., different randomly selected points) and the final k values were determined by averaging the results from each scenario.Figure 6 shows the logarithmic relationship between LAI and FPAR from a set of points extracted from the Evergreen Broadleaf Forest LC type (EBF).The slope of the line is the light extinction coefficient.Results for other LC's are similar with a minimum R 2 value of 0.91.With k estimated, Figure 7 shows a comparison between LAI derived F c using Equation ( 5) and FPAR values for one MCD15A2 tile (h10v05) for all selected time periods.The k values based on the above analysis for all MODIS LC types (Table 3) range from 0.46 to 0.67 and are consistent to values reported in previous studies [64,65].
selected points) and the final k values were determined by averaging the results from each scenario.Figure 6 shows the logarithmic relationship between LAI and FPAR from a set of points extracted from the Evergreen Broadleaf Forest LC type (EBF).The slope of the line is the light extinction coefficient.Results for other LC's are similar with a minimum R 2 value of 0.91.With k estimated, Figure 7 shows a comparison between LAI derived Fc using Equation (5) and FPAR values for one MCD15A2 tile (h10v05) for all selected time periods.The k values based on the above analysis for all MODIS LC types (Table 3) range from 0.46 to 0.67 and are consistent to values reported in previous studies [64,65].

ET Estimations
Daily ET (ETLST,d, ETLST,m) based on daily or mean monthly LST, 16-day albedo, 8-day LAI and annual land cover was estimated for the period of 2000-2006.Note, when monthly temperatures are used, ET is still estimated daily using the same temperatures for each day of the month.In general, daily ET using monthly data is similar for 8-days based on the other MODIS values.However, there are slight 8-day changes in estimated shortwave radiation and LAI.Table 4 shows daily and monthly validations including minimum, maximum, median, and mean ETLST with standard errors, average ETobs and performance metrics over the 39 AmeriFlux sites used in this study (Figure 4).Note that the error measures are conducted based on daily estimates.The spatial resolution of ETLST was 0.05 selected points) and the final k values were determined by averaging the results from each scenario.Figure 6 shows the logarithmic relationship between LAI and FPAR from a set of points extracted from the Evergreen Broadleaf Forest LC type (EBF).The slope of the line is the light extinction coefficient.Results for other LC's are similar with a minimum R 2 value of 0.91.With k estimated, Figure 7 shows a comparison between LAI derived Fc using Equation ( 5) and FPAR values for one MCD15A2 tile (h10v05) for all selected time periods.The k values based on the above analysis for all MODIS LC types (Table 3) range from 0.46 to 0.67 and are consistent to values reported in previous studies [64,65].

ET Estimations
Daily ET (ETLST,d, ETLST,m) based on daily or mean monthly LST, 16-day albedo, 8-day LAI and annual land cover was estimated for the period of 2000-2006.Note, when monthly temperatures are used, ET is still estimated daily using the same temperatures for each day of the month.In general, daily ET using monthly data is similar for 8-days based on the other MODIS values.However, there are slight 8-day changes in estimated shortwave radiation and LAI.Table 4 shows daily and monthly validations including minimum, maximum, median, and mean ETLST with standard errors, average ETobs and performance metrics over the 39 AmeriFlux sites used in this study (Figure 4).Note that the error measures are conducted based on daily estimates.The spatial resolution of ETLST was 0.05

ET Estimations
Daily ET (ET LST,d , ET LST,m ) based on daily or mean monthly LST, 16-day albedo, 8-day LAI and annual land cover was estimated for the period of 2000-2006.Note, when monthly temperatures are used, ET is still estimated daily using the same temperatures for each day of the month.In general, daily ET using monthly data is similar for 8-days based on the other MODIS values.However, there are slight 8-day changes in estimated shortwave radiation and LAI.Table 4 shows daily and monthly validations including minimum, maximum, median, and mean ET LST with standard errors, average ET obs and performance metrics over the 39 AmeriFlux sites used in this study (Figure 4).Note that the error measures are conducted based on daily estimates.The spatial resolution of ET LST was 0.05 and 0.25 degrees, where all required input parameters were spatially averaged from their native resolution to 0.05 or 0.25 degrees.Table 4 shows that our approach using MODIS LST results in an average bias of −0.19 mm day −1 (ET LST,d ) and −0.09 mm day −1 (ET LST,m ), RMSE of 0.81 mm day −1 (ET LST,d ) and 0.77 mm day −1 (ET LST,m ), and MAE of 0.59 mm day −1 (ET LST,d ) and 0.56 mm day −1 (ET LST,m ).Correlations are 0.68 and 0.72 and Taylor Skill Scores are 0.74 and 0.7, for ET LST,d and ET LST,m , respectively.Mu, et al. [19], who used meteorological re-analysis data, reported an average RMSE of 0.87 mm day −1 , MAE of 0.27 mm day −1 , correlation of 0.56 and Taylor score of 0.54, over the same 39 sites.Here, RMSE, R and S metrics are improved but MAE is worse.However, it is important to note that numbers of days with available ET LST estimates at each site are slightly different in this study as compared to those reported in Mu, et al. [19] as discussed above, which may also contribute to differences in performance measures.In addition, Mu, et al. [19] reported MAE values much lower than RMSE.For example, at a specific site (Tonzi Ranch), MAE is reported to be nearly perfect (0.01 mm day −1 ) with a RMSE of 0.67 mm day −1 .While the individual error measure seems reasonable, the ratio of MAE to RMSE is not typical.For example, typical ratios of MAE/RMSE is roughly 0.7 to 0.8 [75].In this study, the ratio of mean MAE/RMSE is 0.73 (for both ET LST,d and ET LST,m ), which is consistent with typically reported ratios.To assess the impacts of using temperature derived energy and VPD method, daily ET was simulated at each site using the in-situ metrological measurements described previously (ET met ).Overall, the mean ET met from all sites for the period 2000-2006 is 1.23 mm day −1 using AmeriFlux tower input data, which is consistent with the mean of the ET obs of 1.25 mm day −1 for the same sites and period (average bias = −0.02mm day −1 ).Results show that RMSE (0.69 mm day −1 ), MAE (0.51 mm day −1 ), R (0.76) and S (0.8) are all improved but remain similar in magnitude to ET LST performance metrics.For comparison, RMSE, MAE, R and S are 0.82 mm day −1 , 0.30 mm day −1 , 0.62 and 0.56, respectively, in Mu, et al. [19].Thus, the revised methods (i.e., not the meteorological inputs) used in this study compares reasonably to MOD16 methods, which uses more inputs.More importantly, the similar comparison between performance metrics for ET LST and ET met suggests that our methods using LST derived terms are reasonable.
To illustrate the variability of individual site/day performance, Figure 8 shows the comparison between daily ET LST and ET obs over the study period for all 39 AmeriFlux sites (Figure 4), where the 29,800 individual daily ET pairs are grouped by land cover type.The Figure shows that model performance varies based on land cover.Although the large number of points represented in each plot make it difficult to assess mean performance, low ET LST rates tend to be consistently overestimated at Evergreen Broadleaf Forest (EBF) sites (i.e., tropical rain forests) while the highest ET rates tend to be underestimated for Crops and Evergreen Needle Forest (ENF) site.However, for the majority of ET LST (<0.94 mm day −1 ), estimates tend to mirror the 1:1 line with the exception of EBF.The range of variability along the 1:1 line also varies with land cover, with Woody Savannas (WL) sites have the greatest spread.
In terms of land cover and mean model performance, when using AmeriFlux site meteorological data, the model performs best for sites characterized as Evergreen Needleleaf Forests (ENF) with a mean bias of −0.01 mm day −1 and an average RMSE of 0.58 mm day −1 , and worst for sites characterized as Evergreen Broadleaf Forests (EBF) with a mean bias of −0.48 mm day −1 and an average RMSE of 0.78 mm day −1 .When using daily MODIS LST (Figure 8), the model performs best for ENF (RMSE = 0.65 mm day −1 ) and Savannas (SV) (RMSE = 0.61 mm day −1 ) but performs worst for Evergreen Broadleaf Forests (EBF) (RMSE = 1.28 mm day −1 ) and Croplands (RMSE = 1.08 mm day −1 ).In both cases, there is only a minimal change in model performance when using daily or mean monthly LST.The poor performance for EBF is consistent with the high LST errors that will be discussed later.The poor performance for cropland is likely due to transpiration resistance terms needing to be crop specific and the potential for soil moisture stress during the growing season.
Remote Sens. 2017, 9, 1138 14 of 24 performance varies based on land cover.Although the large number of points represented in each plot make it difficult to assess mean performance, low ETLST rates tend to be consistently overestimated at Evergreen Broadleaf Forest (EBF) sites (i.e., tropical rain forests) while the highest ET rates tend to be underestimated for Crops and Evergreen Needle Forest (ENF) site.However, for the majority of ETLST (<0.94 mm day −1 ), estimates tend to mirror the 1:1 line with the exception of EBF.The range of variability along the 1:1 line also varies with land cover, with Woody Savannas (WL) sites have the greatest spread.
In terms of land cover and mean model performance, when using AmeriFlux site meteorological data, the model performs best for sites characterized as Evergreen Needleleaf Forests (ENF) with a mean bias of −0.01 mm day −1 and an average RMSE of 0.58 mm day −1 , and worst for sites characterized as Evergreen Broadleaf Forests (EBF) with a mean bias of −0.48 mm day −1 and an average RMSE of 0.78 mm day −1 .When using daily MODIS LST (Figure 8), the model performs best for ENF (RMSE = 0.65 mm day −1 ) and Savannas (SV) (RMSE = 0.61 mm day −1 ) but performs worst for Evergreen Broadleaf Forests (EBF) (RMSE = 1.28 mm day −1 ) and Croplands (RMSE = 1.08 mm day −1 ).In both cases, there is only a minimal change in model performance when using daily or mean monthly LST.The poor performance for EBF is consistent with the high LST errors that will be discussed later.The poor performance for cropland is likely due to transpiration resistance terms needing to be crop specific and the potential for soil moisture stress during the growing season.

Uncertainties in MODIS LST, LAI, Albedo and In-Situ Measurements
The model used in this study estimates daily ET without meteorological data such as humidity, radiation or cloud cover.Model forcings such as land cover (annual), albedo (16-day) and LAI (8-day) are assumed constant between measurements.The only input that varies daily is temperature.Building on the strong links between temperature and ET [76], temperature is used to determine net shortwave radiation (R ns ) by estimating cloud cover from the daily temperature range and is a key variable in calculating the net outgoing longwave radiation (R nl ).Additionally, temperature is used to estimate vegetation and soil resistance terms in the PM equation [19].The primary role of temperature is in determining vapor pressure (humidity and VPD), which is used to determine the wet surface fraction.
We analyzed LST errors (for daily T min and T max ) over 2000 to 2006 (i.e., the period investigated in this study) over all the sites and through different land covers.Other than Ever Green Broadleaf Forests (EBF), the correlation coefficients between remotely sensed temperature and in-situ measurements range from 0.85 to 0.92 for T min and 0.82 to 0.95 for T max .For EBF this value is 0.14 for T min and 0.22 for T max .Mean daily MODIS LST over the studied AmeriFlux sites is on average 0.3 • C lower than the measurements with an RMSE of 4.67 • C and MAE of 3.28 • C. Regardless of these challenges, MODIS LST provides global, consistent temperatures and is used herein.Additional discussion is provided below to highlight potential ET errors that likely results from day time LST uncertainty.
In addition to LST, MODIS LAI and Albedo products uncertainties are also included.Wang, et al. [77] indicated that RMSEs for the MCD43A Albedo product is less than 0.05 over grasslands, croplands and forests during dormant and snow-covered periods.Cescatti,et al. [78] found albedo MAEs of less than 0.02 over different land covers across the global FLUXNET network.For LAI, Yang, et al. [79] suggested a RMSE of 0.66 for the MODIS LAI product.Two recent studies [80,81] reported RMSE of 0.8 for the MCD15A2 V5 product, while Ganguly, et al. [82] suggests a RMSE of 0.5 LAI for different biomes.
To account for uncertainty in LST, LAI, Albedo and F c in ET estimates, we build on the above studies to estimate uncertainties for each product.The ET algorithm is executed 100 times each time-step (daily over 2001-2013) for each pixel (globally), where values for LST, LAI, Albedo and estimated F c are assigned a zero mean Gaussian random error with a standard deviation of 3 • C for LST, 0.7 for LAI, 0.05 for Albedo and ranging from 0.03 to 0.07 based on LC for F c .For each time-step and each pixel the corresponding mean and standard deviation are determined [83].The results of this analysis are provided in Table 4 as the mean ± the standard error of the mean, averaged over the period of 2000-2006 for all AmeriFlux sites used in this study.
Another source of uncertainty is the in-situ measurements.Studies have suggested that ET measurements at Flux towers have errors of about 10-30% [19,84].Also, the scale of our model (i.e., 0.05 or 0.25 degrees) is much coarser than the footprint of the eddy towers that have a horizontal scale of 2-5 m [18,19].To overcome this discrepancy, especially at sites with heterogeneous environments, an upscaling process is recommended as a future work.Here, the uncertainty of the in-situ measurements is assumed to be much smaller in magnitude as compared to the remotely sensed measurements and not included in our analysis.

Sources of Errors in Estimations
Given our use of MODIS derived LST and VPD and their importance, it is likely that ET LST errors are linked to LST and VPD errors, which is consistent with other studies [85,86].For example, Figure 9 shows relationships between daily errors in T max , VPD and ET LST for AmeriFlux sites having grassland (Grass) and evergreen broadleaf forest (EBF) land cover types.As illustrated in Figure 9, our approach tends to under-estimate ET LST for LCs with relatively low vegetation cover and low LAI (i.e., grasslands, croplands and savannas).In these locations, MODIS LST based T max is often higher than the AmeriFlux site measurements, which leads to an over-estimation of VPD (Equation ( 3)).
In low LAI conditions, higher VPD results in lower surface wetness (i.e., smaller F w ) and increased water stress, which all leads to lower ET LST .The lower ET LST results from decreasing the fraction of high efficiency canopy evaporation from the fraction of saturated vegetation [F c F w ] and decreasing transpiration from unsaturated vegetation [F c (1 − F w )] by increasing resistance.In contrast, this is often not the case for the forests that have relatively high vegetation density (high LAI) such as evergreen broadleaf forests (in Amazon), deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests.For these LCs, over-estimation of VPD usually does not lead to an under-estimation of ET LST because the surface wetness and soil moisture do not limit ET (i.e., Energy limited).Thus, higher VPD leads to higher ET.In other words, higher VPD or lower RH makes it easier for water flux to occur because there is less moisture in the air.This is consistent with the argument that in deciduous forests ET is generally energy limited and not water limited.In highly vegetated regions with high water availability, an over-estimation of VPD does not limit ET, and "actual" ET approaches "potential" ET.
Remote Sens. 2017, 9, 1138 16 of 24 low LAI conditions, higher VPD results in lower surface wetness (i.e., smaller Fw) and increased water stress, which all leads to lower ETLST.The lower ETLST results from decreasing the fraction of high efficiency canopy evaporation from the fraction of saturated vegetation [Fc Fw] and decreasing transpiration from unsaturated vegetation [Fc (1 − Fw)] by increasing resistance.In contrast, this is often not the case for the forests that have relatively high vegetation density (high LAI) such as evergreen broadleaf forests (in Amazon), deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests.For these LCs, over-estimation of VPD usually does not lead to an underestimation of ETLST because the surface wetness and soil moisture do not limit ET (i.e., Energy limited).Thus, higher VPD leads to higher ET.In other words, higher VPD or lower RH makes it easier for water flux to occur because there is less moisture in the air.This is consistent with the argument that in deciduous forests ET is generally energy limited and not water limited.In highly vegetated regions with high water availability, an over-estimation of VPD does not limit ET, and "actual" ET approaches "potential" ET.Building on the above for Grasslands, it is more probable that an over (under) estimation of Tmax or VPD leads to under (over) estimation of ETLST.Thus, there is an "inverse" relationship between Tmax and VPD errors with ETLST error.In contrast, for the EBF LC type, there is a "direct" relationship.To further investigate the effects of vegetation density, relationships between the probability of VPD and ETLST positive/negative errors and LAI were investigated.Figure 10 shows these relationships for the period 2000-2006 from sites representing different land covers.In general, the higher the LAI value the more probable the error relationships are direct, which is consistent with our above discussion.For example, when LAI is 4.5, 65% of the time in which VPD is over-estimated, ETLST is also over-estimated, while for an LAI of 0.1 this occurs only 33% of the time (i.e., 67% of the time in which VPD is over-estimated, ETLST is under-estimated).The key point of this discussion is that VPD error is likely a major contributor to ETLST error.Building on the above for Grasslands, it is more probable that an over (under) estimation of T max or VPD leads to under (over) estimation of ET LST .Thus, there is an "inverse" relationship between T max and VPD errors with ET LST error.In contrast, for the EBF LC type, there is a "direct" relationship.To further investigate the effects of vegetation density, relationships between the probability of VPD and ET LST positive/negative errors and LAI were investigated.Figure 10 shows these relationships for the period 2000-2006 from sites representing different land covers.In general, the higher the LAI value the more probable the error relationships are direct, which is consistent with our above discussion.For example, when LAI is 4.5, 65% of the time in which VPD is over-estimated, ET LST is also over-estimated, while for an LAI of 0.1 this occurs only 33% of the time (i.e., 67% of the time in which VPD is over-estimated, ET LST is under-estimated).The key point of this discussion is that VPD error is likely a major contributor to ET LST error.Another source of error is radiation.The method we use to account for cloudiness in estimating shortwave radiation (Rns) links the LST errors to radiation (Equation ( 3)).As expected, our Rns estimation method performs worst in EBF (R = 0) because of high LST errors in this land cover, while in other land covers it performs reasonably well with correlation coefficients range from 0.77 to 0.92.

Global Estimations
Monthly ETLST is estimated globally (based on monthly MODIS LST grids, which have far less missing values than daily LST grids) for the period 2001-2013 both at 0.05 and 0.25 degree resolutions.Figure 11 shows mean monthly ETLST corresponding standard errors and the difference between our approach and MOD16 ET product for the same period.In general, ETLST is similar to MOD16 with an average difference of only 0.2 mm day −1 less than MODIS ET (e.g., mean ETLST for the near-global MODIS extent is 1.2 ± 0.2 mm day −1 as compared to the mean MOD16 ET of 1.4 mm day −1 ).Although the global patterns are similar, there are specific regions near the boundaries of the Congo and Amazon rainforests where ETLST is much higher than MOD16 (e.g., 0.5-2.0mm day −1 ).These overestimations generally result from transpiration over-estimations (i.e., regions on the boundary of tropical rainforests with high temperature ranges and considerable vegetation but lower rainfall than neighboring rainforests).Although the sheltering factor discussed previously, decreased transpiration over-estimations considerably, it cannot over-come the problem completely and future research is needed to resolve the problem (i.e., soil moisture stress function).In contrast, ETLST is lower than MOD16 (e.g., 0.5-2 mm day −1 ) in Amazon (i.e., Energy limited regions).Table 5 shows summary statistics for the ETLST comparison based on LC types.Compared to MOD16, the largest difference occurs in the Evergreen Broadleaf Forests, where estimates are on average 0.97 mm day −1 lower than MODIS.As discussed above, the problem in EBF is likely related to temperature ranges and resulting VPD, where lower VPD leads to lower ETLST.In these regions, energy is typically the limiting factor and VPD should not limit ETLST.Another source of error is radiation.The method we use to account for cloudiness in estimating shortwave radiation (R ns ) links the LST errors to radiation (Equation ( 3)).As expected, our R ns estimation method performs worst in EBF (R = 0) because of high LST errors in this land cover, while in other land covers it performs reasonably well with correlation coefficients range from 0.77 to 0.92.

Global Estimations
Monthly ET LST is estimated globally (based on monthly MODIS LST grids, which have far less missing values than daily LST grids) for the period 2001-2013 both at 0.05 and 0.25 degree resolutions.Figure 11 shows mean monthly ET LST corresponding standard errors and the difference between our approach and MOD16 ET product for the same period.In general, ET LST is similar to MOD16 with an average difference of only 0.2 mm day −1 less than MODIS ET (e.g., mean ET LST for the near-global MODIS extent is 1.2 ± 0.2 mm day −1 as compared to the mean MOD16 ET of 1.4 mm day −1 ).Although the global patterns are similar, there are specific regions near the boundaries of the Congo and Amazon rainforests where ET LST is much higher than MOD16 (e.g., 0.5-2.0mm day −1 ).These over-estimations generally result from transpiration over-estimations (i.e., regions on the boundary of tropical rainforests with high temperature ranges and considerable vegetation but lower rainfall than neighboring rainforests).Although the sheltering factor discussed previously, decreased transpiration over-estimations considerably, it cannot over-come the problem completely and future research is needed to resolve the problem (i.e., soil moisture stress function).In contrast, ET LST is lower than MOD16 (e.g., 0.5-2 mm day −1 ) in Amazon (i.e., Energy limited regions).Table 5 shows summary statistics for the ET LST comparison based on LC types.Compared to MOD16, the largest difference occurs in the Evergreen Broadleaf Forests, where estimates are on average 0.97 mm day −1 lower than MODIS.As discussed above, the problem in EBF is likely related to temperature ranges and resulting VPD, where lower VPD leads to lower ET LST .In these regions, energy is typically the limiting factor and VPD should not limit ET LST .Although not shown here, the effects of spatial resolution were also investigated.ETLST was estimated at both 0.05 and 0.25 degrees where all input datasets where spatially averaged to each respective resolution.Overall, there were minimal differences between results at the AmeriFlux sites based on 0.05 and 0.25 degree resolutions.The average RMSE and MAE difference over the AmeriFlux sites between the two resolutions is less than 0.02 and 0.01 mm day −1 , respectively.The only noticeable differences were along shorelines where the coarse resolution included open water, which results into high ETLST estimates.Although not shown here, the effects of spatial resolution were also investigated.ET LST was estimated at both 0.05 and 0.25 degrees where all input datasets where spatially averaged to each respective resolution.Overall, there were minimal differences between results at the AmeriFlux sites based on 0.05 and 0.25 degree resolutions.The average RMSE and MAE difference over the AmeriFlux sites between the two resolutions is less than 0.02 and 0.01 mm day −1 , respectively.The only noticeable differences were along shorelines where the coarse resolution included open water, which results into high ET LST estimates.

Conclusions
In this study, an approach was developed to simulate global ET using the Penman−Monteith equation and data products derived from only remotely sensed measurements (i.e., no in-situ or reanalysis datasets).Key components of the approach are new methods for estimating the fractions of vegetated and wet surfaces.New relationships are developed to estimate the vegetated fraction of the land areas based on LAI, MODIS land cover types and light extinction constants.A revised relationship for estimating the wet fraction of the land surface based on temperature-derived humidity is also presented.
Comparing ET performance metrics to those reported in Mu, et al. [19] suggests that the presented approach is comparable to MODIS ET while only using remotely sensed measurements.For example the presented approach results in an average RMSE of 0.81 mm day −1 for ET LST compared to a RMSE of 0.87 mm day −1 reported by Mu, et al. [19] over the same AmeriFlux sites.When using the AmeriFlux meteorological datasets as inputs to the model, the presented approach results in an average RMSE of 0.69 mm day −1 for ET met compared to a RMSE of 0.82 mm day −1 reported by Mu,et al. [19].This study also highlights the importance of LST on ET estimates.For example, when maximum LST is over-estimated, our approach over-estimates VPD (i.e., lower humidity) resulting in cascading effects on mass transfer and wet surface contributions and resistance terms that impact ET estimates differently based on vegetation conditions.Over-estimates of maximum LST can lead to under-estimates of ET in sparsely vegetated areas (i.e., water limited regions) due to increased resistance controls, but over-estimates of ET in highly vegetated areas (i.e., energy limited regions) due to increased mass transfer contributions.Additional research is needed to reduce model sensitivity to over-estimations of maximum LST (i.e., humidity and VPD) in both water and energy limited conditions.

Figure 1 .
Figure 1.Variations in total and individual ET components as a function of relative humidity for select conditions (minimum and maximum temperatures of 15 °C and 30 °C, LAI of 4, albedo of 0.19, sea level elevation, 30 degrees latitude, EBF land cover type and time of year is mid-July): MODIS ET wet fraction method (Blue color), and wet fraction determined using Equation (9) with an β of 4 and RHmin of 0 (Red color).

Figure 1 .
Figure 1.Variations in total and individual ET components as a function of relative humidity for select conditions (minimum and maximum temperatures of 15 • C and 30 • C, LAI of 4, albedo of 0.19, sea level elevation, 30 degrees latitude, EBF land cover type and time of year is mid-July): MODIS ET wet fraction method (Blue color), and wet fraction determined using Equation (9) with an β of 4 and RH min of 0 (Red color).

Figure 2 .
Figure 2. 2001 Land Cover (MCD12C1) with shaded squares showing the selected FPAR/LAI tiles used in this study.

Figure 2 .
Figure 2. 2001 Land Cover (MCD12C1) with shaded squares showing the selected FPAR/LAI tiles used in this study.

Figure 3 .
Figure 3. Variations in total ET with LAI with and without the shelter factor (fs) to limit transpiration component of ET for a select location along with potential Priestley-Taylor ET (30 degrees latitude, 100 m elevation, closed shrublands land cover, albedo of 0.19, daily minimum maximum temperatures of 15 °C and 30 °C, time of year is mid-July).

Figure 3 .
Figure 3. Variations in total ET with LAI with and without the shelter factor (f s ) to limit transpiration component of ET for a select location along with potential Priestley-Taylor ET (30 degrees latitude, 100 m elevation, closed shrublands land cover, albedo of 0.19, daily minimum maximum temperatures of 15 • C and 30 • C, time of year is mid-July).

Figure 4 .
Figure 4. Locations of the AmeriFlux Eddy Covariance tower sites using in this study.Figure 4. Locations of the AmeriFlux Eddy Covariance tower sites using in this study.

Figure 4 .
Figure 4. Locations of the AmeriFlux Eddy Covariance tower sites using in this study.Figure 4. Locations of the AmeriFlux Eddy Covariance tower sites using in this study.

Figure 5 .
Figure 5. Relationships between mean error metrics (a-d) over the period of 2000-2006 from all study sites with the constants (RHmin and β) in Equation (4); (e) comparison between mean daily ETmet (shaded surface) and mean daily ETobs (plain surface) over the period of 2000-2006 for all study sites.

Figure 5 .
Figure 5. Relationships between mean error metrics (a-d) over the period of 2000-2006 from all study sites with the constants (RHmin and β) in Equation (4); (e) comparison between mean daily ET met (shaded surface) and mean daily ET obs (plain surface) over the period of 2000-2006 for all study sites.

Figure 6 .
Figure 6.Relationship between LAI and FPAR for points extracted from Evergreen Broadleaf Forests (EBF) land cover type, where the slope of the regression line represents the light extinction coefficient (k).

Figure 6 .
Figure 6.Relationship between LAI and FPAR for points extracted from Evergreen Broadleaf Forests (EBF) land cover type, where the slope of the regression line represents the light extinction coefficient (k).

Figure 6 .
Figure 6.Relationship between LAI and FPAR for points extracted from Evergreen Broadleaf Forests (EBF) land cover type, where the slope of the regression line represents the light extinction coefficient (k).

Figure 7 .
Figure 7.Comparison between LAI derived F c values determined using Equation (5) and FPAR values of a selected MCD15A2 tile (h10v05).

Figure 8 .
Figure 8.Comparison between daily ETobs and ETLST along with mean ETobs and ETLST and error measures for all the AmeriFlux sites used in this study over the period of 2000-2006 categorized into different land covers.

Figure 8 .
Figure 8.Comparison between daily ET obs and ET LST along with mean ET obs and ET LST and error measures for all the AmeriFlux sites used in this study over the period of 2000-2006 categorized into different land covers.

Figure 9 .
Figure 9. Days with different daily Tmax (a,c) and VPD (b,d) errors versus ETLST errors for grasslands (Grass) and Evergreen Broadleaf Forests (EBF) over the period of 2000-2006.Errors are daily estimates minus measurements.

Figure 9 .
Figure 9. Days with different daily T max (a,c) and VPD (b,d) errors versus ET LST errors for grasslands (Grass) and Evergreen Broadleaf Forests (EBF) over the period of 2000-2006.Errors are daily estimates minus measurements.

Figure 10 .
Figure 10.Probability of over-(black lines) or under-estimation (red lines) of ETLST with corresponding over-(solid lines) or under-estimation (dashed lines) of VPD as a function of daily LAI at AmeriFlux sites; note, black solid and red dashed lines represent direct relationships between ETLST and VPD errors while black dashed and red solids lines represent inverse relationships.

Figure 10 .
Figure 10.Probability of over-(black lines) or under-estimation (red lines) of ET LST with corresponding over-(solid lines) or under-estimation (dashed lines) of VPD as a function of daily LAI at AmeriFlux sites; note, black solid and red dashed lines represent direct relationships between ET LST and VPD errors while black dashed and red solids lines represent inverse relationships.

Figure 11 .
Figure 11.Mean annual ETLST estimates over the period of 2001-2013 (top) and corresponding standard error estimates (middle), and difference between top and mean annual MODIS ET MOD16 [19] for the same period (bottom).

Figure 11 .
Figure 11.Mean annual ET LST estimates over the period of 2001-2013 (top) and corresponding standard error estimates (middle), and difference between top and mean annual MODIS ET MOD16 [19] for the same period (bottom).

Table 1 .
Descriptions of key datasets used in this study.

Table 2 .
Summary of ET component equations.

Table 3 .
[19]e specific parameter look-up table for Evergreen Needleleaf Forest (ENF), Evergreen T min_open , T min_close , VPD open , VPD close determine boundary conditions for mT min and mVPD multipliers, which control transpiration, gl_sh and gl_e_wv are leaf conductance terms and Cl is a stomatal conductance term, RBL min and RBL max are boundary conditions that control soil resistance terms.See Mu, et al.[19]for more details.Light extinction coefficient (k) will be explained later.

Table 4 .
Summary of results and statistics at the AmeriFlux sites used in this study, where the min, max, median and mean values correspond to the mean ET (mm day −1 ), RMSE (mm day −1 ), MAE (mm day −1 ), R and S, based on daily ET estimates at each individual site using site measurements or MODIS datasets at a spatial resolution of 0.05 degrees; for MODIS derived values, mean ET LST is the average of the daily means ± the average of the standard error (se) of the daily means (discussed later) from each site.

Table 5 .
Summary statistics for the difference (mm day −1 ) between ET LST estimates and MODIS ET averaged over the period 2001-2013 based on 2001 land cover.