Performance of the METRIC model in estimating evapotranspiration fluxes over an irrigated field in Saudi Arabia using Landsat-8 images

The accurate estimation of evapotranspiration (ET) is essential for hydrological modelling and efficient crop water management in hyper-arid climates. In this study, we applied the METRIC algorithm was applied on Landsat-8 images, 10 acquired from June to October 2013, for the mapping of ET of a 50-ha center pivot irrigated alfalfa field in the Eastern Region of Saudi Arabia. The METRIC estimated energy balance components and ET were evaluated against the data provided by an Eddy Covariance (EC) flux tower installed in the field. Results indicated that the METRIC algorithm provided accurate ET estimates over the study area, with RMSE values of 0.13 mm h and 4.15 mm d. The METRIC algorithm was observed to perform better in full canopy conditions compared to that in partial canopy conditions. On average, the METRIC algorithm 15 overestimated the hourly ET by 6.6 % in comparison to the EC measurements; however, the daily ET was underestimated by 4.2 %.


Introduction
In the Kingdom of Saudi Arabia (KSA), the agricultural sector consumed about 85 % of the total freshwater used in 2008 (Al-Kahtani and Ismaiel, 2010).This share increased to 90 % by 2012 (Elnesr and Alazba, 2013).Hence, efficient use of water for crop production is essential to fulfilling the needs of the increasing population in the KSA (Hussain et al., 2010;Praveen et al., 2012).Various studies (Kassem and Al-Moshileh, 2008;Atta et al., 2011;Al-Ghobari et al., 2013;Mohammad et al., 2013) recommend the development of advanced irrigation systems for improving agricultural wateruse efficiency in the kingdom.As an example, a reduction in the irrigation water of 30-40 % could be attained when sprinkler irrigation is used instead of traditional methods.An additional saving of 10-25 % can be reached with drip irrigation systems (Rizaiza and Al-Osaimy, 1996).Furthermore, implementation of recent innovative precision irrigation technologies, in conjunction with the accurate estimation of crop water requirements through remotely sensed data, could significantly enhance the efficient use of irrigation water in the agricultural sector.
Evapotranspiration (ET) measurements play a crucial role for water management under hyper-arid conditions, particularly in irrigation scheduling, hydrologic modeling and drought monitoring (Bastiaanssen et al., 2000;Allen et al., 2005;Chavez et al., 2005;Senay et al., 2008;Santos et al., 2010;Anderson et al., 2012;Mkhwanazi and Chavez, 2012;Gowda et al., 2013;Lagos et al., 2013;Moorhead et al., 2013;Yee et al., 2014;Madugundu et al., 2017).In semi-arid climates, ET is characterized as one of the most significant components influencing the hydrologic cycle.Hence, accurate determination of ET is considered as one of the crucial and essential factors influencing the optimal management of crop water use through different irrigation systems (Hoedjes et al., 2008).Also, the accurate determination of the energy balance (EB) components such as sensible heat (H ), soil heat (G) and latent heat (LE) fluxes is of great value for wa-Published by Copernicus Publications on behalf of the European Geosciences Union.
ter management practices in these arid and semi-arid regions (Zeweldi et al., 2010).
The ET can be estimated with different direct or indirect methods such as the lysimeter, the water balance, the Bowen ratio, the eddy covariance (EC) and the scintillometer (SC) methods (Allen et al., 2011a;Rana and Katerji, 2000).Lysimeter and EC methods provide direct measurements of ET, while other methods rely on models to estimate heat fluxes using measurable parameters (FAO, 1977;Drexler et al., 2004;Chavez et al., 2009a).The FAO-56 (Allen et al., 2007) modeling for instance, is widely used in numerous studies.This method consists of estimating crop evapotranspiration (ET c ) for a crop canopy using a reference evapotranspiration (ET r ) and a crop coefficient (K c ).The ETr is computed based on the Penman-Monteith method.The FAO-56 computed actual ET can be used to monitor spatially distributed ET developed with remote sensing (RS) approaches (Courault et al., 2005).
Recently, EC systems have gained a lot of popularity in the determination of ET.Several studies have been conducted to investigate the effectiveness of the EC systems in estimating accurate values of H , LE and ET.Yet, the method is often constrained by the lack of surface energy budget closure.For example, Chavez et al. (2009b) reported that the EC measured EB components (H and LE) underestimated about 30 % of these components as inferred from the large weighing lysimeter on an irrigated cotton field in the USA.Similarly, Twine et al. (2000) reported that carbon dioxide fluxes measured over a sorghum crop with four EC systems (equipped with the same models of instruments), underestimated the reference values with the same factor when EB closure was not achieved.However, the EB closure can be improved by adjusting the EC data through the Bowen ratio (Ding et al., 2010;de Teixeira and Bastiaanssen, 2012).Ding et al. (2010) reported that the EC measured ET (ET EC ) compared to the lysimetric data (ET L ) of maize crop showed an EB closure of up to 84 % for the daytime fluxes.However, the forced EB closure with Bowen ratio data improved the EB closure from 79.2 to 95.2 %.Therefore, adjusting ET values measured by the EC system, using the Bowen ratio method improves the accuracy.
In situ point measurements provide accurate ET values; however, these methods are limited to small areas only.Alternative methods, such as RS methods, have been successfully used for assessing the spatial distribution of ET on larger scales, and over different landscapes (Bala et al., 2013;Farah et al., 2004;Colaizzi et al., 2006) and agricultural fields (Kalma et al., 2008;Gowda et al., 2008).One of the most commonly used RS methods for the determination of ET is the Surface Energy Balance Algorithm for Land (SE-BAL) method.This method has been demonstrated to assess ET with an accuracy ranging from 67 to 97 % (Bastiaanssen et al., 2005).Recently, several SEBAL versions have been developed and successfully implemented in different water management practices (Paul et al., 2013).Ex-amples of such methods include the Mapping Evapotranspiration at high Resolution and with Internalized Calibration (METRIC), the Modified SEBAL (M-SEBAL) method, the Simplified Surface Energy Balance (SSEB) method, the Remote Sensing of Evapotranspiration (ReSET) method, the Surface Energy Balance System (SEBS) method and the Surface Energy Balance with Topography Algorithm (SEBTA) method.The METRIC method, which was developed to determine the quantity and spatial distribution of ET over large areas, is considered as one of the most appropriate models for the continuous estimation of ET over crops during the growing season (Allen et al., 2007;Patil et al., 2015).
Over the past decades, several hydrological studies in the KSA have been directed towards accurate RS-based ET estimates on various spatial and temporal scales.As an example, Mahmoud and Alazba (2016) processed a spatial ET dataset over western and southern regions of Saudi Arabia with SE-BAL on MODIS.However, the quantification of ET from satellite data in conjunction with EC data is still limited in the KSA.As the eastern region of the KSA experiences severe dry conditions, long-term EC data in conjunction with Landsat-8 data will enhance the accurate estimates and distribution of ET over various landscapes and agricultural fields.In view of the pressing need to assess the productivity of agricultural fields in the KSA, this study was undertaken in an attempt of apply the METRIC model on Landsat-8 imagery for the determination of spatial and temporal variability in ET aiming at optimizing the quantification of crop water requirement and the formulation of efficient irrigation schedules.

