Actual Evapotranspiration ( Water Use ) Assessment of the Colorado River Basin at the Landsat Resolution Using the Operational Simplified Surface Energy Balance Model

Accurately estimating consumptive water use in the Colorado River Basin (CRB) is important for assessing and managing limited water resources in the basin. Increasing water demand from various sectors may threaten long-term sustainability of the water supply in the arid southwestern United States. We have developed a first-ever basin-wide actual evapotranspiration (ETa) map of the CRB at the Landsat scale for water use assessment at the field level. We used the operational Simplified Surface Energy Balance (SSEBop) model for estimating ETa using 328 cloud-free Landsat images acquired during 2010. Our results show that cropland had the highest ETa among all land cover classes except for water. Validation using eddy covariance measured ETa showed that the SSEBop model nicely captured the variability in annual ETa with an overall R of 0.78 and a mean bias error of about 10%. Comparison with water balance-based ETa showed good agreement (R = 0.85) at the sub-basin level. Though there was good correlation (R = 0.79) between Moderate Resolution Imaging Spectroradiometer (MODIS)-based ETa (1 km spatial resolution) and Landsat-based ETa (30 m spatial resolution), the spatial distribution of OPEN ACCESS Remote Sens. 2014, 6 234 MODIS-based ETa was not suitable for water use assessment at the field level. In contrast, Landsat-based ETa has good potential to be used at the field level for water management. With further validation using multiple years and sites, our methodology can be applied for regular production of ETa maps of larger areas such as the conterminous United States.


Introduction
Water management is becoming more challenging, particularly in the arid western United States because of changes in climate, land use, and demography.One climate change study has indicated that annual runoff in the Colorado River Basin (CRB) will be reduced by 14%, 18%, and 17% for the periods 2010-2039, 2040-2069, and 2070-2098, respectively [1].The reduced availability of water in the basin could also prevent water allocation requirements of the Colorado River Compact from being met [2].Increased demand for water from various sectors has motivated water managers to demand more accurate water budgets at higher resolution so that available water resources may be better allocated.
Remote sensing has great potential but remains underutilized by practicing water resource managers [3,4].Accurately estimating consumptive water use using remotely sensed data helps water managers in planning, distribution, and management of water resources.Agro-meteorological models based on remote sensing are the most suited for estimating crop water use at the field and regional scales [5].Some of the remote sensing based models for estimating actual evapotranspiration (ET a ) include the surface energy balance index (SEBI) [6], two source model (TSM) [7], surface energy balance algorithm for land (SEBAL) [8], simplified surface energy balance index (S-SEBI) [9], surface energy balance system (SEBS) [10], ET mapping algorithm (ETMA) [11], atmosphere-land exchange inverse (ALEXI) [12], mapping evapotranspiration at high resolution with internalized calibration (METRIC) [13], simplified surface energy balance (SSEB) [14], wet METRIC (wMETRIC) [15], and operational simplified surface energy balance (SSEBop) [16].Various publications have reviewed some of these models and other methods for estimating ET a using remotely sensed data [5,[17][18][19].
The US Bureau of Reclamation of the Department of the Interior (DOI) is the nodal agency for managing water in the Colorado River Basin.As per "the 1964 Decree" of the U.S. Supreme Court in Arizona vs. California et al., the Secretary of the DOI must provide complete, detailed, and accurate records of consumptive use of water from the Colorado River.Furthermore, one of the six science strategic directions of the US Geological Survey (USGS) in the decade 2007-2017 is a Water Census (water availability and use assessment) of the United States for quantifying, forecasting, and securing freshwater for America's future [20].The DOI has launched WaterSMART (Sustain and Manage America's Resources for Tomorrow) to achieve a sustainable water strategy to meet the Nation's water needs [21].The CRB is one of the focus areas under WaterSMART initiatives for developing and testing suitable methods to meet the Water Census objectives [22].
Our objective of this study is to produce seamless ET a estimates of the Colorado River Basin at the Landsat scale using the SSEBop model and evaluate the model performance using field measurements, water balance study, and Moderate Resolution Imaging Spectroradiometer (MODIS)-based ET a estimates.

