Evaluation of a MODIS triangle-based evapotranspiration algorithm for semi-arid regions

Abstract The current study investigates a MODIS-based remote-sensing-based algorithm for the assessment of daily latent heat flux (LE), or evapotranspiration (ET), with a focus on semi-arid regions. The approach uses the triangle method, relying on remotely sensed inputs and a previously developed net radiation model. A major difference from previous studies is that we utilize only MODIS products for estimation of ET at the daily timestep for clear and cloudy days. The algorithm is evaluated at four flux tower locations in the San Pedro River basin in Arizona. The mean daily LE varies significantly between the sites, ranging from 144 to 179     W / m 2 in the riparian areas to 36 to 76     W / m 2 at the rangeland sites. Comparison of the flux towers shows good correlation and low root mean square error at the riparian sites (0.81 to 0.82 and 37 to 51     W / m 2 , respectively) with slightly larger errors at the upland sites, where ET is strongly correlated to precipitation events. The model assumption of a linear variation in evaporative fraction across the triangular domain (LST/EVI space) results in more uncertainty under water-stressed conditions such as those found at the upland sites. Overall, the proposed MODIS-based algorithm provides reasonable estimates of riparian and upland plant water use and unique spatial and temporal information.


Introduction
Evapotranspiration (ET) is a key parameter in numerous disciplines, including water resource management, agriculture, ecology, and climate change science. However, ET is arguably one of the most difficult hydrologic components to estimate, given the dependence on a range of climatological parameters (i.e., solar radiation, temperature, wind speed, vapor pressure, etc.) as well as soil physical properties, land cover, and the surrounding environment. Traditional ground-based measurements of land surface heat fluxes (e.g., Bowen ratio system and eddy covariance system) represent a relatively small area or footprint, and spatial interpolation is sometimes used to regionalize these values. 1 A number of ET methods using meteorological and radiosonde observations are available, but the required data are only found typically over small areas of the global land surface. Also, utilization of reanalysis data has limitations due to its coarser scale. 2,3 Alternatively, remotely sensed data have been recognized as a means to provide broader areal coverage, more frequent estimates, and moderate (or high) spatial resolution. [3][4][5] A number of models with varying complexity have been developed to estimate ET based on combinations of remote-sensing observations, ground-based data, and/or high-resolution models. These include surface energy balance algorithm for land (SEBAL), mapping evapotranspiration at high resolution with internalized calibration (METRIC), atmosphere-land exchange inverse (ALEXI), and temperature-vegetation index triangle methods. [5][6][7][8][9][10][11][12][13][14] Although significant progress has been made in the past 30 years on the use of satellite-based products to provide estimates of ET at daily to weekly time scales and higher spatial scales (1 m to 1 km), most efforts have been hindered by two key issues, including: (1) water fluxes and their associated stores must be estimated indirectly using algorithms that relate, for example, leaf area index, vegetation indices, and land surface temperature, to ET, and (2) key parameters, such as wind speed, air temperature, vapor pressure deficit, and vegetation height, are not available at regional scales or from remote-sensing platforms. In addition, deriving accurate spatial representation of ground-based measurements for large heterogeneous regions is problematic, especially in areas with limited observational networks. Residual methods, such as SEBAL and METRIC, minimize the requirements of surface measurements through unique internal calibration methods that retrieve air temperature and convert wind speed at blending height. [6][7][8] However, these procedures (especially air temperature retrieval) also incur some degree of uncertainty. In order to overcome these problems, the temperature-vegetation index triangle method attempts to develop a regional parameterization of ET [evaporative fraction (EF)] through a relatively simple approach that solely uses satellite-derived surface parameters without complex parameterization and the need for extensive ground-based measurements. Further, the need for absolute accuracy from satellite-derived temperatures is avoided. 10,11 Since its inception, several modifications of the triangle (or trapezoid) method have been developed and widely applied to study soil moisture, land use, and drought monitoring, using data from several satellite sensors, including MODIS, AVHRR, GOES, etc. 12,[15][16][17][18][19][20] These studies established their own version of the triangle approach. For example, Jiang and Islam 20 utilized a simplified Priestley-Taylor formulation to build up the triangle method and Jiang et al. 18 introduced a trapezoid shape by taking into account correction parameters for vegetation transpiration that consider water stress in the root zone. 18,20 Wang et al. combined day-night temperature differences by using MODIS Terra and Aqua products with NDVI in their triangle method. 17 However, the idea behind these approaches is similar in that variations of temperatures for a given vegetation index are due to evaporative cooling effects and the proposed triangle methods are able to derive ET without auxiliary data or site-specific relationships. 21,22 In the current study, we investigate a model for estimation of ET utilizing only MODIS satellite data and limited ancillary data. We combine the triangle method developed by Jiang and Islam 12 and Wang et al. 17 with thermal inertia information obtained from the MODIS sensor to estimate a regional EF. Additionally, MODIS remote-sensing data is used for estimating daily net radiation and ground heat flux, which in combination with the triangle-based EF, returns a daily actual evapotranspiration product. Specifically, the objectives of the current study are: (1) to evaluate the modified net radiation model from the Kim and Hogue 23 framework over semiarid regions, (2) to estimate the actual ET (LE) by integrating the Kim and Hogue 23 net radiation model and additional MODIS products, (3) to investigate the ability or restriction of MODISbased data to account for the spatial variation in fluxes for many landscapes, and (4) to evaluate the robustness of this parsimonious approach against ground-based flux tower observations representing various ecosystems in a semi-arid region. 23

