Estimating Riparian and Agricultural Actual Evapotranspiration by Reference Evapotranspiration and MODIS Enhanced Vegetation Index

Abstract

Dryland river basins frequently support both irrigated agriculture and riparian vegetation and remote sensing methods are needed to monitor water use by both crops and natural vegetation in irrigation districts. We developed an algorithm for estimating actual evapotranspiration (ET_a) based on the Enhanced Vegetation Index (EVI) from the Moderate Resolution Imaging Spectrometer (MODIS) sensor on the EOS-1 Terra satellite and locally-derived measurements of reference crop ET (ET_o). The algorithm was calibrated with five years of ET_a data from three eddy covariance flux towers set in riparian plant associations on the upper San Pedro River, Arizona, supplemented with ET_a data for alfalfa and cotton from the literature. The algorithm was based on an equation of the form ET_a = ET_o [a(1 − e^−bEVI) − c], where the term (1 − e^−bEVI) is derived from the Beer-Lambert Law to express light absorption by a canopy, with EVI replacing leaf area index as an estimate of the density of light-absorbing units. The resulting algorithm capably predicted ET_a across riparian plants and crops (r² = 0.73). It was then tested against water balance data for five irrigation districts and flux tower data for two riparian zones for which season-long or multi-year ET_a data were available. Predictions were within 10% of measured results in each case, with a non-significant (P = 0.89) difference between mean measured and modeled ET_a of 5.4% over all validation sites. Validation and calibration data sets were combined to present a final predictive equation for application across crops and riparian plant associations for monitoring individual irrigation districts or for conducting global water use assessments of mixed agricultural and riparian biomes.

1. Introduction

1.1. Need for Evapotranspiration Estimates for River Basin Irrigation Districts

Dryland (arid and semi-arid) agricultural districts are frequently located in riparian floodplains, using water diverted from the river for irrigation of crops. In these districts, crops and natural riparian vegetation often must co-exist on limited water supplies [1]. Riparian corridors provide habitat for fish, mammals and resident and migratory birds, while agriculture is the primary source of income for surrounding communities. Resource managers need accurate estimates of water requirements of both crops and riparian vegetation to construct district water budgets and apportion water. Globally, crop water budgets require robust algorithms that are calibrated for specific agricultural systems [2]. Irrigated lands account for about 18% of total cropland area but 40% of global food production. They also consume 60–90% of all water diverted for human purposes. Hence, the main goal of this study was to test the feasibility of developing a remote sensing algorithm for mixed riparian and crop systems in dryland irrigation districts that can be applied on a global basis.

1.2. Methods to Determine Crop Water Budgets

Several methods exist to estimate consumptive water use by crops. Often, diversions and drainages (return flows) are gauged to estimate each term [3]. Crop consumptive use (i.e., ET_a) can then be estimated as a residual in a water balance equation. Another widely used approach is to census the planted fields within a district and to estimate consumptive use by applying a crop coefficient approach to the fields where:

$${\text{ET}}_{\text{a}}={\text{K}}_{\text{c}}{\text{ET}}_{\text{o}}$$

(1)

where ET_a is actual evapotranspiration, K_c is a crop coefficient and ET_o is reference crop evapotranspiration estimated from meteorological data [3]. K_c is derived from experimental data for each crop and relates ET_a to local meteorological conditions through ET_o. An example of this approach is the Lower Colorado River Accounting System [4]. In well-monitored irrigation districts, ET_o is calculated with data from local micrometeorological stations based on wind speed, relative humidity, temperature and net radiation data according to the FAO-56 formulation of the Penman-Monteith equation [5].

Water balance approaches can provide accurate data for closed basins over long time spans (typically annually or longer) over which changes in surface and groundwater storage are negligible compared to inflows and outflows. However, these conditions are not always met, and water balances are problematic over short time periods or for individual fields. Crop coefficient methods tend to overestimate crop consumptive use because they are typically derived for crops grown under optimal conditions in lysimeters, whereas actual field crops can have uneven stands and are subject to nutrient limitations and periodic water stress between irrigations. Also, neither of these approaches is easily adapted to riparian vegetation, for which crop coefficients might be variable or unknown.

1.3. Remote Sensing Methods for Estimating ET_a

In response to these limitations, numerous methods have been developed to estimate ET_a through remote sensing with satellite imagery [6,7]. These methods fall into two broad categories. Surface energy balance (SEB) methods use data from thermal infrared (TIR) satellite bands to estimate sensible heat flux, then calculate latent heat flux (ET_a) as a residual in the SEB equation. By contrast, vegetation index (VI) methods, used in the present study, combine estimates of green foliage density with meteorological data to estimate evaporation from transpiring vegetation (E_veg). VI-ET_o methods [8–10] substitute a VI for K_c in Equation (1) then calculate ET_a based on locally-measured ET_o:

$${\text{E}}_{\text{veg}}={\text{aET}}_{\text{o}}{\text{VI}}^{\text{n}}$$

(2)

