Landsat Hourly Evapotranspiration Flux Assessment Using Lysimeters for the Texas High Plains

: Evapotranspiration (ET) is one of the biggest data gaps in water management due to limited ET measurements, and further, spatial variability in ET is di ﬃ cult to capture. Satellite-based ET estimation has great potential for water resources planning as it allows estimation of agricultural water use at ﬁeld, landscape, and watershed scales. However, uncertainties with satellite data derived ET are a major concern. This study evaluates hourly satellite-based ET from 2001–2010 for the growing season (May–October) under irrigated and dryland conditions for both tall and short crops. The evaluation was conducted using observed ET from four large weighing lysimeters at the United States Department of Agriculture Agricultural Research Service (USDA-ARS) Conservation and Production Research Laboratory in Bushland, Texas. Hourly ET from satellite data were derived using the Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) model. Performance statistics showed that satellite-based hourly estimates compared to lysimeter measurements provided good performance with an root-mean-square error( RMSE ) of 0.14 mm, Nash–Sutcli ﬀ e e ﬃ ciency ( NSE ) of 0.57, and R 2 of 0.62 for ET for dryland crops, and RMSE of 0.16, NSE of 0.63, and R 2 of 0.65 for irrigated crops. METRIC provided accurate hourly ET estimates that may be useful for irrigation scheduling and other water resources management purposes based on the hourly assessment.


Introduction
Evapotranspiration (ET) is one of the major components of the water budget. It plays a major role in the water cycle and irrigation water management [1,2]. Further, ET is the main consumer of rainfall and irrigation water in irrigated agricultural fields in arid and semiarid regions of the world. ET measurement methods such as pan evaporation, sap flow, and weighing lysimeters are widely used. These measurements are considered point estimates for particular irrigation practices in a homogenous field. Measuring ET on a watershed scale is challenging due to unavailability of direct measurement methods at that scale [3,4]. Other ET methods such as Bowen ratio, scintillometers, and eddy covariance can provide greater spatial coverage but may not adequately represent a watershed scale. Remote sensing models using satellite data can provide large-scale spatial coverage at various spatial resolutions [5,6]. Considering the limitation that multiple weather parameters are required

Site Description
Four fields equipped with lysimeters at the USDA-ARS Conservation and Production Research Laboratory at Bushland, TX (35.19 • N, 102.10 • W) at an elevation of 1170 m above mean sea level were selected for this study. The research field is a square with a surface area of approximately 20 ha, subdivided into four quadrants of approximately 4.7 ha each. A precision, large weighing lysimeter is located at the center of each quadrant; two lysimeters were managed as dryland (NW and SW), and the other two lysimeters were managed as irrigated (NE and SE) during the study years. The observed data were combined for each lysimeter, and statistical performance assessment was performed for the estimated values for each parameter for the two dryland lysimeters (NW and SW) and the two irrigated lysimeters (NE and SE). The soil characteristics for the study field are deep, well-drained Pullman silty clay loam (fine, mixed, superactive, thermic Torrertic Paleustoll) [18]. The local climate is classified as semi-arid, with large daily temperature variations. Cotton, soybean, grain and silage sorghum, sunflower, and cotton were the predominant crops for the research fields.

Weather Parameter Measurements
Dataloggers (CR6, Campbell Scientific, Logan, Utah) were used to collect data using the following sensors: air temperature and relative humidity sensor (HMP155, Vaisala, Helsinki, Finland), four soil water sensors (Acclima 315, Acclima Inc., Meridian, Idaho), six soil heat flux plates (HFT-3, Radiation Energy Balance Systems, Bellevue, WA, USA), an infrared thermometer, and a net radiometer (Q*7.1, Radiation Energy Balance Systems, Bellevue, WA, USA) [19]. These sensors were utilized to measure surface temperature, soil heat flux, and net solar radiation over each lysimeter. The QA/QC protocol from Evett et al. [20] was performed on the collected data. Weather parameter data were recorded every 6 s and summed or averaged for 30-min intervals.
The solar radiation sensor calibration was performed to assess the sensor accuracy through evaluating the observed R n with the sum of measured downwelling and upwelling short-and longwave radiation as well as comparing to the theoretical maximum clear sky solar irradiance. Based on the sensor calibration results, the sensor performance was accurate in solar radiation measurement.
Colaizzi et al. [21] highlighted the soil heat flux calculation and calibration process for 30-min time intervals. The calibration process was performed using the soil water and temperature sensor data to calculate the change in soil heat storage from the surface to the soil heat flux plates.