Study Area
Initial validation of the proposed model is carried out at four flux tower sites in southern Arizona (Table 1; Fig. 1). We utilize two riparian sites and two upland sites to evaluate the performance in important biomes in the region. Specifically, we test our algorithm at a desert grassland (Kendall site), a savanna (Santa Rita site), a riparian mesquite (Charleston site), and a riparian grassland (Lewis Spring site). The climate is classified as semi-arid with average annual precipitation ranging between 354 and 458 mm, and 53% of the total rainfall occurring from July through September during the North American monsoon (NAM). [24][25][26] Precipitation during the NAM is characterized by localized short-duration, high-intensity convective thunderstorms, which closely tie to ecosystem response in this water-limited region. 27,28 For the current study year (2005), the NAM was noted to have begun on July 18 and ended on September 10. 29 One of the flux observation sites (Kendall) is located within the Walnut Gulch experimental watershed. This site is a typical of southwestern rangeland, where cattle graze on gentle hillslopes comprised mainly of C4 grasses with a few scattered shrubs. 30 The Kendall site provides a set of soil temperature and moisture profile measurements, rain gauge, and runoff observations operated by the U.S. Department of Agriculture's Agricultural Research Service (ARS). Another upland site is the Santa Rita site located in the Santa Rita Experimental Range, approximately 45 km south of Tucson, Arizona. Over the past several decades, the rangeland around Santa Rita mesquite site transitioned from a semi-desert grassland into a savanna. 30 Vegetation consists of mesquite growing in a matrix of native and nonnative perennial grasses, subshrubs, and forbs. Soil at the upland site is a coarse-textured sandy loam. 31 Beside the two upland sites (Kendall and Santa Rita), we also evaluate the model performance at two riparian sites (Charleston and Lewis Spring), which reside along the upper San Pedro Basin. The composition and amount of vegetation in the riparian area differ significantly from those in the terrestrial upland areas, reflecting  the influence of water from the adjacent San Pedro River on soil moisture in the riparian zone. Vegetation at the Charleston site is predominantly dense woodland with a maximum vegetation height of 10 m. The Lewis Spring site is a giant sacaton grass. 32