where a is a coefficient determined by regressing measured E_veg against a VI and n is an exponent relating E_veg to VI [11]. Sometimes Equation (2) is assumed to be a linear function, in which case n = 1. In some cases, as in the present study, the calibrating data come from measurements of ET_a (e.g., from lysimeters or micrometeorological flux towers), in which case Equation (2) estimates ET_a rather than E_veg.

Both SEB and VI methods have been reviewed recently, and their strengths and weaknesses compared [10,12,13]. Gonzalez-Dugo et al. [14] found similar levels of accuracy for SEB and VI-ET_o methods compared to eddy covariance data in a mixed-crop, agricultural district in Iowa, subsequently adapting the later method to water-use monitoring in an irrigation district in Spain [15]. In this study, we tested the feasibility of developing a generalized VI-ET_o algorithm that could be used to estimate ET_a of arid-zone irrigation districts, including both crops and riparian vegetation. If successful, this algorithm could be used to characterize water use efficiency and crop and environmental water needs in this type of agricultural biome, as part of an overall strategy to assess global croplands water use [2].

1.4. Applicability of VI-ET_o Methods to Dryland Irrigation Districts

The two main limitations of VI-ET_o methods are: (1) they do not account for direct evaporation from flooded fields or from soil and canopies following rains; and (2) they do not account for short-term (periods of hours to several days) stress effects on plants [6,13]. These problems tend to be minimized in dryland irrigation districts where rainfall is infrequent and 80–90% of applied water is consumed in E_veg, see [5]. Stress of crops is minimized due to the provision of irrigation water, and riparian plants are mainly phreatophytes that have a more-or-less constant water supply from the shallow aquifer under the riparian floodplain [16] (also see [17,18]). Furthermore, an estimate of soil evaporation can be obtained separately with the two-source model for estimating ET_a, which uses drying curves based on soil type to derive evaporation rates following an irrigation or rain event [19]. In districts with significant periods of wet soil due to rain or high-frequency irrigation (e.g., by overhead spray booms), the two-source model can be used [19]. However, in flood-and-furrow irrigated arid irrigation districts, the single-source model, which lumps evaporation and transpiration together in a single K_c, often performs nearly as well as the two-source model [3]. Moisture or nutrient stress can result in short-term decreases in stomatal conductance (G_S), but in the longer term, plants tend to reduce leaf area index (LAI) to match the capacity of the environment to support photosynthesis [20]. Hence, over a crop cycle, stress results in a reduction in LAI which can be detected as a decrease in VI.

1.5. VI-ET_o Methods Applied to Crops

Numerous studies have documented a strong correlation between VIs and crop ET_a [10]. For example, the VI-ET_o concept was tested for a number of crops in California’s Central Valley [21] and forms the basis of a prototype satellite irrigation management support system developed to predict real-time crop water demand in California irrigation districts [22]. Similarly, a recent set of experiments in Spain’s Guadalquivir basin irrigation district found an acceptable level of accuracy between ET_a measured by eddy covariance flux towers and VI-ET_o estimates of ET_a for four test crops [14]. The authors concluded the VI-ET_o method produced valid and robust estimates of ET_a and they, too, developed an operational irrigation management program for farmers in the district based on spatially distributed ET_a estimates from Landsat imagery [15]. Both methods had errors under 10% when compared to ground estimates of ET_a. However, different algorithms were proposed for each irrigation district and no single VI-ET_o method has yet been applied across a wide range of irrigation districts. Furthermore, riparian vegetation has not been included in estimates of irrigation district ET_a. These methods can supplement more conventional water balance methods applied to irrigation districts in which not all the terms are adequately measured [23,24].

1.6. VI Methods Applied to Riparian Vegetation

Several studies have shown that riparian vegetation is also amendable to VI-ET_o or other VI methods for estimation of ET_a. Surface soils are typically dry for most of the year, so direct evaporation of water is low. Nagler et al. [25] developed an algorithm for riparian ET_a based on the Enhanced Vegetation Index (EVI) from the Moderate Resolution Imaging Spectrometer (MODIS) sensor on the Terra satellite that adequately predicted (r² = 0.74) ET_a measured at eight moisture flux towers on three different river systems and over five different plant associations ranging from grasses to shrubs and trees. Scott et al. [26] found that this algorithm reproduced water balance estimates of ET_a on the upper San Pedro River, Arizona, within 10%. Groeneveld et al. [16] found that a single summer Landsat Normalized Difference Vegetation Index (NDVI) image could be used to predict annual ET_a measured at moisture flux towers for a variety of riparian and desert phreatophyte sites in the western US with r² = 0.95.

Riparian plant communities are subject to environmental factors that can reduce ET_a below ET_o. ET_a typically decreases with increasing depth to groundwater below about 2 m [27]. Depth to groundwater frequently increases over the summer as water is withdrawn to support ET_a, hence ET_a/ET_o might not be constant over a growing season [28]. Furthermore, salts can accumulate in soils and aquifers, further limiting ET_a [29]. Finally, soil hydraulic properties can constrain ET_a by limiting the rate at which the capillary fringe from which phreatophytes extract water can be replenished during the day [30]. Nagler et al. [17,18], working at salinized sites on the Lower Colorado River, found that mean ET_a across sites could be predicted by VI-ET_o methods with a standard error of about 20% across sites, but that ET_a at any given site could vary much more due to differences in aquifer and soil properties among sites.

