Comparison of Satellite Driven Surface Energy Balance Models in Estimating Crop Evapotranspiration in Semi-Arid to Arid Inter-Mountain Region

The regional-scale estimation of crop evapotranspiration (ETc) over a heterogeneous surface is an important tool for the decision-makers in managing and allocating water resources. This is especially critical in the arid to semi-arid regions that require supplemental water due to insufficient precipitation, soil moisture, or groundwater. Over the years, various remote sensing-based surface energy balance (SEB) models have been developed to accurately estimate ETc over a regional scale. However, it is important to carry out the SEB model assessment for a particular geographical setting to ensure the suitability of a model. Thus, in this study, four commonly used and contrasting remote sensing models viz. METRIC (mapping evapotranspiration at high resolution with internalized calibration), SEBAL (surface energy balance algorithm for land), S-SEBI (simplified surface energy balance index), and SEBS (surface energy balance system) were compared and used to quantify and map the spatio-temporal variation of ETc in the semi-arid to arid inter-mountain region of Big Horn Basin, Wyoming (Landsat Path/Row: 37/29). Model estimates from 19 cloud-free Landsat 7 and 8 images were compared with the Bowen ratio energy balance system (BREBS) flux stationed in a center pivot irrigated field during 2017 (sugar beet), 2018 (dry bean), and 2019 (barley) growing seasons. The results indicated that all SEB models are effective in capturing the variation of ETc with R2 ranging in between 0.06 to 0.95 and RMSD between 0.07 to 0.15 mm h−1. Pooled data over three vegetative surfaces for three years under irrigated conditions revealed that METRIC (NSE = 0.9) performed better across all land cover types, followed by SEBS (NSE = 0.76), S-SEBI (NSE = 0.73), and SEBAL (NSE = 0.65). In general, all SEB models substantially overestimated ETc and underestimated sensible heat (H) fluxes under dry conditions when only crop residue was available at the surface. A mid-season density plot and absolute difference maps at image scale between the models showed that models involving METRIC, SEBAL, and S-SEBI are close in their estimates of daily crop evapotranspiration (ET24) with pixel-wise RMSD ranged from 0.54 to 0.76 mm d−1 and an average absolute difference across the study area ranged from 0.47 to 0.56 mm d−1. Likewise, all the SEB models underestimated the seasonal ETc, except SEBS.


Introduction
Recent times have seen the unsustainable utilization of water resources, leading to short-and long-term water crises. The degrading soil and water resources coupled with climate change and variability have made it inevitable to scientifically manage agricultural water [1]. In such a scenario, managing the scarce water resources to fulfill increasing demands is a challenge. The accurate quantification of ET c at local and regional scales can aid in water resource-based policy and decision making and help manage our water resources. ET c is an energy-driven process and an important component of water budget [2] and is an essential component of irrigation water requirement quantification, irrigation and elevation [10,13,36,37]. Thus, this research anticipates identifying the best fit model suitable to semi-arid to arid mountainous climatic settings. Likewise, complex models require extensive time and effort to set up and do not necessarily outperform the simpler models [38]. Similarly, the evaluation of different models helps to identify the limitations and uncertainties of a model. Therefore, the specific objectives of this study were to: (i) assess and compare the performance of METRIC, SEBAL, SEBS, and S-SEBI algorithms using Landsat imagery on estimating ET c with respect to measured ET c from the Bowen ratio energy balance system for the different vegetative surfaces in the intermountain region of Wyoming; (ii) quantify, map, and evaluate spatial and temporal distribution (daily, monthly, and seasonal) of ET c over the study region via the METRIC, SEBAL, SEBS, and S-SEBI model.