Satellite Observations (MODIS)
The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, with 36 spectral bands (20 reflective solar and 16 thermal emissive bands), provides unprecedented information regarding vegetation and surface energy, which are critical for developing a remotely sensed ET model. 34 The proposed MODIS-based ET model requires a total of 15 variables obtained from MODIS atmospheric and land surface products for the calculation of net radiation and ground heat flux, as well as vegetation index (VI) and land surface temperature (LST) products for the development of the EF ( Table 2). The required Aqua MODIS data are extracted from the NASA EOS Gateway (https://wist.echo.nasa.gov/api/) and the Land Processes Distributed Active Archive Center (LP DAAC) (http://modis.gsfc.nasa.gov) with a standard hierarchical data format (HDF). The albedo data used in the current study is a Terra and Aqua combined product (MCD43A1). For simplicity, we utilize a solar zenith angle equal to local solar noon and an optical depth of 0.2 as default values in calculating the actual (blue-sky) albedo based on a known black-and-white sky albedo (http://daac.ornl.gov/MODIS/MODIS-menu/MCD43.html). Daily (MYD11A1) and 8-day composite (MYD11A2) LST products were obtained at a 1-km spatial resolution. Both daytime and nighttime products from MYD11A2 are acquired in order to calculate the difference of daytime and nighttime of LST. The VI was obtained using 16-day composite datasets with a 250-m resolution from both Terra (MOD13Q1) and Aqua (MYD13Q1). The MODIS composite technique (constrained-view angle-maximum value composite) includes determination of the two greatest VI values for each pixel per 16-day composite interval. The VI value observed with the nearest to nadir view is then selected for inclusion in the composite product. 35 An 8-day phasing in the production of the 16-day composites between Terra and Aqua allows the generation of a combined 8-day time series of VI data. We note that the compositing process may result in phenological changes going undetected, particularly where localized short-duration, high-intensity convective thunderstorms are characterized. 36 Since enhanced vegetation index (EVI) is optimized to improve the vegetation signal and reduce soil background influence, we utilize EVI instead of normalized difference vegetation index (NDVI) in the current study. 37 We note that the selection of VI may affect the shape of the final EVI-LST triangle; however, the EVI-calculated triangle method has been shown to perform as well as, or better than, than NDVI-calculated triangle method in most cases. 38

Estimation of Available Energy (Rn − G)
To estimate the daily available energy (Rn − G), which in combination with EF returns the daily actual ET, we propose a satellite-based stand-alone methodology.
The net radiation for clear sky conditions is estimated according to where A is the surface albedo, Rs ↓ clear is the downward shortwave radiation, Rl ↓ clear is the downward longwave radiation, and Rl ↑ clear is the upward longwave radiation for clear sky condition. The solar radiation scheme known as the PS parametric model is used to estimate the downward shortwave radiation for clear skies. 39 This scheme accounts for many factors, including gaseous absorption, Rayleigh scattering, and scattering and absorption by ozone and aerosols. Further details of this scheme can be found in Kim and Hogue. 23 Upward longwave radiation is expressed using the Stefan-Boltzmann equation: where ε s is surface emissivity, σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W∕m 2 K 4 ), and T s is the surface temperature. A parameterization scheme derived by Brutsaert 40 is chosen to estimate downwelling longwave radiation. This methodology is chosen based on the comparison among eight longwave parameterization schemes. [40][41][42][43][44][45][46][47] Rl In the Brutsaert scheme, air emissivity (ε a ) is determined by water vapor pressure (e 0 ) and air temperature (T a ). Water vapor pressure is obtained from MYD05, and air temperature is retrieved through Eqs. (6) and (7).
Under cloudy conditions, we use a formulation of downward shortwave radiation proposed by Bisht and Bras. 48 where N is the cloud cover fraction, τ is cloud optical thickness from MYD06, and θ is solar zenith angle from MYD03. In the current study, Rs ↓ clear is the theoretical clear-day shortwave radiation, which is generated by a simple interpolation scheme between adjacent clear days.
The downward longwave radiation for cloudy pixels is estimated as proposed by Bisht and Bras. 48 where ε c is cloud emissivity, and T c is the cloud temperature from MYD06. To retrieve air temperature (T a ), Bisht and Bras 48 proposed a method using the offsets between air temperature from the MYD07 data and surface temperature from the MYD06 data. However, we calculated the interpolated ratio (Ratio int ) between air temperature from the MYD07 data and surface temperature from the MYD06 data using inverse distance weighting [Eq. (6)] and back-calculated air temperature using Eq. (7).
where T s 06 is surface temperature from MYD06. To estimate the upward longwave radiation for cloudy-sky pixels, we obtain the surface temperature (T s 06 ) from MYD06 to substitute for the MYD11 data product, and surface emissivity (ε s A2 ) is obtained from MYD11A2, which is an 8-day composite product.
Using Eqs. (1)-(8), we can now calculate net radiation estimates for all sky conditions at the satellite overpass (instantaneous) time. Finally, we convert the net radiation estimates from the instantaneous time to a daily average net radiation estimate based on a sinusoidal function, assuming that net radiation begins to rise (become positive) at sunrise and declines at sunset. 23 whereRn and Rn i are daily average and instantaneous net radiation in W∕m 2 , respectively. t sunrise , t sunset , and t i are sunrise and sunset times obtained from the U.S. Naval Observatory for the satellite over-passing time. Soil heat flux (G) is estimated from the net radiation estimates and the MODIS EVI, as proposed in the Kim and Hogue 23 methodology. However, due to the relatively poor performance of the ground heat flux scheme noted in Kim and Hogue 23 at the semi-arid study site, we incorporate a modified ground heat flux scheme into our previous algorithm for application in the Arizona region. Given that empirical schemes have strong sensor dependency, we updated the previously formulated Moran and Jackson 2 relationship used in our prior work. We also evaluate the performance of various MODIS-based VI (NDVI, VI, and EVI) in approximating flux tower observations of ground heat flux for the study sites. Huete et al. demonstrated that EVI is optimized to improve the vegetation signal and reduce the soil background influences in the semiarid regions. 37 EVI also provided the best performance in our study, with optimized parameters of 0.22 and −1.4, among the vegetation indices tested. The revised ground heat flux scheme is developed as follows:Ḡ whereRn is the net radiation, andḠ is the ground heat flux at daily average in W∕m 2 .

