A comparison of actual evapotranspiration estimates based on Remote Sensing approaches with a classical climate data driven method

The knowledge of actual evapotranspiration at farm level is a prerequisite for irrigation planning, farm management, to increase production and reduce water consumption. To accomplish this, comprehensive and accurate assessment methods should be applied. In order to evaluate accurately evapotranspiration processes we compared lysimeter evapotranspiration data with MODIS (Aqua and Terra satellites) and LANDSAT (SEBAL algorithm) satellite images as well as with the FAO Penman-Montith method. The findings indicate the low error rate, high correlation (1) and appropriateness of SEBAL in estimating actual evapotranspiration. The error values MAD, MSE and RMSE between lysimeter and the SEBAL algorithm were 0.59, 0.36 and 0.60 respectively. The second best performance was established for the FAO Penman-Montith method. The obtained error values MAD, MSE and RMSE between the lysimeter and FAO-Penman-Montith method are 0.91, 1.29 and 1.13, respectively.


Introduction
Water need is one of the most important parameters in crop cultivation and in terms of planning the irrigation calendar. Water Source deficit estimations are one of the major challenges in dry and semi-arid regions like Iran which is due to the low amount of precipitation (248 mm), high temperatures (average temperature 18 °C, which is 3 °C higher than the global average) long dry season (in some areas up to 8 months), high evaporation, inappropriate cultivation pattern and irregular irrigation methods (Alizadeh 2016; Zakerinejad and Masoudi 2019). In the current condition, the phenomenon of global warming and the occurrence of severe and continuous droughts and desertification aggravate the problem (Zakerinejad and Maerker 2015). Consequently, due to the high impact of evapotranspiration processes in plant and water resources management we assess different models for the Iranian conditions. Especially, we focus in our study on the main crop cultivation in Iran, which is namely wheat.
Evapotranspiration is highly affecting the hydrological cycle and water balance equations. Measuring, calculating and estimating the evapotranspiration volume are essential in water resource management. Different methods exist on the marked such as: direct measurements (Lysimeter), multiple models like Priestley-Taylor, Jensen-Haise, Thornth-Waite, Blaney-Criddle, FAO-Pennman Monteith, Hargreaves-Samani, Turc, Making and Ritchie (Allen et al. 1998;Zare Haghi et al. 2016). In addition, evapotranspiration can be estimatedusing remotely sensed data and respective modelling approaches such as SEBI, SEBAL, S-SEBI, SEBS, METRIC, S-TSEB and P-TSEB (Alizadeh et al. 2016). Moreover, the MODIS sensor also measures evapotranspiration, that is represented in a 8-day composite dataset.
The results of many studies in different countries like China, Poland, Slovakia, Iraq, and Brazil indicate that the SEBAL algorithm is suitable to estimate evapotranspiration even in areas with climate data shortage (Santos et al. 2017;Ndou et al. 2018;Li et al. 2013;Santos 2017;Jaber 2016;Jian 2015;Bezerra 2015;Sun 2011). The MODIS evapotranspiration product provides significant information on variations of evapotranspiration over a wide area. Extensive research has been done in this context by e.g. Yu et al. (2019), Rasmussen et al. (2014), or Sun et al. (2012. In this context, the results obtained through the SEBAL algorithm where compared with experimental methods like Hargrevi-Samani, Blaney-Criddle, FAO Penman-Monteith, Metric, SWAT and lysimeter data in countries like Turkey, India and some cities in Iran (Zanjan, Mazandaran, Neyshabur). Particularly, the SEBAL algorithm has been compared with actual lysimeter data showing small errors, which mainly related to the determination of the cold and hot pixels. The low value of MSE, MAE, MAD, and RMSE obtained in different other studies sustain the findings mentioned above (e.g. Karbasi et al. 2016;Ghorbani et al. 2015;Morshedi et al. 2016;Rezaei Banafsheh et al. 2014;Kamali and Nazari 2018;Atasever and Ozkan 2018;Fu et al. 2018;Rawat et al. 2017). The results obtained by Wagle et al. (2017) on the operation of five of surface energy balance SEB, models of (SEBAL), (METRIC), (SEBS), (S-SEBI (SSEBop)) for evapotranspiration of sorghum prediction indicate that the S-SEBI, SEBAL and SEBS outperform METRIC and SSEBop models with higher accuracy.
In this study we estimate evapotranspiration through: i) the SEBAL algorithm, ii) the FAO-Penman-Montith method and iii) MODIS evapotranspiration products and compare the obtained results con observed lysimeter data. The study area is located in the Shahrekord plain that is characterized by a temperate climate and wheat cropping as dominant agricultural production.

