**179 W/m**in the riparian areas to 36 to

_{2}**76 W/m**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

^{2}**51 W/m**, 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.

^{2 }## 1.

## 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 semi-arid 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 MODIS-based 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}

## 2.

## 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} Using the model of Hsiegh et al., flux tower footprints ranged from 50 to 230 m for upland sites and 25 to 80 m for the riparian sites during daytime conditions.^{33}

## Table 1

Location and characteristics of the four tower sites utilized for the study and single-pixel extraction.

Site | Latitude (deg) | Longitude (deg) | Elevation (m) | Precipitation* (mm/year) | Reference crop evaporation* (mm/year) | Vegetation |
---|---|---|---|---|---|---|

Kendall | 31.73611 | −109.941 | 1083 | 177 | 5.365 | Grass |

Charleston | 31.66333 | −110.178 | 1206 | 236 | 4.560 | Mesquite woodland |

Lewis Spring | 31.56167 | −110.140 | 1242 | 290 | 4.496 | Sacaton tall grass |

Santa Rita | 31.82139 | −110.866 | 1116 | 335 | 5.132 | Desert shrub |

## *

Precipitation and reference crop evaporation (AZMET standard) are the annual averaged value for year 2005.

## 3.

## Data

## 3.1.

### In Situ Measurements

Study flux towers consist of a net radiometer instrument and several soil heat flow transducers with soil temperature/moisture sensors for estimating the available energy. Sensible and latent heat fluxes are measured by eddy covariance (consisting of a sonic anemometer and a fast-response open-path $\text{water vapor}/{\mathrm{CO}}_{2}$ sensor) at all study sites. Soil heat fluxes are measured with soil heat flux plates and are corrected for heat storage above the plate. Thus, relevant components of the surface energy balance, including measurements of net radiation, soil heat flux, sensible and latent heat fluxes, are provided at all study flux towers. In addition, ancillary meteorological observations (air temperature, wind speed, vapor, atmospheric pressure, precipitation, etc.) are available at the four study sites. Measurements at the grassland and shrubland upland sites (i.e., Kendall and Santa Rita) are made at 2 to 3 m above the local terrain. Measurements at the riparian sites (i.e., Charleston and Lewis Spring) are made at 1.5 to 2.5 m above mean canopy height. All sites have data recorded at 30-min intervals. When the daily average is calculated, days with less than 75% of data available are excluded.

## 3.2.

### 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}

## Table 2

Summary of MODIS products used in this study with relevant spatial and temporal information.

Product name | Layer | Pixel (grid) resolution | Temporal resolution | |
---|---|---|---|---|

Net radiation and ground heat flux | MYD03 (Aqua) | Solar zenith angle | 1×1 km | Daily (daytime) |

Geolocation (Lat, Lon) | 1×1 km | |||

MYD04_L2 (Aqua) | Aerosol optical depth | 10×10 km | Daily (daytime) | |

Angstrom exponent | 10×10 km | |||

MYD05_L2 (Aqua) | Water vapor | 1×1 km | Daily (daytime) | |

MYD06_L2 (Aqua) | Cloud fraction | 5×5 km | Daily (daytime) | |

Cloud optical thickness | 1×1 km | |||

Surface temperature | 5×5 km | |||

MYD07_L2 (Aqua) | Total ozone | 5×5 km | Daily (daytime) | |

Air temperature | 5×5 km | |||

Dew-point temperature | 5×5 km | |||

MYD11_L2 (Aqua) | Emissivity | 1×1 km | Daily (daytime) | |

LST | 1×1 km | |||

MOD13Q1 (Terra) | EVI | 250×250 m | 16 days | |

MYD13Q1 (Aqua) | EVI | 250×250 m | 16 days | |

MCD43A3 (Terra+Aqua) | Albedo | 1×1 km | 8 days | |

Evaporative fraction | MYD11A1 (Aqua) | LST | 1×1 km | Daily (daytime) |

MYD11A2 (Aqua) | LST | 1×1 km | 8 days (daytime, nighttime) | |

MOD13Q1 (Terra) | EVI | 250×250 m | 16 days | |

MYD13Q1 (Aqua) | EVI | 250×250 m | 16 days |

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}

## 4.

## Methodology

## 4.1.

### 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

## (1)