Estimation of Evaporative Fraction
The EF is defined as the ratio of LE to available energy. 15 We utilize a model developed by Wang et al., which derives EF through interpretation of the Priestley-Taylor parameter α using the trapezoidal distribution of a ΔT (aqua day-night temperature difference)-NDVI spatial variation. 17 Again, in the current study, we incorporate EVI instead of NDVI in our model. We also chose to utilize the temporal variation of T s , rather than T s itself, to avoid a significant absolute bias from the MODIS LST product as noted in previous studies. 49,51 EF is formulated as where γ is the psychrometric constant (unit: hPa∕K), and Δ is the slope of saturated vapor pressure at the air temperature (unit: hPa∕K), which can be calculated as Due to the instability of air temperature retrieval from the MOD07 product and no significant difference (less than 5%) between the use of air and surface temperatures reported by Wang et al., surface temperature is used to estimate Δ instead of air temperature. 17 At each pixel, ET consists of surface evaporation from a bare soil and transpiration from vegetation. For a given EVI, variations of surface temperature occur primarily as a result of stressed surface soil water content. At ΔT ¼ 0, on the cold edge, bare soil and vegetation will evaporate at the potential rate, implying that ET is not limited by soil moisture. As ΔT increases, evaporation from bare soil decreases and become negligible on the warm edge.
Based on the above interpretation of the (ΔT and EVI) space, we can derive From Jiang and Islam, 12 α max ð¼ 1.26Þ and α min can be calculated as 12 By substituting α i from Eq. (15) into Eq. (12), we obtain EF for any point within the boundary of the triangle domain as Equation (17) can also be rewritten as Since a range of LST and vegetation indices are required to correctly represent the triangle boundary, we employ an 8-day composite LST (MYD11A2; the average LST of all cloud-free data in the compositing window) which minimizes cloud obstacles in determining ΔT max and ΔT min . However, ΔT i is derived from the daily LST (MYD11A1) to represent the LST at the satellite overpass. Thus, Eq. (18) can ultimately be expressed as where ΔT max A2 and ΔT min A2 are the parameters determined by a schematic plot of ΔT [MYD11A2 (8-day composite)]-EVI (8-day composite), and ΔT i A1 is the day and night temperature difference obtained from MYD11A1 (daily).