Study area
The Centre of Shahrekord plain is located at 32°29' to 32°38' N and 50°46' to 50°55' E at 2066 m above sea level (Fig. 1). The annual average precipitation of the plain is 330 mm and the annual temperature average is of 12 °C. The test farm is located at Farrokhshahr Agricultural Meteorological Research Center (AMRC). The farm is equipped with a drainage lysimeter with a diameter of 3 m and area of 7.60 m 2 and a cover consisting of 1200 wheat seeds. This farm is one of the experimental farms that estimats the actual evapotranspiration data through the SEBAL algorithm and we compare the results with that of the nearest wheat farm (Fig. 1).

Method and materials
The method applied in this study is comparative and illustrated in the following.

Method
According to the available databases, the observational data (lysimeter) is used a reference data. We estimate the evapotranspiration processes following the SEBAL algorithm given by Eq. 1-20, and the FAO-Penman-Monteith model reported by Eq. 21. To test the accuracy of these models and select the optimal model the RMSE Eq. 23, MES Eq. 24, MAD Eq. 25 and R Eq. 22 indexes are calculated.

Materials
The data bases consist of three data groups: i) meteorological data collected from Farokhshahr station which include average, minimum and maximum temperature, average, minimum and maximum precipitation, relative humidity (RH), wind speed, sunny hours, ii) evapotranspiration data measured by the lysimeter, iii) Landsat satellite images (2016)(2017), and vi) the MODIS evapotranspiration product (2016-2017) Where λET is the latent heat flux (W/m 2 ), R n is the net radiation flux at the surface (W/m 2 ), G is the soil heat flux (W/m 2 ), and H is the sensible heat flux to the air (W/m 2 ).
The surface energy budget equation is further explained in part 4 of this section.
The net radiation flux at the surface (R n ) represents the actual radiant energy available at the surface. It is computed by subtracting all outgoing radiant fluxes from all incoming radiant fluxes ( Figure 2) as illustrated in the surface radiation balance equation: Where R S� is the incoming shortwave 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 ), and ε o is the surface thermal emissivity (dimensionless) (Waters et al. 2002).

Surface Albedo (α):
The albedo at the top of the atmosphere is compute as follows: Where ρ λ is the reflectivity and ω λ is a weighting coefficient for each band compute as follows:  Where ESUN is elevation of the sun. Albedo is defined as the ratio of the electromagnetic radiation reflected from the surface of the soil and the plant to the incident light emitted by the sun. Surface albedo is computed by correcting the α toa for atmospheric transmissivity: Values for α path_radiance range between 0.025 and 0.04 and for SEBAL we recommend a value of 0.03 based on Bastiaanssen (1998).
τ sw includes the transmissivity of both direct solar beam radiation and diffuse (scattered) radiation to the surface. We calculate τ sw assuming clear sky and relatively dry conditions using an elevation-based relationship from FAO-56: Where z is the elevation above sea level (m).

Incoming Shortwave Radiation (R S� ):
Incoming shortwave radiation is the direct and diffuse solar radiation flux that actually reaches the earth's surface (W/m 2 ). Its value is computed as follows: Where G sc is the solar constant (1367 W/m 2 ), cos θ is the cosine of the solar incidence angle as defined above, d r is the inverse squared relative earth-sun distance, and τ sw is the atmospheric transmissivity. The value R S� is computed for the days specified. Outgoing Long wave Radiation (R L� ): The outgoing long wave radiation is the thermal radiation flux emitted from the earth's surface to the atmosphere (W/m 2 ). It is computed in SEBAL through the following steps:

Computation of vegetation indices of Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), and Leaf Area Index (LAI)
The NDVI is the ratio of the differences in reflectivities for the near-infrared band (5) () and the red band (4) () to their sum: The NDVI is a sensitive indicator of the amount and condition of green vegetation. Values for NDVI range between −1 and +1. Green surfaces have a NDVI between 0 and 1 and water and cloud are usually less than zero.
The SAVI is an index that attempts to "subtract" the effects of background soil from NDVI so that impacts of soil wetness are reduced in the index. It is computed as: Where; L is a constant for SAVI. If L is zero, SAVI becomes equal to NDVI. A value of 0.5 frequently appears in the literature for L.
The LAI is the ratio of the total area of all leaves on a plant to the ground area represented by the plant. It is an indicator of biomass and canopy resistance. LAI is computed for southern Idaho using the following empirical equation: Where; SAVI ID is the SAVI calculated from Equation (9).

Computation of Surface emissivity (ε)
Surface emissivity (ε) is the ratio of the thermal energy radiated by the surface to the thermal energy radiated by a blackbody at the same temperature.

Computation of corrected thermal radiance (Rc)
The corrected thermal radiance (R c ) is the actual radiance emitted from the surface.
Brightness Temperature is obtained from the following relation: Where K 1 and K 2 are constants for Landsat images, L λ (L λ = M L Q cal + A L ) spectral radiance where M L is band-specific multiplicative rescaling factor from metadata, Q cal is quantized and calibrated standard product pixel value and A L is band-specific additive rescaling factor from metadata.in this paper from Brightness Temperature and wavelength of emitted radiance recorded by the sensor (thermal band) is used.

Computation of Outgoing Long wave Radiation (R L �)
This is computed using the Stefan-Boltzmann equation: Where ε ο is the "broad-band" 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).
Choosing the "Hot" and "Cold" Pixels: The "cold" pixel is selected as a wet, well-irrigated crop surface having full ground cover by vegetation. The surface temperature and near-surface air temperature are assumed similar at this pixel. The "hot" pixel is selected as a dry, bare agricultural field where ET is assumed zero.

Incoming Long wave Radiation (R L� ):
The incoming long wave radiation is the downward thermal radiation flux from the atmosphere (W/m 2 ). It is computed using the Stefan-Boltzmann equation: Net surface radiation (R n ) is calculated is computed using Equation (2).

Soil Heat Flux (G):
Soil heat flux is the rate of heat storage into the soil and vegetation due to conduction. Estimates of G/R n for agriculture surfaces is between 0.05-0.15.

Sensible Heat Flux (H):
Sensible heat flux is the rate of heat loss to the air by convection and conduction, due to a temperature difference. It is compute using the following equation for heat transport: Where ρ is air density (kg/m 3 ), c p is air specific heat (1004 J/kg/K), dT (K) is the temperature difference (T 1 − T 2 ) between two heights (z 1 and z 2 ), and r ah is the aerodynamic resistance to heat transport (s/m). z 1 is the height just above the zero plane displacement (d ≅ 0.67 × height of vegetation) for the surface or crop canopy and z 2 is some distance above the zero plane displacement, but below the height of the surface boundary layer. Based on experimental analysis, values of 0.1 meter for z 1 and 2.0 meters for z 2 are assigned. Temperature difference (dt) is given as dT = Tz 1 -Tz 2 . The air temperature at each pixel is unknown, along with explicit values for Tz 1 and Tz 2 . Therefore, only the difference dT is utilized. SEBAL computes dT for each pixel by assuming a linear relationship between dT and T s : dt = b + aT s , where b and a are the correlation coefficients and T S is the land surface temperature. a is obtained by subtracting the dT (dt hot pixel − dt cold pixel) and the LST (LST hot pixel − LST cold pixel). Using Envi software, first the hot and cold pixels are separated according to vegetation and temperature of the pixels and the dt are calculated based on the difference of two hot and cold pixels: a = (dt hot pixel − dt cold pixel) / (LST hot pixel − LST cold pixel). b is obtained by multiplying −a in LST hot pixel and dt hot pixel: b = (−a) × LST(hot) + dt(hot).

