Remote Sensing of Evapotranspiration over the Central Arizona Irrigation and Drainage District, USA

Abstract

Knowledge of baseline water use for irrigated crops in the U.S. Southwest is important for understanding how much water is consumed under normal farm management and to help manage scarce resources. Remote sensing of evapotranspiration (ET) is an effective way to gain that knowledge: multispectral data can provide synoptic and time-repetitive estimates of crop-specific water use, and could be especially useful for this arid region because of dominantly clear skies and minimal precipitation. Although multiple remote sensing ET approaches have been developed and tested, there is not consensus on which of them should be preferred because there are still few intercomparison studies within this environment. To help build the experience needed to gain consensus, a remote sensing study using three ET models was conducted over the Central Arizona Irrigation and Drainage District (CAIDD). Aggregated ET was assessed for 137 wheat plots (winter/spring crop), 183 cotton plots (summer crop), and 225 alfalfa plots (year-round). The employed models were the Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC), the Two Source Energy Balance (TSEB), and Vegetation Index ET for the US Southwest (VISW). Remote sensing data were principally Landsat 5, supplemented by Landsat 7, MODIS Terra, MODIS Aqua, and ASTER. Using district-wide model averages, seasonal use (excluding surface evaporation) was 742 mm for wheat, 983 mm for cotton, and 1427 mm for alfalfa. All three models produced similar daily ET for wheat, with 6–8 mm/day mid-season. Model estimates diverged for cotton and alfalfa sites. Considering ET over cotton, TSEB estimates were 9.5 mm/day, METRIC 6 mm/day, and VISW 8 mm/day. For alfalfa, the ET values from TSEB were 8.0 mm/day, METRIC 5 mm/day, and VISW 6 mm/day. Lack of local validation information unfortunately made it impossible to rank model performance. However, by averaging results from all of them, ET model outliers could be identified. They ranged from −10% to +18%, values that represent expected ET modeling discrepancies. Relative to the model average, standardized ET-estimators—potential ET (ET${}_{\circ }$), FAO-56 ET, and USDA-SW gravimetric-ET— showed still greater deviations, up to 35% of annual crop water use for summer and year-round crops, suggesting that remote sensing of actual ET could lead to significantly improved estimates of crop water use. Results from this study highlight the need for conducting multi-model experiments during summer-months over sites with independent ground validation.

1. Introduction

As one of the world’s largest users of fresh water resources, irrigated agriculture has a profound effect on the water cycle, water availability, distribution of water, and water quality. Irrigated croplands are essential contributors—they comprise 15% of arable land yet provide 36% of food (FAOSTAT Statistical Database, http:\apps.fao.org). However, they are also large consumers of fresh water resources, where 42% of fresh water withdrawals are used for irrigation purposes [1]; amounts that are ∼65–69% when a footprint analysis is considered [2] or when water withdrawals for thermo-electric power generation are excluded [3,4]. Despite high yields and relatively low risk (relative to rainfed lands), the large amounts of water used for irrigation often conflict with competing demands by urban populations [5]. In arid lands, these conflicts are exacerbated by already strained resources, climate change, and population growth. In Arizona, the reality of limits to Colorado Rivers resources have led to the development of a ‘Lower Basin Drought Contingency Plan’ (https://new.azwater.gov/lbdcp), an initiative that potentially has major implications (i.e., reductions in supply) for agricultural water. Addressing these problems requires assessments of irrigation practices to ensure that farms are currently using water efficiently, and the development of tools to facilitate additional water saving strategies at farm to irrigation district scales. These assessments and tools can be realized by using remote sensing data, which are being increasingly used to monitor, manage, and forecast agriculture water use at fine spatial scales. Examples include efforts such as OpenET to make remote tools widely available (https://etdata.org), operational ET mapping in California’s Central Valley [6], Climate Engine, an on-line modeling resource [7], sen-et (esa-sen4et.org), and agricultural monitoring by the European Commission’s Joint Research Centre (https://ec.europa.eu/jrc/en/research-topic/agricultural-monitoring).

Despite the strong interest in remote sensing for ET, the ideal choices for which ET model and under what circumstances they ought to be used, are not resolved. Choices are driven by the kind and quality of remote sensing and ground-based data sets, by the monitoring and frequency requirements, timeliness of results, and by the economic value of the ET maps. In recent years, model inter-comparison studies have been conducted and reported [8,9,10], each showing reasons for adopting one particular approach or another. However, trends in other disciplines such as climate assessment are highlighting that, for forecast accuracy, it is equally important, or possibly more important, to simulate processes with multiple models. Ensemble modeling quantifies uncertainties and leads to better estimates, either because model averaging tends to improve estimates, or because employing multiple models increases the likelihood of adding information to the analysis [11,12,13].

This study was undertaken to evaluate ET over three important crops—wheat, cotton, and alfalfa—with the aim of assessing modeling uncertainties when using differing but established remote sensing paradigms. The research was initiated as part of a joint US-Israeli effort to improve inverse biophysical modeling of ET, for example with approaches such as the Simple Biosphere Model (SiB2, [14]), which can provide a way to use remote sensing to diagnose the minimum crop water requirement. Ideally, many models would be included; however, for practical reasons—available resources, time, and modeling expertise—three were selected: a contextual surface energy balance model, a biophysically focused surface energy balance model, and a strictly empirical vegetation index model. Each have strengths and complementary weaknesses. The contextual model was the Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC; [15]), a thermal infrared approach that provides ET estimates over stressed and un-stressed vegetation and requires minimal ancillary correction, but lacks the capacity to disaggregate transpiration from evaporation. The biophysical model was the Two Source Energy Balance (TSEB; [16]), also a thermal model capable of diagnosing ET under wide ranges of water stress with the capacity to estimate both transpiration and evaporation, a capability that comes at increased model complexity and sometimes difficult-to-parameterize inputs. The third model was the Vegetation Index for the US Southwest (VISW; [17,18]) a visible near infrared (VNIR) empirical approach. VISW benefits are its model simplicity, use of high spatial and high temporal resolution remote sensing data, while its chief drawbacks are needed for crop-specific calibrations and relative insensitivity to rapid changes in vegetation water status. Use and interest in vegetation index (VI) models to estimate are widespread for improving otherwise standardized ET estimates [19,20].

Site selection to test the models was also based on practical matters, namely identifying sites growing familiar crops and managed using practices familiar to the authors. Additionally, it was important that the site be spatially extensive, i.e., encompass areas large enough to include 10’s to 100’s of fields, and that it contribute significant regional economic impact. The chosen site was the Central Arizona Irrigation and Drainage District (CAIDD) a 35,450 ha, agricultural region in South Central Arizona, USA. Typical crops at CAIDD have been cotton, small grains, forage, and vegetables. The useful attributes for this site are long growing seasons, non-deficit irrigation, hot and mostly dry summers, and ample clear skies for remote sensing. CAIDD made water delivery data for 2008 available, hence that year became the assessment period. Crops modeled for ET were durum wheat, cotton, and alfalfa in plots ranging between 1.4 and 313.4 ha.

2. Materials and Methods

2.1. The Central Arizona Irrigation and Drainage District (CAIDD)