Results
We implemented the above algorithm for the entire 2005 year and based on the relevant MODIS remote-sensing products ( Table 2). We evaluate various algorithm products on an instantaneous, daily or 8-day basis. In the first part of the results (Sec. 5.1), we compare MODIS-based net radiation estimates against the net radiation observations at the flux tower sites. Second, we examine the performance of the triangle method to derive the EF. Finally, we evaluate the ET results against the latent heat flux observation from the four flux tower sites.

Validation of Net Radiation
The initial net radiation model developed in the Kim and Hogue 23 methodology incorporates the cloud product from MOD08 (1 × 1 deg) using regional fitting coefficients for downward shortwave radiation under cloudy conditions; however, this study uses the cloud product from MOD06 (1 × 1 km and 5 × 5 km) and does not require regional calibration. Bisht and Bras 48 discussed this method in detail and validated their algorithm against a range of flux towers throughout the United States. Although a semi-arid site (Desert Rock in Nevada) is included in the Bisht and Bras 48 validation, the single semi-arid/dry site was not sufficient to provide strong conclusions on the performance of the net radiation model for semi-arid regions. Thus, before application of the proposed ET model, net radiation estimates derived in the current study are further evaluated against observed net radiation from the flux tower sites.
First, we present results of the comparison between estimated instantaneous downward shortwave radiation for all sky conditions against ground-measured instantaneous solar radiation. For clear sky conditions, the shortwave radiation is estimated using the PS model from Kim and Hogue, 23 but for shortwave radiation under cloudy sky conditions, we calculate shortwave radiation using the Eq. (4) suggested by Bisht and Bras, 48 along with a regional parameterized scheme from the Kim and Hogue 23 methodology. In the comparison during clear sky days, shortwave radiation estimates have bias errors that range from 1 to 5%, and correlations are on the range of 0.94 to 0.98 across all study sites (Fig. 2). Although the Bisht and Bras 48 shortwave radiation scheme has significant advantages in that it is independent from ground-based measurements, the scatter plot between the estimated and observed quantities (Fig. 2) indicates that the result from the Bisht and Bras 48 model tends to be more scattered than the shortwave radiation estimates derived from Kim and Hogue 23 across all study sites. Bias errors in instantaneous shortwave radiation obtained from the Kim and Hogue 23 method range between 12 and 17% of the mean value, while errors from the Bisht and Bras 48 method were slightly higher, between 17 and 22% across the study sites. Root mean square errors (RMSEs) reported in Kim and Hogue 23 were also lower than those reported in Bisht and Bras. 48 Our results are also similar to a previous study, in which GOES-based solar insolation estimates have errors ranging from 5 to 10% and 15 to 30% for clear and all sky conditions, respectively. 53 The downward longwave radiation scheme, 41 as mentioned in Kim and Hogue, 23 showed larger uncertainty in semi-arid regions. 23 Therefore, we compared several existing algorithms that use more readily available meteorological observations such as air temperature and humidity. 54 Although these simpler algorithms may have more uncertainty relative to the more complex methods, they may be more useful for a variety of disciplinary satellite applications. Eight longwave parameterization schemes [40][41][42][43][44][45][46][47] were assessed at the Santa Rita site, which is the only site that provides longwave radiation observations. [40][41][42][43][44][45][46][47] The Brutsaert, 40 Garratt,45 and Keding 46 algorithms show the best overall performance in terms of correlation and RMSE (Fig. 3). Among them, the Brutsaert 40 algorithm showed the best performance with a correlation of 0.94 and an RMSE of 30 W∕m 2 . Due to the limited observed longwave radiation data available, we are not able to perform a more comprehensive assessment of downscaling algorithms across diverse biomes including deserts, semi-arid rangelands, and riparian corridors. However, based on the results noted above, we incorporate the Brutsaert 40 longwave radiation scheme into our net radiation model for our study area.
Comparison between the observed net radiation and the modeled net radiation at satellite overpass time [ Fig. 4(a)] tends to be scattered due to the cumulative effect of bias in both the shortwave and longwave radiation schemes. However, this bias was mostly related to the underestimation in the shortwave radiation. Although not shown here, the MODIS-derived upward longwave radiation is slightly overestimated (up to 52 W∕m 2 for the Santa Rita site) due to the positive bias of the LST 1-km product MODIS (MYD11). Previous studies note that the generalized split-window LST algorithm used in the 1-km MODIS LST retrieval overestimates LST in semi-arid and arid regions. 51,55,56 The overall RMSE of net radiation derived from the MODIS products ranges from 153 to 234 W∕m 2 . The largest bias is observed at the Charleston site, where the modeled net radiation is, on average, 197 W∕m 2 underestimated compared with the observations during the MODIS overpass time. Derived daily scatter plots for each site from the sinusoidal model are shown in Fig. 4(b). Results from the daily net radiation comparison are almost similar or slightly improved relative to the instantaneous net radiation comparisons in terms of %Bias with 32 and 20% at the Lewis Spring and Kendall sites, respectively, in the daily comparison. The RMSE error at Kendall is 69 W∕m 2 and 97 to 116 W∕m 2 for the remaining flux tower sites. Correlation coefficients from 0.65 to 0.69 are found at the daily scale. Due to the lack of detailed spatial information around the flux towers, it is difficult to explicitly identify sources of errors and uncertainties. However, we note that significant uncertainty may result from the obvious scale differences between the MODIS-based value and in situ data. Several model inputs (i.e., MOD06 cloud fraction and surface temperature and MOD11 surface temperature) are relatively coarse (1 to 5 km), which ultimately adds uncertainty and reduces accuracy when evaluating at the flux tower scale. In addition, topographic effects and uncertainties from instrumental measurements and remote-sensing retrieval are also a factor. 57 Relative to results from the semi-arid site (Desert Rock, Nevada) in the Bisht and Bras 48 study, which reported an RMSE of 41 W∕m 2 and a correlation (R) of 0.88, results at our study sites show relatively higher errors. However, a recent study by Tang