Study area
The study was carried out on a 50 ha alfalfa field, which was one of the 48 agricultural fields of Todhia Arable Farm (TAF) located about 250 km southeast of Riyadh, the capital city of Saudi Arabia, at coordinates of 24 • 11 00 E and 48 • 56 14.6 N (Fig. 1).The farm was under an arid climate with hot summers (40 ± 2 • C) and cold to moderate winters (15 ± 3 • C) and a mean air temperature of 35 • C. The annual rainfall was about 90 mm, most of which occurred in the period from November to February.The soil of the study area is sandy loam soil with a mean value of soil pH of 7.58 and soil electrical conductivity (EC) of 2.36 dS m −1 .The field is in a flat terrain with slight undulations in the desert environment with an elevation ranging between 329 and 453 m.The major crops cultivated in the study area were alfalfa (Medicago sativa L.), Rhodes grass (Chloris gayana Kunth) and corn (Zea mays.L).Due to the high crop water demand combined with the highly erratic rainfall, irrigation is a prerequisite for crop growth.It is entirely provided using groundwater delivered by center-pivot irrigation systems.The experimental field was cultivated with an alfalfa crop (Green Master) sown on 6 December 2012 at a seeding rate of 20 kg ha −1 , and was irrigated through a center-pivot system using a groundwater of mean values of EC, pH and sodium absorption ratio of 2.917 dS m −1 , 7.82 and 1.42, respectively.In the study region, the alfalfa crop is usually cultivated for 2 years, with, on average, a cutting each 30-35 days during summer periods and every 45-60 days during winter periods.The mean alfalfa hay yield for each cut is estimated to be 2.0-3.0 kg ha −1 during winter and 4.0-5.0kg ha −1 during spring or summer periods (Kayad et al., 2016).