Landsat Imagery
Landsat Collection-1 was used for Landsat 5 Thematic Mapper (TM) satellite data throughout the growing season (May-October).
The images were obtained through Earth Explorer (https://earthexplorer.usgs.gov/), with two satellite paths (30 and 31) and row 36 from 2001 to 2010. A total of 53 cloud-free images were acquired, as the presence of clouds can cause errors in the estimation process due to aerodynamics and radiometric disturbance [17]. Supplementary Materials includes day of year (DOY) for images used. The Landsat 5 TM consists of six spectral bands (bands 1-5 and band 7) with spatial resolutions of 30 m and one thermal band (band 6) with a spatial resolution of 120 m. Top-of-atmosphere reflectance was used to calculate NDVI Equation (5) using the red and near-infrared bands [22].

Image Analysis
The Bushland Evapotranspiration and Agricultural Remote Sensing (BEARS) software is an image processing and geographic information system software deriving hourly, daily, and seasonal evapotranspiration maps. It also produces other energy exchange outputs between land and atmosphere using Landsat 5, 7, and 8 [23]. The BEARS software was used for this analysis to estimate the hourly ET, T s , R n , and G o for the four lysimeters. The BEARS software is public domain software that offers users the option to select one of five energy-balance-based ET methods: Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC), Surface Energy Balance Algorithm for Land (SEBAL), Surface Energy Balance System (SEBS), Tow Source Model (TSM), and Simplified Surface Energy Balance (SSEB). The METRIC model has been utilized in this study. To extract the average value of each parameter, a grid of 3 by 3 pixels was selected, with each lysimeter located towards the center of this grid.
The Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) model is a satellite image processing model that is used for ET estimation as a residual of the energy balance for various short crops, tall crops, trees, and forest [24][25][26][27][28][29][30][31][32] Equation (1). METRIC has been employed for Landsat image analysis (with a spatial resolution of 30 m) for agricultural water use and other agricultural applications that require field-scale resolution. The energy balance is evaluated internally under two extreme conditions (dry and wet) using surface temperature, vegetation growth, and available weather data. The accuracy of ET estimation depends on user experience, and the estimation accuracy was found to be inversely related with ET levels, where low ET levels had high uncertainty and vice versa [33]. Those two extremes are the strength of METRIC, compared to the Surface Energy Balance System (SEBS) and other satellite-based ET models [25]. An automated statistical calibration method has been developed to better estimate the dry and wet pixels, enhancing ET estimation accuracy [17].
Atmospheric correction is not required for surface temperature (Ts) and reflectance (albedo) measurements using radiative transfer models due to the use of the indexed temperature gradient dT and the internal calibration of the sensible heat computation within METRIC [34]. Another advantage of using the internal calibration is that it reduces the bias with aerodynamic stability corrections and surface roughness.
Net radiation at the surface (R n ) is calculated by subtracting all outgoing radiant fluxes from all incoming radiant fluxes including solar and thermal radiation Equation (2): where R S↓ is the incoming short-wave radiation [W m −2 ], α is the surface albedo (dimensionless), R L↓ is the incoming long-wave radiation [W m −2 ], R L↑ is the outgoing long-wave radiation [ W m −2 ], ε o is the broad-band surface thermal emissivity (dimensionless), and (1 − ε o )R L↓ is the fraction of incoming long-wave radiation reflected from the surface. The surface temperature (T s ) is calculated for Landsat images using the modified Planck equation, based on Markham [35], with both atmospheric and surface emissivity correction, as shown in Equation (3).
where K 1 and K 2 are constants (K 1 = 607.8 and K 2 = 1261 W m −2 sr −1 µm −1 for Landsat 5), ε NB is the narrow-band emissivity corresponding to satellite thermal sensor wavelength band, and R c is the corrected thermal radiance from the surface using spectral radiance L t,thermal from the thermal band of Landsat. The corrected thermal radiance (R c ) [W m −2 sr −1 µm −1 ] is calculated using Equation (4) [36]: where L t,thermal is the spectral radiance of Landsat 5 thermal band [W m −2 sr −1 µm −1 ], R P is the path radiance in the 10.4-12.5 µm band [W m −2 sr −1 µm −1 ], T NB is the narrow-band emissivity of air (10.4-12.5 µm range), and R sky is the narrow-band downward thermal radiation from a clear sky The normalized difference vegetation index (NDVI) is the ratio of the differences in reflectivity for the near-infrared band and the red band divided by the sum Equation (5).
where ρ t,nir is the reflectance of the near-infrared band, and ρ t,red is the reflectance of the red band. The soil heat flux represents the rate of heat storage in the soil and vegetation as a result of conduction. METRIC calculates G as a ratio G R n using Equation (6) [37] to represent values near mid-day.
where T s is the surface temperature (K), and α is the surface albedo. Alternatively, [38] developed the ratio G R n using soil heat flux data collected by USDA-ARS [39] for an irrigated crop near Kimberly, Idaho, as represented in Equations (7) and (8).
where LAI is the leaf area index, and T s is the surface temperature. METRIC sensible heat flux calculation is estimated from the aerodynamic function as listed in Equation (9). H = ρ air C P dT r ah (9) where ρ air is the air density [kg m −3 ], C P is the specific heat of air at constant pressure [J kg −1 K −1 ], r ah is the aerodynamic resistance [s m −1 ] between two near-surface heights, z 1 and z 2 (always 0.1 and 2 m) computed in a particular pixel, and dT is the near-surface temperature difference (K) between z 1 and z 2 .
In METRIC, the r ah calculation uses wind speed extrapolation from blending heights (normally 100 to 200 m) above the ground surface. dT is used in Equation (9) due to complications in estimating surface temperature accurately from satellites, which arise from uncertainties in atmospheric attenuation or contamination and radiometric calibration of the sensor [25]. The surface temperature obtained by the satellite, whether it is a radiometric or kinematic temperature, can also be different by several degrees from the aerodynamic temperature, the main driver for the heat transfer process [3,40,41]. dT can be estimated using Equation (10) [42].
Where a and b are the empirical constants for a given Landsat satellite image, and T s datum is the surface temperature adjusted to a common elevation for each image pixel using a digital elevation model and customized lapse rate.
Determining accurate hot (dry) and cold (wet) pixels is one of the most critical and challenging steps in implementing METRIC to spatially estimate ET [17,33]. A manual selection method was performed to determine the dry and wet pixels using the surface temperature and the NDVI outputs based on the criteria listed in Table 1. The manual selection process was performed using the surface temperature (T s ) and NDVI histogram distribution. Obtaining the T s distribution is necessary to determine the range of the high and low temperature threshold for each image. Several iterations were implemented to obtain the most accurate hot and cold pixels that met the conditions listed in Table 1. The more accurate the dry and wet pixel determination, the better the ET estimates across the satellite scene. Special consideration of the hot and cold pixel selection, such as avoiding the image edges, as well as the hot and cold pixels distribution across the scene instead of centralizing the pixels toward a specific location was used.