Methods and Materials
The operational Simplified Surface Energy Balance (SSEBop) model was used for estimating ET a within the CRB.This model uses pre-defined hot and cold boundary conditions unique to each pixel.Independent weather datasets from the Global Data Assimilation System (GDAS) are used to compute the reference ET (ET o ).ET o is constrained by ET fraction (ET f ) based on Landsat surface temperature and Parameter-Elevation Regressions on Independent Slopes Model (PRISM) air temperature data.A brief overview of the SSEBop model is presented here.Readers may refer to Senay et al. (2013) [16] for more details on model parameterization, hypothesis, and justification.

The Operational Simplified Surface Energy Balance (SSEBop) Model
The SSEBop Model uses a pre-defined hot and cold boundary condition for computing daily evapotranspiration ET a (mm•d −1 ) as where ET o is short grass reference ET (mm•d −1 ), ET f is ET fraction (−), and k is a scaling coefficient (−) based on calibration.In this study we used k as 1.Senay et al. (2013) [16] suggest a calibration process to determine the k value since the calculation of the pre-defined parameters may already incorporate a compensating bias.ET fraction, ET f (−), is computed as where T h is pre-defined idealized reference hot pixel temperature (K), T s is land surface temperature (K) obtained from Landsat images, and T c is pre-defined idealized reference cold pixel temperature (K).
In the SSEBop model, each pixel has a predefined hot and cold boundary values based on maximum air temperature and differential temperature (dT).Pre-defined cold pixel temperature is computed as a fraction of maximum air temperature.Maximum air temperature is not selected over the Landsat scene or the modeling domain but each pixel has its own predefined boundary values.The cold pixel temperature is approximated as being close to the corresponding air temperature based on the assumption that for a given clear-sky day, the land surface will experience an ET a rate equal to the potential rate for healthy and well watered vegetation when its Ts is close to the near-surface air temperature (i.e., little or no sensible heat flux).Calibration of the SSEBop Model using MODIS-based images has shown this coefficient as 0.993 [16].However, our calibration using Landsat images has shown a different value; hence, we computed T c as max 0.985 where T max is the maximum air temperature (K).Pre-defined hot pixel temperature is computed as Differential temperature (dT, K) is computed based on the assumption that latent heat flux and heat storage at the daily time scale for a dry, bare soil will be negligible.Thus based on energy balance, sensible heat flux (H) will be equal to net radiation (R n ) at the dry, bare soil.So we can replace H with R n in the conventional H formulation and compute dT as where r ah is the aerodynamic resistance to heat transfer from a hypothetical dry, bare surface (110 s•m −1 ), ρ a is the density of air (kg•m −3 ), and c p is the specific heat of air at constant pressure (1,004 J•kg −1 •K −1 ).The calculations of net radiation (R n ) under clear-sky condition and air density (ρ a ) are adaptions of Allen et al. (1998) [23] as described in Senay et al. [16].It is important to note that dT is unique for each period and location, but the value does not vary from year to year since it is calculated under a clear-sky condition.

Study Area
The Colorado River originates in the Rocky Mountains of the western United States and flows about 2,300 km through seven states (Wyoming, Colorado, Utah, New Mexico, Arizona, Nevada, and California) before draining into the Gulf of California (Figure 1).The CRB has an area of about 630,000 km 2 , much of it arid.The Colorado River is a major source of water supply to the southwestern United States.The Colorado River supplies water to more than 25 million people and irrigates more than 12,000 km 2 of cropland across the seven basin states [22].The annual flow of the river has ranged from 6.5 billion cubic meters (BCM) to 29.6 BCM during 1906-2000 [1].The combined reservoir storage capacity (74.0 BCM) within the basin is about four times the long-term average annual flow (16.7 BCM).Based on the National Land Cover Database (NLCD, 2006) [24], shrubland (61%) is the most dominant land cover in the study area followed by evergreen forest (19%) (Table 1).Most of the pastureland is in the upper basin, most of the cropland is in the lower basin, and shrubland is distributed throughout the basin (Figure 2).Land use/land cover of the CRB is also used as an input in the Lower Colorado River Accounting System (LCRAS), a water accounting system developed by the USGS and the Bureau of Reclamation [25].Congalton et al. (1998) [26] developed a procedure to accurately map agricultural crops and other land cover in the lower CRB.For additional detail on geo-topographic and climatic conditions of the study area, see Kumar and Duffy (2009) [27].