Study Area, Climate and Satellite Dataset, and Image Processing
The study area is situated in the Rocky Mountain region of the United States covering the majority of the Bighorn Basin, Wyoming (Landsat Path: 37 and Row: 29; Figure 1) [25]. The study area encompasses 7930 Km 2 with an elevation between 1110 m to 3254 m. A total of 19 Landsat 7-Enhanced Thematic Mapper Plus (ETM+) and Landsat 8-Operational Land Imager (OLI) and thermal infrared sensor (TIRS) images (Path: 37, Row: 29) were retrieved for the 2017 (5 images), 2018 (9 images), and 2019 (5 images) growing seasons ( Table 1). The majority of the study area is covered by natural vegetation (83%) comprising evergreen and deciduous forest, grassland, shrubland, and woody and herbaceous wetlands. The growing season is generally short due to a limited number of frost-free days in the year. The major soil type in the study area is sandy loam while major crops are alfalfa, barley, dry bean, sugar beet, and corn.
In general, the study area is dominated by a semi-arid to arid climate, with average annual precipitation (1981-2010) as low as 235 mm [39]. Over the study duration, substantial variations were observed in weather variables during each of the three years. Compared to long-term average values, 2017 had below-normal precipitation and higher ET r . The total seasonal P of 94 mm, 128 mm, and 218 mm was observed in 2017, 2018, and 2019, respectively. Likewise, the cumulative ET r was highest in 2017 (883 mm) followed by 2018 (820 mm), and the 2019 (794 mm) growing season. A detailed description of the study area characteristics and weather conditions were reported in Acharya et al. [25], Sharma et al. [40], and Rai et al. [41]. High quality hourly and daily weather data to estimate ET r and precipitation data were retrieved 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. The quality control of meteorological data was performed by the Water Resources Data System at the University of Wyoming and the data were disseminated via the WACNet website [42]. Quality control involves both automated and manual techniques. The flagged and missing datasets (failed quality control) are estimated via spatial and statistical methods to result in a complete data-set [43]. Based on the quality analyses, all climate data used over the study duration were judged to be of good quality. The meteorological condition at PREC, Wyoming on image acquisition day is provided in Table 1.
Likewise, the BREBS installed in 2017 in a center-pivot-irrigated field in PREC, Wyoming (44.46 • N, 108.45 • W) was used to assess the performance of the SEB model over three irrigated vegetative surfaces, i.e., sugar beet in 2017, dry bean in 2018, and barley in 2019. These crops were selected based on their high acreage, high economic value, and high irrigation demand in the intermountain region of Wyoming. In 2017, sugar beet cultivar 9418RR was planted on 8 May at 0.3 m plating depth with a planting density of 118,600 seeds per hectare and harvested on 13 October. In 2018, dry bean cultivar La Paz was planted on 5 June at 0.05 m soil depth and 0.56 row spacing at a seeding rate of 222,395 seeds per hectare and harvested on 10 September. In 2019, barley was planted on 5 April at a target planting density of 850,000 plants per hectare on 0.19-row spacing Remote Sens. 2021, 13, 1822 5 of 31 and harvested on 14 August. A thorough description of the BREBS station installed at PREC, Wyoming is provided on our precedent paper [25]. The BREBS data were closely supervised and general maintenance was provided on a regular basis (once a week). The meteorological data (e.g., air temperature, relative humidity, incoming shortwave solar radiation, wind speed) collected from BREBS underwent rigorous quality checks with comparisons with nearby weather station data. In BREBS, the accuracy of the calculated LE and H fluxes depends on the accuracy of the Bowen ratio. A negative Bowen ratio during the day and a larger Bowen ratio during sunset and sunrise indicate an error in H and/or LE fluxes. In such instances, H and LE fluxes are recalculated via bulk aerodynamic estimation techniques using wind speed and temperature gradient [44].
ote Sens. 2021, 12, x FOR PEER REVIEW 5 of at a target planting density of 850,000 plants per hectare on 0.19-row spacing and h vested on 14 August. A thorough description of the BREBS station installed at PREC, W oming is provided on our precedent paper [25]. The BREBS data were closely supervis and general maintenance was provided on a regular basis (once a week). The meteorolo ical data (e.g., air temperature, relative humidity, incoming shortwave solar radiati wind speed) collected from BREBS underwent rigorous quality checks with compariso with nearby weather station data. In BREBS, the accuracy of the calculated LE and fluxes depends on the accuracy of the Bowen ratio. A negative Bowen ratio during day and a larger Bowen ratio during sunset and sunrise indicate an error in H and/or fluxes. In such instances, H and LE fluxes are recalculated via bulk aerodynamic estim tion techniques using wind speed and temperature gradient [44].   The SEB models convert the digital numbers (DNs) of each image pixel to comprehensible SEB fluxes. The conversion begins with the computation of top of atmosphere radiance and reflectance from the geo-rectified images using equations as provided in USGS Landsat 7 and 8 handbooks. Moreover, the data loss incurred in Landsat 7 due to scan line correction error was overcome by carrying out natural neighbor interpolation in ArcGIS 10.6. The final output images (instantaneous) were gap filled by finding the closest subsets of missing pixels from the surrounding pixels. Image processing was carried out using the ERDAS IMAGINE 2020 (Leica Geosystems GIS and Mapping, LLC, Atlanta, Georgia, USA) graphical model maker tool. Likewise, ArcMap 10.8 and R-Studio Version 1.4.1103 were employed for the presentation of analyzed images. To account for aspects, slopes, and elevation on T s , the solar incidence angle was computed for each pixel and a lapse rate of 6.5 K km −1 was considered in all the models. METRIC additionally used a second lapse rate of 10 K km −1 in the mountainous terrain to compute elevation-corrected surface temperature (T s_dem ) [13]. A lapse rate change of 2000 m (foot of the mountain) was considered based upon the average elevation of the study area to toggle between the two lapse rates in METRIC. In addition to the above modification, Z om , wind speed, and atmospheric pressure were adjusted [13] as they appeared in the model algorithm. The digital elevation model (DEM) was used to account for slopes, aspects, and elevation in the above adjustments. In this study, a land-use map was used to estimate Z om and extract information for different land cover types during the result presentation. Table 2 provides the source of various datasets and the computation of primary and secondary inputs (intermediate outputs) needed in the SEB models.

Energy Balance Models (SEB)
The SEB models compute ET c as a residual of other energy balance fluxes viz. R n , soil heat (G) and H [13] (Equation (1)) as: The following section provides a brief overview of the SEB models used in this study.