Statistical Analysis
In addition to visual inspection of observed and simulated ET values for pixels covering lysimeters, the satellite-based ET performance was evaluated for the selected day's images that were available. Various statistical parameters were estimated, including Nash-Sutcliffe efficiency (NSE), root-mean-square error (RMSE), and mean bias error (MBE), to evaluate the relationship strength between simulated and measured values [19,43,44]. RMSE and MBE were calculated using Equations (11) and (12), respectively. Linear regression was performed to determine the coefficient of determination (R 2 ), the slope and intercept, and Nash-Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) (Equation (13) [45]. The slope of 1.0, intercept of zero, and R 2 approaching 1.0 indicates a perfect fit. The MBE provides the ability to determine the deviation between the measured and satellite-based estimates, with MBE = 0 indicating no bias in estimation. The NSE values can range from -∞ to +1, with +1 being a perfect agreement between the model and observed data [46].
where y t = the i-th observed value,ŷ t = the i-th simulated value, and n = total number of observations.
where y m is the average measured value, y s is the average simulated value, y m is the measured data on day i, y s is the simulated output on day i, and j represents the rank.

Lysimeter Comparison
The METRIC model [25] was used to estimate ET, and an hourly evaluation was performed for surface temperature, solar radiation, soil heat flux, and evapotranspiration (T s , R n , G o , ET). The evaluation process was conducted for tall crops, including forage sorghum and forage corn; short crops including cotton, sunflowers, and soybeans; and bare soil for the dryland and irrigated lysimeters. Tables 2 and 3 summarize the statistical performance of the METRIC model estimation for T s , R n , G o , and ET versus the observed lysimetric data. The regression line between the observed and simulated values for the dryland and irrigated lysimeters can be found in Figures 1 and 2, respectively.
The observed mean hourly surface temperature value for the dryland lysimeter at the time of satellite overpass was 33.8 • C, and it closely matched the estimated mean value of 33.0 • C. The dryland lysimeter regression line for surface temperature had a good coefficient of determination of 0.76, with a slope and intercept of 0.7 and 9.3 • C. The RMSE was~6.8% (2.3 • C of the mean observed values, with an MBE of 0.8 • C. Figures 1b and 2b represent the R n comparison between dryland and irrigated lysimeters, respectively, compared to the measured data. The dryland model estimated R n value was 521.5 W M −2 , and it closely matched the measured value of 533.5 W M −2 . The irrigated lysimeter estimate was 551.3 W M −2 , and the observed R n was 542.3 W M −2 . The dryland lysimeter model provided good performance for net radiation with an NSE of 0.68 and poor performance for the irrigated lysimeter with an NSE of 0.17. The regression model for dryland lysimeters accounted for 71% of the variability, with a slope of 0.8 and intercept of 98.2 W M −2 ( Table 2). The irrigated regression model captured 22% of the variability, which is lower than for the dryland regression model, with a slope of 0.29 and intercept of 392.3 W M −2 ( Table 3).  Table 2). The irrigated regression model captured 22% of the variability, which is lower than for the dryland regression model, with a slope of 0.29 and intercept of 392.3 W M (Table 3).