Processing of Landsat Images
The CRB is covered by 44 Landsat scenes spread over paths 33-44 and rows 30-38 (Figure 3).First, we downloaded all the Landsat images (Thematic mapper and Enhanced thematic mapper plus) of the study area for 2010 with less than 10% cloud cover from the Earth Explorer site (http://earthexplorer.usgs.gov/).We then selected images that were nearly cloud free or had cloud only along the edges (total 328 images).In general, one Landsat image per month is ideal, but two images are desired for cropland, particularly during the crop growing season.However, based on our Landsat processing experience, reasonable results can be obtained based on 10-12 images per year.Processing of Landsat images for obtaining land surface temperature (T s ) was carried out as follows.First, the digital numbers (DN) of downloaded Landsat images were converted to at-sensor radiance (L) as where , DN is quantized calibrated pixel value (−), Q calmax is quantized calibrated pixel value corresponding to the L max (−), and Q calmin is quantized calibrated pixel value corresponding to the L min (−).The values of L max , L min , Q calmax , and Q calmin for each band are provided in the metafile of unzipped downloaded images.
The at-sensor radiance for the shortwave bands was converted to top of atmosphere reflectance (ρ) as where ρ is planetary top of atmosphere reflectance (−), d r is earth-sun distance parameter (−), ESUN is mean exoatmospheric solar irradiance (W•m −2 •μm −1 ), and θ is solar zenith angle (degree).ESUN values are given in Chander et al. (2009) [28].Top of atmosphere albedo and at surface albedo were calculated as ( ) where α toa is albedo at the top of the atmosphere (−), α is at surface albedo (−), c is weighting coefficient (based on Tasumi et al., 2008) [29], α path is albedo path radiance ranging from 0.025 to 0.04 (−), and τ sw is transmittance as computed below.
where Z is elevation above the mean sea level (m).Surface emissivity and land surface temperature were computed as presented in Allen et al. (2007) [13] without using any atmospheric radiative transfer simulation model such as MODTRAN.In some desert areas where albedo is high (>0.3), the radiometric land surface temperature tends to decrease, potentially due to reduced net radiation [16].We corrected those T s values to avoid erroneous ET fraction as _ 50 ( 0.3) where T s_co is corrected land surface temperature (K).The SSEBop model algorithm was implemented with Landsat images using Model Maker in Erdas Imagine 2011 (version 11.0.4)(Intergraph Corporation, Huntsville, AL, USA).Each individual scene was processed separately for computing ET a on the day of the satellite overpass.All the available daily ET a images were used for upscaling to compute annual ET a as discussed in the Computation of Annual ET section.

Other Supporting Data
Monthly minimum and maximum temperature and precipitation (P) gridded 4 km data for 2010 were downloaded from the Parameter-Elevation Regressions on Independent Slopes Model (PRISM) website.We used these data in their original format without any spatial interpolation to preserve the data integrity.PRISM is a knowledge-based system that uses ground-based meteorological station records, a digital elevation model (DEM), and many other geographic datasets to generate gridded estimates of monthly climatic parameters [30].Monthly runoff (Q) based on Hydrologic Unit Code (HUC) for 2010 was downloaded from the USGS Water Watch website (http://waterwatch.usgs.gov/index.php).The estimates of HUC runoff are generated by the USGS by combining historical flow data from stream gauges, the drainage basins of the stream gauges, and the HUC boundaries [31].Annual ET a based on spatial water balance at the HUC8 level was carried out by subtracting annual HUC-based runoff from the annual PRISM precipitation.Only PRISM-based air temperature (not precipitation) was used in the SSEBop model for computing Landsat-based ET a .Thus ET a based on water balance (residual of precipitation and runoff, i.e., P-Q) should serve as an independent data source for validation of Landsat-based annual ET a .We assume that there is no change in annual storage of water; therefore, P-Q is equivalent to ET a .
Topographic elevation data of the study area were derived from the February 2000 National Aeronautics and Space Administration (NASA) Shuttle Radar Topography Mission (SRTM), a seamless and complete void-filled spatial dataset [32].Data were downloaded from the SRTM data portal.The daily reference evapotranspiration (ET o ) (Senay et al., 2008) [33] was computed from climate parameter data extracted from global data assimilation system (GDAS) based 6-hourly weather datasets [34].The daily ET o data are available at the Famine Early Warning Systems Network website (http://earlywarning.usgs.gov/adds/global/index.php).The global ET o data available at 1° ground resolution were downscaled to 10 km resolution [14].Gridded flux data are a global, spatially and temporally explicit estimates of latent flux upscaled from flux measurements based on remote sensing indices, climate and meteorological data, and information on land use [35].These data were obtained from the Max Planck Institute for Biogeochemistry server.MODIS-based ET a maps using the SSEBop model [16] were also used for result comparison.

Computation of Annual Evapotranspiration
ET a on the day of the satellite overpass was computed using the SSEBop model.Computation of seasonal/annual ET a becomes challenging when daily ET a is not available because of Landsat's 16-day repeat cycle and/or a lack of cloud-free images.We have used upscaling of daily ET a to annual ET a based on the ratio of actual ET to reference ET (ET coefficient).A study by Singh (2009) [36] with multiyear datasets showed that there was no single interpolation method that worked better than other methods under all conditions.Suitability of the method depended on the number of images per season, number of days between two consecutive images, extreme values on any particular day, and duration of the season [37].The fixed method worked well for longer duration (annual).So we used the fixed method for computing annual ET based on ET coefficient: ( ) where ET annual is the total annual ET (mm), ET oi is the reference ET (mm) for period i (days), and ET fi is the representative ET coefficient (−) for period i.
After computing the annual ET a for each path/row, all annual ET a images were mosaicked.Mosaicking of annual ET a for different path/row is done to obtain a nearly seamless ET a map.For this we used average values for the overlapped area.Finally, the study area was extracted from the mosaicked images using a boundary shape file of the river basin.

Field Validation
Modeled ET a was compared with the eddy covariance measurements carried out at the different eddy covariance measurement sites (Table 2).These sires are Flagstaff managed forest [38], Flagstaff unmanaged forest [38], Flagstaff wildfire [38], Santa Rita Creosote [39], Santa Rita Mesquite [40], Kendall grassland [41], and Charleston Mesquite [42].For comparison, we used the ET a values of pixels that have a flux tower without any knowledge of the spatial extent and relative importance of upwind source areas (footprint).Because all of these sites except for the Charleston Mesquite site are part of FLUXNET [43], we have included additional details and discussion on the Charleston Mesquite site.The data used from the FLUXNET sites are Level 2 data, so they have data gaps ranging from about 7% (Flagstaff wildfire site) to 22% (Santa Rita Creosote site).The data collection and methodology for FLUXNET sites were according to Ameriflux protocol described in the references given in Table 2.The Charleston Mesquite site is located at an elevation of 1200 m, about 16 km northeast of Sierra Vista, Arizona, on the east side of the San Pedro River.This site is dense riparian woodland dominated by velvet mesquite (Prosopis velutina), while the understory is composed of sacaton grass (Sporobolus wrightii), greythorn shrubs (Zizyphus obtusifolia), and other summer active annual herbaceous species.The average canopy cover is about 70% with leaf area index ranging from about 1.2 prior to leaf-out to about 2 during most of the growing season.A three-dimensional sonic anemometer (Model CSAT3, Campbell Scientific Inc., Logan, UT, USA) and an open path infrared gas analyzer (IRGA; Model LI-7500, LI-COR Inc., Lincoln, NE, USA) were used at the site.The variables were sampled at 10 Hz by a datalogger (CR5000, Campbell Scientific Inc., Logan, UT, USA) and averaged over 30-min.The average depth to groundwater was about 10 m.Other details of the Charleston Mesquite site and measurement techniques are described elsewhere [42,44].

Annual ET a of Different Land Use/Land Cover
Spatial distribution of the first ever annual ET a of the CRB at the Landsat scale showed different water use patterns within the basin (Figure 4).As expected, open water had the highest mean annual water loss (994 mm) followed by cropland (538 mm) as a result of evapotranspiration (Table 3).There was a wide spatial variation (high standard deviation) in annual ET a of different land cover classes (Table 3).The standard deviation of mean annual ET a ranged from 144 mm (grassland) to 699 mm (water).A wide range of mean annual ET a may be attributed to many variables including thematic accuracy of the NLCD 2006 map, changes in land use from 2006 to 2010, topography, cropping system, cropping intensity, irrigation, local climatic conditions, and Landsat scene availability for ET a mapping.In general, thematic accuracy of the Level II NLCD map is less than 80% [45].This means many land cover classes may be misclassified, resulting in a wide range of annual ET a .However, the classification errors aggregated over a larger area do not result in significant error [46].Cultivation of different crops will also result in a wide variation of actual ET as different crops have different water requirements based on plant physiology and crop duration.Similarly, growing more than one crop per year will result in higher ET a than growing only one crop per year.There will also be a difference in ET a of the same crop under irrigated and non-irrigated conditions because irrigation leads to increased ET a .A study on the effect of irrigation on the water and energy balance of the CRB has shown decreased stream flow and surface temperature [47].The climatic gradient within the basin due to its large size and relief will also cause differences in the water use pattern.
As an example, the histogram of annual ET a for some of the land cover classes for path 37/row 37 (Figure 3) is presented in Figure 5.A majority (76%) of the barren land pixels had an annual ET a within 100-200 mm, and very few pixels (4%) had an annual ET a greater than 500 mm.Similarly, about half of the grassland (47%) had an annual ET a between 200 and 300 mm and less than 5% exceeded 500 mm.In contrast, annual ET a was greater than 800 mm in the majority of cropland.Some of the cropland was fallow based on visual interpretation during 2010 resulting in an annual ET a similar to grassland and barren land.In addition, the effect of misclassification in NLCD 2006 and/or change of land cover may be attributed to low annual ET a of cropland and high annual ET a of barren land and grassland.Westenburg et al. (2006) [48] collected ET a data from three sites in Havasu National Wildlife Refuge during 2002-2004 using Bowen ratio stations.These sites included a large area of medium-to-high density, homogeneous saltcedar (Tamarix spp.); an area of medium density mixed vegetation; and a homogeneous area of low-to-medium density arrowweed.They reported annual ET a for 2003 for saltcedar, mixed vegetation, and arrowweed as 1,076 mm, 728 mm, and 716 mm, respectively.Our results for these sites showed annual ET a for 2010 as 1,291 mm, 1,089 mm, and 1,333 mm, respectively.It appears that ET a increased during these 7 years particularly for mixed vegetation and arrowweed.Doody et al. (2011) [49] listed ET a of some of the cover types based on reviewed literatures.They reported annual ET a of arrowweed as 1,370-1,590 mm, very close to our result for arrowweed.They found a wide range of annual ET a for saltcedar, ranging from 220 to 1,460 mm with a mean value of 765 ± 413 mm.They reported annual ET a for bare soil as 307 mm and for open water as 1,156 mm.Thus, overall our results for these locations and cover types are similar and comparable.2002) with a mean ET a of 851 mm.Our zonal analysis of the same region showed mean annual ET a of 865 mm.Though we did not use exactly the same areal extent because our delineation slightly differed from theirs, the comparable mean ET a shows that our results are reliable.Even though there is a wide range in the standard deviation of different land cover classes, the mean annual ET a values are comparable to some of the measured and reported ET a in the CRB.

Validation with Eddy Covariance Measurement
Figure 6 shows the comparison of annual ET a measured by eddy covariance against Landsat annual ET a for all seven validation sites within the study area.Overall, the SSEBop-estimated ET a was in good agreement (R 2 = 0.78) with the measured ET a for a wide range of elevations (991-2,270 m) and tower heights (4-23 m) (Table 2).The mean difference between measured and estimated annual ET a was about 10%, which is comparable with results from other energy balance based models [5,17].The root mean square error (RMSE) was 106 mm, mainly due to the high discrepancy at the Flagstaff unmanaged forest and Flagstaff wildfire sites (Figure 6).If these two sites are removed, then R 2 improves to 0.95 (Y = 0.997 X − 53.767) and RMSE decreases to 68 mm.This is remarkable considering we had only 6-9 cloud-free Landsat scenes for 2010 (Table 2) and we had not carried out any calibration of our model using eddy covariance measurements.It should be noted that we have compared the measured ET a only with the ET of the pixel that has a flux tower.However, the point distribution in Figure 6 shows that the SSEBop model captured the ET variability well even without footprint analysis.In general, the footprint of a point of observation is affected by thermal stability, surface roughness, observation levels, and wind speed and direction [51].Since the wind direction is variable throughout the year, an exact footprint analysis is simply impossible to carryout [52].The advantage of having a remote sensing based ET map is that we get spatial and temporal distribution of ET, thus not requiring any footprint analysis.For example, Figure 7 shows the ET values for 5 × 5 Landsat pixels centered over the eddy covariance tower.We can see there is variability in ET values even within 150 m × 150 m (5 × 5 pixels), an area smaller than a generally used 100-to-1 fetch-to height ratio for sites (considering the tower height given in Table 2).Based on measured ET values, we can get a good sense of footprint direction and our results are in reasonable agreement.For example, wind direction at the Flagstaff unmanaged forest site was mostly from the southwest, where pixels had lower ET a values (Figure 7).
Though we have not carried out uncertainty analysis in the present study, like any other models, our model result may have some uncertainty.There are also systematic and random errors associated with the eddy covariance measurements [53][54][55].The total annual precipitation and ET a measured at the Charleston Mesquite site during 2010 were 266 mm and 673 mm, respectively.Our SSEBop model-estimated ET a at the site was 730 mm, which is only about 8% higher than the measured value.At this site, there is a lack of energy balance closure by about 11% (on a daily basis) in the eddy covariance measurements, and this is generally the case at most eddy covariance sites.Wilson et al. (2002) [54] evaluated eddy covariance measurements across 22 FLUXNET sites (50 site-years) and reported that the mean imbalance between turbulent energy fluxes (sensible and latent heat fluxes) and available energy (net radiation minus soil heat flux) is on the order of 20%.For forced energy balance closure we scaled the latent and sensible heat fluxes conserving the measured Bowen ratio.The scaled up ET a for energy balance closure at the site was 761 mm.Thus, there was only 4% difference between the scaled up measured ET a and the modeled ET a , giving us good confidence in the annual value at the pixel level.Our result is within the range of general ET a estimation accuracy reported for seasonal ET a using remote sensing models.
Jung et al. (2009) [35] used eddy covariance data from FLUXNET and developed a global, spatially and temporally explicit ET a map (gridded flux data) using model tree ensemble, a machine learning technique.We used this gridded flux data (written communication 2013) to compare our result.The mean annual ET a of the CRB for 2010 using gridded flux data was 292 mm.Our result showed that the mean annual ET a of the CRB is 312 mm.These values are comparable (within 7%) in spite of very coarse (50 km) spatial resolution of gridded flux data against this high resolution (30 m) Landsat data.

Comparison with Annual ET a Based on Water Balance Analysis
The Landsat-based mean annual ET a of HUC8 ranged from 144 mm (Upper Green-Slate in Upper Colorado basin) to 613 mm (Lower Salt in Lower Colorado basin), while mean annual ET a based on water balance ranged from 145 mm (Yuma desert in Lower Colorado basin) to 629 mm (Tonto sub-basin in Lower Colorado basin).Comparison of Landsat-based mean annual ET a with that of water balance-based mean annual ET a at the HUC8 level has shown that some HUC8 had very good correspondence while other HUC8 were scattered (Figure 8).This finding is expected as ET a based on water balance assumes mass balance closure at the annual scale.In reality, basin runoff exhibits strong memory from year to year and water balance may not close even over a year due to time lag.Overall, there was poor correlation (R 2 = 0.16) between annual ET a based on Landsat and annual ET a based on water balance.Most of the mismatch between Landsat and water balance-based ET a is in sub-basins having large irrigated cropland or water bodies or where vegetation has access to groundwater thus invalidating the premise of water balance closure.We used three sub-criteria to remove some HUC8 sub-basins where water balance may not close because (1) mean runoff is more than mean precipitation, i.e., sub-basins where base flow is dominant; (2) mean ET a is more than mean precipitation, i.e., sub-basins where mostly irrigated or open water bodies are dominant; and (3) the runoff ratio (Q/P) is greater than 0.55, indicating large seepage loss/base flow/percolation loss.Out of 144 sub-basins within the CRB, 28 sub-basins were removed from the analysis based on these three criteria.This resulted in significant improvement of R 2 from 0.16 to 0.85, and RMSE decreased from 120 mm to 99 mm.The regression line shows that annual ET a based on water balance was about 12% higher than Landsat-based ET a .This is reasonable considering that water balance-based ET a is computed neglecting the storage term in the water balance equation thus overestimating annual ET a .In general, the accuracy of ET a products at a higher HUC level (HUC8) is lower than the accuracy at a lower HUC level (HUC2) [56].Our results show that even at the HUC8 level, Landsat-based ET a explained about 85% variability in ET a based on water balance.Compared to mean annual ET a , there was wide range in areal extent of sub-basins within the CRB, ranging from 381 km 2 (Cloverdale in Lower Colorado basin) to 11,845 km 2 (Chaco in Upper Colorado basin).A sub-basin having smaller mean ET a but larger area can have larger volume of water lost due to ET a .Thus, spatial distribution of ET a based on remotely sensed data at the sub-basin level can help water managers identify the water surplus and water deficit sub-basins and quantify the water availability.

Comparison with MODIS-Based Annual ET a
We also compared our Landsat-estimated ET with that obtained from MODIS-based annual ET.[56].Spatial distribution of annual ET a using 2010 MODIS images has shown consumptive water use distribution within the CRB (Figure 9).In general, we can see that the ET a distribution pattern using MODIS and Landsat images is similar.Our comparison of Landsat-based annual ET a and MODIS-based annual ET a for all 144 HUCs (level 8) showed that there was good agreement (Figure 10).The coefficient of determination (R 2 = 0.79) indicated that we can obtain a comparable ET a estimate even with the coarse resolution (1 km) of MODIS at the sub-basin (HUC8) level.However, the advantages of Landsat-based annual ET a become obvious when ET a is required at fine resolution (field scale).Figure 11 shows the annual ET a distribution in a small part of Colorado using Landsat and MODIS images.The ET a variation within individual center-pivot irrigation fields is clearly visible in the Landsat-based ET a estimate, but there is no distinction between fields in the MODIS-based ET a estimates.This is of high importance particularly for implementing water management decisions at the local level.For example, water conservation measures can be adopted in the fields that have higher ET a than surrounding fields.Besides improving efficiency, this can save energy and water that can be used elsewhere.As tackling water challenges is one of the highest priority initiatives under the DOI WaterSMART Program, ET a estimates at the field scale will help in water conservation, availability, and water conflict resolution.As modelers, we recognize that model development is a continuous process.As with any other models, our model has continued evolving since its first publication [14].With the availability of better data sources, our model output will also improve as a result of using better input data.For example, the use of Daymet (Thornton et al., 1997) [57], daily climatological data is expected to give better results than monthly PRISM air temperature data.Similarly, North American Land Data Assimilation System (NLDAS) (Mitchell et al., 2004) [58], which is a quality-controlled and spatially and temporally consistent land surface model dataset (1/8th degree spatial resolution), may improve ET o and thus ET a estimates.

Conclusions
Accurate information on water availability and water use is necessary for planning sustainable use of water, particularly in an arid region like the southwestern United States.Though detailed information on water availability in the Colorado River Basin is available, there is a need for reliable estimation of spatially distributed water use.We have developed a basin-wide ET a map of the Colorado River Basin using images from Landsat, the longest existing polar orbiting satellite.The first ever ET a map of the Colorado River Basin at Landsat scale has shown the distribution pattern of ET a , which is based on the scientifically robust and operationally viable SSEBop model.We have demonstrated that, in spite of the complexities of the CRB, seamless ET a based on remote sensing can be estimated at the Landsat scale.Though cropland covers less than 1% of the CRB, it has the highest ET a after open water bodies.There was a wide variation in ET a within most land use/land cover classes in the basin due to hydro-climatic variation within the basin and uncertainty in the land cover classification.Validation of Landsat-based estimated ET a with ET a measured using eddy covariance has shown that there was good correspondence (R 2 = 0.78) and the model captured the variability of ET a over a wide range of elevation and tower height.The R 2 improved to 0.95 with a slope of 0.997 when two sites (FUF and FWF) were removed from validation.The SSEBop model-estimated ET a was within 4% of the measured ET a at two sites (SRC and CMS).Comparison with ET a based on water balance at the HUC8 level showed good agreement (R 2 = 0.85) for sub-basins most likely having water balance closure.ET a maps using Landsat and MODIS images showed a similar water use pattern at the regional scale.However, the advantages of the Landsat images were clearly observed when the ET a map was zoomed in for field-level ET a distribution.In contrast with MODIS-based ET a maps, individual fields are clearly distinguishable in Landsat-based ET a maps.Thus, Landsat-based ET a mapping is better suited for planning and management of water resources at the field scale.
ET a based on remote sensing may have biases, but these biases can be easily identified and removed based on validation using reliable field data.Once the bias is removed, consistent ET a maps can be generated for regular ET a mapping and water resources planning and management.A regular production of Landsat-based ET a maps using a uniform method and data source at a seasonal and annual scale has the potential for resolving inconsistencies arising from using different methods and variable data sources.With further validation using multiple years and sites, our methodology can be applied for larger areas such as the conterminous United States.

Figure 1 .
Figure 1.Location of the study area showing the Colorado River Basin.

Figure 2 .
Figure 2. Land use/land cover distribution in the Colorado River Basin based on National Land Cover Database 2006 [24].

Figure 3 .
Figure 3. Coverage of Landsat TM/ETM+ images (path/row) in the Colorado River Basin based on the worldwide reference system (WRS-2).

Figure 4 .
Figure 4. Spatial distribution of evapotranspiration in the Colorado River Basin for 2010 using Landsat images.

Figure 5 .
Figure 5. Histogram of the annual evapotranspiration for selected land cover classes (a) Barren land; (b) Grassland; and (c) Cropland within the Landsat path 37/row37.
Senay et al. (2013) [16] used the SSEBop model with the MODIS images for 2005, and their monthly ET validation with 45 Ameriflux data sites for the conterminous United States resulted in an R 2 of 0.64 with RMSE of 27 mm.With more years (2000-2007), MODIS-based annual ET a using the SSEBop model has resulted in a high skill score across all the climate zones in the western United States

Figure 9 .
Figure 9. Spatial distribution of evapotranspiration in the Colorado River Basin for 2010 using MODIS images.

Figure 10 .
Figure 10.Comparison of mean annual evapotranspiration of sub-basins (HUC8) using Landsat and MODIS images.

Figure 11 .
Figure 11.Example showing the spatial distribution of annual evapotranspiration map of the same area within the Colorado River Basin with (a) zoomed in view of the Landsat-based evapotranspiration map for the area in the black square; and (b) zoomed in view of the MODIS-based evapotranspiration map in the black square area.Legend and scale shown at the bottom are for images Figure 11a,b.

Table 1 .
[24] use/land cover distribution within the study area based on National Land Cover Database[24].

Table 2 .
Details of eddy covariance tower validation sites used in this study.

Table 3 .
[24]al evapotranspiration (mm) during 2010 for different National Land Cover Database[24]land use/land cover within the study area.