The CAIDD district is located in South Central Arizona and is a triangularly shaped region ∼35 km North–South and ∼26 km East–West at its southern boundary Figure 1. CAIDD lies 90 km to the southeast of Phoenix, and 80 km northwest of Tucson. CAIDD water deliveries were ∼339,207,000 m${}^3$/year (ed4.biz). Irrigation water is supplied in approximately equal portions from ground water wells (350), and Colorado River water via the Central Arizona Project (CAP). Water deliveries within CAIDD are conveyed by three main canals: Santa Rosa, Central Main, and South Main. On-farm deliveries are provided by canal turnouts; in 2008, there were 327. Per agreement with CAIDD directors, water volumes and release dates at each turn-out during 2008 were provided to the authors, but with the condition that proprietary water use information not be disclosed. The nearest weather station data for CAIDD was provided by the Coolidge AZMET [21] station, approximately 50 km to the North. Using a CAIDD-provided map and imagery from Google Earth™ (mention of trademarked or commercial products in this manuscript is for reader convenience, USDA provides no endorsement or recommendations for these products), CAIDD canals and turn-out locations were digitized and incorporated into an in-house GIS database. Using Google Earth™ imagery and Landsat 5 (TM5) scenes, boundaries of 754 active cropping areas were digitized as shape files within the ArcGIS editor environment. Considering the objective to assess average crop evapotranspiration across the district, a small sample (∼30) would likely be sufficient for analysis; however, a larger sample size was needed to offset errors in a multi-cropping environment where planted areas change shape and size seasonally. Thus, this study identified defined crop areas for most of the 2008 season at CAIDD. Using the USDA Cropland map for Arizona 2008 [22], crop mask files for wheat, cotton, and alfalfa were generated at its source 30 m resolution. Estimated crop classification accuracy exceeded 80%. Dominant crops at CAIDD include wheat, cotton, alfalfa, corn, and sorghum, with several combinations occurring within a year. Using the total area of digitized cropping areas, ∼22,000 ha (54,305 acres), and provided water delivery for 2008 of 374,506,700 m${}^3$ (303,618 ac-ft), average irrigation water depth was 1700 mm (5.6 feet). Crop varieties grown at CAIDD in 2008 were not available.

2.2. Remote Sensing Data

This study evaluated ET over the CAIDD considering five main sources of remote sensing data collected in 2008: Landsat TM5, Landsat TM7, ASTER, MODIS Terra and Aqua. Additionally, the study implicitly used Resourcesat-1 (IRS-P6) and AWiFS via USDA Cropland crop cover maps.