1.7. Experimental Approach

Our overall goal was to develop a robust algorithm that could be used to estimate annual rates of water consumption by crops and riparian vegetation over entire irrigation districts. Our intention is to apply this algorithm, if successful, to assessing water use efficiency and pin-pointing opportunities for water savings in global irrigation districts [2]. However, ground data for ET_a is often available only as short-term estimates over limited spatial areas. Furthermore, we wanted to develop an algorithm that could be applied over mixed vegetation districts, with crops and riparian plants of different roughness lengths and other biophysical properties. In this study, we developed an algorithm for ET_a of mixed crop and riparian areas based on ET_o and the EVI from the MODIS sensor on the NASA Terra satellite. MODIS images are acquired near-daily with 250 m resolution, hence are able to closely track the phenology of crops and riparian plants over a growing season. The algorithm was calibrated with eddy covariance flux tower data for riparian plant communities on the San Pedro River, Arizona, supplemented with data from cotton and alfalfa fields in southwestern irrigation districts. The algorithm was then tested against season-long results from seven irrigation districts or riparian zones in dryland areas in the US, Spain and Australia, determined by water balance or flux tower measurements as reported in the literature. We used seasonal or longer data sets for validation because the purpose of the algorithm is to predict season-long ET_a, not just to reproduce short-term data as in the calibration set. Examples of the calibration and validation sites are shown in Figure 1.

2. Materials and Methods

2.1. Form of the Algorithm for ET_a

The algorithm for ET_a took the form used by Nagler et al. [25] and Guerschman et al. [31]:

$${\text{ET}}_{\text{a}}={\text{ET}}_{\text{o}}[\text{a}(1-{\text{e}}^{(-\text{bEVI})})-\text{c}]$$

(3)

where a, b and c are fitting coefficients and (1 − e^(−bEVI)) is derived from the Beer-Lambert Law modified to predict absorption of light by a canopy [32]. As formulated for plant canopies, LAI normally replaces N in the Beer-Lambert equation as an estimate of density of light absorbing particles, and in our application EVI replaces LAI, assuming a linear relationship between EVI and LAI over LAI values from 0 to ca. 4 as demonstrated in [17,33]. The coefficient c accounts for the fact that EVI is not zero at zero ET_a since bare soil has a low but positive EVI. Note that we are predicting ET_a rather than E_veg with this equation because all our ground data sources for ET_a included both E_veg and direct evaporation from soils. Ground data for ET_a were used to determine the coefficients a, b and c by regression analyses.

2.2. Moisture Flux Tower Data

Riparian ET_a data came from three eddy covariance moisture flux towers set in different plant communities on the upper San Pedro River, Arizona. The Charleston Mesquite (CM) tower (Figure 1A) was in a dense, mesquite woodland; the Lewis Springs Mesquite (LSM) tower was in a less-dense mesquite shrubland site; and the Lewis Spring Sacaton (LSS) tower was in a big sacaton grass (Sporobolus wrightii) site. Site details and procedures for measuring ET_a and meteorological variables have been described in several previous publications listed below.

Scott et al. [34] described the woodland site in detail and documented the water and carbon dioxide fluxes measured in 2001 and 2002. Scott et al. [35] compared all three sites for water and carbon dioxide measurements made in 2003. This study includes these previously reported 2003 measurements plus an additional four years of data from each of these sites collected in 2004–2007. Scott et al. [34,35] gave details about the eddy covariance instrumentation and methods. Studies incorporating eddy covariance instrumentation commonly use the standard energy balance closure method to evaluate the accuracy and efficacy of their measurements [36]. The average daily closure ratio ([lE+H]/[Rn_G]) for the sites used in this study ranged annually from 0.89 to 0.92 for the woodland, 0.81 to 0.83 for the shrubland, and 0.76 to 0.78 for the grassland. These values are consistent with numerous other studies [36]. For our analysis we chose to follow Twine et al. [37] who suggested that forcing closure was justified when available energy was known and errors in its measurement modest. Consequently, we scaled our latent and sensible heat fluxes to force daily closure while conserving the measured Bowen ratio. We computed ET_o as the evaporation which would occur from a short, well-watered grass with a fixed-height of 0.12 m following the FAO-56 method [5].

2.3. ET_a and ET_o of Selected Crops