Eddy covariance (EC) flux data
The EC system was installed on 27 May 2013 over the selected alfalfa field.The EC flux tower was powered by solar panels with rechargeable batteries.The meteorological and gas exchange (flux) measurements were made at a height of 3.67 m.The EC was equipped with response sensors: including an open-path infrared gas analyzer, IRGA (LI-7500); a 3-axis ultrasonic anemometer (3D Master Pro), Gill Instruments; a net radiometer (CNR-4, Kipp and Zonen); and a quantum sensor.Soil heat flux plates (HFP01, Hukseflux) were installed at various depths (5, 10, 15, 20, 25 and 30 cm).The sonic anemometer and IRGA were set to record flux data logged at 10 Hz (CR3000) and stored as 30 min files (*.ghg).Data on crop condition, growth and phenological parameters were recorded during the frequent field visits.
The recorded EC data for the period from June to October 2013 was used for this study.The collected 30 min raw datasets ( * .ghgfiles) were processed with the "advanced mode" of an automated software program EddyPro (ver.5.0) developed by LI-COR Biosciences, USA.With the EddyPro, wind speed measurement offsets are processed, axis rotation for tilt is corrected, detrending is made for turbulent fluctuations, the time lag is optimized, spectral corrections are made and footprint estimations are made following the procedures described in the EddyPro software instruction manual (LI-COR).As the sonic anemometer was tilted 243 • towards the north, an angle-of-attack correction for wind com-ponents was also performed.During the filtering process, EC data were excluded whenever there was a condensation, a covariance with missing values (−9999) or periods with incorrect sonic temperatures (i.e., > 50 • C).The missing data were filled with the standard methods (basic interpolation) as described in the manual (LI-COR).The missing data were linearly interpolated when gaps were no longer than 1 hour, whereas gaps due to a malfunction of the EC system for longer than 24 h were not removed.Subsequently, the processed EC data were used for the computation of the ET (i.e., ET EC ).

Footprint analysis (FTP)
As described in Nappo et al. (1982), the footprint (FTP) is "the extent to which a set of measurements was taken in a given space-time domain".In this study, an effort was made to understand the spatial representation area or FTP of flux measurements from the EC system.The FTP mainly depends on the function of surface roughness with respect to the wind speed, direction and wind shear.Thus, the FTP for a single 30 min data record will be unique and may vary for the next segment as atmospheric conditions change.Moreover, the shape and length of the FTP may vary with upwind direction, as well as the relative weights (magnitude of flux contributions); in each area inside the FTP, the weighted FTP function is used to attribute the measured flux to a weighted areal estimate (Gockede et al., 2006).The areas very close to the EC tower may contribute less to the total flux sensed by the instrument, whereas areas away (upwind) from the EC tower contribute significantly (up to a point where a peak is reached).Thereafter, the contribution may decrease rapidly (Zhao et al., 2014).
Depending on the complexity, numerous FTP estimation models have been presented in the literature.The widely used Flux Source Area Model (FSAM), described by Schmid (1994), was used in this study for the determination of FTP.The FSAM was computed based on contributions (up to 90 %) of the sensed fluxes by the EC tower.Subse-  quently, the obtained FTP was used to compare the EC tower and Landsat-8 estimates of energy fluxes.As described by Schmid (1994), the EC contributions starting 1.5 h after sunrise and ending 1.5 h before sunset were considered and the FTP was computed for every 30 min corresponding to the averaging period used for the EC tower analysis.The FSAM variables such as observation point height, Obukhov length, the standard deviation of lateral wind speed fluctuations and friction velocity were taken from the EC measurements.Surface roughness length was estimated by measuring the crop (alfalfa) height at regular site visits using linear interpolation between observations.Crop height was transformed into roughness length as described by Allen et al. (1998).
For the FTP analysis, the utilized wind parameters such as wind speed, wind direction and yaw angle at the time of the Landsat-8 overpass are provided in Table 1.The produced oval-shaped "fetch area" was oriented to the mean wind direction corresponding to that half hour time-period.A weighting matrix for the FTP surrounding the EC tower was developed for each day.Subsequently, a 30 m grid matching to the Landsat-8 pixels was created for each 30 min interval, with a binary coding.Grid cell pixels which fell in the "fetch area" were coded as 1, and the pixels outside the fetch area were given the value of 0. The daily FTP fetch calculation was computed by scoring the grid cells coded with 1 divided by the total number of observations recorded during the day.

Energy balance (EB)
The EC system allows accurate measurements of the EB components.Hence, using the measured EB components, the surface energy budget can be estimated (Eq. 1) along with the EB closure.The four flux components (R n , H , G and LE) measured by the EC system were assessed for the EB closure (Kustas et al., 2005).The EB closure rate for the EC system ranged between 71 and 99 %, with a mean value of 87 % (Table 2 and Fig. 2).The EB closure was significantly influenced by the amount of unaccounted energy lost (i.e., lack of closure in the EB equation due to advection and measurement errors in the field) in the partition of the net radiation into latent, sensible heat and soil heat fluxes.It may be due to the portion of energy that is missing in the budget.It was accounted for between 5 and 20 % of the total R n .This was attributed to the inherent errors in measuring the EB components by the ground systems, which was reported to be 15-20, 5-10 and 20-30 % according to Weaver (1990), Field et al. (1994) and Twine et al. (2000), respectively.