Sensitivity of Evaporative Fraction Model
Before validation of the final ET product, we evaluate the sensitivity of the triangle relationship (i.e., LST and EVI responses). Since the trapezoidal domain may be different during the growing (mid-May to early September) or senescent seasons, we separate results and explore triangle development during the two distinct periods. Figure 5(a) shows the relationship between the EF and the difference between daytime and nighttime LST from the Aqua platform, and Fig. 5(b) shows the correlation between the EF and EVI at individual sites. The observed EF is calculated from observed (flux tower) net radiation, ground heat flux, and latent heat flux through the relationship EF ¼ LE∕ðRn − GÞ. First, correlations between EF versus LST and EF versus EVI are very different during the growing and nongrowing seasons. For all sites, the correlations are significantly higher during the growing season than during the nongrowing season, due to the fact that EF is closely related to root-zone soil moisture during conditions of extreme soil water deficit (nongrowing season), which is difficult to represent only through either LST or EVI variables. But, in general, LST appears better correlated with EF than EVI during the nongrowing season, because the influence of soil moisture is more dominant when vegetation is small. However, during the growing season, the correlation of both LST and EVI with EF are reasonably high, ranging in magnitude from 0.63 to 0.86 and 0.50 to 0.86, respectively, indicating that LST and EVI can be used to parameterize EF. More importantly, it is clear that the slopes in the scatter plots (Fig. 5) are different for the upland grass (Kendall), savanna (Santa Rita), and riparian sites (i.e., Charleston and Lewis Spring), indicating that the influence of LST and EVI relative to EF varies from site to site. Given these results, we note that a single-triangular LST-EVI relationship may not be satisfactory for all of our study sites. Since riparian vegetation is distinctly different from the upland vegetation (density and species and better access to water compared with upland vegetation), the dominant factors affecting ET patterns are different. Ultimately, precipitation may be a more appropriate variable for predicting ET at these dry upland sites, as noted in Nagler et al. 58 We evaluate the EF derived from the triangle method against observed (flux tower) EF [ Fig. 6(a)]. The EF comparison demonstrates that the performance of the proposed triangle model in riparian areas is reasonably good, but the comparison is not as good at the semi-arid rangelands or upland sites. Particularly, the lowest correlation occurs for the Kendall site, which has the smallest vegetation cover and limited water supply, which ultimately results in extremely low EF even during growing season. Comparison at the Savanna mesquite (Santa Rita) site shows a relatively good correlation (0.65) during the growing season. At the riparian sites, the correlation varies from 0.62 to 0.80 and mean absolute error (MAE) varies from 0.10 to 0.15 (unitless) during the growing season. However, the triangle-derived EF during the nongrowing season is generally overestimated at all sites. Similarly, results reported by Tang et al. showed that the triangle method significantly overestimates EF when compared with observed EF from in situ measurements from May to September during 2004 to 2007 at two flux tower sites in Arizona. 57 Another study by Wang et al. showed that the correlation coefficients of EF ranged from −0.081 to 0.89 at a range of sites over the southern Great Plains with data collected from May to October in 2004. 17 They also presented that EF is nonuniform regardless of the vegetation uniformity having low correlations when the mean soil moisture content is low.