Discussion
The hourly comparison for the lysimeter site was performed for surface temperature, evapotranspiration, net radiation, and soil heat flux. The observed surface temperature for the dryland lysimeter was higher than the irrigated lysimeter (Figures 1c and 2c), as expected. The hourly evaluation of the METRIC model illustrated good performance for both water management conditions; however, the irrigated lysimeter overall hourly comparison provided slightly better statistics than the dryland lysimeter. The dryland lysimeter observed and estimated temperatures are higher than for the irrigated lysimeter, and this is due to the irrigation cooling effect.
The observed mean surface temperature value for the irrigated lysimeter was 31.5 ℃, and it closely matches the estimated mean value of 31.2 ℃. The irrigated lysimeter regression line for surface temperature has better estimates than for the dryland lysimeter, with a coefficient of determination of 0.8, a slope and an intercept of 0.63 and 11.2 ℃. The RMSE was 14.6% (4.

Discussion
The hourly comparison for the lysimeter site was performed for surface temperature, evapotranspiration, net radiation, and soil heat flux. The observed surface temperature for the dryland lysimeter was higher than the irrigated lysimeter (Figures 1c and 2c), as expected. The hourly evaluation of the METRIC model illustrated good performance for both water management conditions; however, the irrigated lysimeter overall hourly comparison provided slightly better statistics than the dryland lysimeter. The dryland lysimeter observed and estimated temperatures are higher than for the irrigated lysimeter, and this is due to the irrigation cooling effect.
The observed mean surface temperature value for the irrigated lysimeter was 31.5 • C, and it closely matches the estimated mean value of 31.2 • C. The irrigated lysimeter regression line for surface temperature has better estimates than for the dryland lysimeter, with a coefficient of determination of 0.8, a slope and an intercept of 0.63 and 11.2 • C. The RMSE was 14.6% (4.6 • C) of the mean observed values with MBE of 0.55 • C. These predicted errors are close to those reported in the literature [43,[47][48][49][50]. Moorhead and Gowda have similar unpublished reports for METRIC hourly comparisons. The NSE values of 0.65 and 0.76 for the dryland and irrigated surface temperatures, respectively, are viewed as a good match [27,31,43].
The RMSE for the solar radiation for dryland lysimeters was lower than for irrigated lysimeters, with values of 8.8% (46.9 W M −2 ) and 15.6% (84.8 W M −2 ), for the dryland and irrigated lysimeters, respectively. These errors may have resulted from estimation errors in air emissivity and surface albedo as well as the water content in the air above the irrigated field being higher compared to the dryland lysimeter. This water content may distort the electromagnetic reflectance more than for dry fields, resulting in poor solar radiation estimates. These errors ranged within reported values in the literature [31,43,[51][52][53][54].
As a result of soil heat flux underestimation, the NSE values were 0.17 and −0.47 for dryland and irrigated lysimeters, respectively, which indicated poor performance for the irrigated field compared to the dryland lysimeter. The RMSE and MBE for the dryland were 55% (21.2 W M −2 ) and 4.2 W M −2 , respectively; the RMSE and MBE for the irrigated lysimeter were 89.6 and 0.8 W M −2 , respectively. The dryland regression model accounted for 20% of the variability in measured data with a slope of 0.2 and intercept of 25.7 W M −2 ( Table 2). The irrigated regression model accounted for only 8% of variability, with a slope of 0.3 and intercept of 25.5 W M −2 ( Table 3). The range of satellite-estimated soil heat fluxes is much lower than that with the observed data, indicating METRIC is underestimating soil heat flux, possibly due to the time difference between the physical process of conduction and the instantaneous acquisition of the satellite image. The wet conditions of the irrigated fields may have an impact on the estimation accuracy of R n and G 0 due to the humidity above the soil surface. Accuracy in the estimation of R n and G 0 may be affected by increased humidity near the top of the canopy on the irrigated lysimeter, which may disturb the electromagnetic radiation reflectance compared to the dryland lysimeter. These results agree with [27,31,43].
The METRIC ET estimation provided a good correlation for the dryland model; however, METRIC ET estimation performed better for the irrigated conditions. The hourly ET regression model explained 62% of the ET variability in the observed data for dryland lysimeters, with a slope of 0.68 and an intercept of 0.1 mm h −1 . The regression model for irrigated conditions estimated 65% of the variation in the observed data, with a slope of 0.64 and an intercept of 0.19 mm h −1 .
The dryland RMSE for the hourly ET was 0.14 mm, which was 50% higher than that of the mean observed ET; the MBE was −0.05 mm. The irrigated modeled hourly ET RMSE was 0.16 mm h −1 , about 37.2% higher than that of the mean observed hourly ET; the MBE was −0.04 mm. The irrigated lysimeters error values were less than those of the dryland lysimeter (Tables 2 and 3), possibly due to increased overall ET and less sensible heat flux, which may have caused attenuation in the reflected radiation from the dryland field as captured by the satellite. Overall, these results agreed with [23,45], illustrating that the accuracy level of METRIC hourly T s , R n , and G o for dryland and T s and R n for irrigated lysimeter estimation are considered good for the Texas High Plains, with potential applications for other geographic locations. These results are within the range reported in the literature [12,27,[29][30][31]44,[54][55][56].

Conclusions
The METRIC model was tested against lysimetric hourly surface temperature (T s ), net radiation (R n ), soil heat flux (G o ), and evapotranspiration (ET) in the Texas High Plains. For this purpose, 53 cloud-free images from Landsat 5 TM acquired during 2001-2010 were used. The hourly values were extracted and compared against observed T s , R n , G o , and ET for different (tall and short) crops as well as bare soil conditions managed under dryland and irrigated conditions. The METRIC model performance was good in estimating T s , R n , and ET under dryland conditions and T s and ET for irrigated conditions. The wet conditions in the irrigated lysimeter fields affected electromagnetic reflectance.
The overall accuracy level of the Landsat estimates of the evaluated parameters (surface temperature, net radiation, soil heat flux, and evapotranspiration) indicate that the METRIC model may be suitable for deriving hourly input at a watershed scale with fairly good accuracy.