Landsat-8 images and preprocessing
Eight cloud-free Landsat-8 TIRS (  (ver.5.1).Subsequently, spatially distributed ET fluxes were weighted and integrated using an FTP for comparison with ET measured with the EC system.

Land surface temperature (LST) -split-window algorithm
Landsat-8 TIRS bands (10 and 11) and OLI bands (2 to 5) were utilized in the estimation of LST through the splitwindow algorithm executed in ENVI (ver.5.1).During the process of LST estimation, the algorithm used brightness temperatures (TB) of two TIRS bands, and mean and the difference in land surface emissivity (LSE) of an area under observation.The split-window algorithm for LST determination is given in Eq. ( 1) as described by Skokovic et al. (2014).
where LST is the land surface temperature (K), C 0 to C 6 are the values of the split-window coefficients (Table 4), TB 10 and TB 11 are brightness temperatures of Landsat-8 bands 10 and 11 (K), ε is the mean LSE of TIRS bands, W is the atmospheric water vapor content, and ε is the difference in LSE.
The brightness temperature (TB), a microwave radiation radiance traveling upward from the top of Earth's atmo- sphere, was calculated using Eq. ( 2).
where K 1 and K 2 are thermal conversion constants of the TIRS bands (Table 5) and L λ is the Top-of-Atmospheric spectral radiance.This Top-of-Atmosphere spectral radiance was determined by multiplying multiplicative rescaling factors (Table 5) from Eq. ( 3).
where M L is the band-specific multiplicative rescaling factor (radiance_mult_band_10/11), Q cal is the band 10 or 11 image, A L is the band-specific additive rescaling factor (ra-diance_add_band_10/11). Subsequently, LSE was calculated using Eq. ( 4).The ε s and ε v are soil and vegetation emissivity values of the corresponding bands (Table 5).The fractional vegetation cover (FVC) was estimated based on the Normalized Difference Vegetation Index (NDVI) obtained over the experimental area (Eq.5).
where NDVI s and NDVI v are the reclassified NDVI for soil and vegetation areas, respectively.OLI bands 2, 3, 4 and 5 were layer stacked and NDVI was calculated using ENVI (ver.5.1) software.The output value of NDVI ranged between −1 and 0.59.To get NDVI s and NDVI v , the NDVI image was reclassified into soil and vegetation.After generating LSE for both bands of the TIRS, the mean and difference LSE were obtained according to Eqs. ( 6) and ( 7).

METRIC algorithm and ET estimation
The METRIC algorithm has been developed exclusively for the estimation of ET from Landsat data (Allen et al., 2005).A total of eight Landsat-8 images were processed and used for the estimation of ET (ET METRIC ) over the experimental alfalfa field.During the ET METRIC computation, surface characteristics such as surface albedo, vegetation indices, surface emissivity and surface temperature were estimated as intermediate products.Anchor pixels (hot and cold) were selected, and the energy components such as the net radiation (R n ), the soil heat flux (G) and the sensible heat flux (H ) were estimated as well.Finally, the latent heat flux (LE) was predicted as a residual of the land surface balance (Allen et al., 2005(Allen et al., , 2007)), see Eq. ( 8).Consequently, the instantaneous ET (ET inst ) for each pixel was calculated.Moreover, the obtained leaf area index (LAI G ) and the reference ET (ETr) over the alfalfa field were used as inputs to the METRICbased Landsat-8 ET prediction (Table 6).The ground-based LAI (LAI G ) of alfalfa was recorded on the day of satellite overpass using a plant canopy analyzer (LAI-2200).
The first step in the METRIC model was to compute the net radiation (R n ) using the surface radiation balance (Eq.9).The R n estimation was accomplished in a series of steps by summing up the net shortwave radiation and net longwave radiation (Hipps, 1989;Brunsell and Gillies, 2002;Allen et al., 2007).
where R S↓ is the incoming shortwave radiation (W m −2 ), α is the broadband surface albedo (dimensionless), and R L↓ and R L↑ are the incoming and outgoing longwave radiation (W m −2 ), respectively.ε o is the broadband surface thermal emissivity (dimensionless).The (1-ε o ) R L↓ term represents the fraction of incoming longwave radiation reflected from the surface.
The incoming broadband and shortwave radiation (R S↓ ), which represents the principal energy source for ET, is calculated for the Landsat-8 image with time as a constant for the whole image assuming clear sky conditions as per Eq.(10).
where G SC is the solar constant (1367 W m −2 ), θ ref is the solar incidence angle, τ sw is the broadband atmospheric transmissivity and d 2 is the square of relative Earth-Sun distance.The τ sw is calculated using Eq. ( 11) drafted in the ASCE-EWRI (2005).
where P is the atmospheric pressure (kPa), W is the amount of water present in the atmosphere (mm), Z is the solar zenith angle (extracted from the image metadata) and K t is the air turbidity coefficient (K t = 1.0 for clean air and 0.5 for extremely turbid or polluted air.K t = 1.0 was used in this study).P and W are calculated using the measured or estimated near-surface vapor pressure, as per Eqs.( 12) and ( 13 where the constant 293 is the standard air temperature (K), z is the elevation above the sea level (m) and the e a is the nearsurface vapor pressure (kPa).The parameter d 2 was computed, from Eq. ( 14), as a function of the day of the year (DOY) as described in Duffie and Beckman (2013).
where α toa is the planetary albedo of each pixel; α atm atmospheric albedo and τ sw is obtained from the Eq. ( 11) following da Silva et al. (2016).
Outgoing longwave radiation (R L↑ ) emitted from the surface is driven by surface temperature and surface emissivity.The R L↑ is computed using the Stefan-Boltzmann Eq. ( 16).
where ε o is the broadband surface emissivity (dimensionless), σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ) and T s is the surface temperature (K).In this study, T s was accounted for as LST and obtained from Eq. ( 1).The surface emissivity was computed using an empirical Eq. ( 17) after Tasumi et al. (2008) based on soil and vegetative thermal emissivities.The LAI is computed as per Eq. ( 18) proposed by Bastiaanssen et al. (1998).During the correction, the typical bare soil and the fully vegetated surface values were set as 0.93 and 0.98, respectively.The soil-adjusted vegetation index (SAVI) was calculated based on TOA reflectance of bands 4 and 5 (Huete, 1988).
The incoming longwave radiation (R L↓ ), a downward thermal radiation flux from the atmosphere (W m −2 ), was estimated using the Stefan-Boltzmann Eq. ( 19) described in Allen et al. (2007).
where τ sw is the broadband atmospheric transmissivity calculated from Eq. ( 11).
where T s is the surface temperature (K) and α is the surface albedo.Subsequently, the G METRIC was obtained by multiplying G/R n with R n .
The sensible heat (H METRIC ) was estimated from an aerodynamic function as expressed in Eq. ( 22).In the calculation of r ah , wind speed measurements were used.
where ρ is the air density (kg m −3 ), C p is the specific heat capacity of the air (J kg −1 K); T is the near-surface air temperature gradient and r ah is the aerodynamic resistance for heat transfer (S m −1 ) between two near-surface heights (i.e., at alfalfa canopy height and the EC measurement height).For further details of the H METRIC algorithm, we refer to Allen et al. (2007).

Calculation of ET
Following the establishment of R n , G and H from the Landsat-8 processing, LE was calculated as a residual of the EB equation.The obtained LE is equivalent to the ET inst at the time of Landsat-8 overpass (Eq.23).
where ET inst is the instantaneous ET (mm h −1 ), 3600 converts from seconds to hours, ρ ω is the density of water (∼ 1000 kg m −3 ) and λ is the latent heat of vaporization (J kg −1 ) representing the heat absorbed when a kilogram of water evaporates.The λ component was computed as per Eq. ( 1).
Finally, as presented in Eq. ( 25), the reference ET fraction (ET r F) was calculated as the ratio of the computed ET inst from each pixel to the reference ET (ET r ) computed from the weather data.
where ET inst is from Eq. ( 23) and ET r is for the standardized 0.5 m tall alfalfa at the time of the image.The EC systemrecorded weather parameters were used to calculate ET r as described in ASCE-EWRI (2005).The obtained ET r F was subsequently extrapolated to daily values.In the processes, ET 24 was computed by assuming that the instantaneous ET r F computed at satellite overpass is the same as the average ET r F over the 24 h average (Allen et al., 2007;see Eq. 26).
where C rad is a correction term used to correct for variation in 24 h versus instantaneous energy availability.

Data analysis
The spatially-estimated energy flux components, derived from the METRIC algorithm and the EC flux tower, were subjected to heat flux correction, energy balance and FTP analyses, as described in Schmid (1994).Different statistical performance indicators (RMSE, mean bias error (MBE) and Nash-Sutcliff coefficient).As the samples were limited, the Mann-Whitney U test and/or the Kruskal-Wallis H test (Gisondi et al., 2004;McCune and Grace, 2002) were performed for the assessment of the METRIC performance in estimating ET against the EC system.These statistical indicators are often used when small sample sizes are considered.
3 Results and discussion

Footprint analysis (FTP)
The simple arithmetic averages of weighted and integrated heat fluxes over the fetch areas are the FTPs (Chavez et al., 2005).FTPs are widely used in validating the spatially distributed fluxes obtained from remote sensing (RS) with EC measured fluxes.The fetch area of the FTP was classified into six classes, based on the cumulative contribution of the eddies, as provided in Fig. 3. Based on the FTP analysis, about 90 % of the EC system observes eddies originating from the 122 to 144 m in the upwind direction.The peak (vertical) fetch is within 42 to 51 m.It is also observed that 10 % of the fluxes are registered between 6 and 17.5 m from the flux tower.About 30, 50 and 70 % of the fluxes are recorded between 35 to 63 m, 52 to 97 m and 81 to 136 m, respectively.Energy fluxes are analyzed as shown by Chavez et al. (2005) by integrating Landsat-8 image acquisition (i.e., overpass).This allows the EC footprint to overlap the METRIC footprint by 90 % (Kustas et al., 2005).

Surface temperature (T )
The linear regression analysis of Landsat-8-derived surface temperature (T LST ) against the EC flux tower measured temperature (T EC ), obtained from the upwind longwave component measured by the CNR-4, shows a good correlation with an R 2 value of 0.71 (p = 0.0084; Fig. 4).However, the T LST is underestimated compared to the T EC as evidenced by the RMSE and MBE performance indicators of The recorded errors (4.2 to 19.6 %) in the T LST are slightly higher than the values reported in previous studies.For example, compared to the ground-level infrared thermometer, Chavez et al. (2009a) reported METRIC errors of 11.1 and 1.9 % in estimating surface temperatures for corn and sorghum fields, respectively.The LAI significantly affects the accuracy of T LST .Vegetation with higher LAI records lower temperatures because the amount of heat stored is reduced through transpiration (Omran, 2012).In this study, the Landsat-8-derived T LST is underestimated by about 20 % at full foliage cover of alfalfa crop (i.e., LAI ∼ 6) compared to T EC .High LAI surfaces can trigger large coherent eddies that are efficient in heat convection, whereas low LAI surfaces are less efficient in generating energetic eddies (Voogt and Grimmond, 2000).Thus, the temporal variation in LAI would re-  sult in oscillations of the land surface temperature.As an example, in the case of low LAI (< 1.23) conditions, the T LST was about 40.99 • C, whereas at higher LAI (5.92) the T LST was 32.6 • C.However, there is a discrepancy at higher temperatures when the crop density of alfalfa was low (i.e., bare soil was visible to the radiometer) as evidenced by the low LAI values over the footprint and the recorded T EC and T LST (Chavez et al., 2005).Hence, the cooling effect of vegetation on the T LST accelerated the error by 1 to 20 % on both the low LAI and high LAI conditions.

Energy balance (EB) components
The mean values of EB components obtained from both the EC system and the METRIC algorithm are provided in Table 7.The performance indicators (RMSE, MBE, Nash-Sutcliff coefficient, Mann-Whitney U test and Kruskal-Wallis H test) are given in Table 8.An illustration (Fig. 5) of the comparison of EC and METRIC estimates of EB components are provided.The temporal trend of EC (CNR-4) measured R n (R nEC ) was analyzed for discrepancies.About 9 to 16 % of the collected data exhibit infrequent errors.Hence, the abnormal data have been discarded, and a gap-filling process was performed using the EddyPro software program (version 5.0).The scatter plot (Fig. 5) represents the relationship between R nEC and METRIC-estimated R n (R nRS ).The moderate linear correlation (R 2 = 0.54 and p > F = 0.0381), low RMSE (18.32 W m −2 ; 3.8 %) and very low MBE (8.66 W m −2 ; 1.76 %) indicate that the METRIC model accurately (96 %) estimate the R n .The obtained RMSE and MBE values are in agreement with those reported by earlier studies of Mkhwanazi and Chavez (2012), where the RMSE value was 4.1 % and MBE value was 3.3 %, and Chavez et al. (2007), where the RMSE value was 9.8 %.In addition, the performance of METRIC was significant with a Kruskal-Wallis H tests (χ 2 = 0.893).However, the null hypothesis cannot be rejected with the Mann-Whitney U test (U and Z value are 23 and −0.892, respectively).These variations in R n between METRIC and CNR-4 estimates are likely due to variations in the computation of LAI along with LST and soil moisture conditions (Zhang et al., 2013).

Soil heat flux (G)
Figure 5 illustrates a scatter plot between the EC measured (G EC ) and the METRIC-estimated (G RS ) soil heat flux values.The relationship was relatively fair with an R 2 value of 0.67 and a Nash-Sutcliffe coefficient of 0.59.The RMSE and MBE, however, were determined at 28.46 W m −2 (37.33 %) and 12.42 W m −2 (16.29 %), respectively.These recorded errors were relatively higher than the values reported in similar previous studies.Mkhwanazi and Chavez (2012) reported an RMSE value of 14.2 W m −2 (27.6 %) and an MBE value = 0.397) and the null hypothesis was not rejected.Furthermore, the ratio of R n to the G (G/R n ) derived from the METRIC model against the LAI inferred from the Canopy analyzer (PCA 2200, LI-COR, USA) is presented in Fig. 6.The G/R n is one of the essential components in the analysis of the accuracy of Bowen ratio (Allen et al., 2011b).Correlating the LAI with the G/R n produced a moderate polynomial (3rd order) relationship for both the METRIC (R 2 = 0.288 and p > F = 0.23) and the EC (R 2 = 0.313 and p > F = 0.15) methods.However, the results suggested that beyond a certain value of LAI (up to 4.2), the relationship between G/R n and LAI was decreased.The scatter in the G/R n , where the LAI increased to 4, showed an increasing trend as the values of G/R n for full canopy cover are ranging between 0.05 and 0.15 as shown in Waters et al. (2002).