$$R{n}_{\text{clear}}=(1-A)Rs{\downarrow}_{\text{clear}}+Rl{\downarrow}_{\text{clear}}-Rl{\uparrow}_{\text{clear}},$$^{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 ${\epsilon}_{\mathrm{s}}$ is surface emissivity, $\sigma $ is the Stefan–Boltzmann constant ($5.67\times {10}^{-8}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}\text{\hspace{0.17em}}{\mathrm{K}}^{4}$), and ${T}_{\mathrm{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}

In the Brutsaert scheme, air emissivity (${\epsilon}_{\mathrm{a}}$) is determined by water vapor pressure (${e}_{0}$) and air temperature (${T}_{\mathrm{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}

## (4)

$$Rs{\downarrow}_{\text{cloudy}}=Rs{\downarrow}_{\text{clear}}[(1-N)+N{e}^{-\tau /\mathrm{cos}\text{\hspace{0.17em}}\theta}],$$The downward longwave radiation for cloudy pixels is estimated as proposed by Bisht and Bras.^{48}

## (5)

$$Rl{\downarrow}_{\text{cloudy}}={\epsilon}_{\mathrm{a}}\sigma {T}_{\mathrm{a}}^{4}+(1-{\epsilon}_{\mathrm{a}}){\epsilon}_{\mathrm{c}}\sigma {T}_{\mathrm{c}}^{4},$$^{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 (${\text{Ratio}}_{\mathrm{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}_{\mathrm{s}\_06}$ is surface temperature from MYD06.

To estimate the upward longwave radiation for cloudy-sky pixels, we obtain the surface temperature (${T}_{\mathrm{s}\_06}$) from MYD06 to substitute for the MYD11 data product, and surface emissivity (${\epsilon}_{\mathrm{s}\_\mathrm{A}2}$) is obtained from MYD11A2, which is an 8-day composite product.

## (8)

$$Rl{\uparrow}_{\text{cloudy}}={\epsilon}_{\mathrm{s}\_\mathrm{A}2}\sigma {T}_{\mathrm{s}\_06}^{4}.$$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}^{,}^{49}^{,}^{50}

## (9)

$$\overline{Rn}=R{n}_{\mathrm{i}}\times \frac{2}{\pi \text{\hspace{0.17em}}\mathrm{sin}\left[\right(\frac{{t}_{\mathrm{i}}-{t}_{\text{sunrise}}}{{t}_{\text{sunset}}-{t}_{\text{sunrise}}}\left)\pi \right]},$$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 semi-arid 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:

## 4.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 $\alpha $ using the trapezoidal distribution of a $\mathrm{\Delta}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}_{\mathrm{s}}$, rather than ${T}_{\mathrm{s}}$ itself, to avoid a significant absolute bias from the MODIS LST product as noted in previous studies.^{49}^{,}^{51} EF is formulated as

## (13)

$$\mathrm{\Delta}=\frac{26297.77}{{({T}_{\mathrm{a}}-29.65)}^{2}}\mathrm{exp}\left[\frac{17.67({T}_{\mathrm{a}}-273.15)}{{T}_{\mathrm{a}}-29.65}\right].$$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 $\mathrm{\Delta}$ 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 $\mathrm{\Delta}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 $\mathrm{\Delta}T$ increases, evaporation from bare soil decreases and become negligible on the warm edge.

Based on the above interpretation of the ($\mathrm{\Delta}T$ and EVI) space, we can derive

## (14)

$$\frac{\mathrm{\Delta}{T}_{\mathrm{i}}-\mathrm{\Delta}{T}^{\mathrm{min}}}{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}^{\mathrm{min}}}=\frac{({\alpha}^{\mathrm{max}}-{\alpha}^{\mathrm{min}})-({\alpha}_{i}-{\alpha}^{\mathrm{min}})}{{\alpha}^{\mathrm{max}}-{\alpha}^{\mathrm{min}}}.$$## (15)

$${\alpha}_{\mathrm{i}}=\frac{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}_{\mathrm{i}}}{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}^{\mathrm{min}}}({\alpha}^{\mathrm{max}}-{\alpha}^{\mathrm{min}})+{\alpha}^{\mathrm{min}}.$$From Jiang and Islam,^{12} ${\alpha}^{\mathrm{max}}(=1.26)$ and ${\alpha}^{\mathrm{min}}$ can be calculated as^{12}

## (16)

$${\alpha}^{\mathrm{min}}=1.26{f}_{\mathrm{veg}}=1.26\frac{{\mathrm{EVI}}_{\mathrm{i}}-{\mathrm{EVI}}^{\mathrm{min}}}{{\mathrm{EVI}}^{\mathrm{max}}-{\mathrm{EVI}}^{\mathrm{min}}}.$$By substituting ${\alpha}_{\mathrm{i}}$ from Eq. (15) into Eq. (12), we obtain EF for any point within the boundary of the triangle domain as

## (17)

$$\mathrm{EF}=\frac{\mathrm{\Delta}}{\mathrm{\Delta}+\gamma}[\frac{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}_{\mathrm{i}}}{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}^{\mathrm{min}}}({\alpha}^{\mathrm{max}}-{\alpha}^{\mathrm{min}})+{\alpha}^{\mathrm{min}}].$$Equation (17) can also be rewritten as