Validation of Evapotranspiration
Based on the net radiation estimates (Sec. 5.1) and the triangle-derived EF (Sec. 5.2), we calculate the actual ET (LE) and validate against the ground-based LE at the four study sites [ Fig. 6(b)]. Evaluating our MODIS-derived net radiation results against observed net radiation from flux towers for all sites, we hypothesize that the MODIS-derived net radiation contributes significant uncertainty or bias in the ET model. Compared with the result from the EF [ Fig. 6(a)], the estimated daily LE is slightly improved for all sites, likely due to the increased number of days used in the comparison. Particularly, Lewis Spring has the highest correlation (0.82) and low RMSE (36.67 W∕m 2 ). Results from Charleston also show similar values (0.81, 51 W∕m 2 and 55 W∕m 2 for correlation, RMSE, and MAE, respectively). Overall, the LE rates for the riparian areas are consistently higher and show lower variability, reflecting the typical features of riparian LE and the influence of a consistent supply of soil water. However, we note that underestimation of LE in riparian sites during the growing season may be related to the inclusion of LE from dry, nonriparian portions within the MODIS pixels. Unlike the riparian sites, the LE at the upland sites tends to be more strongly tied to precipitation events. During the summer monsoon season, the variability of LE is especially dependent on precipitation. We also observe that the LE is moderately correlated with precipitation (correlations of 0.58 and 0.63 for Kendall and Santa Rita sites, respectively). In addition, both of the upland sites (Santa Rita and Kendall) are sparsely vegetated and contain more bare soil exposure. Hence, the patchy distribution of vegetation at the upland sites likely causes a higher level of noise due to the scale difference between the satellite-derived VI (∼250 m) and the flux tower (footprint).
It is important to note the additional degree of uncertainty in this analysis due to the accuracy of eddy covariance systems and our comparison to these estimates. From previous work in this region, 59 LE measurements from eddy covariance tower are noted to be systematically underestimated due to the lack of energy balance closure. 30 In the current study, we applied a Bowen ratio correction method (BR method), in which the available energy is repartitioned into H and LE by conserving their ratio to enforce closure in the energy components. 57 We observed that the correct LE is consistently greater than the observed LE (MAE: 6 to 7 W∕m 2 for Kendall and Santa Rita). These differences were much larger at both riparian sites showing about RMSE of 30 W∕m 2 and MAE of 19 to 22 W∕m 2 , which is about 18 to 26% of the mean. Our results are similar to previous studies, 30 noting that LE is systematically underestimated by 17 to 27% from the lack of energy balance closure at these sites (i.e., Kendall and Santa Rita in Arizona). 59 Also, due to the long, narrow strips of vegetation along the riparian corridor, flux measurements at the riparian sites may not meet the required fetch (i.e., upwind area), which can be a significant source of error. 60 Performance of the proposed triangle model for an 8-day time series is highlighted in Fig. 7. The results are significantly improved over the daily results. The correlation is 0.72 on average for all four sites, and the RMSE is also decreased. The multiple relationships between ET, LST, and EVI are complex in semi-arid regions. However, the various relationships may be categorized depending on the land surface type (upland versus riparian). Water use efficiency is also variable depending on plant species, as seen at the grass-dominated Kendall site and the Mesquite Santa Rita site.