Sensible heat flux (H )
After the adjustments of EB closure, the error in recorded H EC flux ranged from 10.89 W m −2 (low LAI conditions) to 38.91 W m −2 (full canopy condition).It may be due to advection, which varies very strongly in the hyper-arid environment.The H METRIC estimated with an MBE value of 15.72 W m −2 (13.87 %) compared to H EC .The high RMSE value of 72.01 W m −2 (63.54 %) for the H METRIC might be due to the selection of anchor pixels especially during 22 August, where the Landsat-8 image was viewed at the time of irrigation.Hence, the advection and variability in the wetness across the fetch or FTP area and canopy reflectance affected the calculations in RS estimation of H METRIC (Brutsaert and Stricker, 1979;Yang et al., 2014).During the earlier (June) and late summer (August) periods, the observed  sharp humidity variations are linked to changes in wind direction.However, during the alfalfa post-harvest practices, the fields are under fallow condition.Hence, the LE flux was always less as most of the available energy partitioned as H rather than LE flux.This was evident in the linear regression analysis (Fig. 5), where a good correlation between the H RS and H EC (R 2 = 0.61) was observed; however, it was not significant (p > F = 0.021), and it was also confirmed with the MBE of 13.87 %.The obtained results resemble the reported values on H METRIC performance (MBE of 10 %) by Carrasco-Benavides et al. (2013).