The riparian vegetation in this study only covered the low end of possible ET_a values based on ET_o. The mesquite and grass plots in this study were not fully vegetated, with bare soil or sparse vegetation occurring between trees, and LAI tended to be lower than for agricultural crops. For example, typical woodland riparian zones on three river systems had ET_a of only 40% of ET_o [25]. We included literature values for alfalfa and cotton crops to cover the higher end of ET_a values. One source of alfalfa ET_a data was from Bowen ratio flux tower measurements made in a commercial, flood-irrigated, 30.1 ha alfalfa field in the Palo Verde Irrigation District (PVID) on the Lower Colorado River in 2007–2008 (Figure 1B) [38]. Alfalfa exhibited a saw-toothed pattern of ET_a as it was cut approximately every 30 days in summer in the PVID. A single MODIS EVI pixel selected over the tower location was acquired for 16-day collection periods matched to ET_a peaks in 2007 and 2008 (3 peaks per year). ET_o was obtained from the Parker, AZ AZMET station in the same irrigation district [39].

A second source of alfalfa ET_a data was from the Imperial Irrigation District (IID) in California, using ET_a data collected at five commercial alfalfa farms over three years (2006–2008). ET_a was measured by eddy covariance towers as reported in Hanson et al. [40]. ET_o data was from the California Irrigation Management System (CIMIS) [41]. MODIS EVI pixels were acquired for the period 2006–2008, and mean annual ET_a across sites was correlated with mean annual EVI across sites. The exact location of the commercial fields was not given in [40] so we chose five fields identified as alfalfa fields by having year-round high EVI exhibiting the characteristic sawtooth patterns due to cutting. Hence, the fields we chose were not necessarily the same fields where flux towers were located but both data sets represented commercial fields within the IID.

Cotton ET_a values were obtained from the USDA-ARS Conservation & Production Research Laboratory, Bushland, Texas [42]. ET_a data were collected in a 10 ha field planted to cotton in 2008 as part of an experiment to compare lysimeter ET_a values with those collected by neutron hydroprobe soil water balance and eddy covariance flux towers (Figure 1C). Two precision weighing lysimeters (9 m² each) in the field were surrounded by cotton planted at the same density and irrigated at the same frequency as the lysimeter crops. We obtained the MODIS EVI pixel centered on the middle of the field for the dates of maximum ET_a and LAI (DOY 219-247) using neutron probe results as they were distributed over the whole field. ET_o data was from an on-site micrometeorological station and was calculated by the FAO-56 method [5].

2.4. Sources of Validation Data

We tested the ET_a algorithm against annual data collected for five dryland irrigation districts and two riparian areas based on independent estimates of ET_a not included in the calibration data sets (

Table 1. Comparison of MODIS EVI estimates of 2008 ET_a with ground measurements. PVID = Palo Verde Irrigation District, CA; IID = Imperial Irrigation District, CA; VAU = Victoria, Australia; BRID = Bear River Irrigation District, UT; LVSP = La Violada Irrigation District, Spain; CNWR = Cibola National Wildlife Refuge, CA; PR = Platte River, NE.

Location	MODIS EVI (mm·yr−1)	Ground Measurement (mm·yr−1)	No. Years Data	ETa/ETo	Difference between Means (%)	Ground Method	Ref.
PVID Crops	1,369	1,268	1	0.72	7.7	WB	[43]
IID Crops	1,144	1,194	1	0.65	4.3	WB	[44]
VAU Maize2	726	666	2	0.72	8.6	WB	[46]
BRIP3	903	877	1	0.75	2.9	WB	[47]
LVSP	676	616	8	0.66	9.3	WB	[23]
CNWR	886	817	1	0.46	8.1	BR	[38]
Riparian
PR Riparian1	431	370	2	0.63	7.7	BR	[48]
Mean	876	830	-	-	5.4

¹May–August growing season;

²November–June growing season;

³April–October growing season.

). Crop ET_a estimates were based on water balance estimates over one or more growing seasons, in which ET_a was reported as a residual of inflows minus outflows. Taghvaeian and Neale [43] constructed a water balance for the PVID for 2008 (Figure 1D), using gauged inflow and outflows and estimates of canal losses and seepage generated by the irrigation district. Cultivation is year-round and the main crops are alfalfa, cotton, small grains and winter vegetables. Their estimate for ET_a was 1,268 mm·yr⁻¹. We obtained MODIS EVI pixels for the dates covered in the water balance study but created a 500 m buffer zone to exclude MODIS pixels that were mixed with adjacent desert habitat. ET_a estimated by MODIS was then compared to ET_a by the water balance study.

Allen et al. [3] conducted a detailed water balance for IID for the period 1990–1996 based on gauged inflow and outflow data and a census of planted areas. Most of the IID has a subsurface drainage system through which the drainage fraction is collected and conveyed to the Salton Sea for disposal; hence all outflows can be quantified. Reference evapotranspiration (ET_o) is monitored at three stations in the district [41]. Over 40 different crops are grown in IID including alfalfa, cotton, sugar beets, corn, wheat, Sudan grass, Bermuda grass and winter vegetables. They reported ET_a of 1,243 mm·yr⁻¹. We updated their estimate to 2001, the first complete year for which MODIS imagery was available, using consumptive water use data from the Lower Colorado River Accounting System [44]. They reported diversions minus returns of 3.61 × 10⁹ m³·yr⁻¹ for 2001, equal to 1,763 mm·yr⁻¹ over the 2.16 × 105 ha cropping area as reported in [3]. However, returns did not include flows to the Salton Sea from subsurface drainage in IID, estimated at 1.23 × 10⁹ m³·yr⁻¹ [45] and yielding an ET_a estimate of 1,194 mm·yr⁻¹, close to the more detailed estimates in Allen et al. [3] for 1990–1996. We compared this estimate to a MODIS ET_a estimate for 2001 obtained by acquiring pixels incorporating the IID minus a buffer zone around the perimeter of the district.