## (18)

$$\mathrm{EF}=1.26\frac{\mathrm{\Delta}}{\mathrm{\Delta}+\gamma}[\frac{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}_{\mathrm{i}}}{\mathrm{\Delta}{T}^{\mathrm{max}}-\mathrm{\Delta}{T}^{\mathrm{min}}}(1-{f}_{\mathrm{veg}})+{f}_{\mathrm{veg}}].$$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 $\mathrm{\Delta}{T}^{\mathrm{max}}$ and $\mathrm{\Delta}{T}^{\mathrm{min}}$. However, $\mathrm{\Delta}{T}_{\mathrm{i}}$ is derived from the daily LST (MYD11A1) to represent the LST at the satellite overpass. Thus, Eq. (18) can ultimately be expressed as

## (19)

$$\mathrm{EF}=1.26\frac{\mathrm{\Delta}}{\mathrm{\Delta}+\gamma}[\frac{\mathrm{\Delta}{T}_{\mathrm{A}2}^{\mathrm{max}}-\mathrm{\Delta}{T}_{\mathrm{A}1}^{\mathrm{i}}}{\mathrm{\Delta}{T}_{\mathrm{A}2}^{\mathrm{max}}-\mathrm{\Delta}{T}_{\mathrm{A}2}^{\mathrm{min}}}(1-{f}_{\mathrm{veg}})+{f}_{\mathrm{veg}}],$$Finally, daily LE is derived from Eq. (20), given estimated net radiation ($Rn$) [Eq. (9)], ground heat flux ($G$) [Eq. (10)] at a daily basis and evaporative fraction [Eq. (19)] with the assumption of constant EF during a day.^{52}

## 5.

## 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.

## 5.1.

### Validation of Net Radiation

The initial net radiation model developed in the Kim and Hogue^{23} methodology incorporates the cloud product from MOD08 ($1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$) using regional fitting coefficients for downward shortwave radiation under cloudy conditions; however, this study uses the cloud product from MOD06 ($1\times 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{km}$ and $5\times 5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$. The largest bias is observed at the Charleston site, where the modeled net radiation is, on average, $197\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ and 97 to $116\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ and a correlation ($R$) of 0.88, results at our study sites show relatively higher errors. However, a recent study by Tang et al. noted an RMSE of 57 and $84\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ and a correlation of 0.50 and 0.34 for MODIS-based net radiation model over the Kendall and Audubon Research Ranch flux tower sites during satellite overpass time in clear days.^{57}

## 5.2.

### 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 $\mathrm{EF}=\mathrm{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.

## 5.3.

### 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$). Results from Charleston also show similar values (0.81, $51\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ and $55\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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 ($\sim 250\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ for Kendall and Santa Rita). These differences were much larger at both riparian sites showing about RMSE of $30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ and MAE of 19 to $22\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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.

## 6.

## 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 ground-based 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{m}}^{2}$ of LE, while stands of shrubs on Santa Rita and grass on Kendall have much lower ET rates of 36 to $76\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}/{\mathrm{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.

## Acknowledgments

Financial support for this work was partially supported by a NASA Earth System Science (ESS) Fellowship (#NNX07AO53H) as well as a NSF CAREER award (#EAR0846662).

## References

*in situ*measurements at two Ameriflux sites in a semiarid region,” J. Geophys. Res. 116, D04106 (2011), http://dx.doi.org/10.1029/2010JD014543.JGREA20148-0227 Google Scholar

## Biography

**Jongyoun Kim** received her MS degree in atmospheric and oceanic sciences and her MS and PhD degrees in civil and environmental engineering from the University of California, Los Angeles (UCLA). She is currently a project scientist with the Department of Botany and Plant Sciences, UC Riverside. She specializes in the development and application of remote sensing products within hydrologic modeling systems and in the evaluation of land–atmosphere interactions at a range of temporal and spatial scales.

**Terri S. Hogue** is an associate professor in the Department of Civil and Environmental Engineering at the Colorado School of Mines. She received her BS from the University of Wisconsin-Eau Claire, and MS and PhD degrees from the Department of Hydrology and Water Resources at the University of Arizona. Her research centers on understanding hydrologic and land surface processes, with much of her work focused in semi-arid regions. Projects include wildfire impacts, urbanization and ecosystem dynamics, and hydrologic response to climate change.