Latent heat flux (LE)
A scatter plot was established between the METRICestimated LE RS and the EC-determined LE EC values (Fig. 5).The correlation was relatively good (R 2 = 0.66); however, it was not statistically significant (p > F = 0.14).This is attributed to the fact that the LE RS is calculated as a residual component of the EB equation that mainly depends on the calibration accuracy of the H RS .The value of the latter is purely based on the quality of the selected anchor pixels (Weaver, 1990;Field et al., 1994).
Although the regression analysis showed a non-significant correlation, the performance indicators such as the MBE (2.45 W m −2 ; 0.73 %) and the RMSE (115.04W m −2 ; 34.33 %) are in accordance with the earlier reported values indicating that the METRIC algorithm is rather accurate for determining LE over large areas.Carrasco-Benavides et al. ( 2013) stated that the METRIC algorithm overestimated the LE by 14 % (RMSE).The obtained average absolute error of the LE RS for the eight images is 35 % in our study.The obtained inaccuracy may be due to the energy partition at peak growth stages, where the LE accounted for 80-85 % of the net radiation.In addition, the H component accounted only for 5 to 8 % of the R n during the growing season on the daily scale.For the fallow periods, the recorded partitions of LE and H were in the range of 6.9 to 8.2 % and 77 to 84 %, respectively.On the monthly time scale, H and LE fluxes varied throughout the cropping periods.The METRIC algorithm is therefore observed to perform relatively better in full canopy conditions compared to that in partial canopy conditions.Due to the cooling effect of vegetation, the selection of anchor pixels for H calculation is challenging in continuous irrigation regimes.
The LE was high during the full coverage of LAI and irrespective of alfalfa growth stages.Three environmental factors including atmospheric water demand, humidity and wind speed during the growing season have a high impact on the LE measurements, which in-turn triggered the computation of anchor pixels especially at the time of irrigation.Moreover, in the fallow period, H can be accurately determined with METRIC by selecting the targeted pixels after the harvest of alfalfa.Advection in dry months is common in hyperarid regions, which affect the H computation, which in-turn affects the LE and ET estimates on a regional and long-term basis.