Greenwood et al. [46] conducted a water balance study of maize under center-pivot irrigation in a 30.4 ha field in Northern Victoria, Australia (VAU). They measured crop water use over two growing seasons (2003–2005) based on measured application rates and soil moisture depletion rates measured weekly by neutron hydroprobe. Mean ET_a was 666 mm·yr⁻¹ over both growing seasons (ca. 150 days each). We acquired data for the single MODIS EVI pixel located over the center of the field for the dates over which crop ET_a was estimated by Greenwood et al. [46].

Lecina et al. [47] conducted a water balance study of the Bear River Irrigation Project (BRIP), Utah, for the 2008 growing season over a 1,213 ha surface irrigated area in which on-farm measurements of irrigation applied and infiltration were made at 13 sites over the season. Main crops are alfalfa and other forages, corn and winter/spring cereals. Mean ET_a was estimated to be 877 mm for the April to October irrigation season. MODIS EVI pixels covering the dates of the water balance study were acquired for the study area minus a buffer zone. Barros et al. [23,24] conducted a water balance study in La Violada Irrigation District (LVSP) (4,000 ha) (Figure 1 E) in the Ebro River Basin in Spain for the years 1995–2008. The balance was based on measured diversions and returns and estimates of ET_a based on soil moisture depletion rates following irrigation events. Main crops are corn, alfalfa and winter cereals. ET_a for 2001–2008 (the period for which MODIS imagery was available) was 616 mm·yr⁻¹. We obtained the MODIS EVI pixels encompassing the area of LVSP minus a buffer area around the perimeter for model evaluation.

The final validation included data for two riparian areas within irrigation districts: The Cibola National Wildlife Refuge (CNWR), a saltcedar (Tamarix ramosissima) dominated terrace on the Lower Colorado River below PVID (Figure 1F), and a site on the Platte River dominated by common reed (Phragmites australis). At the first site, ET_a was measured at three Bowen ratio moisture flux towers set in stands of saltcedar chosen to represent the range of conditions on the terrace in 2007 and 2008 [38]. Mean ET_a over sites and years was 817 mm·yr⁻¹. Corresponding MODIS EVI pixels were acquired to encompass the CNWR minus a buffer area around the perimeter. Irmak et al. [48] measured ET_a with a Bowen ratio moisture flux tower in 2009 and 2010 at a site on the Platte River dominated by common reed (Phragmites australis) and surrounded by rain-fed and irrigated crops in Merrick County, Nebraska (PR). Mean ET_a of the riparian vegetation over the two years was 379 mm·yr⁻¹. We acquired the single MODIS EVI pixel centered on the tower site for the dates covered by the tower ET_a measurements.

2.5. Acquisition of MODIS EVI Pixels