K1 K2
Landsat 8     By replacing the unknowns in the dt equation, the temperature difference of z 1 and z 2 is obtained. Latent Heat Flux (λET), Instantaneous ET (ET inst ), and Reference ET Fraction (ET r F) computation Latent heat flux is the rate of latent heat loss from the surface due to evapotranspiration. It can be computed for each pixel using the following Equation: λET = R n − G − H Where λET is an instantaneous value for the time of the satellite overpass (W/m 2 ).
An instantaneous value of ET in equivalent evaporation depth is computed as: Where ET inst is the instantaneous ET (mm/hr), 3600 is the time conversion from seconds to hours, and λ is the latent heat of vaporization or the heat absorbed when a kilogram of water evaporates (J/kg) is computed as: The Reference ET Fraction (ET r F) is defined as the ratio of the computed instantaneous ET (ET inst ) for each pixel to the reference ET (ET r ) computed from weather data: Daily values of ET (ET 24 ) are often more useful than instantaneous ET.
Where ET r-24 is the cumulative 24-hour ET r for the day of the image. This is calculated by adding the hourly ET r values over the day of the image (Waters et al. 2002).
Where: ET o is the reference evapotranspiration (mm/day), R n is the pure radiation entering the surface of the plant (MJ/m 2 /day), G is the soil heat flux (MJ/m 2 /day), T is the mean daily air temperature at 2 meters (°C), u 2 is the average daily wind speed at 2 m height (m/s), e s is the saturation vapor pressure (kPa), e a is the real vapor pressure (kPa), e s − e a is the lack of saturation vapor pressure (kPa), Δ is the slope of the vapor pressure curve (kPa/°C), γ is the constant coefficient.

MODIS evapotranspiration product
Evapotranspiration product is an 8-day composite dataset produced at 500-meter (m) pixel resolution. The algorithm used for the Evapotranspiration data product collection is based on the logic of the Penman-Monteith equation, which includes inputs of daily meteorological reanalysis data along with Moderate Resolution Imaging Spectroradiometer (MODIS) remotely sensed data products such as vegetation property dynamics, albedo, and land cover. The download data set is the MODIS evapotranspiration product for Aqua Satellite MYD16A2 and for Terra Satellite MOD16A2. The MOD16A2 and MYD16A2 layers provide the following products: i) The composited Evapotranspiration (ET), ii) Latent Heat Flux (LE), iii) Potential ET (PET), iv) Potential LE (PLE) along with a quality control layer.
The pixel values for the two Evapotranspiration layers (ET and PET) are the sum of all eight days within the composite period and the pixel values for the two Latent Heat layers (LE and PLE) are the average of all eight days within the composite period.

Model evaluation
To assess model performance a comparison is run between observations obtained from applying the three models and the SEBAL algorithm. These indicators include R (Eq. 18), RMSE (Eq. 19), MSE (Eq. 20) and MAD (Eq. 21), as follows: In all these error detection indexes is the modeled data and is the observational data and N is the data count.

Results
In total the Chaharmahal and Bakhtiari province is covered by circa 248,000 hectares of agricultural land, 79,854 hectares of these are located in the Shahrekord plain, of which, 58,553 hectares are under irrigation cultivation and 21,301 hectares are under rain fed cultivation. In total (in Shahrekord city), there are 29,917 hectares under cultivation of crops and 21,759 hectares dedicated to vegetable and garden farming.
According to the available statistics, the average productivity of agricultural water use of the province is 0.87 kg/m 3 while it should be increased to 1.51 kg/m 3 in the Iranian 5-years program and to 1.9 kg/m 3 in the Iranian 10-years planning. The water resources of this province decrease by 46 million m 3 in average, annually (Ministry of Agriculture 2016).
Due to global warming, water demand and related tensions regarding water supply increased. Being aware of the fact that the knowledge about the actual evapotranspiration is essential in water supply and management in this research, we assess the productivity and determine the appropriate pattern of proportional water resources under the present day climatic and hydrological conditions.