Evapotranspiration (ET)
The correlation between the ET values obtained from both the METRIC algorithm and the EC system is found to be highly significant (R 2 = 0.93 and p > F = 0.0001) as illustrated in Fig. 7. Further assessment of the accuracy of the METRIC algorithm in estimating the ET is performed using the RMSE and MBE indicators.The performance of the METRIC model in estimating ET on hourly and daily intervals was significant with the Mann-Whitney U test (U and Z values are 32 and 0.052, respectively).Similarly, the null hypothesis is not be rejected with the Kruskal-Wallis H test (χ 2 = 0.916).Comparing the hourly ET calculated using the METRIC algorithm to that using the EC system results in an RMSE value of 0.13 mm h −1 (25.91 %) and an MBE value of 0.04 mm h −1 (6.6 %).For the daily mean ET, the RMSE and MBE values were 4.15 mm d −1 (34.33 %) and 0.38 mm d −1 (4.2 %), respectively.These results are in agreement with the results reported by Mkhwanazi and Chavez (2012), where the performance errors in estimating the ET using the MET-RIC model were determined at 0.14 mm h −1 (17.6 %) and −0.08 mm h −1 (−10.3 %) for RMSE and MBE, respectively.Similarly, Chavez et al. (2007) reported a METRIC error of 0.7 mm d −1 (7.4 %) in predicting the ET compared to lysimeter data.On average, the present study found that the MET-RIC algorithm overestimated the hourly ET by 6.6 % in comwww.hydrol-earth-syst-sci.net/21/6135/2017/ Hydrol.Earth Syst.Sci., 21, 6135-6151, 2017 parison to the EC measurements.The daily ET is underestimated by 4.2 % (Fig. 8).There are fluctuations in EC and Landsat-estimated latent heat flux (Fig. 9).This might be due to the advection process in case of hourly ET assessment and the variability in the canopy density with respect to the studied footprint.Advection may vary strongly in hyper-arid environments.