MODIS MOD 13 gridded 250 m VI products are supplied to end-users by NASA as atmospherically and radiometrically corrected 16-day composite images (for details of image acquisition and processing see: http://modis.gsfc.nasa.gov/); MODIS EVI pixels were obtained from the Oak Ridge National Laboratory [49]. The pixel selection tool displays the approximate MODIS pixel footprint on a high-resolution Google Earth image; hence it is possible to positively co-locate ground sites with pixel acquisition sites if images or maps are available for the ground site locations. This condition was met for all sites except the IID alfalfa sites in this study. The choice of EVI was based on its superior performance in predicting ET_a noted in previous studies [6,25]. EVI is calculated from band reflectance values as:

$$\text{EVI}=2.5\times (\text{NIR}-\text{Red})/(1+\text{NIR}+(6\times \text{Red}-7.5\times \text{Blue}))$$

(4)

where the coefficient “1” accounts for canopy background scattering and the blue and red coefficients, 6 and 7.5, minimize residual aerosol variations [6].

2.6. Statistical Analyses and Other Methods

Statistical and graphical analyses were conducted with Systat, Inc., software (SigmPlot 12.5), Chicago, IL Exponential curves and linear regression equations were fit to data points by the least squares method, and goodness of fit is reported as the coefficient of determination (r²), standard error of means (SEM) and P-value for the regression coefficient (slope) of each equation. ET_a estimates by MODIS were compared to ground estimates by the t-test for paired values. Since 6 of 7 ET_a estimates were lower by MODIS compared to ground methods, potential bias was also assessed by the Wilcoxon Signed-Rank test using an on-line calculator ( http://www.vassarstats.net/wilcoxon.html).

3. Results and Discussion

3.1 EVI and ET_a Data for San Pedro Riparian Vegetation

Sixteen-day time-series of EVI are shown in Figure 2A and ET_a at each tower are compared to ET_o in Figure 2B. Peak summer EVI and ET_a values tended to overlap among sites. Peak ET_a values were about 60% of ET_o, and were only 42% of ET_o on an annual basis, as the vegetation was dormant during winter. ET_a lagged behind ET_o, with peaks in July vs. June, respectively. ET_a was strongly correlated with EVI across sites (Figure 3), and slopes of ET_avs. EVI were not significantly different at the 95% confidence level (P > 0.05) among sites.

3.2. EVI vs. ET_a/ET_o for Combined Riparian and Crop Sites

Since ET_o differs among irrigation districts due to local meteorological conditions, ET_a was normalized as ET_a/ET_o across combined sites (Figure 4). ET_a/ET_o increased with EVI but was distinctly nonlinear. An algorithm in the form of Equation (3) fit the data with r² = 0.73 and a standard error of the mean (SEM) of 0.136 of ET_a/ET_o. Maximum ET_a/ET_o projected for EVI = 1.0 is 1.29, within the range expected for freely transpiring alfalfa crops compared to grass reference crop ET_o [3,5]. On the other hand, a linear model would greatly overestimate ET_a/ET_o if based on riparian data, and would underestimate riparian ET_a/ET_o if based on crop data. Therefore, Equation (5), which has a biophysical justification based on the Beer-Lambert Law, succeeds in unifying riparian and crop data into a single algorithm:

$${\text{ET}}_{\text{a}}={\text{ET}}_{\text{o}}\hspace{0.17em}[1.73(1-{\text{e}}^{-2.25\hspace{0.17em}\text{EVI}})-0.220]$$

(5)

3.3. Validation Data Sets

Estimates of ET_a by MODIS EVI using Equation (5) are compared to ground estimates of ET_a in Table 1. Differences between estimates for individual sites ranged from 2.9 to 9.7% and the difference between means across all sites was 5.4%. Individual site differences were not significant by paired t-test (P = 0.894). However, ET_a at 6 of 7 sites was higher by MODIS than by ground estimates, a marginally significant difference (P = 0.051) by the Wilcoxon Signed-Rank test (W = 24), indicating a possible small positive bias for the MODIS estimates compared to the ground estimates.

3.4. Final Equation for ET_a

Since ET_a across sites did not differ significantly between methods for the validation data, they were combined with the calibration data to produce a final equation of best fit (Figure 5). The r² (0.77) and SEM (0.124) improved slightly compared to calibration data (Figure 4). The final equation of best fit was:

$${\text{ET}}_{\text{a}}={\text{ET}}_{\text{o}}\hspace{0.17em}[1.65(1-{\text{e}}^{-2.25\text{EVI}})-0.169]$$

(6)

The projected ET_a/ET_o at EVI = 1.0 was 1.29. Hence, ET_a was constrained between 0 at EVI = 0.05 (e.g., bare soil or dormant vegetation) and 1.28 ET_o at EVI = 0.973, the highest EVI value observed in this study. Nagler et al. [18] developed a preliminary algorithm for ET_a for agricultural and riparian associations on the Lower Colorado River, using a limited data set and a linear regression model. Equation (6) should provide more accurate predictions over a wider range of ET_a than the previous algorithm.

ET_a maps of the Palo Verde Irrigation District agricultural fields and the Cibola National Wildlife Refuge riparian zone are displayed in Figures 6 and 7, respectively. These maps are based on three MODIS collection periods in June and July 2007, during the period of ground ET_a measurements by Bowen ratio flux towers. The maps show the general features of ET_a within each area, but due to the low resolution of MODIS imagery, they could not capture the distribution of ET_a within individual fields and could not be used for tasks such as irrigation scheduling on specific farms.

3.5. Limitations and Sources of Error and Uncertainty in ET_o and ET_a Estimates

Despite the apparent success of Equations (5) and (6) in predicting ET_a across crops and riparian vegetation, some caveats must be considered. First, the ground data used to calibrate and validate remote sensing methods are themselves subject to error and uncertainty. For example, the source of ET_o data can produce error among estimates [50]. While all of the present results were based on ET_o calculated for a hypothetical grass reference crop by the FAO-56 method [5], the required data are not available for all irrigation districts on a global basis and are rarely available for riparian corridors away from irrigation districts. Simpler methods of calculating ET_o, for example from temperature or solar radiation data, are available, but they can differ from the FAO-56 method by 10% or more and can exhibit different seasonal patterns [50]. The FAO-56 method is specifically designed for short crops with radiation as the main constraint on ET_o and might not be the best scalar for tall crops or natural vegetation units, e.g., [18]. Furthermore, ET_o is subject to measurement error. Allen et al. [3] pointed out that the CIMIS system replaced measurements of net radiation with estimates based on solar radiation in 1989, resulting in a 13% lowering of ET_o values at IID after the change in methods, a discrepancy which has not been corrected. ET_o estimates by CIMIS and AZMET stations in the same irrigation districts along the Arizona-California border also differ by about 11% (higher by AZMET) [51].

Ground methods for measuring ET_a are also subject to error. Lysimetry is perhaps the most accurate way to measure ET_a, with errors under 5% [52] but extrapolating lysimeter results to field conditions can produce errors in the order of 10–20% even under controlled conditions [42]. Water budget methods for estimating ET_a typically have error rates of 5–10% even under ideal conditions, and the error rates increase when some of the terms have to be estimated rather than measured. Allen et al. [52] estimated that Bowen ratio and eddy covariance moisture flux towers have typical errors in the order of 15–30%, which can perhaps be reduced to 10–15% under ideal conditions. In the present study, the SEM for the calibrating data set was about 23% of the mean value of ET_a/ET_o, and the SEM for the combined validation and calibration data sets was 21%, within the error bounds of the likely ground measurement errors. However, the average error for predicting ET_a/ET_o by Equation (5) was less than 10% for each individual validation site and was only 5.4% across sites. This suggests that errors are likely to be random rather than systematic, so the SEM decreases as more sites and longer measurement periods are added [53].

3.6. Comparison with Other VI Methods for ET_a

Several remote sensing methods have been developed recently to estimate ET_a with VIs or other optical band combinations (i.e., surface reflectance in visible + NIR bands, in contrast to SEB methods that use TIR bands). VI methods cannot be used to measure instantaneous rates of ET_a as can thermal band methods, but for long-term monitoring of ET_a, this is an advantage because VI values tend to be stable over days to weeks, while thermal bands provide only a snapshot of the SEB at the time over satellite overpass. We briefly review three remote sensing methods below to show the basic similarities in approach and levels of accuracy of these types of methods.

Guerschman et al. [31] developed a procedure for estimating ET_a at the continental scale in Australia. It was based on MODIS EVI in the same exponential-decay, rise-to-a-maximum form as in Equation (3), but it also used the Global Vegetation Moisture Index (GVMI) as a means to differentiate between open water and bare soil when EVI is low. The GVMI detects water by absorption in the short wave bands on MODIS; a low EVI and high GVMI indicate an open water area. The model also included a precipitation interception term (high EVI and high GVMI). Hence, the algorithms used in [31] can be applied over a wider range of landscape surfaces than the present algorithm. They calibrated it with flux tower data from two forests, two open savannas, a grassland, a floodplain and a lake, and validated it with water balance data for 227 unimpaired catchments across Australia. The model has a Root Mean Square Error (RMSE) of 22% of the mean ET_a for calibration sites, and RMSE of 19% for the validation sites, similar to the present results.

Mu et al. [54] developed a MODIS-based ET_a product meant to be applied across all global biomes, from rainforest to arctic tundra. The model was based on the MODIS fPAR/LAI product, which uses the Red and NIR bands and is linearly related to NDVI. It applies values for fPAR into the Penman-Monteith equation for ET_o to solve for ET_a. It has additional features such as corrections for view angles and a biome look-up table that accounts for differences among plant communities in the relationship between fPAR and LAI. The algorithm also adjusts leaf-level G_S as a function of biome type, and accounts for nighttime transpiration and direct evaporation from wet surfaces, among other features. When compared with ET_a measurements from 46 moisture flux tower sites in diverse biomes, the remote sensing algorithm differed from tower results by an average of about 24%. This model has been proposed as an operational MODIS ET_a product.

Fisher et al. [55] developed a global model for terrestrial ET_a that used the normalized difference vegetation index (NDVI) and the soil adjusted vegetation index (SAVI) as VI inputs. As in the present study and that of Guerschman et al. [31], transpiration by the plant canopy was assumed to be non-linearly related to LAI through the Beer-Lambert Law. ET_o was estimated by the Priestly-Taylor equation, and requires only net radiation and surface temperature, which can be derived from remote sensing inputs. The algorithm also included estimates of bare soil evaporation and direct evaporation of water from leaves following rain events, based on daytime estimates of relative humidity, assuming periods of high relative humidity are reflective of rainfall events. The model was tested against 16 flux towers representing a wide diversity of habitat types, and produced r² values of 0.55–0.96 among sites, and an overall r² of 0.90, with a SEM of 13% of mean annual ET_a across sites.

Yebra et al. [53] compared ET_a from 16 Fluxnet tower sites over six different land cover types with six different vegetation measures provided by MODIS. None of the VI or other MODIS products was uniformly superior, although EVI performed best over most of the land cover types. An ensemble of the three best-performing VI products combined with local Penman-Monteith parameters derived from the flux towers had r² = 0.82 and a RMSE of about 20% of the mean ET_a across towers. Similarly good results were obtained by scaling ET_o with EVI, as used in the present study. Both the Yebra et al. [53] and our approach benefit from including locally-measured meteorological data. Hence, an important part of the present method is availability of high-quality micrometeorological data within the irrigation districts used for calibration and validation.

3.7. Comparison of VI and Thermal Band Methods to Estimate ET_a

Thermal band and VI remote sensing methods can both give unbiased estimates of ET_a with errors or uncertainties of 10–30% under ideal conditions [6,7]. When compared head-to-head in agricultural fields equipped with flux towers, both methods can give similar results [14]. However, their operating assumptions and limitations are different, and the choice of a remote sensing method to estimate ET_a should consider their differences in light of the objectives of the measurement program.

VI methods either assume that G_S is constant for a given crop or biome type, or they adjust G_S based on remotely-sensed estimates of surface soil moisture or vapor pressure deficit [54,55]. These are reasonable assumptions for crops, which usually grow under non-stress conditions [8,10]. However, riparian plants are subject to reductions in G_S and ET_a due to factors such as salinity and depth to groundwater that cannot be detected by remote sensing, so the ratio of ET_a to LAI is not constant or predictable without knowledge of the physiological constraints on ET_a at a given site [17,18,56]. Nagler et al. [17,18] found that G_S varied by four-fold among stands of saltcedar at Cibola National Wildlife Refuge, whereas LAI only varied by two-fold. By contrast to VI methods, thermal band methods can account for reduced G_S, which is detected as an increase in surface temperature [7]. Hence, if the goal is to detect short-term water stress or to schedule irrigations, thermal band methods would be preferred over VI methods.

On the other hand, VIs provide robust indications of canopy conditions that are stable over days or weeks, and are directly related to green LAI and fractional plant cover. Hence, they are well suited for monitoring biophysical properties of plants, including ET_a and photosynthesis, over a growing season, as in the present application. By contrast, thermal band methods provide only a snap-shot of ET_a at the time of satellite overpass, which typically is from 9:30–10:30 am for Landsat and MODIS imagery in the Northern Hemisphere. ET_a is then projected over a full day by assuming that the Evaporative Fraction (EF, ratio of ET_a to net radiation or ET_o) is constant during the daylight hours [7]. This is often the case for crops [7,13]. However, Nagler et al. [17,18] found that saltcedar under salt stress showed marked midday depression of ET_a, which would lead to an overestimate of ET_a based on the assumption of constant daytime EF. Glenn et al. [56] reported that both thermal band and VI methods over-estimated ET_a by 50% or greater at a high-salinity site in Cibola National Wildlife Refuge due to violations in the assumptions underlying each method. Hence, remote sensing estimates of ET_a should consider the environmental limitations on ET_a for each application.

4. Conclusion

A significant outcome of the study is the finding that ET_a can be estimated across different crops and riparian plant communities in dryland climates using MODIS EVI and locally-available calculations of ET_o. This study provides further evidence that ET_a monitoring can be accomplished using reflectance-based remote sensing information from satellites, combined with ground data to constrain ET_o. In areas without complete meteorological information to calculate ET_o by the Penman-Monteith equation, temperature-based estimates can also be used, e.g., [17,18], and the algorithm can be adjusted to use NDVI instead of EVI. In fact, the Blaney-Criddle method for calculating ET_o, which requires only mean monthly temperature and hours of daylight from a table, had a higher correlation with ground measurements of ET_a than the Penman-Monteith methods for riparian crops [18]. This was because riparian plants in that study were deciduous—dropping leaves in winter—and were better matched to the annual temperature curve than the radiation curve, which drives the FAO-56 grass reference crop ET_o [5]. Hence, the algorithm developed here can be applied (with modification) even in data-poor regions.

The simple algorithm used in this study was successful because in dryland irrigation districts many of the confounding effects on ET_a that are unrelated to plant cover are minimized. Hence, this algorithm can provide a tool for monitoring dryland irrigation districts and riparian areas but cannot necessarily be extended to other biomes. As with other remote sensing methods for estimating ET_a, accuracy improves with repeated measurements over a growing season. For example, Allen et al. [57], found that monthly thermal-band METRIC ET_a estimates over mixed crops in Idaho irrigation districts deviated from lysimeter values by 15%, while annual estimates varied by only 1–4%. Similarly, King et al. [58] reported that monthly error rates of >20% were reduced to <10% for annual values, similar to the results reported in this paper. For random variations in ET_a due to measurement error or short-term natural fluctuations, the standard error of the annual mean value is expected to decrease by the square root of the number of replicates (months of data or number of flux towers). Hence, we conclude with Allen et al. [57] and King et al. [58] that error rates under 10% can be achieved for annual estimates of ET_a if monthly or more frequent estimates are available through the growing season.

The present algorithm is most suitable for developing annual, district-level water budgets for agricultural and riparian areas, because the resolution of a MODIS pixel is too low to evaluate individual fields (see Figures 6 and 7). The sharp field boundaries visible in Figure 1 become less distinct in MODIS images. Combining this algorithm with estimates of district water application rates can be a tool for estimating district irrigation efficiency, an important topic in arid regions of the world and where there is intense competition for water on a global basis.