The core data for the study were 22 scenes from TM5 (

Table 1. TM5 data used for the 2008 CAIDD study. Quality, representing scanline integrity, clouds, sky cover, solar elevation, and solar azimuth were provided by USGS. Weather data: 2 m air temperature (Tair), relative humidity (RH), wind speed (U), incoming solar radiation (Rsolar), and relative solar radiation (Rsolar/Clear Sky Solar Radiation; [23]), were derived from the Coolidge AZMET station https://cals.arizona.edu/azmet/az-data.htm.

Index	Date	Organized by Day of Year (DOY)	Quality	Clouds	Sky Cover (%)	Solar Elev. (${}^{\circ }$)	Solar Azi. (${}^{\circ }$)	Tair (${}^{\circ }$C)	RH	U (m/s)	Rsolar (MJ/m${}^2$)	Rel. Rsolar
1	21 January 2008 10:49	21	9	-	51	31.20	150.60	11.03	38.97	1.12	1.567	0.795
2	6 February 2008 10:48	37	9	-	1	34.60	147.50	6.50	76.36	0.70	1.738	0.813
3	22 February 2008 10:48	53	9	-	94	39.40	144.20	13.98	62.22	1.84	1.208	0.509
4	9 March 2008 10:48	69	9	-	2	45.00	140.80	16.54	28.72	0.76	2.084	0.793
5	25 March 2008 10:48	85	9	-	0	50.80	137.00	22.84	27.04	0.56	2.438	0.849
6	10 April 2008 10:46	101	9	-	0	56.30	132.20	18.70	27.19	4.23	2.845	0.929
7	26 April 2008 10:46	117	9	-	1	61.00	126.30	24.60	17.35	0.42	3.047	0.952
8	12 May 2008 10:46	133	9	-	0	64.30	119.30	27.41	18.28	3.62	2.735	0.835
9	28 May 2008 10:46	149	9	-	0	65.90	112.70	27.46	19.69	1.00	3.043	0.923
10	13 June 2008 10:48	165	NA	-	-	-	-	33.52	13.36	1.52	3.182	0.964
11	29 June 2008 10:46	181	9	-	7	65.40	107.30	36.92	16.31	1.71	2.971	0.909
12	15 July 2008 10:45	197	9	-	10	63.90	110.10	33.48	38.58	1.50	2.723	0.844
13	31 July 2008 10:43	213	9	-	2	61.80	115.90	36.09	21.95	0.43	2.693	0.849
14	16 August 2008 10:43	229	9	-	2	59.10	123.50	31.72	49.84	0.98	2.711	0.873
15	1 September 2008 10:43	245	9	-	19	55.70	131.70	30.30	51.79	0.96	2.663	0.887
16	17 September 2008 10:43	261	9	-	0	51.50	139.40	30.99	33.60	1.93	2.553	0.893
17	3 October 2008 10:43	277	9	-	65	46.80	145.80	30.16	31.58	2.52	2.169	0.811
18	19 October 2008 10:41	293	9	-	0	41.80	150.60	27.66	24.53	1.24	2.196	0.895
19	4 November 2008 10:41	309	9	-	49	37.00	153.60	21.82	40.40	2.83	1.946	0.871
20	20 November 2008 10:40	325	9	-	0	32.90	155.10	18.00	25.20	1.11	1.818	0.895
21	6 December 2008 10:40	341	9	-	0	30.00	155.10	12.89	54.38	1.10	1.648	0.874
22	22 December 2008 10:41	357	9	-	51	28.50	153.70	8.93	69.78	0.84	1.476	0.810

). Data were downloaded from the United States Geological Survey (USGS) Explorer website (https://earthexplorer.usgs.gov) as GeoTIFF files; geo-positioning errors relative to well-defined ground-reference structures such as irrigation canals were small and usually less than the nominal 30 m pixel size. In addition, 23 TM7 were also assessed, but, due to missing scan lines and inconsistent reflectance relative to TM5, their use was limited to assessment of land surface temperatures.

Normalized Difference Vegetation index (NDVI) maps were derived from 30 m red and near infrared (bands 3 and 4, respectively) obtained from the USGS LEDAPS data collection [24]. LEDAPS data are atmospherically corrected land surface reflectance, obviating the need to independently perform atmospheric corrections in the visible and near infrared bands. Land surface temperature (LST) maps were derived from 120 m TM5 band 6 and 60 m TM7 band 6 and resampled to 30 m by USGS (https://landsat.gsfc.nasa.gov/landsat-4-7-thermal-data-to-be-resampled-to-30-meters/). Note that the spatial resampling to 30 m greatly simplifies the modeling process but does not address the inability of TM to resolve temperature variations at scales less than 120 m. The impact on ET estimates for large fields would be minimal, but for small fields ET effects are significant due to mixed pixels. Bare soil would be over-estimated, while full-cover ET would be under-estimated. Considering irrigated settings, the only reliable way to address this problem is to provide higher resolution sensors, but data fusion of thermal infrared TIR with VNIR data provides a meaningful way of downscaling coarse resolution LST data in the absence of such sensors [25]. In this CAIDD study, the majority of fields exceeded 100 m extents and data fusion was not used. LST values were computed by converting TM top-of-atmosphere radiances ($L_{sensor,\lambda }$) to land surface radiances ($L_{surf,\lambda }$) via Equation (1):

$$L_{surf,\lambda }=\frac{L_{sensor,\lambda }{-L\uparrow }_{\lambda }}{{\tau }_{\lambda }}-(1-{ϵ}_{\lambda }){L\downarrow }_{\lambda }.$$

(1)

The spectrally dependent atmospheric correction terms are ${\tau }_{\lambda }$ (transmissivity), ${L\uparrow }_{\lambda }$ (upwelling atmospheric radiance), ${L\downarrow }_{\lambda }$ (downwelling radiance), ${ϵ}_{\lambda }$ (surface emissivity, 0.97 assumed), all obtained by using radiosonde atmospheric profile data from the National Oceanic and Atmospheric Administration (NOAA) FSL Tucson site [26] and the radiative transfer modeling package MODTRANv5 [27] (Spectral Sciences, Inc., Burlington, MA, USA). To convert $L_{surf}$ into LST, a temperature-spectral radiance lookup table was used. This table was created by solving the Planck blackbody function for 0.2 $\mathsf{\mu }$m steps, each weighted by the corresponding TM5 spectral response function F:

$$\overline{L}=\frac{{\sum }_{{\lambda }_{ll}}^{{\lambda }_{ul}}{L}_{\lambda }{F}_{\lambda }{R}_{\lambda }\Delta \lambda }{{\sum }_{{\lambda }_{ll}}^{{\lambda }_{ul}}{F}_{\lambda }{R}_{\lambda }\Delta \lambda }.$$

(2)

Both TM5 and TM7 were downloaded, merged and geopositioned. Although TM data are nominally georegistered, on occasion, large positioning errors occurred and had to be corrected. Use of both TM5 and TM7 would provide 8-day samples under clear sky conditions. Unfortunately, two factors made use of TM7 difficult: inconsistent reflection patterns with respect to TM5, and the well-known gaps due to the broken scan line corrector. While the gaps could be accommodated when developing plot averages, reflection inconsistencies were troublesome. For relatively slowly developing crops such as cotton, reflection anomalies could be readily identified, but, for highly dynamic alfalfa cropping systems, where cuttings occur at intervals of ∼6 weeks, uncertainties were high as to whether the TM7 variations were real or artifacts.

To help validate LST data obtained from TM sensors, 7 clear to mostly clear, L1B ASTER [28] scenes covering all or parts of CAIDD were downloaded (

Table 2. ASTER observations over CAIDD in 2008.

Index	Date	DOY	Notes
1	21 February 2008	51	northwest CAIDD
2	1 March 2008	60	few clouds, all of CAIDD
3	20 May 2008	140	small coverage of east CAIDD
4	28 June 2008	179	patchy clouds
5	24 August 2008	236	west 1/2 of CAIDD
6	25 September 2008	268	west 1/2 of CAIDD
7	2 October 2008	275	NW CAIDD

). Data from its 5 thermal infrared bands were processed in a similar way to that shown in Equations (1) and (2), but with additional steps to refine the surface temperature and emissivity estimates via an approach known as ‘normalized emissivity’ (NEM, [29]). By considering temperatures observed in different parts of the 8–11 $\mathsf{\mu }$m thermal infrared window, one is able to improve LST accuracy over estimates provided by a single thermal band such as for TM5 and TM7 [30]. ASTER is mounted on the sun-synchronous Terra platform and has similar overpass times of day as for TM5 ( 18:20 vs. 17:46 UTC), though with different overpass days. This means that, while direct LST validations were not possible, it was reasonable to expect temperature differences to be on the order of a few ${}^{\circ }$C.

Lastly, 250 m red/near infrared MODIS from both the morning overpass platform Terra, and the afternoon overpass platform Aqua were downloaded and incorporated into the ET analysis routines. The relatively coarse spatial resolution of MODIS would normally preclude its consideration for this study–typical field sizes at CAIDD are ∼180 m—these data sets have superior temporal resolution relative to TM5 and may provide improved time precision of cropping activities such as planting and harvesting. 47 MOD09Q1 and 46 MYD09Q1, weekly composited, red/NIR surface reflectance data were downloaded and included in the CAIDD data set. Map tiles, sourced in sinusoidal grid, were re-projected to UTM zone 12 coordinates as GeoTIFF files using the MODIS Reprojection Tool (MRT, [31]).

2.3. Evapotranspiration Models

To assess the efficacy of different remote sensing approaches for estimating spatially distributed ET, three different models were implemented: the Satellite-Based Energy Balance for Mapping Evapotranspiration with Internalized Calibration (METRIC [32]), the Two-Source Energy Balance (TSEB [16]), and Evapotranspiration Mapping with Vegetation Index ([17,18]). Summaries of and deviations from their standards are described below; complete and authoritative details for METRIC and TSEB can be found in the source citations. Implementation by the authors of METRIC and TSEB has been reported previously [10], while implementation and comparisons with VISW are new. To provide context for model estimates, three standardized, non-remote sensing ET estimates were included: Penman–Monteith reference evapotranspiration, ET0 [33,34]; FAO-56 ET ([15] and referenced below in Section 2.3.3); and USDA-SW ET [35], crop specific consumptive water use estimates derived from gravimetric sampling at 1-foot intervals over multiple seasons within the Salt River valley of Central Arizona. The models obtain ET using distinctly different approaches:

• METRIC: primarily relies upon spatial contrast within LST images, solves the surface energy balance, and resolves rapid changes in crop water status.
• TSEB: uses explicit LST, air temperature data, and distinguishes between soil and canopy biophysical properties, solves the surface energy balance, and resolves rapid changes.
• VISW: uses VI and local weather observations to estimate crop coefficients, potential ET ($E{T}_{\circ }$), and ET itself, without solving the surface energy balance, and only resolves changes over several days.

2.3.1. METRIC

ET estimation with METRIC is based on the assumption of surface energy balance:

$$LE=R_n-G-H,$$

(3)

where ET is represented in its energy flux form, latent heat ($LE$), and equal to net radiation ($R_n$), minus soil heat flux (G) and sensible heat flux (H). Units for all components are W/m${}^2$; sign convention is positive $R_n$ for incoming radiation, and positive G, H, $LE$ for fluxes away from the soil/canopy surface.

METRIC estimates $R_n$ using satellite-based reflectance and thermal infrared data [32,36]:

$$R_n{=R\downarrow }_s{-\alpha R\downarrow }_s{+R\downarrow }_l-R{{\uparrow }_l-(1-{ϵ}_{surf})R\downarrow }_l,$$

(4)

where the radiation components—all estimated using remotely sensed observations—are downwelling shortwave radiation (${R\downarrow }_s$), surface albedo ($\alpha$), downwelling longwave radiation (${R\downarrow }_l$), upwelling longwave radiation (${R\uparrow }_l$), and land surface emissivity (${ϵ}_{surf}$). METRIC estimates of G were derived from a SEBAL empirical relationship with surface albedo and NDVI [37]:

$$G=R_n\left[T_{surf}(0.0038+0.0074\alpha )(1-0.98NDV{I}^4)\right],$$

(5)

where $T_{surf}$ is surface temperature in Celsius. The CAIDD study obtained $R_n$ and G using Equations (4) and (5) and parameters as described in [32].

The remaining term on the right-hand side of Equation (3), H, is solved using METRIC’s distinctive feature of representing vertical temperature gradients between the surface soil/vegetation canopy and the overlying atmosphere in terms of horizontal temperature gradients. The usual resistance equation for H is described:

$$H={\rho }_a{c}_p\frac{\Delta T}{r_a},$$

(6)

where ${\rho }_a$ is density of moist air, $c_p$ is specific heat of the moist air, $\Delta T$ is the air temperature gradient between the surface and overlying air, and $r_a$ is a resistance term representing the effectiveness of heat transport between the surface and overlying air. The crux of remotely sensed surface energy balance models is how they resolve $\Delta T$ and $r_a$. In some models, such as TSEB (described below), $\Delta T$ is determined from radiometric surface temperature and ambient air temperature. The $r_a$ term is estimated from models of surface roughness, wind speed, and atmospheric stability. Following the SEBAL method, METRIC solves for these terms indirectly by selection of so-called ‘cold’ and ‘hot’ pixels that are presumed to represent reference H end-members. Knowing H for the extreme conditions allows estimation of H fluxes at all other pixels. These reference pixels are respectively denoted $H_{cold}$ and $H_{hot}$, and are used together to compute the coefficients needed to convert LST observations into an apparent $\Delta T$:

$$\Delta T=a+b{T}_s,$$

(7)

where the parameters a and b are solved used the selected hot and cold pixels:

$$b=\frac{\Delta T_{hot}-\Delta T_{cold}}{T_{s,hot}-T_{s,cold}},$$

(8)

$$a=\Delta T_{hot}-b{T}_{s,hot}.$$

(9)

$\Delta T_{hot}$ and $\Delta T_{cold}$ are temperature gradients derived from solutions of Equation (3) using reference conditions and observations. This approach simplifies and potentially increases accuracy in ET computations because biases due to errors in land surface temperature (LST) estimation are canceled by the hot/cold differencing. Thus, a reference cold pixel is assumed to represent ET at a standardized rate,$1.05*E{T}_r$, a rate typically corresponding to either a short (0.12 m) or tall (0.5 m) crop ET as defined in [23]; the additional 5% above reference ET accommodates conditions with evaporation at the soil surface. The surface energy balance equation for cold pixels becomes:

$$H_{cold}={(R_n-G)}_{cold}-L{E}_{cold},$$

(10)

where $L{E}_{cold}$ is the energy flux equivalent of 1.05 $E{T}_r$. At the ‘hot’ pixel, METRIC allows residual evaporation based on antecedent surface moisture, and follows the pattern in Equation (10). For this study, however, antecedent moisture was considered unknown and difficult to accurately model. Thus, the traditional SEBAL approach was used, where $LE$ was assumed to be zero:

$$H_{hot}=(R_n-G).$$

(11)

Having solved H, $R_n$, and G, Equation (3) could then be solved, resulting in an estimate of instantaneous $LE$ flux. Options to transform this estimate into daily ET include the use of diurnal incoming solar radiation [38] or use of a reference quantity, for example by assuming a constant evaporative fraction [39], i.e., presumption that the ratio $\frac{LE}{R_n-G}$ is constant over most of the day. These options have been reviewed and discussed in [40]. The METRIC approach, and in this study, also used for TSEB, is to compute a reference ET fraction ($E{T}_{r}F$):

$$E{T}_{r}F=\frac{E{T}_{inst}}{E{T}_r},$$

(12)

where $E{T}_{inst}$ in instantaneous ET (mm/s) computed:

$$E{T}_{inst}=\frac{LE}{{\lambda }_v{\rho }_w},$$

(13)

where ${\lambda }_v$ is latent heat of vaporization of water (J/kg) and ${\rho }_w$ is density of liquid water (kg/m${}^3$). Daily ET ($E{T}_d$) is then computed by multiplying $E{T}_{r}F$ by the cumulative $E{T}_{\circ }$:

$$E{T}_d=E{T}_{r}F\times \sum _{h=1}^{24}E{T}_{\circ ,h}.$$

(14)

A significant obstacle for wide deployment of METRIC is the stated need for trained experts knowledgeable in surface energy balance. This requirement not only limits the number of practitioners, but also introduces subjectivity and non-reproducibility of ET results [41]. This study assessed the need by implementing an alternative objective and reproducible approach based on the hypothesis that quantile selection of hot and cold pixels can return meaningful ET reference locations. The implementation is similar to that previously described in [41] who demonstrated the approach over almonds in a California study. Given a 30 m resolution TM scene, choose as cold reference pixels all those within 95–99% Normalized Difference Vegetation Index (NDVI) quantiles and also with 1–10% LST quantiles as illustrated in Figure 2. While the quantile selection could be fully automated, subjectivity when choosing quantiles was retained in order to balance the need for large (>1000) sample sizes while minimizing outliers. Next, for each selected cold pixel, neighboring pixels within a fixed search radius were scanned for maximal LST. The radius was set to 300 m, a value to approximate the typical CAIDD field size × 2 to improve the probability of finding bare soil pixels. Third, using all previously selected cold/hot value pairs, median hot and cold pixel values were computed and used to estimate the a and b coefficients in Equations (8) and (9).

2.3.2. TSEB

The TSEB model is a thermal-based remote sensing ET approach that is distinguished by its partitioning of surface fluxes into separate transpiration and evaporation components. This partitioning allows for differing transport resistances between soil, canopy, and atmosphere and accommodates use of radiometric instead of aerodynamic temperature. The separate components are obtained by remote sensing of fractional cover, where canopy and soil fractions are respectively f and $1-f$. In this study, f was obtained from NDVI observations following [42,43]:

$$f=1-{\left[\frac{NDV{I}_{max}-NDVI}{NDV{I}_{max}-NDV{I}_{min}}\right]}^{1/\xi },$$

(15)

where $NDV{I}_{max}$ is the NDVI value corresponding to full canopy cover, and $NDV{I}_{min}$ the value for bare soil. $\xi$ is a leaf angle distribution function parameter, ranging between 0.8 for planophile to >1.4 for erectophile canopy geometry. For the wheat and alfalfa plots, $\xi$ was set to 1.4, and, for cotton, 1.0 was used. TSEB modeling of the canopy requires estimation of leaf area index (LAI) to obtain net radiation values using [42]:

$$LAI=\frac{ln(1-f)}{-\beta },$$

(16)

where $\beta$ is a leaf orientation parameter. Due to uncertain calibration of Equation (16), this study used $\beta$ = 0.67, which models randomly oriented leaves. For row crops, LAI estimates need to be reduced to adjust clumped vegetation at sub-pixels scales [44,45,46]. Lacking details for cotton plots at CAIDD, a nominal clumping factor of 0.85 was used. No clumping factor was applied to wheat or alfalfa plots.

Having solved f and $LAI$, the implementation of TSEB used weather data from the AZMET Coolidge station (with the exception for air temperature as noted below) and computed soil and canopy $R_n$ values. In an iterative solution—required to solve atmospheric stability equations—the soil heat flux ($H_{soil}$) was computed:

$$H_{soil}=\rho c_p\frac{\Delta ,T}{r_a+r_{soil}},$$

(17)

where $\Delta ,T$ is the soil-air temperature difference, and a single resistance term (used by METRIC) is replaced by two: one for the canopy-air transport, another for soil-air transport. A third resistance term ($r_x$) to regulate fluxes between soil and canopy is added to create a series resistance formulation (see [16] for details). Thus, soil evaporative heat flux ($L{E}_{soil}$) is computed by residual within the soil fraction:

$$L{E}_{soil}=R_{n,soil}-G-H_{soil}.$$

(18)

For the canopy flux fraction, TSEB models transpiration at close to potential rates as estimated by the Priestley–Taylor parameter $\alpha$:

$$L{E}_{canopy}=\alpha f_G{R}_{n,canopy}\left[\frac{\Delta }{\Delta +\gamma }\right],$$

(19)

where $f_G$ is the fraction of vegetation that is green (assumed to be 1.0 in this study), and $R_{n,canopy}$ is net radiation at the canopy. Complementary to the soil flux fraction, TSEB estimates canopy heat flux ($H_{canopy}$) by residual:

$$H_{canopy}=R_{n,canopy}-L{E}_{canopy},$$

(20)

when solutions of TSEB lead to negative $L{E}_{canopy}$ values, an indicator of condensation, the $\alpha$ term is adjusted until $L{E}_{canopy}$ equals or exceeds zero. Having obtained component latent fluxes, the total surface $LE$ is:

$$LE=L{E}_{soil}+L{E}_{canopy}.$$

(21)

To obtain ET, $LE$ fluxes from the soil and canopy sources are summed and then converted as noted in Equation (13).

For the CAIDD study, the series network TSEB was implemented using methods in [16] with two exceptions. The first was to adopt median cold pixels (as computed for METRIC) for the near surface air temperature. This approach has been used in other contexts [47] and helps accommodate the fact that the nearest air temperature observations lay 50 km away from CAIDD and might have been significantly different. The second amendment was to set plant heights for wheat, cotton, and alfalfa to nominal heights as indicated in

Table 3. Nominal plant heights (m) for TSEB.

Crop	Planting	Mid-Season	Height at Mid-Season
Wheat	1 January	15 March	0.5
Cotton	15 April	15 July	0.5
Alfalfa	-	-	0.5 (before cutting)

and surface roughness over bare soil of 0.02 m. These simplifications were imposed because heights were unknown but approximations are needed to run TSEB.

2.3.3. VISW

The Vegetation Index for the South Western US (VISW) model is an empirical ET estimation approach based on experimental data collected in Arizona [17,18,48,49,50], where wheat, cotton, and alfalfa crops were grown under controlled irrigation regimes and frequently monitored with proximal sensing of NDVI. The approach aims to use remotely sensed vegetation indices as a proxy for the crop coefficient $k_{cb}$, thereby replacing standardized values with real-time observations. Use of empirically constrained vegetation indices to derive crop coefficients is a concept pursued previously by several researchers, including [14,51,52,53,54,55]. The original concept for a crop coefficient is a single scaling factor:

$$ET=k_{c}E{T}_{\circ },$$

(22)

where $E{T}_{\circ }$ is potential evapotranspiration. In the FAO-56 [15] approach, the crop coefficient is re-partitioned into transpiration and evaporation components:

$$ET=(k_{cb}+k_e)E{T}_{\circ },$$

(23)

where $k_{cb}$, a basal crop coefficient, represents transpiration, while $k_e$, represents surface evaporation component.

The simplicity and efficacy for using crop coefficients to model, forecast, and manage crop water requirements have been noted for many years, notably documented in [15,51,56]. VISW in this study did not include soil evaporation ($k_e$) because doing so would require knowing irrigation dates at plot scales, information unavailable. Consequences for the simplification depends upon irrigation practices, crop type, crop stage, and time of year. The greatest impact could be low total ET estimates of 10–15% on weekly time steps (Hunsaker personal communication). However, the underestimation is less during winter or where crop canopy is full. Examples of $k_{cb}$ coefficients are illustrated for wheat, cotton, and a hypothetical alfalfa crop with six cuttings:

A chief shortcoming of the crop coefficient methodology is its need for local calibration [15,56], which has led to suggestions by [53,54,57,58] that remote sensing observation of vegetation, via an index such as NDVI, could greatly alleviate this constraint. Adoption of VI-based ET estimation is becoming well-known and to a limited extent, commercialized or facilitated by governmental agencies at local (Crop Circle (Holland Scientific, Lincoln, NE); Greenseeker (Trimble, Sunnyvale, CA)) and regional scales (CIMIS California TOPS Satellite Irrigation Management Support (SIMS), [6], https://ecocast.arc.nasa.gov/simsi). Vegetation indices are inherently slow reactors to changing surface conditions such drought or delayed irrigation and hence are not ideal for monitoring, forecasting or managing crop water use under water deficits. Thus, while future water shortfalls may change operations, we believe few farms in the US Southwest practice deficit irrigation. This means that the VI approaches calibrated under standard conditions could be a viable water management tool.

In this study, $E{T}_{\circ }$ is derived from the tall-crop (0.5 m) standardized ET equation as presented in [23]:

$$E{T}_{\circ }=\frac{\left[0.408\Delta (R_n-G)\right]+\left[\gamma \frac{C_n}{T+273}{u}_2(e_s-e_a)\right]}{\Delta +\gamma (1+C_d{u}_2)},$$

(24)

where the notation is changed for convenience in this manuscript from ‘$E{T}_{SZ}$’ in [23] to the potential ET notation ‘$E{T}_{\circ }$’. The new terms are $\gamma$ for the psychrometric ‘constant’, $e_s$ for saturation vapor pressure, $e_a$ for ambient vapor pressure, $u_2$ for 2-m wind speed, and numerator/denominator constants $C_n$ and $C_d$ accommodating different computational time steps and crop heights (citation Table 1 is partially reproduced here as

Table 4. Constants for ASCE standardized ET computation. Units for $C_n$ in hourly time steps K mm s${}^3$ Mg${}^{-1}$h${}^{-1}$, and in daily time steps K mm s${}^3$ Mg${}^{-1}$day${}^{-1}$. Units for $C_d$ are s m${}^{-1}$.

Time Step	Short Crop		Tall Crop
	${\mathit{C}}_{\mathit{n}}$	${\mathit{C}}_{\mathit{d}}$	${\mathit{C}}_{\mathit{n}}$	${\mathit{C}}_{\mathit{d}}$
Daily	900	0.34	1600	0.38
Hourly, daytime	37	0.24	66	0.25
Hourly, nighttime	37	0.96	66	1.7

). Computation of $E{T}_{\circ }$ in Equation (24) is done by following the procedures in [23] for hourly time steps.

To provide a basis for comparison of ET results between models, the FAO-56 parameters for wheat, cotton, and alfalfa (

Table 5. FAO-56 crop parameters.

Crop		DOP	Init	Dev	Mid	Late	End
Wheat	Development Days	0	25	48	62	17	-
	$k_{cb}$	-	0.15	-	1.17	-	0.22
Cotton	Development Days	106	31	51	50	60	-
	$k_{cb}$	-	0.15	-	1.20	-	0.52

) were defined based on documentation and common planting date practices in Central Arizona. The FAO-56 methodology is illustrated for wheat and cotton in Figure 3, but not for alfalfa: it is a multi-year crop with multiple and difficult-to-predict cutting times.

Transformation of vegetation indices to $k_{cb}$ is done empirically based on replicated field experiments of different crops and reflectance measurements with accurately calibrated radiometers or imagers. Due to variations in sensor spatial resolution, spectral band placement, sensor calibration, atmospheric clarity, and soil moisture conditions, raw NDVI values were not consistent between experiments without a normalization procedure. For the calibrations used in this study, normalization re-scales original NDVI values to lie between 0 and 1:

$$N{D}^{*}=[NDVI-NDV{I}_l]/[NDV{I}_u-NDV{I}_l],$$

(25)

where TM5 NDVI values were re-scaled to new values ($N{D}^{*}$) that lie between empirically chosen minimum ($L_l$) and maximum ($L_u$) limits, roughly defined as bare-soil and full-cover.

For the three crops considered in this study, the transformation equations used are polynomials as shown in Figure 4.

For wheat, a cubic equation [18,50] is used to accommodate significantly increased ET while NDVI approaches saturation:

$$k_{cb}\left[\mathrm{wheat}\right]=\min \left\{\begin{array}{c}0.15\hfill \\ 1.75913+1.325341N{D}^{*}-1.46574N{D}^{*,2}+1.146003N{D}^{*,3}.\hfill \end{array}\right.$$

(26)

For cotton, a quadratic was adequate for experimental results [17]:

$$k_{cb}\left[\mathrm{cotton}\right]=\max \left\{\begin{array}{c}0.15\hfill \\ 0.2471+1.2012N{D}^{*}-0.2183N{D}^{*,2}.\hfill \end{array}\right.$$

(27)

For alfalfa, a quadratic was also used (Hunsaker, personal communication, based on experiments reported in [48]):

$$k_{cb}\left[\mathrm{alfalfa}\right]=\left\{\begin{array}{cc}0.3\hfill & \mathrm{if}N{D}^{*}<0.25\hfill \\ min\left\{\begin{array}{c}1.25\hfill \\ 0.40925-0.86315N{D}^{*}+1.97011N{D}^{*,2}.\hfill \end{array}\right.\hfill \end{array}\right.$$

(28)

NDVI normalization is needed to ensure that the full range of crop coefficients can be represented while viewing a range of different background soil colors and moisture [18]. To avoid biasing results based on local conditions, this study adopted a quantile selection approach, where all plot-mean NDVI values were aggregated by crop type—Wheat, Cotton, Alfalfa—and their resulting distributions were modeled with Beta functions. The Beta function is specified by two shape parameters, by convention identified as $\alpha$ and $\beta$. Because the Beta distribution domain spans 0–1, NDVI values were transformed from a nominal range found for the NDVI 2008 CAIDD data set, 0.2 to 0.85, to 0.0 to 1.0. This transformation step simplifies the function fitting process. Using empirical histograms for each crop, Beta shape parameters were determined in multiple steps. Beta shape parameters were computed by first determining sample means and variances, and then estimating $\alpha$ and $\beta$ terms via the method of moments (Equation (29)):

$$\begin{array}{c}\alpha =\overline{X}\left[(\overline{X}(1-\overline{X})/X_{var})-1\right],\\ \beta =(1-\overline{X}\left[(\overline{X}(1-\overline{X})/X_{var})-1\right].\end{array}$$

(29)

These estimates initialized their final determination with the maximum likelihood method, implemented via the R ‘fitdistr’ function found in the MASS package. Quantiles at the 10% and 90% probability levels were then computed and transformed back to source NDVI equivalents. Selection of these levels was subjective but not arbitrarily chosen: probabilities spanning a larger range of 5% and 95% were found to produce anomalously low crop coefficient values for wheat and alfalfa, while probabilities spanning a lesser range fail to reproduce representative crop coefficients for sparse cover. The normalization procedure is illustrated for wheat crops (Figure 5). Using instances with early January planting (green lines, Figure 5A), an empirical histogram was generated (Figure 5B), bins were re-scaled to lie between 0 and 1, a Beta distribution function was fit, quantiles at 10% and 90% probabilities identified, then indices transformed back to original NDVI space by inverting Equation (25). After the inversion, the recomputed quantiles are respectively the values for $NDV{I}_l$ and $NDV{I}_u$. Resulting NDVI normalization parameters are shown in

Table 6. Crop normalization parameters.

Crop	Beta Dist.		NDVI Limits		${\mathit{k}}_{\mathit{cb}}$ Limits
	$\mathit{\alpha }$	$\mathit{\beta }$	Lower	Upper	Lower	Upper
Wheat	0.58	1.26	0.06	0.65	0.15	1.25
Cotton	0.42	0.83	0.10	0.74	0.15	1.25
Alfalfa	2.20	1.89	0.27	0.72	0.3	1.25

.

3. Results

In order to evaluate ET results, data quality of input weather and remote sensing needed to be assessed as they have a direct impact on the potential accuracy and variability of ET.

3.1. CAIDD Air Temperature

Accurate near surface air temperature is required by all three ET models, either as a constraint defining $E{T}_{\circ }$ or as a boundary condition for estimating sensible heat flux values. In 2008, no weather station was deployed within CAIDD, and the nearest station was at Coolidge [59]. To address this shortcoming, two approaches were adopted. For the thermally driven models (METRIC and TSEB), median air temperature maps were created from averaged LST values extracted from cool, well-watered vegetated canopy pixel values as identified from NDVI images. The aim was to achieve air temperature estimates within ±2 ${}^{\circ }$C for every modeled satellite overpass time. Considering all available TM5 and TM7 data (respective spatial resolutions of 120 m and 60 m) and median values of automatically selected cold pixels (Figure 6), a linear model showed close correspondence between cool canopy temperatures obtained on clear sky days and AZMET Coolidge (423.0 m MSL) air temperatures adiabatically adjusted (${\widehat{T}}_{2m}$) for the mean CAIDD elevation (479.0 m MSL). Where correspondence was poor, as observed for four days- DOY 77, 285, 333, and 357-ET modeling was not done. These days were confirmed to be partly cloudy by USGS data (Table 1). Considering the high correlation between temperatures (Root Mean Squared Error (RMSE) ∼2.4 ${}^{\circ }$C), the cold-pixel air temperature proxy was adopted for use with inputs needed for TSEB and METRIC, though not for VISW, where adiabatically adjusted AZMET air temperatures were used for $E{T}_{\circ }$ estimation.

3.2. Land Surface Temperature Assessment

Land surface temperature (LST) data, critical to performance of the TSEB and METRIC models, were atmospherically corrected using the radiative transfer model, MODTRAN (Spectral Sciences, Inc., Burlington, MA, USA), and atmospheric temperature and humidity profiles from Tucson radiosonde data. Considering 44 overpass times for TM5 and TM7 in 2008, results showed ∼4.5 ${}^{\circ }$C temperature differences between the land surface and at-sensor values. Exceptions occurred on the same four overpass times as previously noted for estimating near surface air temperatures. Three relevant facts about atmospherically correction thermal infrared data are illustrated in Figure 7A. First, MODTRAN results show that requirements for accurate profiles are not demanding when LST values were ∼22 ${}^{\circ }$C, which means that for these conditions standard profiles could be used to estimate corrections. On the other hand, accurate estimates of atmospheric conditions become increasingly important as LST values increase beyond ∼35 ${}^{\circ }$C, in such cases correction errors could approach 8–10 ${}^{\circ }$C when using inaccurate atmospheric profile input data. These Tucson model results show that at-sensor temperatures are almost always cooler than LST and that the range of observed temperatures is substantially less than actual surface conditions. This latter point is illustrated in Figure 7B, where the impact of temperature span reduction is seasonal. For 2008, during monsoonal conditions, satellite observations of LST are reduced by ∼6 ${}^{\circ }$C.

Having performed atmospheric corrections to the TIR data, relative accuracy of TM5 data were assessed by comparing them with TM7 and ASTER data. Figure 8 shows the seasonal pattern for temperatures in 754 CAIDD plots: symbols represent mean values, bars represent ±1 standard deviation. None of the three sensors had coincidental overpass times so direct comparisons were not possible. However, a distinct LST bias was observed between TM5 and the other two sensors. Considering ASTER to be the most accurate, TM5 bias for 2008 was −3.8 ${}^{\circ }$C (

Table 7. Comparing LST (${}^{\circ }$) estimates from ASTER, TM7, and TM5 over CAIDD in 2008.

Sensor	DOY								Bias
	51	60	140	179	220	236	268	275
ASTER	23.9	31.0	42.5	45.7	42.6	45.5	39.6	41.2	-
TM7	29.4	29.4	42.8	41.1	39.2	42.8	40.4	37.8	−1.1
TM5	20.1	23.5	38.2	40.8	41.3	38.1	40.5	39.5	−3.8

). This offset was applied to all TM5 images prior to running TSEB, though not for METRIC since one of its features is insensitivity to LST bias.

3.3. Evapotranspiration Results: Wheat

Remotely sensed ET estimates for durum wheat plots were assessed for 137 plots across CAIDD with planting dates spanning late 2007 to early 2008. Figure 9 shows daily estimates made from the three models, METRIC (red), TSEB (blue), and VISW (green), linearly modeled between from eight overpass days (days indicated by solid squares) to daily time steps. Had precipitation or other large weather events occurred between overpasses, a physically based infill method would need to be used for greater accuracy. VISW values were smoothed with Savitzky–Golay polynomials to approximate weekly ET patterns. The average of the models, using equal weighting, is shown as a heavy purple line. Standardized ET estimates are shown as symbols: FAO-56 (solid triangles), and USDA-SW ET values (open squares). Rainfall (×0.1 mm) is shown in light blue along the bottom.

Results show good agreement over the ∼150 growing days with peak ET of 8 mm/day reached by all models and supported by FAO-56 estimates. METRIC and TSEB estimates generally closely tracked within 1 mm, indicating no discrepancies between using contextual LST values vs. using best estimates of LST and near surface air temperatures. Discrepancies existed between thermal models and the VISW between late January and mid-February, where the latter produced low daily ET values, ∼1 mm/day, while METRIC and TSEB estimates ranged 2–3 mm/day, this difference could be due to the ability of the LST approaches to detect surface evaporation while VISW was modeled solely for transpiration. Considering consistency of estimates and their average, the USDA-SW model appears to be the outlier.

When daily ET values are accumulated (Figure 10), modeling consistency is confirmed, seasonal transpiration amounts were estimated to be 726.0, 842.8, 657.1 mm respectively for METRIC, TSEB, and VISW (

Table 8. Cumulative ET results for three crops, six models, and model average. Seasonal crop water use in mm. Deviation from model average indicated in parentheses.

Crop	ET Model
	${\mathit{ET}}_{\circ }$	FAO-56	USDA-SW	METRIC	TSEB	VISW	Model Average
Wheat	770.6 (28.6)	566.8 (−175.2)	655.6 (−86.4)	726.0 (−16.0)	842.8 (100.8)	657.1 (−84.9)	742.0
Cotton	1234.6 (252.0)	972.9 (−9.7)	1046.5 (63.9)	803.6 (−179.0)	1142.8 (160.2)	1001.3 (18.7)	982.6
Alfalfa	1733.7 (306.7)	-	1887.7 (460.7)	1317.4 (−109.6)	1746.1 (319.1)	1217.3 (−209.7)	1427.0

). These values roughly agree with general estimates for durum wheat crop water requirements in Central Arizona to be 610 mm plus another 300–450 mm of water to accommodate irrigation inefficiencies [60]. If soil evaporation were included, VISW estimates would likely exceed 750 mm. Considering high model variabilities (∼85 mm), differences between models are likely not significant, except for FAO-56, which is relatively low by 175 mm (Table 8).

3.4. Evapotranspiration Results: Cotton

In contrast to winter-time estimates, results for ET modeling over cotton showed broad agreement with seasonal growth but significant inter-model biases ranging between 1 and 4 mm/day (Figure 11). For days shortly after planting METRIC, VISW indicate 2 mm/day, while TSEB returned 4–5 mm/day, an outcome indicating that TSEB estimates high surface evaporation in April. Mid-season ET estimates were unexpectedly low, ∼6 mm/day, and much less than the 8–10 mm/day estimates provided by other models. Seasonal total ET results (Figure 12) confirm estimation biases, with transpiration values of 803.6, 1142.8, 1001.3, respectively, for METRIC, TSEB, and VISW. These values fall within cotton ET obtained under detailed controlled experiments 2014 and 2015 in Maricopa, AZ, USA, where irrigation amounts ranged between 757–939 mm [61]. As for wheat, plot-to-plot variability of ET over cotton for METRIC was about twice that of TSEB. Plot variability for the VISW, however, was about the same as for TSEB, possibly indicative of close coupling of NDVI with LST values.

Because a greater number of cloudy days occurred for summer months, a background test using weekly MODIS Terra and Aqua data was conducted. Although the spatial resolution, 250 m for red and NIR bands, is questionable for fields typically less than 180 m in extent at CAIDD, the test could help answer the question about potential ET modeling improvement with more frequent satellite overpasses. Using the MOD09Q1 and MYD09Q1 climatology products, VISW was implemented, but without daily interpolation steps. Resulting ET estimates (Figure 13) closely match higher spatial resolution TM5 VISW data, suggesting that frequent overpasses could be a key step for improving remote sensing modeling of ET.

3.5. Evapotranspiration Results: Alfalfa

ET estimates over alfalfa (Figure 14) reflect a combination of the daily patterns observed over wheat and cotton, namely that there was close agreement between models in the cooler months, and significant biases in the summertime. Alfalfa is a year-round crop in Central Arizona with cuttings occurring approximately every two months, diagramatically illustrated in the hypothetical FAO-56 crop coefficient chart (Figure 4). Unfortunately for this study, temporal sampling by satellites was not sufficient to detect cutting dates with daily-weekly precision, and FAO-56 models were not run for alfalfa. Early and late year ET values ranged between 2 and 5 mm/day, while mid-summer values ranged 5–8 mm/day, less than seen in a canopy resistance study, 8–9 mm/day [62], and compatible with values obtained between September and November in 1983: a range of 4.2–7.4 mm/day reported from a Coolidge, AZ, USA study [63]. The USDA-SW model indicating peak ET values over 9 mm/day were not confirmed by remote sensing estimates. The USDA-SW model also predicts a late September ET depression from 9 to 7.5 mm/day on ∼DOY 270, a feature possibly detected by METRIC, but not confirmed by any of the other models. In general, remote sensing data overpasses were insufficiently frequent to resolve rapid ET fluctuations, meaning that the fluctuations present in the display are representations of local weather. With frequent overpasses and high spatial resolution—as observed by the Ven$\mathsf{\mu }$s satellite (https://www.theia-land.fr/en/products/venus), Sentinel 2 (https://sentinel.esa.int/web/sentinel/missions/sentinel-2) and Landsat (https://landsat.usgs.gov/landsat-science-data-products)—peak cover and harvest days could be detected and crop coefficients accurately mapped over time. However, TM periodicity was insufficient to resolve these patterns, which could mean a low-bias remote sensing of ET from all VI-based models.

Cumulative ET estimates for alfalfa (Figure 15) were the highest of the three examined crop types: 1317.4, 1746.1, and 1217.3 mm/day for METRIC, TSEB, and VISW respectively, in broad agreement with lysimeter data collected in Phoenix, AZ [64] and greater than observed at higher elevations in New Mexico, ∼1190 mm [65]. As for wheat results, plot-to-plot variability of ET, as estimated by METRIC and VISW were approximately twice that of TSEB estimates.

4. Discussion

Results from assessing multiple remote sensing data sets demonstrated strategies for handling sparse air temperatures, biased TIR satellite observation (critical for TSEB), automated LST end-member selection for the METRIC approach, and quantile selection of NDVI extremes to help standardize implementation of VI-based models. The study found that LST data over well-watered, cool canopies were reasonable proxies for near surface air temperatures, easing ET modeling problems when no local weather station is available; this addresses the problem of unavailable nearby air temperature data needed for modeling the surface energy balance. Eventually consistent, accurate, and finely gridded weather data may make this approach unnecessary. The study implemented a cross-comparison analysis of sensors—TM5, TM7, and ASTER—to improve radiometric accuracy. This improvement is similar to findings reported in [25], where multiple sensors were utilized for data fusion. By analyzing LST trends over all of 2008 observed by TM5 with respect to TM7 and ASTER, we were able to justify a 3.7 ${}^{\circ }$C LST correction to TM5, a change that significantly reduced TSEB’s ET difference with respect to the model average. Using atmospherically corrected TM5 data, we used an LST/NDVI quantile selection method to identify so-called ‘cold’ pixels, and then used those anchor points to search the neighborhood of each for ‘hot’ pixels. Very large (tens of thousands) samples were generated, greatly exceeding capabilities of manual end-member picking. Lastly, a related quantile selection approach for vegetation index normalization was implemented and tested. The approach is initially subjective in setting quantiles, but is objective and can be replicated for subsequent NDVI normalization. Unresolved in this study is an understanding of the impact of spatial resolution between sensors used to calibrate $k_{cb}$ estimates and sensors used to employ them. Data underpinning calibrations in this study are based on footprints less than 2 m in diameter with NDVI values exceeding 0.9, while satellite-based footprints were 30 m where NDVI values uncommonly exceeded 0.65. As new satellites with finer spatial and temporal sampling, as for example provided by the Ven$\mathsf{\mu }$s and ECOSTRESS missions, the scale divergence problem will be greatly reduced.

Results from this remote sensing study at CAIDD showed ET differences and uncertainties were affected by model and crop. As shown in Table 8, relative to model averaging results for cumulative ET values for all crops, METRIC was low ∼100 mm, TSEB high ∼190 mm, and VISW low ∼90 mm. These respectively represent apparent differences −10%, +18%, and −9%. For wheat, a winter crop at CAIDD, the METRIC model had superior performance with a −16 mm difference, while TSEB returned +100 mm and VISW −85 mm. For cotton, VISW came closest to ensemble average with a 19 mm bias, while METRIC and TSEB returned differences of −180 mm and +160 mm, respectively. For alfalfa, all models diverge from the ensemble average of 1427 mm by over 100 mm with differences consistent with those observed for the other two crops, i.e., METRIC and VISW were relatively low ET estimators while TSEB was relatively high. For METRIC in particular, this result is at odds with findings reported by [66], who report anomalously high ET values for a sorghum crop. Absent ground validation data, the significance of the model differences is unclear, but may be indicative of problems underlying model assumptions and sensitivities to inaccurate or infrequent inputs. In the case of METRIC, having seasonal dependency was not expected, suggesting that the horizontal temperature differencing methodology be checked for LST dependency. For TSEB, efforts to reduce TM5 LST bias effects may have been insufficient, possibly meaning that bias corrections should be adjusted on a scene-by-scene basis: bias LST corrections of 1–2 ${}^{\circ }$C could readily remove any ET bias between TSEB and model averages.

Considering model differences exhibited in this study were larger than desired—one could seek differences in the range of 10%– future implementation strategies that could be considered are adopting ensemble modeling as a standard practice and creating shared standardized data sets to be used for model verification. Inputs, outputs, parameter settings, and intermediate results from runs against such data sets could be provided as a part of the publication and documentation. Currently, ET models results are reported based on the author’s own implementation with no apparent means to verify consistency with other implementations of the same model namesake. Tracking down inconsistencies due to implementation differences could significantly improve community ET modeling results. In this study, a simple to implement, equally weighted model averaging approach, using three models, showed a benefit: outlier estimates were reduced and total ET estimates appear to track ET values known from experience. An additional benefit for averaging is a de-emphasis of model-intercomparisons, with a focus on accepting that each model has contributed complementary information. There are many documented remote sensing ET models that could be added to an ensemble framework, and doing so would add credibility to average ET estimates. Such a framework might not be practical on single-user platforms, but with rapidly improving IT technology and cloud-based computing, as being proposed by the OpenET organization, model sample sizes could be readily increased to 10 or more while remaining accessible for multiple users. The other strategy to address ET output differences—incorporation of model verification—needs community support, cooperation, and extensive sharing of common and well-documented data sets.

5. Conclusions

An ET study using remote sensing models was performed over irrigated farms in South Central Arizona and quantified daily and cumulative water use for wheat, cotton, and alfalfa. Multiple tests to evaluate input data quality were performed. Actions were taken to adjust biases due to atmospheric and sensor effects, while standardized procedures were established for LST end-member selection and normalizing vegetation indices to remove arbitrariness. Three models were implemented, two surface energy balance approaches, TSEB and METRIC, and a vegetation index approach, VISW. Differences had a seasonal dependence, where all models agreed within 18% for winter durum wheat plots, while differed up to 35% during summertime cotton and alfalfa trials. Considering that the three remote sensing models modeled ET in distinctly different ways, model averaging was found to be an objective and simple way to summarize divergent ET estimates. Future experiments will incorporate ground validation of ET to provide an independent assessment of model accuracy.