Quantification and Mapping of Satellite Driven Surface Energy Balance Fluxes in Semi-Arid to Arid Inter-Mountain Region

Crop evapotranspiration (ETc) estimates, on a regional scale, hold enormous potential in managing surface and groundwater resources. This is particularly important for the headwater state of Wyoming, which provides water to found major river basins of the US. In this study, METRIC (Mapping evapotranspiration at high resolution with internalized calibration), a satellite-based image processing model, was used to map and quantify daily, monthly, and seasonal ETc and other energy balance fluxes, i.e., net radiation (Rn), sensible heat (H), and soil heat flux (G) dynamics for different land-use classes. Monthly and seasonal ETc estimated were further used to approximate regional water consumption patterns for different land-use types for nine irrigation districts in semi-arid to arid intermountain region of Big Horn Basin (BHB), Wyoming. The validation of METRIC retrievals against Bowen ratio energy balance system (BREBS) fluxes measured over three vegetative surfaces, viz. sugar beet in 2017, dry bean in 2018, and barley in 2019, indicated high accuracy. The pooled correlation observed between estimated (pooled) and measured instantaneous fluxes had R2 values of 0.91 (RMSE = 0.08 mm h−1, NSE = 0.91), 0.81 (RMSE = 49.6 Wm−2, NSE = 0.67), 0.53 (RMSE = 27.1 Wm−2, NSE = 0.53), and 0.86 (RMSE = 59.2 Wm−2, NSE = 0.84) for ETc, Rn, G, and H, respectively. The biggest discrepancy between measured and estimated monthly ETc values was observed during times when BREBS flux tower footprint was devoid of any crops or the crops at footprint were not actively transpiring. Validation results improved when comparisons were made on monthly scales with METRIC underestimating growing season ETc in the range between 3.2% to 6.0%. Seasonal ETc by land-use type showed significant variation over the study area where crop ETc was 52% higher than natural vegetation ETc. Furthermore, it was found that, in the arid to semi-arid intermountain region of Wyoming, the contribution of irrigation to total seasonal ETc varied in the range of 73–81% in nine irrigation districts that fall within the study area. The high relative contribution of irrigation highlights the importance of identifying and quantifying ETc for improved management in irrigation system design and water allocation.