Lysimeter
Evapotranspiration obtained for the initial growth period is 1.2 mm/day. The evapotranspiration volume increases with the advance of the growth period and the highest evapotranspiration (4.09) is related to the active growth period of the plants. Although the highest temperature is recorded at the end of the growth period, the evapotranspiration volume of this stage is lower (Table 8).

FAO Penman-Monteith method
The FAO Penman-Monteith estimates, before applying the crop factor, the highest evapotranspiration volume (7.43 mm/day) in the harvest time (2017/07/25). Because at this time of the growth period, the plant is completely ripe and no irrigation takes place. The largest volume of the plant consists of seed and chaff.
Consequently, to estimate the most accurate evapotranspiration volume, the volumes obtained from the FAO Penman-Monteith model are corrected through crop coefficients as shown in Tables 7 and 8. After applying the crop coefficient, the maximum evapotranspiration volume is obtained for active plant growth time. The minimum evapotranspiration volume is recorded in the initial period of growth. Fig. 13 Evapotranspiration of initial stage of wheat growth.

MODIS evapotranspiration product
The results of the study of variations of wheat evapotranspiration on the sample plots through MODIS images show that the maximum evapotranspiration of wheat in the field occurred in the middle growth stage. The average evapotranspiration at this stage is 1.13 mm/day based on Aqua satellite and 1.09 mm/day on the Terra satellite (Table 8).
The final stage of wheat growth with 0.59 mm evapotranspiration per day for Aqua satellites and 0.58 mm evapotranspiration per day for Terra satellites had the highest evapotranspiration after the middle wheat growth stage. The lowest evapotranspiration occurs in the initial stages of growth (Figures 13-14).

SEBAL algorithm
The obtained SEBAL parameters are illustrated in Table 9. As observed in figures (4-5, 8-9 and 16-17), the NDVI, R n and ET are expressed in hot and cold pixels indicting one of the critical functions of this algorithm.

Data evaluation
The results of the error comparison and the correlation coefficient of the above-cited methods are tabulated in Table 10. We show that the SEBAL algorithm together with satellite imagery and the FAO Penman-Monteith method, after applying the vegetation  coefficient, are considered the best methods for situations where sufficient instruments are not provided in obtaining evapotranspiration through the installation of the Lysimeter. Although the correlation between the evapotranspiration measured through MODIS images and Lysimeter is high (0.98), the error between their different outputs is also quite high. The error of evapotranspiration measured by the Aqua Satellite in comparison to the Lysimeter data is less than for the Terra satellite products.

Conclusion
Agriculture development and food security are threatened by a decrease in rainfall, rising temperatures, droughts, a decrease in water table level and an increase in evaporation. Evapotranspiration is an effective parameter in providing water balance and food security because of its contribution in determining plant water need. Accurate estimations of the water needs and water supply for plants, especially for wheat, as a strategic product in Iran, are of major interest.
Evapotranspiration estimations by experimental methods and different algorithms as used in this paper is an important step forward especially in data scarce areas. However, a careful validation of the different methods should be carried out using observed data e.g. from Lysimeter. The evapotranspiration estimations obtained through the experimental models and the SEBAL algorithm revealed that the SEBAL algorithm have highest correlation and the least error with the observed Lysimeter data. The error values of MSE is 0.36, of is MAD is 0.59, and of the RMSE is 0.60 in terms of the SEBAL algorithm. Comparable good values we found for the FAO penman Monteith method with MSE equal to 1.29, an MAD of 0.91, and a RMSE of 1.13 compared to the lysimeter output.
In general, the results of this study indicate that applying remote sensing and satellite images, allows estimating evapotranspiration volume in areas with data deficits. The results reported in this study reveal particularly, a high correlation between the SEBAL algorithm and the lysimeter data. Consequently, we suggest using data derived with the SEBAL algorithm in areas with similar environmental conditions where no data are available and as input information for a further assessment of hydrological process dynamics. We show that the results of this study can be applied in studies of water resources management and appropriate irrigation management on farm level.