Summary and Conclusions
In the current study, a remote-sensing based scheme, based only on MODIS products, is developed for the assessment of daily ET distribution for all sky conditions at regional scales with a focus on performance in semi-arid regions. The approach is based on the triangle method, which relies on two remotely sensed inputs (EVI and LST) and a previously developed net radiation (Rn) model which uses a range of MODIS products. A major improvement over previous studies is that we attempt to utilize only MODIS satellite products for a simple and direct estimate of actual ET using limited ancillary data. Ground-based data is only used for regional coefficients developed in the previous PET model from Kim and Hogue. 12,17,23 The proposed model shows the potential for operational monitoring of ET, since the algorithm is not associated with groundbased meteorological data, and all required datasets can be obtained from freely available MODIS platforms. Therefore, we advocate that the model is useful for a range of applications, especially where in situ data are not readily available.
The proposed approach was initially tested in southern Arizona. Flux tower data was available in the riparian corridor of the upper San Pedro River Basin (i.e., Charleston and Lewis Spring Sacaton) in the rangeland Walnut Gulch watershed (Kendall) and in a natural upland region (Santa Rita). The net radiation model used in this study was shown to systematically underestimate the surface net radiation, resulting in underestimation as large as 122 W∕m 2 during satellite overpass. These distinct discrepancies might result from the relatively unsatisfactory performance of shortwave radiation scheme using Bisht and Bras. 48 The observed LE rates are significantly different between the riparian and rangeland upland sites with Charleston and Lewis Spring having approximately 144 to 179 W∕m 2 of LE, while stands of shrubs on Santa Rita and grass on Kendall have much lower ET rates of 36 to 76 W∕m 2 during the growing season (mid May to early September). Although the proposed ET model provides reasonable estimates of LE in the riparian sites, it is obvious that heterogeneity within the MODIS pixel (including both riparian and nonriparian) contributes to the underestimation of ET during the growing season. At the rangeland sites, LE estimation is slightly worse. Our initial results demonstrate some limitations of the proposed triangle method. Since the model presumes a linear variation in EF across the triangular domain of LST/EVI space, the assumption may result in larger errors under water-stressed conditions. Because deeper root zone soil moisture is a controlling factor in severely water stress areas, the triangular relationship between LST and EVI may not represent this soil moisture that is available to plants. Further, the relationship between EVI and root-zone soil moisture is dependent on the vegetation species and local climatic patterns. Alternatively, a modified model by Stisen et al., where decomposition is regarded as nonlinear, may provide improved estimation of EF by allowing some degree of water stress and heterogeneity within the pixel (i.e., mixture of vegetation and bare soil). 15 We note that there are still significant challenges for large-scale, regional ET modeling in semi-arid regions, and ongoing work is aimed at reducing this uncertainty. However, results from the current study indicate that MODIS-based algorithms can be utilized to track temporal and spatial changes in LE across heterogeneous, semi-arid domains. The independence from ancillary data and near real-time applicability makes the proposed MODIS-based ET model well suited as a near real-time operational model, which if we advocate can lead to improved water management in water-stressed regions.