Conclusions
The METRIC algorithm was applied on Landsat-8 images for mapping ET.Its performance was evaluated against the EC flux tower measured EB components and ET at hourly and daily intervals over an irrigated alfalfa field in the eastern region of Saudi Arabia.The following are the specific conclusions of this study: -The METRIC algorithm was successful for estimating the ET with an average error of 6.6 % (hourly ET) and 4.2 % (daily ET).The performance of the METRIC algorithm was found to be more accurate in estimating the hourly ET (R 2 = 0.81) than the daily ET (R 2 = 0.66) compared to the ET EC .
-Compared with the EC data, the estimated R n component from remote sensing was found to be highly accurate (an accuracy of more than 95 %) was obtained.
In addition, the H obtained from satellite data was associated with errors ranging between 10.89 W m −2 and 38.91 W m −2 at low (< 1.2) and high (> 5.0) LAI values, respectively.The corresponding LE estimates have a low mean value of MBE of 2.45 W m −2 (0.74 %).
-The METRIC algorithm was observed to perform relatively better in full canopy conditions compared to partial canopy conditions.Due to the cooling effect of vegetation, the selection of anchor pixels for H calculation is challenging in continuous irrigation regimes.

Figure 1 .
Figure 1.Location map of the study area.

Figure 2 .
Figure 2. Flux Source Area Model (FSAM) footprint, as percentage of the fetch area, of the eddy covariance (EC) system overlaid on a LAI map of the alfalfa crop (date of Landsat-8 overpass -21 July 2013).

Figure 3 .
Figure 3. Energy balance (EB) closure of corrected eddy covariance (EC) data for the studied Landsat-8 date of overpass.

Figure 4 .
Figure 4. Comparison of EC measured (T EC ) and remote sensing derived (T LST ) surface temperature.

Figure 5 .
Figure 5.The relationship between EB components (W m −2 ) measured by the METRIC algorithm and the EC system.

Figure 6 .
Figure 6.LAI and G/R n relationship of the study crop for remote sensing (RS) and eddy covariance (EC).

Figure 8 .
Figure 8.The relationship between daily estimates of ET RS and ET EC .

Figure 9 .
Figure 9. Comparative assessment of eddy covariance (EC) flux tower measured against the METRIC-estimated latent heat flux (LE, W m −2 ) over the study period.(The values on the x axis are hour of the day in local time.)data for the studied Landsat-8 date of overpass.

Table 1 .
Details of wind parameters at the time of satellite overpass.

Table 2 .
Energy balance (EB) analysis of EC measured heat fluxes.

Table 3 .
Details of Landsat-8 data used in the study.

Table 5 .
Thermal conversion constants and rescaling factors of Landsat-8 TIRS bands.Landsat-8 TIRS Thermal constants a Rescaling factor a Emissivity values b

Table 6 .
Input parameters used for the METRIC model for the alfalfa crop.DOY T EC ( • C) T LST ( • C) LAI G ET r u (m s −1 ) h c (m) Z m (m)EC -temperature from flux tower; T LST -estimated temperature from Landsat-8 images; LAI G -ground measured leaf area index; ET r -reference ET; u -average horizontal wind velocity; h c -alfalfa canopy height; Z m -height of wind speed measurement.

Table 7 .
EC measured and METRIC-estimated Energy Balance (EB) components over an alfalfa field.

Table 8 .
Performance indicators results for the METRIC algorithm estimated EB components.
Z values are reported in the parenthesis; * not-significant; * * significant; χ 2 -is the Chi-square #