Introduction
The availability of freshwater resources for agroecosystems has been an important issue for the sustainability of agricultural production in the U.S. and around the world. This is particularly important to the arid to semi-arid regions of the western U.S.A, including Wyoming, where climatic multispectral data from visible and thermal sensors onboard different satellites to identify the differences in land surface characteristics and to estimate ET c .
In this study, METRIC (mapping evapotranspiration at high resolution using internalized calibration), a satellite-based surface energy balance algorithm formulated by Allen et al. [19,20], based on the (surface energy balance algorithm for land) SEBAL model [15,[21][22][23][24], was used to evaluate ET c and other energy balance flux, viz. Rn, H, and G. Allen et al. [19] report the accuracy of METRIC to be heavily dependent upon choosing appropriate hot and cold pixel points (anchor pixels) within a given image during the calibration process, also know as calibration using inverse modeling at extreme conditions (CIMES). In this process, it uses an internal calibration procedure to estimate the near-surface temperature gradient (dT function) as an indexed function of radiometric surface temperature (Ts) for estimating the H component of the surface energy balance. This procedure further eliminates the need for atmospheric correction of Ts and albedo. Generally, a cold pixel is a well-irrigated field with full crop cover, transpiring actively and a dry pixel is a dry field with approximately zero ET c . The selection of these extreme pixels is extremely important and requires careful consideration of maps of Ts, normalized difference vegetative index (NDVI), leaf ares index (LAI), surface albedo, land use, and image in true and false colors. Moreover, this procedure uses the hourly reference evapotranspiration (ET r ) to auto-calibrate the H calculation for each satellite image, which makes ET c estimates more robust [20,25].
Since its introduction, METRIC has been applied worldwide to evaluate and estimate spatial and temporal variability of ET c under different vegetative and climatic conditions at different spatial and temporal scales. In the US, METRIC has been extensively used to manage water use in cropland, riparian and forest vegetation [19,[26][27][28][29][30], regulate water allocation in a river basin and water rights [9], and regulate groundwater withdrawal [31]. Outside the U.S., METRIC has been applied to several regions of the world [26,29,32,33]. For example, Allen et.al. [19] applied METRIC with Landsat images in Idaho, US to predict monthly and season ETc and reported the difference between METRIC estimated and lysimeter measured growing season ET c to be less than 5%. They further used METRIC monthly and seasonal regional ETc estimates to quantify the net ground-water pumpage and total depletion of water from Bear River systems and Snake River aquifer systems to regulate water rights [9]. Likewise, Singh & Senay [34] tested the METRIC model using Landsat images over three Ameriflux cropland sites in the Mid-western United States. The daily ET c from three flux sites had an average R 2 of 0.92 and an average RMSE of 0.93 mm d −1 . Singh et al. [35] used the METRIC model in conjunction with the wet METRIC (a modified version of METRIC) model for analyzing Landsat images in the mid-western United States. They reported that the model-estimated daily ET c was in good agreement with eddy covariance measurements with an average R 2 of 0.91 and a standard error of 0.6 mm day −1 . Furthermore, Mokhtari et al. [36] reported that monthly ET c estimated using the METRIC model falls within ±10% of the actual ET c . Tang et al. [37] concluded that METRIC estimated instantaneous and daily ET c values fall within 10-15% of the actual ET c . They also mentioned that METRIC estimated seasonal ET c was overestimated as compared to that of flux measurements. The METRIC model also showed a promising result in the study conducted on apple and olive orchards [38][39][40].
Unlike its applications over flat cropland areas, only limited studies have been focused on its application to mountainous terrain where there is significant relief and a wide range of slopes and aspects. Thus, a local modification of the algorithms is required. This includes the adjustment of the solar and thermal radiation estimates for slope, aspect, and elevation. Moreover, wind speed and aerodynamic roughness are adjusted to account for the impacts of drag due to undulating topography [19,33,41]. Therefore, the main goals of this study were: to (i) assess the performance of the METRIC algorithm to estimate ET c and other energy balance fluxes using Landsat imagery with respect to measured surface energy balance variables from Bowen ratio energy balance system for the different vegetative surface in the intermountain region of Wyoming; (ii) quantify, map, and evaluate spatial and temporal distribution (daily, monthly and seasonal) of ET c and other energy balance fluxes over the study region; and (iii) quantify the seasonal crop water use and percent contribution of irrigation toward ET c for different vegetative surfaces and different irrigation districts within the study area.

Study Area
The study area covers the arid to semi-arid region of the Big Horn Basin (BHB), Wyoming which lies in the coverage of Landsat Path: 37 and Row: 29 centered at 44.40 • N, 108.52 • W (Figure 1). Situated in the Rocky Mountain region of the Western United States, the BHB is dominated by mountain ranges and high elevation plains and covers the highly developed irrigation-based agricultural system within nine irrigation districts which are fed by three major rivers in the regions, namely the Shoshone River, Big Horn River, and Greybull River (Figure 1). The surface elevation of the study region varies from 1110 m to 3254 m above mean sea level. The major soil type in the study area is sandy loam, with a total study area of 7930 km 2 . Annual average P over the study region can vary from over 639 mm in the high elevation mountainous regions to as low as 142 mm in the lower elevation agricultural plains, with an average of 235 mm compared to the state and national average (Contiguous US) of 406 and 762 mm, respectively [42]. In the inter-mountain region of Wyoming, the growing season is generally short due to a limited number of frost-free days in the year. Despite the limited season, the region provides an ideal environment to grow high-yielding and high-quality crops. Natural vegetation covers about 83% of our study area whereas the cropland area is limited to about 14%. Alfalfa is the major crop grown in terms of total acreage, followed by barley, sugar beet, dry bean, and maize, primarily under furrow irrigation with many growers switching to center pivot irrigation [43]. The growing season in the study area starts with the sowing of barley beginning mid-March, sugar beet in Mid April, maize and dry bean in mid to late May, and ends with harvesting that usually starts from August for barley, early September for dry bean, mid to late September for sugar beet until mid-November for Maize. Because of relatively low growing season precipitation (P) and low relative humidity (RH), disease incidence is minimal, and cool winter temperatures suppress insect and disease cycles and much acreage of the region is dedicated to producing certified seed, including dry beans. The natural vegetation comprises evergreen and deciduous forest, grassland, shrubland, and woody and herbaceous wetlands (Figure 1).
High-quality hourly and daily weather data including maximum and minimum air temperature (T max and T min ), maximum and minimum relative humidity (RH max and RH min ), wind speed at 2-m height (u 2 ), incoming solar radiation (Rs), and precipitation (P) were collected from the Wyoming Agricultural Climate Network (WACNet) weather station at the University of Wyoming Powell Research and Extension Center (PREC), located within the Landsat scene footprint [43]. The measured weather data were used to estimate alfalfa reference ET (ET r ), which is utilized in METRIC algorithm to internally calibrate (hourly ET r ) and upscale (daily ET r ) instantaneous ET c to daily and periodic ET c . Daily precipitation data are helpful to estimate bare soil evaporation and evaluate prior Landsat image-acquisition date precipitation. The ground-based windspeed values are used to calculate wind speed at blending height (200 m) needed during H calculation. Air temperature and RH values are utilized to calculate dew point temperature (T dew ) and vapor pressure which is then used to compute atmospheric transmissivity. Atmospheric transmissivity signifies the ability of Earth's atmosphere to allow/pass incoming radiation from the sun and is required during calculation of R n . Table 1 outlines the average meteorological conditions at PREC on the Landsat image acquisition date. Table 1. Image acquisition date, spacecraft ID, and daily weather parameters measured at BREBS station, Powell Research and Extension Center, Wyoming. Daily weather parameters include minimum air temperature (T min ), maximum air temperature (T max ), minimum relative humidity (RH min ), maximum relative humidity (RH max ), wind speed (u 2 ), solar radiation (R s ), precipitation (P), and ASCE PM alfalfa-reference evapotranspiration (ETr). Assessing the uncertainties in satellite estimated ETc using a ground-based measured ETc is an important and challenging task. While the objective of this study is to assess the performance of METRIC estimated ETc for the complex mountainous terrain, however, no ground-measured fluxes were available for those areas. In this study, the accuracy of METRIC estimated surface energy balance fluxes were assessed using the field measured dataset for different vegetation from Bowen ratio energy balance system (BREBS) located at the center of the study area at PREC (44.46 • N, 108.45 • W) in a large center pivot irrigation production field (14 ha) with a fetch distance of 180 m, 230 m, 223 m, and 160 m in the north, east, south, and west directions, respectively ( Figure 2). The BREBS station was installed on a center-pivot field in June 2017 and continuously sampled the energy flux components and other microclimate variables every 60 s and averaged over 15-mi scales during both crop growing and non-growing seasons using a CR10X datalogger and AM416 relay multiplexer (Campbell Scientific, Inc., Logan, Utah). The BREBS measures ET as latent heat (LE; also known as ET c ) based on the energy balance equation as: where β a is the Bowen ratio, which is defined as the ratio of sensible heat flux to latent heat flux [44]. The Bowen ratio is calculated from the differences in air temperature (T air ) and vapor pressure measured at two heights above the crop canopy. The systems measure T air and RH using two platinum resistance thermometers and monolithic capacitive humidity sensors (REBS models THP04015 and THP04016, respectively) with resolutions of 0.0055 • C for temperature and 0.033% for relative humidity. The T air and RH sensors are located at the center of two coaxial anodized aluminum radiation shield at two heights. It is important to note that a small bias in temperature and relative humidity measurements can cause incorrect estimation of gradients which can cause significant error in β and LE. Therefore, to minimize the impact of any bias, the system automatically exchanges the temperature and humidity sensors every 15 min at two heights above the canopy. Additionally, aluminum housing units for T air and RH sensors face the north direction to avoid any impact of direct sunlight on T air and RH measurements. The system is designed to measure all the components of radiation balance, i.e., incoming and outgoing shortwave and longwave radiation, from which R n is calculated. The R n is measured using a REBS Q*7.1 net radiometer that is installed 3 m above the crop canopy using the tripod assembly. To eliminate the effect of reflection of heat and radiation other instruments installed on the tripod and to capture the R n at the vegetation surface all the R n measurements are taken away from the tripod station. Other radiation balance components, i.e., incoming and outgoing shortwave and longwave radiation are measured using the REBS models THRDS7.1 double-sided total hemispherical radiometer. The P is recorded using a model TR-525 tipping bucket rain gauge with 0.1 mm resolution. Wind speed and direction above the canopy at 3 m height are measured using a model 034B cup anemometer. The G is measured using three REBS HFT-3.1 heat flux plates and three soil thermocouples installed in close proximity, 0.06 m below the soil surface under the net radiometer assembly. One soil temperature and one SMP1R surface soil moisture sensor were also installed nearby of heat flux plated to G correction. Detailed information on BREBS functioning and instrumentation is provided by Irmak [5]. In this study, the estimated fluxes from METRIC model were    (Table 1). Landsat satellite imagery was preferred over other contemporary satellites because it collects thermal data and has a temporal and spatial resolution of 16 days and 30 m by 30 m, respectively making it easier to map individual agricultural fields. The scan line correction for Landsat 7 was carried out on final instantaneous images using the natural neighborhood interpolation tool of ArcGIS 10.6.1. This tool finds the closest subsets of missing pixels from the pixels surrounding the query point to interpolate a value. It is important to note that the pixels corresponding to the BREBS station location did not have any scan line error. The land use map/cropland data layer (CDL) of the study area for 2017, 2018, and 2019 were obtained from the USDA National Agriculture Statistical Service (NASS). The CDL is a georeferenced raster layer created using moderate resolution satellite imagery and extensive agricultural ground truth. METRIC optionally uses a land-use map to improve the parameterization and estimation of the roughness length for momentum transfer [19,20]. Roughness length for momentum transfer is the height above the plant canopy at which windspeed extrapolates to zero because of form drag of air as well as skin layer friction in between air molecules and surface terrain. The CDL layers were also used to extract water use information of different vegetative surfaces for different irrigation districts in the study regions. Irrigation district boundaries shapefiles were collected from the Wyoming Water Resources Data System (WRDS) and State Climate Office (SCO). The digital elevation model (DEM) provides data on the relief of the earth's surface in the form of points of elevation. METRIC requires DEM in mountainous terrain where parameters like temperature, pressure, and solar elevation angle vary by slope and aspect. These data were acquired from USDA geospatial data gateway with spatial resolution same as that of Landsat image.

Satellite Datasets and Image Processing
ERDAS imagine (Leica Geosystems Geospatial Imaging, LLC, Englewood, CO, USA) graphical model maker tool was used for geospatial processing of Landsat images. Firstly, the digital numbers (DN) of the images were converted to top of atmosphere radiance and reflectance values using the procedure outline in Landsat 7 and 8 data users' handbook. Normalized difference vegetation index (NDVI) and soil adjusted vegetation index (SAVI) were computed using reflectance from red and near-infrared bands [45,46]. Leaf area index (LAI) was computed using the empirical equation provided by Bastiaanssen [15], which utilizes SAVI value. Likewise, the solar incidence angle was computed for each pixel of the image using Duffie and Beckman [47] model that considers slope and aspect for each pixel. The slope and aspect of each pixel were computed using DEM. Surface albedo was calculated by integrating the at-surface band reflectance using weighting coefficients [48][49][50]. Surface temperature (Ts) was computed using a modified Planck equation with atmospheric and surface emissivity correction [20]. A neutral-stability-condition atmosphere lapse rate at 6.5 K km −1 in flat terrain and 10 K km −1 in mountainous terrain was considered while calculating DEM corrected Surface temperature (T s-dem ). Likewise, lapse rate change of 2000 m was considered keeping in mind the elevation data of the study area. Moreover, momentum roughness (Z om ), wind speed, and atmospheric pressure were adjusted as mentioned by Allen et al. [20] for use in the mountain model of METRIC.

METRIC Model
METRIC utilizes multispectral raster images to estimate the energy balance at each pixel of an image and is primarily based upon the principle of energy conservation, i.e., energy arriving at the surface must equal the energy leaving the surface for the same period.
where net radiation (R n ) is the total energy available at the earth surface, soil heat flux (G) is the rate of conductive heat exchange between soil and vegetation, sensible heat flux (H) is the rate of convective and conductive heat exchange between plant canopy and the surrounding air due to temperature difference, and Latent heat (LE) is the instantaneous residual energy available for evapotranspiration which is computed by subtracting G and H from R n . Figure 3 shows the flowchart of energy balance components and other intermediate parameters estimation by METRIC.

Net Radiation (R n )
R n is calculated by subtracting incoming and outgoing radiant fluxes at the top-of-the-atmosphere as: where R s ↓ is the incoming shortwave radiation (W m −2 ), R L ↓ is the incoming longwave radiation (W m −2 ), R L ↑ is the emitted outgoing longwave radiation (W m −2 ), α is surface albedo (dimensionless) and ε o is the surface thermal emissivity (dimensionless). The term ((1 − ε o ) R L ↓) in Equation (3) represents a fraction of incoming longwave radiation that is reflected from the earth's surface. These intermediate parameters of long and short-wave radiation are computed as per Allen et al. [20] and are a function of satellite images, albedo, atmospheric transmissivity, DEM, ground-based weather data, and satellite estimated surface temperature. Surface emissivity (ε o ) was computed using LAI in an empirical equation formulated by Tasumi [51].

Soil Heat Flux (G)
G was empirically calculated as a G/R n fraction using vegetation indices and surface temperature [51].
where T s is the surface temperature in Kelvin and LAI (dimensionless) is the leaf area index.

Sensible Heat Flux (H)
H is a function of aerodynamic observations such as u 2 , vegetation type and roughness, and surface to air temperature differences (T s − T a ). H is determined using an one-dimensional aerodynamic function within METRIC [15,20] as: where ρ is the air density (kg m −3 ), C p is the specific heat of the air (1004 J kg −1 K −1 ), dT is the near-surface air temperature difference (T s − T a ) between two reference heights Z 1 and Z 2 , T s is the surface temperature (K), T a is the air temperature (K), and r ah is the aerodynamic resistance to heat transfer (s m −1 ) over the vertical distance. METRIC utilizes the internal calibration procedure at extreme conditions following surface energy balance algorithms for land (SEBAL) to estimate the dT function. The radiometric temperature available in Landsat thermal band is generally different from the aerodynamic temperature of the surface because of the difference in thermal sensor view angle, land use type and cover fraction, and the surface roughness length of heat and momentum [19]. Therefore, to compensate for uncertainties in the estimation of Ts from satellite, dT function is established for two extreme conditions present in the image viz: cold and hot anchor pixels and is estimated assuming a linear relationship between dT and T s-dem [15,19,20,24].
where T s-dem is DEM adjusted surface temperature (T s ), a is the slope, and b is the intercept of the plot between dT and T s-dem for a given satellite image. Generally, a cold pixel is a well-irrigated field with full crop cover transpiring at potential evapotranspiration rate and a dry pixel is a dry field with approximately zero ETc. Hot and cold pixels are the evaporative extremes of the Landsat scene being analyzed and selection of these extreme pixels is extremely important in the METRIC processing. Therefore, to select the appropriate hot and cold anchor pixels, the maps of NDVI, surface albedo, LAI, Ts, land use, and maps of the image in true and false colors were taken into consideration. Figure 4 provides a general overview of the hot and cold pixel selection procedure in the study. Cold pixel candidates were selected from the well-irrigated densely vegetated area with NDVI between 0.76 and 0.84, mid season LAI greater than 3 or 4 m 2 m −2 , surface albedo between 0.18 to 0.24, and a lower surface temperature. Similarly, the hot pixel candidates were selected from the bare agricultural field with NDVI less than 0.20, surface albedo between 0.17 to 0.23, and relatively higher surface temperature. Final selection of the anchor hot and cold pixel for each image was made using the candidate pixels already selected. Tasumi [51] and Tasumi et al. [52] estimated cold pixels to have ET c rates of about 1.05 times the reference evapotranspiration of alfalfa crop (ETr) except during early and late growing season. To account for any residual evaporation and wetting events before image acquisition, METRIC approximates LE at the hot pixel to be 0.1 [20]. A bare soil evaporation model on a daily time step was used to adjust LE hot-pixel h.igher than 0.1 in situations with residual soil evaporation from antecedent rainfall. Moreover, the cold and hot pixels were selected within a 20-25 km radius of the meteorological station to counter heterogeneity in weather variables over the Landsat image. Sensible heat is estimated for both the selected hot and cold pixel using the energy balance equation in respective pixels.
where λ in H cold computation refers to latent heat of vaporization (JKg −1 ). The temperature gradient (dT) is computed for both cold and hot pixels before the final H is calculated.
At hot pixel, dT hot = H hot × r ah−hot ρ hot × C p (11) The aerodynamic resistance to the transfer of the heat from the evaporating surface to the air above the canopy is calculated as given by Bastiaanssen et al. [15].
where z 1 and z 2 are lower and upper reference heights (m) above zero plane displacement (d) of the vegetation and are assigned to be 0.1 m and 2 m respectively [15]. u * is the friction velocity (m s −1 ) computed during the first iteration of H and k is the Von Karman's constant (0.41). The u * at the first iteration, is defined by the relation: where U 200 is the wind speed at the blending height (200 m) above the weather station, z om (m) is the momentum roughness length calculated using a land-use map for each image pixel. z om is calculated as a function of LAI [51].
The sensible heat flux (H) is then calculated for each pixel using Equation (6). However, this event is the first estimation of H and is calculated assuming neutral atmospheric conditions. To counter unstable atmospheric conditions, H is computed through the number of iterations until r ah is stabilized. At each iteration, values for dT cold , dT hot , air temperature, air density, r ah , U * , and H is corrected and is continued until successive values for dT hot and r ah at the hot pixel have stabilized. Usually, stability is achieved after 4-5 iterations. The Monin-Obukhov length (L) is used to stabilize the atmosphere (boundary layer) where turbulence is generated by the buoyancy effect of surface heating. The value of L is zero for neutral, positive for stable, and negative for unstable atmospheric conditions.
where g is the gravitational constant (9.8 ms −2 ).

Estimation of Instantaneous, Daily, Monthly and Seasonal Crop Evapotranspiration (ET c )
Instantaneous ET (ET inst ) at the time of Landsat overpass is computed as residual of the surface energy balance by subtracting instantaneous G and H are from R n as: where ρ w is the density of water (~1000 kgm −3 ). METRIC utilizes reference ET fraction (ET r F) value to compute daily and periodic ET c [20]. ET r F (dimensionless) is calculated for each pixel as a ratio of ET inst to the reference ET values at the time of image acquisition.
where ETr is the alfalfa reference evapotranspiration calculated using the standardized ASCE Penman-Monteith equation [53] on an hourly basis. It is assumed that ET r F (equivalent to crop coefficient) values remain constant throughout the day and the use of it minimizes the impact of advective heat on ETc [54]. The daily evapotranspiration (ET 24 ) rate (mm d −1 ) for each image pixel is then computed as: where ET r-24 is 24-h ET r and is calculated by summing hourly ET r values over the day of image.
To estimate periodic ET c values, both linear and cubic spline methods of interpolation were employed. Linear interpolation [55] was done by interpolating ET r F between successive image dates. The interpolated values were then multiplied by ET r-24 to obtain ET 24, which is then summed up on monthly and seasonal baisis. On the other hand, cubic spline [19] is reported to get a smooth and curvilinear bend as of the typical crop coefficient curve of annual crops [56]. Rather than taking all the data points at once to fit one polynomial, spline interpolation fits a polynomial function between two adjacent data points. For example, 4 images (data points) will have three polynomials fitted in between. A third-degree polynomial (cubic polynomial) is used to get a smooth fit between each spline. Thus, cubic spline interpolation requires two images on either side of the month and at least one ET r F image per month in a growing season. As such, it requires at least 5 images to interpolate for a month. Likewise, synthetic points (ET r F synthetic ) were computed for the beginning and end of the growing season (ET c is dominated by evaporation) to be able to start and end the interpolation. Daily actual evaporation (E a ), estimated using the soil water balance method [10], and ET r were used to calculate the synthetic points.
ET r F synthetic = where ET month (mm) is the monthly cumulative ET c from day 1st to last day of the month, ET r F i is the interpolated ET r F for a month and ET r24i is the 24-h ET r for a month. Seasonal estimates of ET c were then obtained by summing up the monthly ET c values. The METRIC estimated instantaneous and periodic fluxes were compared with BREBS measured fluxes. Root mean square error (RMSE) and Nash-Sutcliffe's efficiency (NSE) were used as the main criterion to judge the accuracy and reliability of the METRIC model for estimating ET c and other surface energy fluxes.
where O and P is the observed and METRIC predicted surface energy balance fluxes, O mean is the mean of observed energy balance fluxes, and n is the number of observations.

Comparison of METRIC Estimated and Measured Surface Energy Balance Fluxes
Several uncertainties may be associated with the satellite retrievals of energy balance systems including sensor pixel scale, sensor accuracy, etc. Thus, a comparison of satellite retrievals with ground-based measurements becomes inevitable to garner confidence on the final output being accurately estimated. Figure 6 Table 2 provides a statistical comparison between METRIC estimated and BREBS measured energy balance components. The R 2 , RMSE, and NSE were used to assess the error associated with each component. The following section discusses the comparison of METRIC estimated instantaneous energy balance fluxes (R n , G, H, and ET c ) of the pixel corresponding to the geographic location of the BREBS towers.  Statistical comparison between METRIC estimated and BREBS measured crop evapotranspiration (ET c , mm h −1 ), net radiation (R n , Wm −2 ), soil heat (G, Wm −2 ), and sensible heat (H, Wm −2 ) for 2017, 2018, 2019, and pooled data from all the growing seasons (N = number of data points, R 2 = coefficient of determination, SD = standard deviation, S = slope, RMSE = root mean square error, NSE = Nash-Sutcliffe's efficiency).

Flux
Year  Figure 6a indicates a good correlation between METRIC estimated and BREBS measured ET inst with R 2 values ranging between 0.21 to 0.95 and slopes between 0.2 to 1.13. For pooled data points, the linear regression model explains 91% of all the variance in the estimated ET inst values. The slope of the regression coefficient of pooled data (0.91) close to unity (p-value < 0.05) indicates a strong fit and low systematic error between modeled and measured ET inst . The average ET inst for all data points over three vegetative surfaces was 0.52 mm h −1 (355 W m −2 ) as compared to the BREBS measurement of 0.49 mm h −1 (336 W m −2 ), with pooled RMSE of 0.08 mm h −1 and NSE of 0.90. Table 2 represents a comparison of METRIC estimate ET inst with BREBS measurement for individual surfaces. Overall, a good correlation was observed between measured and estimated ET c for all three vegetative surfaces, with R 2 ranging from 0.21 for sugar beet in 2017 to 0.95 for dry bean in 2018. A relatively lower R 2 (0.21), slope (0.20) and NSE (0.21) for 2017 ET inst indicate low performance of the model ( Table 2). The regression analysis showed that the correlation is not significant and slope of the regression line was significantly different from unity (p-value > 0.05). The BREBS tower was installed on June 2017 which resulted only handful of images (5) (Table 2). Higher overestimation in 2018 is due to the fact that three out of the nine images were analyzed before planting and after harvest when only crop residue was available at the surface. For example, the combined overestimation of 41% was observed for the three images when only crop residue was available at the surface compared to 9.8% when the image was analyzed within the growing season. Variation in crop residue at the surface can alter the surface temperature which may result in more variation in ET c available at the surface [13,57]. A similar improvement in results was reported by Allen et al. [19]. They compared METRIC estimated ET c with lysimeter measured ET c using eight Landsat images acquired from April to September and reported 30% averaged absolute differences for sugar beet crop (R 2 = 0.82). When they omitted one image date of drying bare soil following P, the average absolute difference was only 14%. The maximum difference in ET inst of 0.13 mm h −1 was observed in 2019 on 21 July, which could attribute to the continuous P events from 14 July to 17 July which resulted in increased soil surface moisture and decreased surface temperature, thus resulting in higher bias from the measured data [58].

Net Radiation (R n )
METRIC estimated R n had an R 2 value of 0.83 to 0.96 and a slope of 0.62 to 0.97 when compared with corresponding BREBS fluxes (Figure 6b). The average METRIC estimated R n for the pooled dataset was 507 W m −2 (SD = 88 W m −2 ) compared to BREBS measured R n of 537 W m −2 (SD = 67 W m −2 ). The RMSE value for R n was found to be ranging from 36.2 to 55.3 W m −2 and NSE between 0.70 to 0.95 (Table 2). Pooled data points had R 2 of 0.81, slope of 0.68, RMSE value of 49.6 W m −2 , and NSE of 0.67 (Figure 6b). Table 2 indicates the average R n to be overestimated by 10.5%, 3%, 6.5%, and 5.8% for 2017, 2018, 2019, and pooled data points, respectively. Similar observations were made by Singh et al. [35], who used METRIC model in conjunction with wet-METRIC (a modified version of METRIC) model for analyzing Landsat images in the three sites of Mid-western United States. They reported model-estimated net radiation (R n ) to be within −12%, −11.4%, and −6.3% of the eddy covariance measurement for three different study sites with R 2 of 0.95. R n is a very complex variable and varies with the complexity of the underlying surface heterogeneity. For example, a heterogeneous surface influences the surface radiation properties, e.g., surface albedo and emissivity, and thereby impacts the fraction of photosynthetically active radiation absorbed crop transpiration rate and surface energy budgets. For example, as observed from Figure 6b, the maximum variation in R n was observed on 12 September and 22 October 2018, when only crop residue was available at the surface.
Horton et al. [57] reported that crop residue that alters surface evaporation and soil moisture content can affect both the shortwave albedo and longwave emissivity, and hence R n . Figure 6c and Table 2 represent the comparison of METRIC derived H with BREBS measured H along with statistical analysis. Overall good correlation was observed between METRIC estimated H and BREBS measured instantaneous H values with R 2 between 0.21 to 0.9 while the value for slope was between 0.19 to 0.75 (Figure 6c). When pooled data were taken into consideration, R 2 value was 0.86 with a slope of 0.71. Lower RMSE value was observed in 2018 (46.8 W m −2 ) as compared to 2017 (72.5 W m −2 ) and 2019 (64 W m −2 ) image dates (Table 2). However, when pooled data points were considered, the average METRIC estimated H was 121 W m −2 which is 12.6% higher than field measurements (average H = 107 W m −2 ) with RMSE 59.2 W m −2 and NSE of 0.82 (Table 2) (Figure 6c). This indicates the movement of energy from the air to the plant canopy and additional energy being used by the crop to meet the high ET demand on that day. This may be due to the already harvested barley field surrounding the BREBS sugar beet field that acted as a source for advective heat. In addition, natural vegetation covers about 83% of our study area surrounding the cropland which also acted as a source of advective heat. However, compared to 2017, in 2018 and 2019, the ratio of LE/(R n -G) was observed to be less than 1 on most image dates; this may be due to higher P and lower T air in 2018 and 2019, where the average P was 19 mm and 81 mm higher than 2017 and average T air was 0.50 • C and 0.74 • C lower than 2017, which reduce the effect of advective heat from surrounding areas. On average, H was better estimated during 2018 image dates where average METRIC estimated H differed from BREBS H by −9.1% as compared to 2017 (−86.4%) and 2019 (50.8%) image dates (Table 2). Poor estimation of H in 2017 and 2019 can be associated with the irrigation and P events that occurred before Landsat overpass time. This created difficulty in choosing a suitable hot pixel candidate having LE approximately zero [35,59]. Similar observations were made by Singh et al. [35] who reported estimated and measured sensible heat had predictive error and R 2 value of 48.3% (R 2 = 0.73), 15.2% (R 2 = 0.71), and −2.4% (R 2 = 0.03) for three different sites in the Mid-western United States. Kilic et al. [60] reported RMSE of 46.5 W m −2 and 31.5 W m −2 in H during the 2005 and 2006 growing seasons in Southeast Nebraska, similar to the RSME in this study in 2017.

Soil Heat Flux (G)
G is one of the key components in the surface energy budget used to estimate the available energy to be partitioned into H and LE. A moderate correlation was observed between METRIC estimated G and BREBS measured G with R 2 in between 0.23 to 0.72 and slope in between 0.24 to 0.46 (Figure 6d). Pooled G data points had R 2 of 0.53 and a slope of 0.5 (Figure 6d). METRIC estimated average G differed by 23.9%, 0.75%, −16.1%, and −0.8% for 2017, 2018, 2019, and pooled data, respectively, when compared with BREBS fluxes ( Table 2). The RMSE value for the growing seasons ranged between 15.7 to 30.8 W m −2 while pooled data points had RMSE of 27.1 W m −2 (Table 2). Likewise, NSE value for the growing seasons considered ranged between 0.1 and 0.46 while pooled data points had NSE of 0.53 (Table 2). In general, G is a relatively small component in energy balance, and its estimation via model is usually difficult because of changes occurring in soil properties such as in conductivity, temperature, and heat capacity with a change in moisture and vegetation of soil. Singh et al. [35], in their study in mid-west USA, reported estimated and measured soil heat flux (G) differing by 16.5%, −6.6%, and −30.3% for three different sites. The same study showed a low R 2 value of 0.03, 0.36, and 0.33 for soil heat flux. Similar results were found by Irmak et al. [61], who reported a higher RMSD of 48.6 W m −2 in 2005 on irrigated continuous maize at three sites near Mead, Nebraska. They used a different remote sensing model to predict G at the Mead and Clay Center locations, and the reported difference in G ranged from 36.3 to 62.6 W m −2 due to differences in soil properties, measurement of G (sensor depth and distance between sensors), and surface conditions at the different locations.
It is important to note that there are many sources of uncertainties associated with the comparison of BREBS measured and modeled fluxes including model algorithm (model assumptions), tower observations (systematic and unsystematic bias), scaling issues (flux footprint), and management practices. For example, BREBS assumes the eddy diffusivities of heat and water vapor to be equal. Studies [62,63] have revealed that these diffusitives may not be equal in some cases resulting a force closure of surface energy budget. Likewise, Barr et al. [62] reported that BREBS favoured prediction of LE as compared to H. The potential effect of condensation in the surface energy balance fluxes is small and has not been considered in this research. Condensation is more pronounced during the night time when RH is at 91-99% and the leaves are cooler than the surrounding air and approaches dew point temperature [64]. Careful analysis of average hourly RH data from BREBS indicated RH barely goes above 90% in the actice growing season where most of the staelite images were collected (data not shown). A similar oberevation was made by Sharma [65] in the semi-arid grass land region of New South Wales, Australia where it was reported that dew fall amounted to 8.7% of ET c in winter, 4.7% of ET c in autumn, 1.6% of ET c in spring season, and as low as 1% of actual or pan evaporation in warmer months.

Comparison of METRIC Estimated and BREBS Measured Monthly ET c
METRIC estimated monthly ET c values during 2017 (3 data points), 2018 (6 data points), and 2019 (4 data points) growing seasons were obtained employing linear and cubic spline interpolation method. Cubic spline interpolation observed a better correlation between estimated and measured monthly values with high R 2 ranging from 0.82 to 0.9, slope between 0.83 to 1.32, and RMSE in between 16.8 to 20.8 mm (Figure 7b). On the other hand, linear interpolation had moderate R 2 in between 0.7 to 0.85, slope between 0.68 to 1.1, and RMSE in between 17 to 28.5 mm (Figure 7a). Pooled data (13 data points) for cubic spline had R 2 and slope at 0.9 and RMSE at 18.2 mm, whereas for linear interpolation the values were at 0.84, 0.85 and 23.6 mm for R 2 , slope and RMSE respectively (Figure 7a,b). Cubic spline interpolation underestimated seasonal ET c by 6%, 5.4%, and 3.2%, while linear interpolation underestimated seasonal ET c by 4%, 12%, and 2% during 2017, 2018, and 2019 growing seasons ( Table 3). The biggest discrepancy between measured and estimated monthly ET c values was observed during times when the BREBS flux tower footprint was devoid of any crops or the crops at footprint were not actively transpiring (Table 3). During the 2018 growing season, METRIC estimation of monthly ET c values ranged between −49% (underestimation) in June (dry bean at BREBS footprint planted on 5 June) to 75% (overestimation) in October (dry bean at BREBS footprint harvested on 10 September) ( Table 3). Allen et al. [19], in a study regarding sugar beet, observed estimated periodic ET c differing by −28% on 7 July to 107% on 4 May as compared to lysimeter measured ET c . However, this difference in data decreased to less than 1% when cumulative ET c values for the growing season were compared. Similarly, Allen et al. [19], considering an irrigated meadow, observed METRIC estimated monthly ET c to be ±16% of the lysimeter measurement. However, in the same study, the difference was reduced to 4% when METRIC estimated and lysimeter measured growing season ET c (July to October) was considered. This decrease in the difference between estimated and measured growing season ET c values can be linked with the error associated being randomly distributed (monthly ET c values being both over and underestimated) tending to cancel while computing a cumulative growing season value [19]. However, both interpolation methods performed well, with cubic interpolation performing slightly better than linear interpolation, owing to its lower and consistent RMSE and % error. Therefore, in this study, cubic interpolation estimated ETc values were used to quantify and map the seasonal ETc and for other analysis. Further quantification and mapping of periodic ET c in this study was done using cubic spline interpolated images.

Mapping Spatio-Temporal Variation of Surface Energy Balance Components
Satellite images, unlike traditional and conventional methods of measurement, provide an excellent means to determine and map the spatial and temporal variation of surface energy balance fluxes and vegetative indices. The spatial (space) and temporal (time) variability on a regional scale is a function of topography, weather parameters, climate, soil characteristics, land-use type, and management practices. Figure 8 depicts the early, mid, and late-season Spatio-temporal variability of Rn, H, and G for 2018 crop growing season. In general, the intensity of R s and thus the R n increases gradually after the March equinox and starts decreasing when the September equinox is reached [66]. A similar pattern is observed in our case where a higher R n value is observed from May to August. Clouded parts in the images over estimates Rn value as evidenced in the 8 June image where the Rn climbed up to 980 W m −2 . As more and more energy gets partitioned into ET c , the H value decreases (occasionally negative) as the growing season moves forward. Negative sensible heat estimation occurred occasionally when the temperature and vapor pressure of surrounding air were higher as compared to plant surface temperature and air vapor pressure surrounding the plant, resulting in a flow of heat from surrounding air to the plant. Clouded parts in the image can cause a drastic reduction in the value of H, obvious in 8 June image in Figure 8f, where the H value hovered around −125 W m −2 . Similarly, lower soil heat value was estimated for cropland during mid-season because of higher LAI and NDVI during this time. High ground cover decreases the conduction of R n to the soil surface. The opposite happens during the early season when crops have less effective ground cover. Similarly, a lower soil heat value at the end of the growing season indicated less availability of net radiant energy. The overestimation of G value can be apparent in clouded parts of the images. A higher G value reaching up to 500 W m −2 is from the open water area of the image. Unlike land surface, G is a larger and more important component for water bodies [67].

Spatio-Temporal Variation of Evapotranspiration (ET c )
In METRIC, ET c is estimated as a residual of the surface energy balance fluxes and occurs through energy exchange at vegetation surface via H or R n . In general, ET c is a function of weather parameters, crop characteristics, soil properties, and crop cultivation and management practices. Figure 9 shows the spatial distribution of daily ET c for Landsat overpasses at scattered times in the 2017 and 2018 growing seasons. Over the production crop areas, the ET c rate was low during the beginning of the crop season and increased as the season moved forward. During the late season, as plants approached physiological maturity, the ET c rate declined. On the contrary, higher ET c was observed over the natural vegetation early in the season (May and June) and strats decrasing as the season progressed. During the early season, evaporation dominates the ET c while as the season progresses, transpiration forms the major part of ET c . Most of the variability in ET c was due to the diverse cropping systems and agronomic practices across the study area. Clouded parts of the image, depending upon the cloud thickness, will experience an overestimated ET values, reaching up to 14.8 mm (5 June 2017) in this study.  , respectively. Natural vegetation suffers from water scarcity during the majority of the growing season which leads to a lower ET c rate. On the other hand, P is accompanied by frequent irrigation in the cropland which causes crops to transpire more. A Higher ET c rate was observed in 2017 because of the even distribution of the total P amount (144 mm from March-September) during the growing season. On the other hand, the 2018 growing season experienced uneven distribution of total P amount (158 mm from March to September) where most of the P was observed during May (49 mm) and June (49 mm). A higher cumulative ET c is evident on the conifer forested part (top right and bottom left) of the image. Brown patches (a higher ET c estimation as compared to surrounding vegetation) seen in the naturally vegetated area are due to cloud masses present in some images. However, cloud masking was performed to estimate seasonal ET c for natural vegetation.

Comparison of METRIC Estimated Mean Monthly Evapotranspiration for Different Land Cover Type
Monthly ET c is the summation of the amount of water lost from a cropped surface during each day of the month. Monthly and seasonal ET c estimates are valuable in determining the irrigation, storage, and conveyance capacities of an irrigation system. Figure 11 shows the METRIC estimated mean monthly ET c rate for sixteen different land cover type viz. alfalfa, dry bean, sugar beet, barley, Maize, Oats, other hay/non-alfalfa, sod/grass seed, spring wheat, fallow/idle crop, barren, open water, shrubland, grassland, woody wetlands, and herbaceous wetlands during 2018 growing season (May-September). While estimating periodic ET c rate for each landcover type, the clouded part seen on some of the images were cropped out. The monthly ET c estimates for crops exhibited a bell-shaped ET c curve. This implies monthly ET c rate will be lower during the beginning and end of the growing season and higher during mid-season when plants have good canopy cover to transpire at their potential rates. Early planting date (starts mid-March) for barley means it reached its peak growth during the period from June to July resulting in a higher ETc rate during those months. The monthly ET c values from the barley field start decreasing by August as the crop reaches physiological maturity. The harvesting of barley starts from early August and can last until Mid-August. Monthly ET c values observed from barley fields after harvesting accounts to surface residue ET c and soil evaporation. Dry bean, sugar beet, and maize had the highest ET c rate during July and August which coincides with the peak growth and flowering stage of these crops. Dry bean fields experienced lower monthly ET c values as compared to fields of Maize and sugar beet ( Figure 11). Dry bean tends to have less dense vegetation and a shorter growing span, resulting in a lower rate of cumulative ET c . Unlike dry bean, sugar beet and maize field had considerably higher ET c values until September. Maize and sugar beet have fair vegetation cover until September and subsequently senesce prior to harvest in October. In Wyoming, alfalfa is grown primarily for hay under irrigation. Alfalfa fields have higher ET c values in May as compared to other crops as it is planted by early-August of the previous year and is already established by winter. Likewise, alfalfa has a longer growing season, a deep root system, and a dense canopy of vegetation that results in higher cumulative seasonal ET c values as compared to other crops. Alfalfa is normally harvested 2-3 times a year with the most active harvesting period between June to September that results in fluctuation of monthly ET c values during this period. As expected, the naturally vegetated areas had lower average monthly ET c values as compared to those of crop fields. Similarly, Figure 11 indicates open water evaporation (E water ) to be greater than wetland ET (ET wetland ) in the study area. E water − ET wetland was less during early (May) and late growing season (October). However, the difference increased considerably during the mid-season summer and fall images (June, July, August, and September). This increase in the difference in mid-season can be associated with vegetation senescence (resulting in low ET wetland ) in the wetland region. On the other hand, water has high specific heat and heat capacity that induces it to absorb and store thermal energy. The storage of thermal energy in the water body is greater in summer and fall as more radiant energy is available. The greater release of stored heat in summer and fall from open water increases the net energy available for evaporation. The E water − ET wetland value was much higher for herbaceous wetlands (E water − ET herbaceous-wetland = 33.7%) as compared to woody wetlands (E water − ET woody-wetland = 15%). Stannard et al. [68] used eddy covariance to measure ET wetland rates from two wetland sites and BREBS to measure E water rate from two open water sites for three years and reported overall E water to be 20% greater than overall ET wetland . Their study also reported E water − ET wetland to have a lower difference from late June to early August, while the difference increased during late summer.
The majority of the study area in this study is under natural vegetation (about 83%). A comparison was done to assess METRIC estimated surface energy balance fluxes between cropland and natural vegetation ( Figure 12). For that, average flux across the study area on five different image acquisition dates comprising early (15 May), mid (8 June, 18 July, and 11 August), and late-season (12 September) during the 2018 growing season were considered. Natural vegetation had comparatively lower R n values (up to 15% lower) as compared to cropland as the growing season progressed. The difference in the net radiation values can be associated with the variation in surface albedo value. As the growing season progresses, the bare soil in the naturally vegetated area is dry (because of water stress) and has higher albedo values (dry soil albedo was reported to be around 0.31 by Idso et al. [69]) causing a decline in R n . However, during the early season, the soil in the natural vegetation area is wet (because of possible precipitation and snow cover) with lower albedo values (wet soil albedo was reported to be be around 0.14 by Idso et al. [69]) causing the R n to be similar or sometimes greater than that of cropland, as evidenced in May in our study area. Natural vegetation can also have lower albedo value due to the effect of shading that might happen in taller vegetation like forest and woody wetlands. As compared to cropland, naturally vegetated areas had higher H value (up to 94% higher in mid-season). A huge proportion of R n energy is used up by crops in the ET c process, thus lowering the H and G estimates in cropland. Natural vegetation, on the other hand, have a lower ET c rate, and the net available energy is partitioned more into H and G. Lower H values in cropland during mid-season indicate crops transpiring at a higher rate whereas higher sensible heat values for crops at the end and beginning of growing season indicate lowering ET c rate. Higher ET c rate in the early season lowered the H value whereas lower ET c rate during the mid-season resulted in a higher H. Lower H values at the end of the growing season is caused by the decline in net available radiant energy [58]. G has a reciprocal relationship with high ground cover (shading or insulation effect) and leaf area index and a direct proportional relationship with albedo. These factors cause a significant difference between crop land and natural vegetation G. The decreasing G value at the end of the growing season in the natural vegetation area reflects less availability of R n energy. During the early season, the ET c rate from the cropland can be similar or in some cases less than that of natural vegetation as evidenced during the month of May in this study. The precipitation received early in the growing season helps natural vegetation induce new growth and thus transpire more. However, as the season progresses, the natural vegetation could suffer from water stress, resulting in a lower ET c rate as compared to that of cropland, where P gets accompanied by irrigation.

Irrigation District Average Water Consumption
One of the main objectives of this study is to use satellite estimated ET c to quantify the average crop water use and crop water use for different cropping systems for different irrigation districts within the study area. Table 4 presents the average ET c , P, and percent of irrigation-based ET c for 2017 and 2018 crop growing seasons over nine irrigation districts in the study region. For this analysis, monthly P values for 2017 and 2018 for each irrigation district were extracted from gridded 4-km monthly parameter-elevation regression on independent slope model (PRISM) dataset [70] obtained from the PRISM Climate Group website. Within the study region, the area of each irrigation district varies from 51 Km 2 for the Hunt and Godfrey irrigation district to 1284 Km 2 for the Greybull Valley irrigation district. No significant difference was observed in seasonal crop ET c among the irrigation districts, with average ET c (average of two years) ranging from 598 mm in Greybull Irrigation District to 689 mm in Heart Mountain irrigation district. The seasonal ET c for 2017 was approximately 34 mm greater than 2018, except for the Will Wood irrigation district where greater ET c was observed for 2018. This may be due to higher precipitation P and lower T air in 2018 where average precipitation P was 19 mm higher than 2017 and average T air was 0.50 • C lower than 2017. To understand the amount of ET contributed by irrigation, irrigation volumes for each irrigation district were calculated by subtracting the average seasonal P from ET c (ET c − P) for each irrigation district, assuming other variables such as streamflow and change in soil moisture are negligible. This assumption is valid considering the semi-arid to arid intermountain region of Wyoming, where seasonal precipitation is generally less than 200 mm. The results indicated that, in 2017, an average of 81% (maximum of 83% in Hunt and Godfrey irrigation district and minimum of 79% in Greybull irrigation district) and, in 2018, an average of 73% (maximum of 77% and minimum of 69%) of ET c is contributed by irrigation only. To further understand the dynamics of seasonal ETc from different vegetative surfaces, ET c and percent irrigation contribution to ET c for four major crop types within each irrigation district were analyzed (refer Supplementary Materials, Tables S1-S9). For both years, higher ET c was observed for alfalfa followed by sugar beet, and the lowest ET c was observed for dry bean in all irrigation districts. The percent irrigation contribution to ET c (ET c -P) varies from 60% for dry beans in the Heart Mountain irrigation district in 2018 to a maximum of 83% for alfalfa in Lovell and Deaver irrigation districts in 2017. This information could be very useful for irrigation district personnel, water managers, decision-makers, policymakers, and state and federal water regulatory agencies for quantification/re-evaluation of the overall distribution of evaporative losses for a given irrigation district relative to available water resources for evaluating various management practices.

Summary and Conclusions
Estimates of ET c represent a key component of water management in croplands, naturally vegetated areas, open water sources, wetlands, and riparian basins. Today, determining short term and long term ET c rates to guide water resource-based policy and decision making has become a necessity. In the face of population increment and water scarcity around the globe, the importance of ET c has grown many folds. A total of 22 cloud-free Landsat 7-ETM+ and Landsat 8-OLI and TIRS images were analyzed using a satellite-based energy balance model, namely mapping evapotranspiration with internalized calibration (METRIC). The study was carried out for 2017, 2018, and 2019 growing season period in the semi-arid to arid intermountain region of Wyoming (Path/Row: 37/29). The accuracy check of METRIC estimates was done by comparing estimated instantaneous and periodic fluxes with respective Bowen ratio energy balance system (BREBS) measured fluxes. The BREBS station was installed on a center pivot irrigation field and had sugar beet in 2017, dry beans in 2018, and barley in 2019 growing season. The correlation observed between estimated and measured instantaneous energy balance fluxes (pooled data points) had R 2 values of 0.91 (RMSE = 0.08 mm h −1 , NSE = 0.91), 0.81 (RMSE = 49.6 W m −2 , NSE = 0.67), 0.53 (RMSE = 27.1 W m −2 , NSE = 0.53), and 0.86 (RMSE = 59.2 W m −2 , NSE = 0.84) for ET c , R n , G, and H respectively. Likewise, estimated instantaneous energy balance fluxes (pooled data points considered) differed with BREBS fluxes by 5.7%, 5.8%, −0.8%, and 12.6% for ET c , R n , G, and H, respectively. Cubic spline interpolation performed a bit better as compared to linear interpolation, and a comparison of estimated and measured monthly ET c values (pooled data) for cubic spline had R 2 and slope of 0.9 and RMSE of 18.3 mm with growing season ET c underestimated in the range between 3.2 to 6.02%. A spatio-temporal distribution maps of instantaneous energy balance fluxes and seasonal ET c showed an anticipated variation across the study area as the plant growth stage and surface condition changed. In addition, mean monthly ET c values estimated across the study area for different land cover types showed a substantial variation throughout the growing season as climate and crop growth stages differed. Furthermore, it was found that, in the arid to semi-arid intermountain region of Wyoming, the contribution of irrigation varied in the range of 73-81% of the total seasonal ET c in nine irrigation districts that fall within the study area. The high relative contribution of irrigation highlights the importance of identifying and quantifying ET c for improved management in irrigation system design and water allocation.
In this study, METRIC algorithms performed reasonably well, particularly in estimating cumulative ET c for the entire growing season. However, the model performance was poor during the early and late growing season when fields had little active leaf area, with METRIC estimated ET c values were within 10% of the BREBS measured values. A companion paper to the research compared the performance of four surface energy balance models viz. METRIC, SEBAL, SEBS, and S-SEBI in estimating the ETc in the intermountain region of Wyoming. The METRIC model can be an important tool in quantifying and mapping surface energy balance components especially in places where quality hourly weather data are available. Accurate hot and cold pixel selection during sensible heat calculation is a defining moment of the METRIC algorithm. A selection error can result in an undesired outcome. The selection of hot pixels was an arduous task, especially when P events occurred prior to Landsat overpass time, raising the LE value. Similarly, the estimation of negative reference ET fraction (ET r F) in some naturally vegetated and deserted areas was observed altering ET c estimation for those surfaces. Thus, further research work on securing accurate hot pixels during prior P events, controlling wayward estimations during the early and late growing season, and minimizing the estimation of negative reference ET fraction in naturally vegetated and deserted areas can provide important support for a more accurate estimation of instantaneous and periodic fluxes. Furthermore, future research will also focus on multi-sensor data fusion to increase the temporal variation and thereby reduce the uncertainties in ET c estimations.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/12/24/4019/s1, Table S1: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Cody Canal Irrigation District for 2017 and 2018. Table S2: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Deaver Irrigation District for 2017 and 2018. Table S3: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Greybull Valley Irrigation District for 2017 and 2018. Table S4: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for our major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Heart Mountain Irrigation District for 2017 and 2018. Table S5: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Hunt & Godfrey Irrigation District for 2017 and 2018. Table S6: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Lovell Irrigation District for 2017 and 2018. Table S7: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Shoshone Irrigation District for 2017 and 2018. Table S8: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Sidon Irrigation District for 2017 and 2018. Table S9: Average seasonal crop evapotranspiration ET c , seasonal precipitation (P), and Irrigation contribution to ET c (ET c -P) for four major crops i.e., Sugar beet, Dry bean, Alfalfa, and Maize within Will Wood Irrigation District for 2017 and 2018.

Conflicts of Interest:
The authors declare no conflict of interest.