METRIC Model and SEBAL Model
The net radiation (R n ) is derived for each pixel by deducting reflected short-(αR s ↓) and longwave ((1 − ε o ) R L ↓) as well as emitted long-wave (R L ↑) from incident short-(R s ↓) and longwave (R L ↓) [13].
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 ), and ε o is the surface thermal emissivity (dimensionless). For METRIC, G was empirically calculated as a G/R n fraction using vegetation indices and surface temperature [52].
where T s is the surface temperature in Kelvin and LAI (dimensionless) is the leaf area index. For SEBAL, an empirical relation put forward by Bastiaanssen et al. [14] was used to compute G as: The main feature of the METRIC and SEBAL model is the assumption of linearity between T s and the near-surface air temperature difference [10,[13][14][15]. In fact, the METRIC model was developed on the back of the SEBAL model [10,13], emphasized to accurately predict energy balance fluxes, particularly in advective conditions and undulating ter-Remote Sens. 2021, 13, 1822 8 of 31 rain [10,13]. The major difference between these two models lies in how anchor pixels (hot and cold pixel) are selected during H calculation. METRIC assumes a wet, well-irrigated crop surface with full cover as a cold pixel candidate, whereas a dry, bare agricultural field is its hot pixel candidate [10,13]. On the other hand, local water bodies and dry, bare agricultural fields are cold and hot pixel candidates in SEBAL, respectively [10,[13][14][15]. NDVI, α, LAI, T s , and land use map was taken into consideration while selecting the anchor pixels for both the model. In this study, a manual selection of hot and cold pixels was performed for both METRIC and SEBAL models. Both the anchor pixels were selected within a 20 km distance from the meteorological station. For METRIC, cold pixel candidates were selected from the well-irrigated densely vegetated area with the NDVI between 0.76 and 0.84, mid-season LAI higher than 3 or 4 m −2 m −2 , surface albedo between 0.18 to 0.24, and comparatively lower T s . However, for SEBAL, cold pixel candidates were picked from the nearest water body to the meteorological station. The hot pixel candidate for both the models had an NDVI less than 0.2, surface albedo between 0.17 to 0.23, and comparatively higher T s . Thus, hot pixels were the same for METRIC and SEBAL models in all the images considered in the study. A detailed description of how the selection of hot and cold pixels is performed is presented in our companion paper [25]. For both METRIC and SEBAL, H is calculated as a function of aerodynamic observations such as wind speed at 2-m height (u2), vegetation type and roughness, and surface to air temperature differences (T s − T a ) as [13,14]: 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. Our companion paper, Acharya et al. [25], provides a detailed description of the METRIC algorithm as well as for H computation of the SEBAL model. Unlike SEBAL, METRIC uses ground-measured instantaneous ET r to internally calibrate the cold pixel during H calculation [10,13]. SEB models compute ET inst as residual of the surface energy balance components (Equation (7)).
where ρ w is the density of water (~1000 kg m −3 ) and λ refers to the latent heat of vaporization (J Kg −1 ). METRIC uses the reference ET fraction (ET r F) value to upscale ET inst to ET 24 and periodic ET c [13]. ET r F (dimensionless) is computed as a ratio of ET inst to the ET r .
where ET r is computed using the standardized ASCE Penman-Monteith equation [53] on an hourly basis. METRIC utilizes ET r F to minimize the impact of advective heat on ET c . Likewise, ET r F is presumed to be equivalent to the crop coefficient and its values remain constant throughout the day [54]. The ET 24 rate (mm d −1 ) is then calculated by multiplying ET r F values with hourly ET r values summed over 24 h (ET r-24 ) (Equation (9)).
ET 24 = ET r F × ET r−24 (9) For SEBAL ET 24 is computed using Λ and daily average net radiation. Λ is the ratio of latent heat (R n − G − H) to the energy available at the surface (H + LE) and denotes the fraction of energy that gets partitioned towards LE. SEB fluxes (R n , G, H, and ET c ) have considerable diurnal variation. However, the ratio of these fluxes (Λ) is reported to be relatively constant during the day [55]. Thus, Λ is being used in various SEB models to upscale ET inst to ET 24 .
where R n24 is the daily average net radiation (W m −2 ) and G 24 (assumed to be zero) is the daily average G flux (W m −2 ). R n24 was computed using daily extraterrestrial solar radiation (R a ; W m −2 ) and daily atmospheric transmittance (τ sw ) [3]. This R n24 is the same as that used in ET r computation using the ASCE Penman-Monteith equation. Likewise, cubic spline interpolation was favored over linear interpolation [25] to interpolate ET 24 to monthly and seasonal ET c . A detailed description of cubic spline interpolation [13] is also provided in our companion paper [25].
where ET cumulative (mm) represents the summation of ET 24 from the first day to the last day of the month, K c−i is the interpolated K c for a month and ET r 24 i is the daily ET r summed over a month. Figure 2 shows the flowchart of SEB components and other intermediate parameters estimated by SEBS. For the SEBS model, R n is computed as in Equation (2). In this study, an empirical relationship between vegetation indices and T s as provided by Tasumi [52] was used to calculate G (Equations (3) and (4)). Unlike METRIC and SEBAL, where anchor pixels delineate the energy balance boundary conditions, the selection of hot and cold pixels is not required in the SEBS model. In SEBS, H is calculated using an iterative procedure that solves the relationships for the layers of the friction velocity, the difference (∆θ, K) between the near-surface potential air temperature (T a ) and T s , and Monin-Obukhov length (L) which is expressed as:

SEBS Model
where u* = (τ o / ρ) 1/2 is the friction velocity, τ o is the surface shear stress, ρ is the density of air, k is von Karman's constant, d o is the zero-plane displacement height, z is the height above the surface, z om is the roughness height for momentum transfer, z oh is the scalar roughness height for heat transfer, Ψ m and Ψ h are stability correction functions for momentum and H transfer, C p is the specific heat of air at constant pressure, g is the acceleration due to gravity, θ v is the virtual temperature near the surface. The Paulson [56] and Webb [57] method was used to determine the roughness height for heat and momentum transfer and the stability correction function. An empirical relationship provided by Brutsaert [58] was used to compute wind speed at blending height (z) and the zero-plane displacement height (d o ). Figure 2. Flowchart illustrating the computation of surface energy balance fluxes using SEBS algorithm. Figure 3 is the computational flowchart of SEB components by S-SEBI. The peculiar feature of the SSEBI model is that it does not require any meteorological data [19] and the model algorithm is also comparatively simpler. The S-SEBI model is based on the correlation observed between Ts and α ( Figure 4). The concave Ts vs. α relationship is divisible into the evaporation-controlled and radiation-controlled zone. At the evaporation-controlled zone, Ts is first more or less constant and then starts to increase with increasing α value. Beyond a certain threshold value of α, Ts decreases with increasing α which is the radiative zone. This correlation is used to calculate Λ (Equations (20)- (22)). Further, the SEBS algorithm utilizes relative evaporation (Λ r ) derived from H at dry (H dry ) and wet limits (H wet ) to compute Λ (Equations (17) and (18)), ET inst (Equation (7)), and ET 24 (Equation (11)) [16]. At H wet , evaporation occurs at a potential rate and H is at its minimum value. It is calculated using the Penman-Monteith parameterization equation [59], assuming the bulk internal resistance to be zero. Likewise, H is at its maximum value (LE = 0) at H dry . It corresponds to fields where soil moisture is limited for evaporation to occur. Periodic ET c in SEBS is computed as in Equation (13). Figure 3 is the computational flowchart of SEB components by S-SEBI. The peculiar feature of the SSEBI model is that it does not require any meteorological data [19] and the model algorithm is also comparatively simpler. The S-SEBI model is based on the correlation observed between T s and α ( Figure 4). The concave T s vs. α relationship is divisible into the evaporation-controlled and radiation-controlled zone. At the evaporationcontrolled zone, T s is first more or less constant and then starts to increase with increasing α value. Beyond a certain threshold value of α, T s decreases with increasing α which is the radiative zone. This correlation is used to calculate Λ (Equations (20)- (22)).

S-SEBI Model
where, T H is the temperature of a dry pixel, and T λE is the temperature of a wet pixel, r o is threshold α, a h and b h are the regression coefficients for the dry boundary, and a λE and b λE are the regression coefficients for the wet boundary. The T s vs. α regression coefficients were computed by excluding α's below the threshold α value at the wet boundary [19]. Threshold α (r o ) distinguishes the evaporative and radiative regions of the scatterplot and corresponds to the maximum temperature of T s vs. α relationship. Λ is then used to compute H (Equation (23)) and ET inst (Equation (7)). ET 24 and periodic ET c is computed using Equations (11) and (12), respectively.
where, TH is the temperature of a dry pixel, and TλE is the temperature of a wet pixel, ro is threshold α, ah and bh are the regression coefficients for the dry boundary, and aλE and bλE are the regression coefficients for the wet boundary. The Ts vs. α regression coefficients were computed by excluding α's below the threshold α value at the wet boundary [19]. Threshold α (ro) distinguishes the evaporative and radiative regions of the scatterplot and corresponds to the maximum temperature of Ts vs. α relationship. Λ is then used to compute H (Equation (23)) and ETinst (Equation (7)). ET24 and periodic ETc is computed using Equations (11) and (12), respectively.

Models Validation
Model accuracy was assessed using the standard regression statistics, i.e., mean, standard deviation, and coefficient of determination (R 2 ) along with three error-index sta-y = −20.92x + 324.09 y = 19.06x + 294.74

Models Validation
Model accuracy was assessed using the standard regression statistics, i.e., mean, standard deviation, and coefficient of determination (R 2 ) along with three error-index statistics, i.e., root mean square difference (RMSD), percent bias error (PBE), and Nash-Sutcliffe's efficiency (NSE) for estimating ET c as compared to corresponding BREBS flux as: Besides, intercomparison between the models was carried out using density plots and an absolute difference map.

Comparison between SEB-Estimated and BREBS-Measured Instantaneous Fluxes
To evaluate the performance of SEB models, model-estimated instantaneous fluxes are compared with corresponding BREBS fluxes measured at the BREBS flux tower footprint in the Powell Research and Extension Center, Powell, Wyoming (Figures 5 and 6) using the approach described by Acharya et al. [25]. BREBS values were retracted at satellite overpass time (11:00 a.m. MST). A total of 19 different data points acquired for the three vegetative surfaces (sugar beet in 2017, dry bean 2018, and barley in 2019) along with the pooled data were used in the analyses. Table 3 provides a statistical comparison between the model-estimated and BREBS-measured ET inst . Our results indicated that all SEB models are effective in capturing the variation of ET inst with R 2 ranging in between 0.06 to 0.95 and RMSD from 0.07 to 0.15 mm h −1 . Pooled data over three vegetative surfaces for three years under irrigated conditions revealed that METRIC (NSE = 0.9) performed better across all land cover types, followed by SEBS (NSE = 0.76), S-SEBI (NSE = 0.73), and SEBAL (NSE = 0.65). Although no significant difference was observed between R 2 values of METRIC and SEBS-estimated vs. BREBS-measured ET inst (0.91 vs. 0.87), the RMSD, NSE, and PBE of SEBS model was 27% larger, 18% lower, and 54% higher compared to the METRIC model.
Since the performance of SEB models varies with the underlying surface, this study also assessed the performance of SEB models individually for each surface. All SEB models performed poorly in 2017 over sugar beet, except the SEBS model. The scatter plot between the model-estimated vs. BREBS-measured ET inst revealed an R 2 of 0. 21  To further understand the variation in SEB model-derived ETinst, estimated H was compared with the BREBS-measured H ( Figure 6). Overall, all models are effective in capturing the variation of H with R 2 ranged from 0.69 to 0.86 and RMSD from 59 to 88 W m −2 . The average SEB models estimate of H was 13%, 7%, and 5% higher for METRIC, SEBAL, and S-SEBI and 14% lower for the SEBS model as compared to ground observations (mean: 108 W m −2 ). However, the biggest difference in modeled and measured H was observed for SEBS (RMSD = 88 W m −2 , PBE = −14%), where most of the underestimation was observed in the dry bean (35%) in 2018 and barley (10%) in the 2019 growing season after harvest. The fact that the 2018 growing season had nine images spread across the growing season can lead to disproportionally more accuracy for the dry beans as compared to sugar beet (2017 with five images) and barley (2019 with five images) (Table 1). An analysis of temporal biases towards the end of the 2018 growing season was carried out to better understand the model performance. For that, 2018 images from mid-July to September    To further understand the variation in SEB model-derived ET inst , estimated H was compared with the BREBS-measured H ( Figure 6). Overall, all models are effective in capturing the variation of H with R 2 ranged from 0.69 to 0.86 and RMSD from 59 to 88 W m −2 . The average SEB models estimate of H was 13%, 7%, and 5% higher for METRIC, SEBAL, and S-SEBI and 14% lower for the SEBS model as compared to ground observations (mean: 108 W m −2 ). However, the biggest difference in modeled and measured H was observed for SEBS (RMSD = 88 W m −2 , PBE = −14%), where most of the underestimation was observed in the dry bean (35%) in 2018 and barley (10%) in the 2019 growing season after harvest.

Comparison of Estimated and Measured Monthly Crop Evapotranspiration (ET c )
In general, monthly ET c is more desirable for hydrological applications (e.g., seasonal irrigation requirement, conveyance capacities of irrigation systems, etc.) as compared to ET inst and ET 24 . In this study, the monthly ET c was calculated from ET 24 using cubic interpolation [25]. Figure 7 depicts the performance of the SEB models on estimating the monthly ET c . For that, the model predicted monthly ET c for the 2018 growing season (May to September) were plotted against BREBS-measured corresponding flux. For all the models, moderate to high correlation was observed between measured and estimated ET c with R 2 ranged from 0.70 to 0.95 and RMSD ranged from 18 mm for METRIC to 34 mm for S-SEBI. In general, all the models underestimated the seasonal ET c , except SEBS. However, METRIC stood out to be best-performing model owing to its lower RMSD (17.6 mm) and lower percentage error (−7.9%) followed by SEBS (RMSD = 25.8 mm, % error = +6.6), SEBAL (RMSD = 33 mm, % error = −29) and S-SEBI (RMSD = 34.3 mm, % error = −31). METRIC monthly ET c ranged between -38% (underestimation) in June to 18% (overestimation) in October. The over and underestimation of the estimated monthly ET c was higher when the soil surface was devoid of active leaf area cover. The dry bean at the BREBS footprint was planted on 5 June and harvested on 10 September. Similarly, overestimation (46%) was observed for SEBS-derived ET c in October after the dry bean harvest. Contrary to this, consistence underestimation ranging between 18% to 55% was observed for SEBAL and S-SEBI.
The over and underestimation of the estimated monthly ETc was higher when th surface was devoid of active leaf area cover. The dry bean at the BREBS footprin planted on 5 June and harvested on 10 September. Similarly, overestimation (46% observed for SEBS-derived ETc in October after the dry bean harvest. Contrary t consistence underestimation ranging between 18% to 55% was observed for SEBA S-SEBI.

Intercomparison of SEB Model Estimated Daily Evapotranspiration (ET 24 )
An intercomparison of the SEB model-derived ET 24 was carried out using density plots (Figures 8 and 9) and absolute difference maps ( Figure 10). Figure 8 reflects a concentrated range of estimated ET 24 on 2018 image acquisition dates. The shape of the curve fluctuated (unimodal, bimodal, and multimodal) throughout the growing season as surface conditions changed. The study area comprises vast swathes of natural vegetation (83% of the total study area) whereas the cultivated area (14%) is limited. Bimodal peaks (one of the peaks relate to cropland and the other for natural vegetation) is observed in SEB models (especially for METRIC) in mid and late season (18 July, 11 August, 4 September, and 12 September) when crops are at full growth and the naturally vegetated area is under water stress. As expected, the distribution of density plots for all SEB models in mid-season (2 July, 18 July, 26 July, and 11 August) is more skewed toward higher ET c . Likewise, to understand the behavior of SEB-modeled fluxes at a spatial scale larger than a single tower footprint, the gridded SEB model ET 24 output over the study domain was compared with each other using the density plots ( Figure 9). The R 2 and RMSD values were calculated from approximately 400,000 pixels randomly selected using ordinary linear regression. The mid-season density plots involving METRIC, SEBAL, and S-SEBI for 2 July 2018 showed a good degree of correlation (R 2 > 0.81) (Figure 9a,b,d). METRIC and SEBAL models were comparatively much closer on their estimates of ET 24 with R 2 of 0.93 and RMSD of 0.54 mm d −1 (Figure 10a). The highest difference was between the S-SEBI vs. SEBS and SEBAL vs. SEBS models with RMSD of 1. 21    Similarly, to further investigate the difference in ET24 estimated using different SEB models over the heterogeneous surface, the image-scale pixel-by-pixel absolute difference map on 2 July 2018 was constructed ( Figure 10). The absolute difference map provides the inter-comparison of SEB models based on the distance from zero on the number line. On an image scale, the average difference ranged between 0.47 mm d −1 for METRIC-SEBAL to 1.23 mm d −1 for S-SEBI-SEBS. It was found that the METRIC-SEBAL model had a comparatively lower absolute difference followed by METRIC-S-SEBI (0.53 mm d −1 ). On the other hand, maps involving SEBS had greater absolute differences (1 to 1.23 mm d −1 ) which is consistent with the greater spread-out distribution of ET24 (Figure 10c, 10e, and  10f). Similarly, Figure 11 also reflects the absolute difference being lower in the cropped region as compared to the naturally vegetated area. The average absolute difference between the models in the cropped area ranged between 0.33 to 0.78 mm d −1 , while for natural vegetation it ranged between 0.45 to 1.3 mm d −1 . Similarly, to further investigate the difference in ET 24 estimated using different SEB models over the heterogeneous surface, the image-scale pixel-by-pixel absolute difference map on 2 July 2018 was constructed ( Figure 10). The absolute difference map provides the inter-comparison of SEB models based on the distance from zero on the number line. On an image scale, the average difference ranged between 0.47 mm d −1 for METRIC-SEBAL to 1.23 mm d −1 for S-SEBI-SEBS. It was found that the METRIC-SEBAL model had a comparatively lower absolute difference followed by METRIC-S-SEBI (0.53 mm d −1 ). On the other hand, maps involving SEBS had greater absolute differences (1 to 1.23 mm d −1 ) which is consistent with the greater spread-out distribution of ET 24 (Figure 10c, 10e, and  10f). Similarly, Figure 11 also reflects the absolute difference being lower in the cropped region as compared to the naturally vegetated area. The average absolute difference between the models in the cropped area ranged between 0.33 to 0.78 mm d −1 , while for natural vegetation it ranged between 0.45 to 1.3 mm d −1 .

Mapping Spatial Variation of SEB Model-Estimated Seasonal Evapotranspiration (ETc)
To investigate the difference in SEB models for different land use, spatial-interco parison of the SEB models were also carried out using spatial variation of daily and s sonal ETc during the 2018 growing season. For daily ETc, one image on 11 August 2 from the middle of the growing season was selected. Spatio-temporal variability throu out the growing season were also analyzed and added as Supplementary Materials (F ures S1-S4). For seasonal maps, images from May to September were considered ( Figu  11 and 12). On average, METRIC (4.7 mm d −1 ) estimated a higher ET24 for cropland lowed by SEBS (4.1 mm d −1 ), S-SEBI (3.8 mm d −1 ), and SEBAL (3.6 mm d −1 ) (Figure 8). natural vegetation (forest, shrubland, grassland, and wetlands), the average ET24 for

Mapping Spatial Variation of SEB Model-Estimated Seasonal Evapotranspiration (ET c )
To investigate the difference in SEB models for different land use, spatial-intercomparison of the SEB models were also carried out using spatial variation of daily and seasonal ET c during the 2018 growing season. For daily ET c , one image on 11 August 2018 from the middle of the growing season was selected. Spatio-temporal variability throughout the growing season were also analyzed and added as Supplementary  Materials (Figures S1-S4). For seasonal maps, images from May to September were considered (Figures 11 and 12). On average, METRIC (4.7 mm d −1 ) estimated a higher ET 24 for cropland followed by SEBS (4.1 mm d −1 ), S-SEBI (3.8 mm d −1 ), and SEBAL (3.6 mm d −1 ) (Figure 8). For natural vegetation (forest, shrubland, grassland, and wetlands), the average ET 24 for all the models was less than 1 mm d −1 ,except for S-SEBI, which had an average ET 24 of 1.6 mm d −1 . A comparison between and within SEB model(s) throughout the growing season (spatio-temporal variability) was also performed using the ET 24 maps (refer to Supplementary Materials, Figures S1-S4). At the seasonal scale, as expected, cropland observed higher ET c as compared to the naturally vegetated area. The study area received an inconsistent rainfall distribution during the 2018 growing season, where the majority of the rainfall was observed during May (49 mm) and June (49 mm) with a total rainfall of 158 mm between March to September. Likewise, coniferous forests (top right and bottom left of the image) accounted for a higher seasonal ET c as compared to other land cover types. Intermittent cloud masses (brown patches in the image) were captured above naturally vegetated areas, resulting in a higher ET c estimation as compared to surrounding vegetation. However, to extract the ET c rate (Table 4) cloud masking was performed. Comparison among SEB models revealed that SEBS had a higher estimation of seasonal ET c for all the major crops. SEBS seasonal ET c estimates were 4, 21, and 22% higher compared to METRIC, S-SEBI, and SEBAL, respectively. However, for natural vegetation, S-SEBI (291 mm) predicted higher average seasonal ET c , followed by METRIC (250 mm), SEBAL (186 mm), and SEBS (284 mm). ever, to extract the ETc rate (Table 4) cloud masking was performed. Comparison among SEB models revealed that SEBS had a higher estimation of seasonal ETc for all the major crops. SEBS seasonal ETc estimates were 4, 21, and 22% higher compared to METRIC, S-SEBI, and SEBAL, respectively. However, for natural vegetation, S-SEBI (291 mm) predicted higher average seasonal ETc, followed by METRIC (250 mm), SEBAL (186 mm), and SEBS (284 mm).

SEB Models and Their Strengths and Limitations
This study analyzed 19 Landsat images captured during the 2017, 2018, and 2019 growing season to assess and compare the performance of the METRIC, SEBAL, SEBS, and S-SEBI algorithms in the semi-arid to arid intermountain region of Wyoming using four SEB models. The study anticipates putting forward a best suited SEB model for the region. Without a proper model comparison study, it would be unwise to choose a model and apply it for real-world application by policy and decision-makers. It is a matter of fact that SEB models can respond differently in different climatic and geographic settings. A model performing better in a humid condition might not be able to produce the same result in rather arid conditions [26,30].
Our result indicates that the SEB model-estimated fluxes (instantaneous and periodic) have a substantial differences between each other. The fact that the SEB models differ in the model structure, input data requirements, and assumptions made in the algorithm has contributed to the variability in model output. For example, METRIC and SEBAL models use CIMEC (calibration using inverse modeling of extreme conditions) calibration approach during H calculation to remove the systematic biases in the estimation of surface temperature and surface reflectance [13]. In METRIC and SEBAL, two extreme conditions (dry and wet pixels) are chosen to internally calibrate the SEB using hourly ET r . Therefore, the final ET c estimates in these models can be accurate even if other SEB components face uncertainties. However, models such as SEBS and S-SEBI do not undergo such rigorous calibration, and delineation using anchor pixels and ground-based hourly ET r . The SEBS algorithm theoretically determines the wet and dry limit for H and Λ. Likewise, S-SEBI utilizes a T s vs. α relationship to determine the wet edge and dry edge.
In this study, the METRIC estimate of ET c was found to have a better correlation with the corresponding BREBS flux. This performance can be directly linked with the internal calibration procedure (CIMEC) performed in METRIC. SEBAL ET c estimates were lower compared to METRIC. This may be due to the difference: (i) in the selection of cold pixel between METRIC and SEBAL; (ii) use of ET r F compared to Λ to upscale ET inst to ET 24 [26]. Likewise, all the SEB models overestimated ET inst , with percent biases ranging between 2.2% for SEBAL to 12.3% for SEBS. Similar overestimation in SEB-estimated ET c ranging between 8 to 32% was observed by Wagle et al. [30], who compared the performance of five SEB models viz METRIC, SEBAL, SEBS, S-SEBI, and SSEBop. They further ranked the SEBS models based on four statistical measures and concluded that S-SEBI and SEBAL performed better, followed by SEBS and METRIC. They reported that the poor performance of METRIC is due to overestimation of METRIC-estimated ET c because of drier and rainfed conditions, which is contrary to the irrigated setting used in this study. Similar to this study, the higher performance of METRIC and SEBS were observed by Singh and Senay [26] who compared METRIC, SEBAL, SEBS, and SSEBop-estimated ET c using Landsat 5 and 7 images over center-pivot irrigated continue corn and center-pivot irrigated maize-soybean rotation in mid-western U.S. The performance of METRIC was also found to be better compared to SEBS by Liaqat and Choi [33] in Northeast Asia over four different vegetative surfaces. They reported that internal calibration in METRIC helped reduce the biases in atmospheric correction and other input parameters. On the other hand, the inconsistent estimation of G and H in SEBS resulted in high ET c .
A consistent overestimation with PBE ranged from 11.4% for SEBAL to 38.3% for SEBS in 2018 is due to the fact that three images out of nine were captured when the soil surface was devoid of active leaf area (only crop residue was available at the surface). The three images resulted in a combined overestimation of 41%, 15%, 71%, and 10% for METRIC, SEBAL, SEBS, and S-SEBI compared to less than 10% overestimation when all the images within the growing season were analyzed. It has been reported [17,60] that the presence or absence of crop residue at the surface can impact the T s , α, and emissivity, which can ultimately impact ET c . Allen et al. [10] compared METRIC-estimated ET c with lysimetermeasured ET c using eight Landsat images acquired from April to September. They reported 14% averaged absolute differences for sugar beet crop when one image date with dry bare soil was omitted. However, the difference increased to 30% when all the image dates were considered. Sharma et al. [17] reported higher overestimation in SEBS-estimated ET 24 (21% in 2009-2010 in winter wheat crop season and 21% in 2011 in maize crop growing season) after harvest when only crop residue was available at the surface. Other studies [30,61] have also reported that SEB models overestimate ET c over crop residue and vegetation under dry conditions. Wagle et al. [30] found that all SEB models, i.e., METRIC, SEBAL, SEBS, S-SEBI, and SSEBop overestimate ET c when soil moisture is less than 10 percentiles and their performance improved with increasing soil moisture. It is worth pointing out that the SEB model outputs (ET inst ) in this study exhibited considerable temporal biases, particularly when the model estimates from different times of the growing season were considered. For example, the BREBS tower considered in this study was installed in June 2017, which resulted in only a handful of images (five) in the mid to late growing season for analysis in 2017 (sugar beet growing season), compared to nine images in 2018 (dry bean growing season). This resulted in the skewness of data and errors in one direction, causing lower NSE and R 2 values in 2017 [25].
As compared to the BREBS-measured H, a slight overestimation from METRIC and SEBAL ( Figure 6) can be accounted to the internal calibration process (CIMEC), which accumulates all the biases from other variables into H estimation [13]. Similarly, slight overestimation in S-SEBI-derived H suggests that the assumption of zero evaporation over the dry edge ( Figure 4) may not be valid. The fact that SEBS-derived H had a higher difference with measured H can be linked to the sensitivity of SEBS models toward temperature gradient and vegetation properties [33]. In general, SEBS does not require extreme anchor pixels in the H calculation and uses a pixel-by-pixel calculation of H, using the iterative procedure by solving the relationships for the profiles of the friction velocity (roughness parameterization) and the difference between the near-surface potential air temperature and potential surface temperature, which also questions its reliability because of inaccuracies in the temperature gradient and aerodynamics resistance length [33,62]. This also tends to under and overestimate ET inst for dense and sparsely vegetative conditions, respectively, which is also observed in this study, where SEBS under and overestimate ET inst and H, respectively, over a dense sugar beet surface in 2017, and over and underestimate ET inst and H, respectively, when only crop residue is available in 2018 [63]. These results were also supported by Liaqat and Choi [33]: they observed that a −5K difference between the absolute surface and radiometric temperature can result in 107% overestimation in ET inst by SEBS compared to 3% overestimation in ET inst using METRIC. Gokme et al. [64] associate this over and underestimation of H with the fact that most SEB models do not consider the soil moisture dependency and assume the variation of Ts and NDVI as a surrogate for soil moisture, which causes uncertainty in the estimated H. Figure 6 indicates negative values of H (observed in 2017) for both SEB models and BREBS station, indicating the movement of energy from the air to the plant canopy. The already harvested barley field, as well as vast swathes of natural vegetation surrounding the BREBS station, can be a potential source of advective heat, causing negative H values and in turn fulfilling the high ET c demand. However, in case of 2018 and 2019 image dates, a higher P and lower T air might have reduced the effect of advective heat from surrounding areas to the BREBS station. In 2018 and 2019, the average P was 19 mm and 81 mm higher than 2017 and the average T air was 0.50 • C and 0.74 • C lower than 2017, respectively. The ratio of LE/(R n -G) was also observed to be less than 1 on most 2018 and 2019 image dates [25]. A similar variation in modeled H was observed by Wagle et al. [30]: they suggested the importance of adjustment in SEB model for accurate energy partitioning.
For all the SEB models, a higher discrepancy between models was observed on their monthly estimates (Figure 7) as compared to instantaneous estimates (Table 3). This difference can directly be correlated to the difference between using ET r F and Λ for interpolating instantaneous estimates to monthly values. The METRIC model uses ET r F while the rest of the SEB models use Λ. The ET r F is tied down with ground-based ET r and the use of it is expected to result in better monthly estimates when compared with corresponding BREBS values (METRIC RMSD = 18 mm; Figure 7). A comparison of ET r F vs. Λ SEBAL at the BREBS footprint showed an RMSE of 0.18 and percentage bias ranging from −37% (underestimation) to 72%. A similar comparison between ET r F vs. Λ S-SEBI had RMSE of 0.2 and percentage bias between −12% to 67%. Likewise, ET r F vs. Λ SEBS had RMSE of 0.13 and percentage bias between −57% to 9% (data not shown).
It is important to note that various model assumptions, systematic, and unsystematic measured flux bias, scaling issues, user technical skills, and management practices can lead to many uncertainties and inaccuracies in the SEB model comparison. For example, the manual selection of anchor pixels for METRIC and SEBAL can use variation in T s of hot and cold pixels which can result in significant bias in ET inst and H outcomes [65]. Similarly, the sensitivity of the SEBS model to aerodynamic and surface roughness estimation [62] and the assumption of a linear relationship between T s and albedo to define the hot and cold edge in S-SEBI can cause significant bias in the final estimation of H and ET inst . However, the bias in SEB-derived ET inst and H is not only resulted from the uncertainties in SEB model parameters but also by errors in flux measurements [25]. For example, BREBS assumes the eddy diffusivities of heat and water vapor to be equal. Studies [66,67] have revealed that these diffusivities may not be equal in some cases, resulting in a force closure of the surface energy budget. Likewise, Barr et al. [66] reported that BREBS favored the prediction of LE as compared to H.
In general, the spatio-temporal patterns of ET c can be highly variable due to the heterogeneity of land surface and environmental factors that control the SEB of the land-atmosphere system [17]. Since there is no way to validate the models on a pixelby-pixel basis over a heterogeneous land cover, in this study, model performance was evaluated based on the spatio-temporal maps (Figures 8 and 12; Figures S1-S4), density plots (Figures 9 and 10), and absolute difference maps ( Figure 11). The diverse cropping systems and agronomic practices across the study area are responsible for a significant fluctuation in ET c . However, the early-season ET c rate from the cropland can be com-parable to or in some cases less than that of natural vegetation, as evidenced by the unimodal curves (Figure 9) during the early season (15 May and 8 June). The rainfall events early in the growing season induce new growth in natural vegetation and thus help transpire more. However, as the growing season advances, the naturally vegetated area normally suffers from water stress due to scant rainfall, resulting in a lower ET c rate as compared to that of cropland, where scarce rainfall is replenished by irrigation. On the other hand, the higher air temperature and ample amount of solar radiation during mid-season led to the greater availability of energy for crop evapotranspiration in a well-watered crop surface, resulting in higher ET c . Likewise, the image on 26 July ( Figure 9) coincided with the harvesting of barley and alfalfa from cropland, resulting in a multimodal density plot (METRIC) due to heterogeneous surface characteristics. As compared to other models, the METRIC-estimated ET 24 curve showed consistent and anticipated fluctuation over the growing season ( Figure 9). Figure 10 indicates some degree of linearity between METRIC, SEBAL, and S-SEBI models. In SEBS, because of the heterogeneity of land surface and the rather arid climatic condition in the intermountain terrain of Wyoming, the boundary delimitation can undergo some error, resulting in a bit wayward estimation of ET c [17]. Strong linearity between METRIC, SEBAL, and SEBS-derived ET 24 was observed by Singh and Senay [26] over cropland and grassland in the mid-western United States. They reported that high linearity between models is because of the use of thermal data as the main driving factor in estimating ET using the energy balance approach. A higher estimation of seasonal ET c for natural vegetation by S-SEBI ( Figure 12) implies that the T s vs. α relationship utilized in the S-SEBI model to determine the wet edge and dry edge may not be representative of naturally vegetated areas. A similar overestimation of S-SEBI-estimated (16%) ET c was observed by Wagle et al. [30] in the dry 2012 growing season in Oklahoma.

SEB Models Implication
The model selection primarily depends upon its performance on a particular geographical area, data availability, and the expertise of the user [26]. Although we found several discrepancies between models and their comparison with ground-based flux measurements, this study demonstrated the usefulness of SEB models to estimate and map surface energy balance fluxes over heterogeneous surfaces in the semi-arid to arid intermountain region of Wyoming. To further understand the variation on monthly ET c over the intermountain region of Wyoming, Table 4 provides the mean monthly, as well as seasonal ET c , from all the SEB models for sixteen different land cover types during the 2018 growing season (May-September). These monthly and seasonal data can directly contribute to predict and regulate irrigation diversions from a river and aquifer system [68], water allocation in a river and aquifer system [68], groundwater consumption [3,68], determine irrigation reservoir, storage, and conveyance system capacities [3,68], and establishing and regulating water rights [68,69]. Our companion paper [25] utilized the METRIC model (chosen based upon better performance in this study) in quantifying average seasonal water consumption for different cropping systems and each irrigation district within the study area. The paper also quantified the percent of irrigation contributing to seasonal ET c for each cropping system and irrigation district in semi-arid to arid regions of the intermountain region of Wyoming.
The daily, monthly, and seasonal ET c estimates over the cropland, natural vegetation from different models can also be useful in validating the performance of various physically based hydrological models. In general, to achieve the improved performance, hydrological models are often calibrated against the ground measurements, e.g., streamflow measurements. However, in most cases, ground observations are insufficient and of poor-quality. Under such circumstances, satellite-based estimations of water fluxes such as ET c provide valuable information over large geographical areas and diverse land-use types (as presented in Table 4) at regular intervals and with sufficient record length. The selection of an appropriate remote sensing model becomes critical under such scenarios, considering the availability of different remote sensing models differ in model performance, input data requirement, and suitability in a particular geographical area. Over the years, many studies used remote sensing-based modeled ET c ; for example, Uniyal et al., [70] evaluated the Soil and Water Assessment Tool (SWAT) for the upper 300 mm of a soil profile with the indirect measurement of soil moisture estimates from Landsat images in 2016. They used the NDVI, the thermal vegetation difference index (TDVI) and brightness temperature from Landsat images to evaluate the spatio-temporal variation of soil moisture and compared with SWAT output. Similar analysis was performed by Parajuli et al. [71] in northwestern Mississippi, where they used monthly ET c estimates derived from the SEBAL model to evaluate the performance of the SWAT model. Jiang et al. [72] provided the detailed review of studies that integrate physically based process models with remote sensing models (METRIC, SEBAL, MOD16, etc.).

Summary and Conclusions
The SEB model selection is primarily carried out keeping in mind the difference in input data requirements, level of complexity, various assumptions made in the models, and time required to set up a model and produce the final output. This study was conducted to secure a suitable SEB model that is a better fit for the agro-climatic and elevated landscape setting of Wyoming. For that, a total of 19 cloud-free Landsat 7-ETM+ and Landsat 8-OLI and TIRS images were analyzed for the 2017, 2018, and 2019 growing season using four satellite-based energy balance models-METRIC, SEBAL, SEBS, and S-SEBI. The correlation observed between METRIC-estimated and BREBS-measured ET inst had R 2 between 0.21-0.95 and RMSD between 0.07-0.09 mm h −1 for the three growing seasons considered. METRIC vs. BREBS statistics for pooled data points were R 2 = 0.91, RMSD = 0.08 mm h −1 , PBE = 5.7%, and NSE = 0.9. SEBAL had a comparatively lower correlation with BREBS with R 2 falling in between 0.06-0.75 and RMSD between 0.13-0.17 mm h −1 for the three growing seasons. Similarly, SEBAL pooled data points had R 2 of 0.69, RMSD of 0.14 mm h −1 , PBE of 2.2%, and NSE = 0.7. S-SEBI had R 2 ranging between 0.21-0.81 and RMSD between 0.11-0.15 mm h −1 . S-SEBI pooled data had R 2 of 0.76, RMSD of 0.13 mm h −1 , PBE of 3.9%, and NSE of 0.70. Likewise, SEBS, when compared with BREBS flux, produced good correlation (R 2 = 0.72 −0.9, RMSD = 0.09-0.14 mm h −1 ). Pooled data points statistics for SEBS was R 2 = 0.87, RMSD = 0.11 mm h −1 , PBE = 12.3%, and NSE = 0.80. No significant difference was observed between R 2 values of METRIC and SEBS-estimated vs. BREBS ET inst (0.91 vs. 0.87), the RMSD, NSE, and PBE of the SEBS model were 27% larger, 18% lower, and 53% higher compared to the METRIC model. These findings can explain the role of the internal calibration procedure (CIMEC) in METRIC to estimate ET inst more accurately. All the SEB models overestimate ET inst with percent biases ranged from 2.2% for SEBAL to 12.3% for SEBS. Overall, METRIC proved to be a better model for estimating ET c as its RMSD values were lower and consistent for three consecutive growing seasons considered in the research. Comparison over individual vegetative surfaces indicated that the largest discrepancies were observed under drier conditions when only crop residue was available at the surface. On comparing the seasonal outputs, METRIC again was a standout model with relatively low RMSD of 17.6 mm and percentage error of 7.9, followed by SEBS (RMSD = 25.8 mm, % error = 6.6), SEBAL (RMSD = 33 mm, % error = 29) and S-SEBI (RMSD = 34.3 mm, % error = 31). Likewise, a mid-season density plot and absolute difference map showing the model intercomparison revealed the METRIC and SEBAL model were close on their estimates of ET 24 with pixel-wise RMSD of 0.54 mm d −1 and overall absolute difference across the study area of 0.47 mm d −1 . The highest difference was observed between S-SEBI and SEBS with density plot RMSD of 1.21 mm d −1 and absolute difference of 1.23 mm d −1 . Likewise, quantification and mapping of the model-estimated ET 24 and seasonal ET c reflected an anticipated variation across the study area as the growing season progressed, with overall estimates of METRIC being relatively higher as compared to other models. This study was carried out to identify a best-fit model for the intramountainous terrain of Wyoming as well as outline some of the limitations and uncertainties associated with the SEB models on estimating ET c . The results indicated that the METRIC model performed comparatively better in this geographical and agroclimatic setting when model estimates were compared with corresponding BREBS fluxes. However, a pixel-wise density plot and an absolute difference map depicted closeness between the models involving METRIC, SEBAL, and S-SEBI estimates of ET 24 .