Evaluation of SEBS, METRIC-EEFlux, and QWaterModel Actual Evapotranspiration for a Mediterranean Cropping System in Southern Italy

Abstract

Remote sensing-based evapotranspiration (ET) models with various levels of sophistication have emerged recently with the possibilities of user-defined model calibrations. Their application for water resources management and climate studies from regional to global scale has been rapidly increasing, which makes it important to validate field scale ET in a complex crop assemblage before operational use. Based on in situ flux-tower measurements by the eddy-covariance (EC) system, this study tested three single-source energy balance models for estimating daily ET from fennel/maize/ryegrass-clover cropland rotations in a Mediterranean context in southern Italy. The sensitivity of three user-friendly ET models (SEBS, QWaterModel, and METRIC-EEFlux) with reference to the EC system over a center pivot irrigated cropland is discussed in detail. Results in terms of statistical indicators revealed that SEBS and METRIC-EEFlux showed reasonable agreements with measured ET (r² = 0.59_SEBS, RMSE = 0.71 mm day⁻¹; r² = 0.65_METRIC, RMSE = 1.13 mm day⁻¹) in terms of trends and magnitudes. At 30 m spatial resolution, both models were able to capture the in-field variations only during the maize development stage. The presence of spurious scan lines due to sensor defects in Landsat L7 ETM+ can contribute to the qualities of the METRIC-Efflux’s ET product. In our observation, the QWaterModel did not perform well and showed the weakest congruency (r² = 0.08_QWaterModel) with ground-based ET estimates. In a nutshell, the study evaluated these automated remote sensing-based ET estimations and suggested improvements in the context of a generic approach used in their underlying algorithm for robust ET retrievals in rotational cropland ecosystems.

1. Introduction

Evapotranspiration (ET) is a most important component flux in the annual water balance of semi-arid and arid regions [1]. Since the early eighties, ET estimation through modeling and remote sensing techniques has gained much attention in the context of enhancing water productivity in agriculture, also in the context of enhanced water scarcity and climate change [2]. Global climate projections postulate that not just the arid and semi-arid regions, but even irrigated agroecosystems in the Mediterranean regions are highly vulnerable to water scarcity [3]. For climate change preparedness, accurate water accounting from irrigated fields is important for sustainable water use and enhanced water productivity, especially in semi-arid regions where water availability is limited and irregular. ET, as a versatile biophysical indicator, can be considered as a proxy for crop water demand and helps in estimating the optimal water need for crops at the right time and also helps in the pre-determination of rotations in the crop field under stress conditions [4].

Being a complex eco-hydrological process that is governed by the biophysics operating in the soil–plant–atmosphere [5,6], ET estimation depends on various soil hydrological, plant-physiological, and meteorological variables that also bring uncertainties in the measurements and estimations [7,8,9]. It has been observed in various ET estimation approaches such as hydrological (e.g., soil water balance depletion method, weighing lysimeters) [10], micrometeorological (energy balance and Bowen ratio, aerodynamic method, eddy covariance) [11,12], and plant physiological (sap flow method, Chambers system or Penman–Monteith model) or empirical (methods based on crop coefficient) [13,14], soil water balance modeling) that various complexities operate and scale issues are inherently rampant [15,16,17]. There are also methods of field-based modeling approaches (i.e., ETo calculator [18], CROPWAT8.0, or MY-SIRR [19,20]). Mostly, these are point-based methods, tedious and data demanding, and spatially insufficient for large-scale needs. In contrast, remote sensing-based ET estimates are a cost-effective option that can estimate ET at larger scales with various degrees of granularities. As the satellite sensors are not capable of measuring near-surface vapor content, ET fluxes are generally determined by the energy balance approach on pixels in the spatial domain [21,22,23].

In the energy balance technique, models such as SEBS, TS-SEBS, SSEBop, METRIC, SEBAL, etc. use auxiliary surface parameters (albedo, emissivity, elevation, vegetation, and leaf area index, etc.) to predict turbulent fluxes on the Earth’s surface. These models differ in terms of computing sensible heat flux (H) and estimating fraction of reference ET (ETrF) [24,25,26,27,28,29,30,31,32]. Some models use land surface temperature (LST) and vegetation index (VI) relations to derive actual ET (ET_a) and crop coefficients (K_c) or either deploy the normalized difference vegetation index (NDVI) and basal crop coefficient (K_cb) to estimate crop water and transpiration [26]. This spatial ET estimation has great importance in hydrology, water management, and agriculture. Based on these concepts, real-time ET products have been developed at multiple temporal and spatial resolutions in recent years using cloud-based computing opportunities by adopting user-friendly Graphical User Interface (GUI). Some examples are WAPOR (developed by the Food and Agriculture Organization (FAO) to monitor water productivity), EEFlux (Google Earth Engine Evaporation flux developed by a consortium of the University of Nebraska-Lincoln, Desert Research Institute (Desert Research Institute) and University of Idaho), SSEBop (the operational simplified surface energy balance model, provided by USGS Famine Early Warning Systems Network), designed for agri-water-food assessment from local to continental scales [27,28].

In the literature, various single-source models (SEBS, METRIC and SEBAL) exist with low model complexity and simple parameterization scheme which can simulate turbulent heat fluxes. To improve the accuracy of turbulent heat fluxes, the more advanced dual-source energy balance models (TSEB and ALEXI) was developed with complex calibration and parameterization scheme which also require large number of input datasets which thereby compromise the computational efficiency of those models [29,30,31]. In order to increase the computational efficiency in dual-source energy balance models as well as simulate regional to global fluxes with higher accuracy, an advanced and a more robust dual-source energy balance model (TS-SEBS) has been recently developed which require fewer number of input datasets [32] to yield fluxes from regional to global scale. Although all of the above-mentioned models can simulate energy balance components with acceptable accuracies based on the manual calibration, complex parameterization, field-level and programming-based expertise. However, due to the unavailability of a user-friendly GUI, those models are difficult to handle by non-technical water managers. Nevertheless, in the recent era, the complex calibrations and parameterizations of energy balance model-based calculations are being substituted with ET applications consuming simple user-friendly GUI. These models require limited climate and surface information for flux calculations that can even detect ET dynamics in response to ground energy fluxes. Despite the landscape heterogeneity and water attenuations, these simple models can even predict crop water stress. This study was particularly designed to evaluate the predictive ability of open-source ET estimation tools developed in three different platforms, like METRIC-EEFlux on Google Earth Engine, SEBS model on PCRaster (programming software for environmental modeling), and QWaterModel as the QGIS plugin. These models are comparable due to their single-source energy balance scheme that uses lumped resistance and assume plant and soil as one surface layer, which detects the flux turbulence via calculating aerodynamic and surface resistance at the top of the canopy [33,34,35].

This study particularly tested the existing models’ robustness in capturing the field-level ET estimates for a rotational Mediterranean maize cropping system that was irrigated by central pivot system. It is important to mention that the study site had irrigation impairments in 2009, which were taken into consideration. In this context, the study aimed to: (1) assess the predictive ability of freely available single-source energy balance models in estimating actual ET in comparison to ground-based flux tower measurements and (2) the EC system’s ability to sense field-based variations in a rotational maize agroecosystem. For comparison purposes, three energy balance applications (METRIC-EEFLUX, SEBS, and QWaterModel) were tested to investigate how well these models performed in a rotational irrigated field (maize, fennel, ryegrass-clover), which was also an experimental site of a managed agroecosystem under ten years of close surveillance (2004–2014).

2. Site Features

The Borgo Cioffi Farm, (40.52° N, 14.95° E, 15 m a.s.l.) is located in the Sele Plain, Southern Italy at a 5 km distance from the Tyrrhenian Sea (Figure 1). This 15-hectare Buffalo Farm is an experimental field for an agroecosystem where maize and forage crops are grown. It is watered by a central pivot system while neighboring fields are on the mobile irrigated system. It has sandy-clay soil (clay: 52%; silt: 28%; sand: 20%), mostly alluvial due to the nearby Sele River. The higher clay content in the soil means a larger soil water-retaining capacity. Its field capacity and wilting points are 30.9 ± 3.8% and 17.2 ± 1.7%, respectively, with 2.5 ± 0.3% organic matter content.

The mean annual air temperature is 18 °C, with a mean maximum up to 30 °C in the summer, and the mean annual precipitation is 600 mm. The main cultivated species are maize and alfalfa, along with some winter grass crops (i.e., fennel, ryegrass, and clover). The same crop rotations were in place before the field was administered as a flux site. However, maize was only grown until 2010, later, replaced by alfalfa until 2015. The terrain is mostly flat with an approximately 2% gentle slope toward the south. The eddy covariance and the met station lie approximately in the center of the field. This tower monitors carbon, water, and energy fluxes between the agroecosystem and atmosphere. Its footprint is 182 m along the prevalent wind direction (NE–SW) [36]. Further details on the installed instruments are summarized in

Table 1. Specific characteristics of the instruments installed at the Borgo Cioffi farm in Italy.

Instrument	Field Measurements
Eppley Laboratory Inc., Newport, RI, USA	Incident and reflected global radiation (4 m)
Rebs Q7.1, Seattle, WA, USA	Net radiation (4 m)
LI-190SA, LI-COR Inc., Lincoln, NE, USA, Rotronic MP, Campbell Sci. Ltd., Shepshed, UK	Incident and reflected photo-synthetically active radiation, air temperature and relative humidity (4 m)
Sonic anemometer R3 (Gill Instruments Ltd., Lymington, UK)	Wind speed and direction (4 m)
Rain Collector II, Davis Instruments, Hayward, CA, USA	Precipitation and irrigation (ground level)
TCAV thermocouple probes	Soil temperature (0.3 m depth)
CS water content reflectometer (Campbell Scientific, Ltd., Shepshed, UK)	Volumetric soil water content (SWC)—(0.3 m depth)

.

3. Materials and Methods

3.1. Ground-Based Actual Evapotranspiration (ET_a) Estimates

Level 3 daily observations of the IT_BCi EC flux site were downloaded from the ICOS-Carbon portal, a network of flux sites measuring. The eddy covariance data were used to derive ET_a by following the Bowen ratio method, applied as [37,38];

$${\mathrm{ET}}_{\mathrm{a}\_\mathrm{Bowen}}=\frac{1}{\lambda {\rho }_w}\left[\frac{R_n-G}{1+\left(\gamma \frac{\Delta T}{\Delta e}\right)}\right]$$

(1)

where ‘R_n’ is the net radiation; ‘G’ is the ground heat flux; ‘γ’ is the psychometric constant; ‘∆T’ and ‘∆e’ are the temperature and humidity difference between surface and height; ‘ρ_w’ is the water density; and ‘λ’ is the latent heat of vaporization.

3.2. Ground-Based Reference ET (ET_{O_FAO56})

ET_O was computed from daily observations using the standard FAO56 Penman–Monteith model that refers to the water consumed by a well-watered grass reference crop; the model is fully described in [39].

$${\mathrm{ET}}_{\mathrm{O}\_\mathrm{FAO}56}=\frac{0.408\Delta \left(R_n-G\right)+\gamma \left[\frac{900}{T+273}\right]u_2\left(es-ea\right)}{\Delta +\gamma \left(1+0.34u_2\right)}$$

(2)

where ‘ET_O’ is the reference ET (mm day⁻¹); ‘R_n’ is the net radiation at the crop surface (MJ m⁻² day⁻¹); ‘G’ is the soil heat flux density (MJ m⁻² day⁻¹); ‘T’ is the mean daily air temperature (°C) at a height of 4 m; ‘u₂′ is the wind speed(ms⁻¹) at a height of 4 m; ‘es’ is the saturation vapor pressure (kPa); ‘ea’ is the actual vapor pressure (kPa), ‘Δ’ is the slope of the vapor pressure curve (kPa °C⁻¹); and ‘γ’ is the psychometric constant (kPa °C⁻¹). The soil heat flux is ignored (G = 0) in daily applications.

These methods were applied on EC data to determine weekly ET trends of ET_a and ET_o across a ten-years span (2004–2014), since some data on net radiation were missing for January–April 2004 and September 2011–April 2013. Therefore, a temporal subset of 2008–2009, related to the rotational pattern was used in study because later, the crop pattern was replaced with alfalfa only, with no maize rotations.

3.3. Remote Sensing Data and Retrievals

The study used 36 near cloud-free Landsat 5 Thematic Mapper and 7 Enhanced Thematic Mapper Plus (ETM+) images (path 189/row 32) retrieved through Google Earth Engine (GEE) for the selected time range (2008–2009) covering maize–fennel–ryegrass–clover crops (Table S1). Landsat was chosen due to its medium spatial resolution and reasonable seasonal coverage for historical analysis (2008–2009). Image pre-processing and variable retrievals were also done on the GEE Platform. Processing involved steps like cloud masking, at-sensor digital number to radiance conversion, and conversion to Top of Atmosphere (TOA) or at-sensor reflectance [40] by applying algorithms, such as:

• ee. Algorithms.Landsat. simpleCloudScore
• ee. Algorithms.Landsat. calibratedRadiance
• ee. Algorithms.Landsat.TOA

Furthermore, gap-filling of Landsat ETM+ images was done through the USGS L7 Phase-2 Gap filling protocol using the single kernel method [41,42]. Additionally, SEBS dependent variables (albedo, NDVI, emissivity, and surface temperature) were derived using the equations given in

Table 2. List of equations used for deriving input variables from Landsat data for the Surface Energy Balance System (SEBS) model processing.

Variables	Equations
Albedo	$\frac{\mathsf{\alpha }A+\mathsf{\beta }B+\mathsf{\gamma }C+\mathsf{\delta }D+\mathsf{\epsilon }E-\mathsf{\theta })}{\mathsf{\mu }}$
NDVI	$\frac{NIR-Red}{NIR+Red}$
Fractional Cover	$fc=\frac{{(NDVI-NDV{I}_{min})}^2}{NDV{I}_{max}-NDV{I}_{min}}$
Emissivity	$fc\times \left(a+b\right)$

where A is Blue, B is Red, C is NIR, and D as SWIR1 and E as SWIR2 are the top-of-the canopy reflectances of the Landsat bands, respectively, and the respective coefficients are α = 0.356; β = 0.13; γ = 0.37; δ = 0.085, ε = 0.072, θ = 0.018, and µ = 1.016 [43,44]. For emissivity, fractional vegetation cover was computed first from the minimum and maximum NDVI values and coefficients a = 0.004 and b = 0.986 [45].

.

Six environmental variables (NDVI, emissivity, and 30 m-SRTM DEM [46], LST, fractional cover, and albedo) were incorporated in the SEBS model. The same variables were used in the automated METRIC-EEFlux application, except for the National Land Cover Database (NLCD) because SEBS uses fractional vegetation cover derived from LAI for the surface roughness and aerodynamic stability estimations. The study did not modify the model scheme to test their applicability for the site in the present state.

After flux estimations, ET_a was computed from images by solving the energy balance equation for each pixel. ET_a average was extracted from daily ET_a maps (36 images) within the 182 m circular footprint. Later, the averaged modeled ET_a was compared with ground-based ET_a for assessment.

3.4. Crop Field Related Measurements

In field, daily soil water content (SWC) was measured at a depth of about 0.03 m via two soil moisture sensors (Theta Probe, Delta-T devices Ltd., Cambridge, UK) while gross primary production (GPP) was calculated with The Eddymeas and Eddysoft packages through the partitioning of net ecosystem exchange (NEE) into ecosystem respiration (Reco), and gross primary production (GPP) on the CarboEuropeIP database [47]. Leaf area index (LAI) was estimated using an empirical relation between LAI and the weighted difference of vegetation index (WDVI) derived from the SPOT image. Crop yield was determined as the product of GPP accumulated over the growing season of the respective crop and a crop-specific harvest index, and processing detail was also mentioned in [35]. In this study, agrometeorological variables were submitted for descriptive statistical analysis: range, maximum, average for maize and fennel crop cycle (i.e., GPP, above ground biomass (AGB), LAI, canopy height, yield (

Table 3. Sum of rain, irrigation and ET_a for the period 2008–2009.

Crop	ETa (mm)	Irrigation (mm)	Rain (mm)	Rainy Days	Duration
Maize 2008	305.41	34.21	68.59	18.00	1 May–22 August 2008
Fennel 2008–2009	179.61	0.00	924.67	88.00	15 September 2008–5 February 2009
Maize 2009	275.90	291.60	24.90	17.00	11 June–8 September 2009
All crops	1336.88	399.58	2966.35	309.00	1 January 2008–31 December 2009

and

Table 4. Yield and crop productivity of maize and fennel.

Crop	GPP (g cm−2)	AGB (m2m−2)	Max Yield (kg m−2)	WUE (kgha−1mm−1)	Canopy Height (m)	Max LAI (m2m−2)
Maize 2008	9.55	2.44	1171.30	5.20–30.00	0.31–2.95	4.98
Fennel 2008–2009	2.41	0.87	406.39	6.80–36.40	0.10–0.80	3.00
Maize 2009	8.20	1.92	923.50	3.60–34.00	0.26–2.64	4.80

)). Moreover, crop productivity data was used to calculate water use efficiency (WUE) as follows [48]:

$$\mathrm{WUE}=\frac{\mathrm{Yield}}{{\mathrm{ET}}_{\mathrm{a}}}$$

(3)

where ‘Yield’ means the yield of maize or fennel (kg·ha⁻¹) in each year and ‘ET_a’ is the total measured ET_a (mm) in the corresponding year throughout the respective crop cycle.

3.5. Energy Balance Models

Energy balance models require spatially distributed visible, near-infrared, and thermal infrared inputs from satellite imageries to use their spectral information for deriving surface fluxes. Some of these models are available in the public domain.

With simple parameterization of land cover micrometeorology, these models calculate energy balance through net radiation (R_n), soil heat flux (G), and sensible heat flux (H) at each remote sensing pixel and estimate the latent heat (λE) as the residual term as follows:

$$\mathsf{\lambda }\mathrm{E}={\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\mathrm{H}$$

(4)

3.5.1. Mapping Evapotranspiration at High Resolution with Internalized Calibration- (METRIC)-EEFlux

The Earth engine evapotranspiration flux (EEFlux) model operates on the Google Earth Engine platform to derive RS-based ET estimates using the METRIC algorithm [49]. The application deploys Landsat’s thermal and shortwave bands to derive actual ET (ET_a) from albedo, vegetation, and other surface parameters [50,51]. Images in METRIC-EEFlux were internally calibrated from the alfalfa reference ET (ET_o) using gridded weather data to retrieve the fraction of ET_o (ET_oF). This fraction was used to extrapolate instantaneous ET (ET_ins) according to the given equations.

The first term, ET_ins belongs to ET at each pixel at the time of satellite overpass, derived from latent heat flux as;

$${\mathrm{ET}}_{\mathrm{ins}}=3600\times \frac{\mathrm{LE}}{\lambda \times {\rho }_w}$$

(5)

where ‘ET_ins’ is instantaneous ET (mmh⁻¹), ‘λ’ is latent heat of vaporization (Jkg⁻¹), and ‘ρ_w’ is density of water (kgm⁻³).

The fraction of reference ET is computed as;

$${\mathrm{ET}}_{\mathrm{O}}\mathrm{F}=\frac{{\mathrm{ET}}_{\mathrm{inst}}}{{\mathrm{ET}}_{\mathrm{O}}}$$

(6)

In the end, the ET_a at each pixel is derived as;

$${\mathrm{ET}}_{\mathrm{a}}={\mathrm{ET}}_{\mathrm{O}}\mathrm{F}\times {\mathrm{ET}}_{\mathrm{O}}$$

(7)

The EEFlux approach is practically congruent with the traditional FAO method for actual ET (reference ET_O × K_c), where ET_OF is similar to crop coefficient (K_c). In this automated process, hot and cold pixels estimated from ET_O are used in modeling sensible heat flux (H). This automation has reduced the intermediate estimation errors and biases inherent to Rn, G, and subcomponents of H.

In this study, the Google EEFlux application (version 0.20.2) was used. The EEFlux portal also provides ready-to-use intermediate products (i.e., NDVI, LST, DEM, albedo, land cover). The temporal analysis of ET_a maps showed that some image pixels had null ET_a values at the field location due to a scan line issue in ETM+. Therefore, cloud masking and de-stripping of the ETM+ images were done in data processing for SEBS, model explained below. About six out of 36 images had the scan lines issue over the farm location, which were then improved and are shown in the final SEBS-ET product.

3.5.2. Surface Energy Balance System (SEBS)

The SEBS model [52] implemented in PCRaster Python Library [53] can be downloaded from GitHub repository. Previously, it was tested on Advanced Spaceborne Thermal Emission and Reflectance (ASTER) images for homogenous and heterogeneous land cover, the details of which are in [54]. In summary, the model contains three submodules: (1) for energy balance terms derivation; (2) roughness length derivation for heat transfer; and (3) for stability parameters [52]. Consequently, it gives the energy balance terms (LE, G, and H), relative evaporation, evaporative fraction (EF), and evapotranspiration. Rn and G computations are almost similar in most SEB models and they do not require hot and cold pixel selection, like in METRIC. Instead, it determines sensible heat flux (H) by solving the similarity relationships for the wind speed profile (u; ms⁻¹) and the difference between potential temperature (K) at the surface and a reference height. While in METRIC, the near-surface temperature difference (dT) is approximated from surface temperature (Ts) using a simple linear relationship dT = a + bTs (where a and b are correlation coefficients).

3.5.3. QWaterModel

QWaterModel is the QGIS plugin for calculating ET_a from LST maps developed on the energy balance model DATTUTDUT (deriving atmosphere turbulent transport useful to dummies using temperature). This model uses LST as a key parameter while assuming that ET is low as a result of hot pixels and ET is high for cold pixels [55]. In the LST map, the maximum temperature (T_max [K]) is taken as the hottest pixel, whereas the minimum temperature (T_min [K]) is taken as the 0.5% lowest temperature in the image to avoid extreme conditions (open water). Previously, the QWaterModel was applied in a tropical environment using high resolution drone-based thermal images for predicting ET_a in an Indonesian oil palm plantation against eddy covariance measurements. Landsat images were not tested in the study due to high cloud cover in their specific time span [56], however, they were used in the original DATTUTDUT model study for its comparison with other conventional energy balance models (SEBAL, SEBS) [55]. The original DATTUTDUT model estimates the radiation budget from a set of parameters and Sun–Earth geometrics, applicable to cloud-free images only. In this study, model parameters were configured according to in situ field measurements, and the application was tested for clear sky conditions to explore its usability and generality for Mediterranean crops. The performance was tested against other open-source applications (i.e., SEBS and METRIC). In the application, ET_a was derived from latent heat fluxes by calculating the latent heat of vaporization using in situ air temperature and the lowest pixel temperatures from the LST maps for the plugin estimates following the methodology explained in [56].

3.6. Validation of Model-Based ET_a

Modeled ET_a estimates were validated against EC-based measured ET_a. Validation was done only for one EC site located at the Borgo Cioffi Farm. Considering the height of the tower and the direction of prevailing winds, a circular footprint of 182 m was used. For statistical comparison, the determination of coefficients (r²), root mean square error (RMSE), mean bias error (MBE), and mean percentage error (MPE) were calculated to understand the agreement between the model-based estimates and EC measurements. The detailed mathematical formulations are as follows.

$$r^2=\frac{{{\displaystyle \sum }}_{i=1}^N{\left(E{T}_{estimated,i}-{\overline{ET}}_{measured}\right)}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(E{T}_{measured,i}-{\overline{ET}}_{measured}\right)}^2}$$

(8)

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\lfloor E{T}_{estimated,i}-E{T}_{measured,i}\rfloor }^2}$$

(9)

$$MBE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N\left(E{T}_{estimated,i}-E{T}_{measured,i}\right)$$

(10)

$$MPE=\frac{1}{N}{\displaystyle \sum }_{i=1}^N(\frac{{\lfloor E{T}_{estimated,i}-E{T}_{measured,i}\rfloor }^2}{E{T}_{measured,i}})\times 100$$

(11)

where ‘N’ is the number of data (days); ‘ET_measured,i’ is the EC-based ET_a each day at the ith position; and ‘ET_estimated,i’ is the ET_a of the corresponding day at ith computed by one of the models considered herein.

4. Results and Discussion

4.1. Averaged Ground Based Evapotranspiration (ET) Trend and Temporal Selection

The time series analysis as presented in Figure 2, demonstrated the weekly averaged ET calculated for ten years (2004–2014), later, used as climate normal. The applied method is explained in Section 3.1 and Section 3.2. According to the magnitude, ET_o was higher than ET_a every year except 2012. Unparalleled ET_o trends in the 2012 and 2013 series were due to the absence of measured net radiation (Rn) during January–April 2004 and September 2011–April 2013, which is a principal energy input in the ET estimation. To assess the accuracy of the energy balance models specifically, those unparalleled times series were excluded from the analysis, and only a subset of 2008–2009 was chosen wherein ET_o was temporally aligned with ET_a.

4.2. Irrigated Water Management and Crop Performance

Regarding weather conditions in the selected period, the field showed average maximum temperature during summer between 25–30 °C and a minimum temperature between 10–15 °C.

Average relative humidity (RH) was comparatively higher in winter, ranging between 60 and 80% (Figure 3a). Rainfall frequency was also higher in winter with a maximum weekly sum equivalent to 140 mm, the field experienced higher rainfall at the end of 2008 with a significant precipitation decline in the next summers, whereas soil water was also reduced during the no-crop season and subsequent ryegrass-clover coverage in winter 2009.

The graph also explained crop management practice, when soil water content dropped to 5% in 2009, it was irrigated for maize crops during the dry period, as being irrigated in 2008 too. It was observed that soil water content significantly dropped in 2009 compared to 2008, possibly due to fennel cultivation, which consumed sufficient soil water. Through 2008 irrigation routine, it was observed that four rounds of irrigation were not enough for soil, as obvious from high seasonal ET_a, this moisture decline get stabilized after 11 rounds of irrigation in 2009, at a maximum of 60 mm and a minimum 20 mm (Figure 3b). This field behavior was further examined through measured and estimated ET_a comparisons during the crop cycles.

Field behavior was further examined by descriptive statistics on field productivity and water consumption (Table 3 and Table 4) including all crops and fennel–maize rotation, as illustrated in Section 3.4. Overall, plant density was lower in 2009 than in 2008 due to a fault in the irrigation system that occurred at the early vegetative stages of maize in 2009, by which seed establishment and seedling survival were impaired; this aspect justifies the low biomass productivity in 2009, equivalent to 1.92 m², AGB (Figure 4b,d), whereas the total water influx [precipitation and irrigation] was lower in 2008, compared to ET_a, indicating water stress conditions in 2008. The higher ET_a in 2008 means that significant water was consumed in maize senescence period and evaporated during the dry period; in contrast, ET_a was relatively lower in 2009 during maize growth, with slightly lower gross primary productivity (GPP = 8.2) and its estimated WUE was higher than last year, though with lower canopy height. In comparison to maize, the fennel showed greater WUE under Mediterranean conditions. Nevertheless, due to lower yield and biomass by fennel (Figure 4c), maize is still important for the field productivity (Figure 4b,d). Whether this crop arrangement was optimum or not for the site, it was difficult to answer with two years of agronomic data. However, it can be explored by long-term assessment of WUE and water productivity of this experimental site. This study, however, focused more on the energy balance models’ capacity to address farm-level complexities at a 30 meters’ resolution.

To retrieve model estimations about ET_a, Landsat 5 and 7 images were processed for 2008–2009. Their temporal correspondence with ET_a is shown in Figure 5, along with a comparison of the study year’s ET_a with the ET_a climate normal (2004–2014). The analysis showed high evaporative demand during summers, which was also the maize growing season in both years. This temporal behavior is in agreement with the ET_a climate normal in the field, indicating that ET_a is more sensitive to seasonal variations than crop type, because of the fact that this field grew maize up to 2010. Later, it was replaced with alfalfa up to 2014, but no significant signature regarding change in crop rotation was found in the average trend of climate normal.

This aspect also needs further research as to which crop, in particular, can deviate seasonal ET_a from the ET_a climate normal conditions, which can be analyzed with different crop ET_a comparisons. The present trend, however, justifies the use of general energy balance models that estimate ET, regardless of crop type.

4.3. Evaluation of Model-Based Evapotranspiration

To confidently use ET tools in hydro-meteorological studies, it is important to assess their independent accuracies. Therefore, we analyzed, assessed, and applied three models (METRIC-EEFlux, SEBS, and QWaterModel) against EC-based ET_a measurements. The results of our scatterplots and temporal analysis are shown in Figure 6 and Figure 7. According to the results, both SEBS and METRIC-EEFlux captured seasonal ET_a estimations in correspondence with flux tower ET_a, except for the QWaterModel. This model overestimated the ET_a magnitude during the maize period in 2009. Among all three models, METRIC-EEFlux performed better as indicated by the error matrix (r² = 0.65, RMSE = 1.13 mm day⁻¹, MPE = 10.7%, MBE = 0.34) in Figure 7. This showed correspondence with measured ET_a trends throughout crop cycles and was deemed reasonable in the case of short crops (fennel, rye-grass, clover). The performance of the SEBS was better than METRIC-EEFlux in terms of RMSE (r² = 0.58, RMSE = 0.71 mm day⁻¹, MBE = −0.05), however, was similar in MPE = 10.3% with the lowest MBE (−0.05). Though SEBS and METRIC-EEFlux had a comparable statistical evaluation, the results showed that METRIC-EEFlux overestimated ET_a in summer 2008, whereas SEBS slightly underestimated, as illustrated in Figure 6. The 1:1 line in Figure 7 also indicated that METRIC-EEFlux overestimated the daily ET_a of rotational crops for ET_a values higher than 1.13 mm day⁻¹, however, this was lower by SEBS in comparison (RMSE = 0.71 mm day⁻¹).

In contrast, the QWaterModel showed significant dissimilarity with the measured ET_a (r² = 0.09 with a difference of RMSE = 2.65 mm day⁻¹ and highest MPE = 14%, indicating the least fitting model for crop arrangements under study. While same meteorological parameters were used for simulations in both SEBS and QWaterModel, while METRIC-EEFlux used the automated input of gridded weather data. QWaterModel showed tremendous difference with 0.4–0.8 mm day⁻¹ values in June, lower than the expected range (3–4 mm day⁻¹), whereas neither it could give satisfactory results during the maize growing season. The reason could be its sensitivity to the in situ measurements of net radiation and shortwave irradiance, as mentioned in [55,56] and also observed during analysis. This application is perhaps more effective to thermal camera images, as tested previously on oil palm plantations since no good results were obtained for the Landsat images in this study. QWaterModel developers also suggested excluding other artifacts from the thermal image as a necessary step and manually defining minimal and maximal temperatures for further estimations, since the model solely depends on the LST image. It is also recommended to first calibrate satellite images from ET by planes, drones, and handheld or fixed thermal cameras (e.g., using QWaterModel) to be used for temporal and spatial analysis [56]. This research also confirms that the use of Landsat images in the QWaterModel with ET calibration is essential.

Overall, the lumped results suggest that METRIC-EEFLUX and SEBS-ET products can be equally used considering their better performance and correspondence in trend with the ground measurements. However, it would be an advantage if SEBS could also be available on a GEE platform like METRIC-EEFlux.

Figure 8 showed only those images where scanlines were observed over the flux site with the METRIC-EEFlux and QWaterModel based ET product. It illustrates that the use of intermediary products of METRIC-EEFlux also has no advantage if the scan line issue persists. Thus, gap-filling of METRIC-EEFlux end products is equally important to get a better ET_a output over the rotational crop field to avoid different spatial ET_a ranges by each model. The map also highlights the indifferent response by the QWaterModel-based ET_a in the L7 SLC-off based product.

At a 30 m resolution, SEBS and METRIC-EEFlux could show the in-field ET_a variations during the maize development period only, as obvious in the cases of May 2008 and July 2008, while QWaterModel could not retrieve them. However, variations were not visible in September 2008 and November 2008 during the fennel growth period, where field ET_a was very low, ranging between 0.8–1.8 mm day⁻¹. This indicates that the spatial ET estimations for short crops could be particularly challenging for these applications at a 30 m resolution. It also implies the use of high resolution ET_a for short crops water assessment in agrometeorological studies. Moreover, current applications showed reasonable agreement for tall crops like maize; the same was observed in the maize crop study conducted in Brazil using the METRIC-EEFlux tool [28] and tested on MODIS images for the Nile Delta River in Egypt [57]. Since it was the first ET based study conducted over the Borgo Cioffi experimental site, which also contained short crops, further work on crop WUE and productivity is needed to better understand irrigated field performance in water productivity and for the selection of the optimum crop pattern considering water stress conditions in the summer.

Overall, the newly built Google Earth Engine and other open-source applications made ET_a estimations convenient, though poor product quality was observed in the generalized energy balance schemes. Such minor differences in ET_a magnitude may cast a significant impact on large-scale crop water need assessments as realized in [58]. Improvements in online energy balance applications can make ET_a assessment more precise and accessible globally, however, cross-validations with in situ data at the farm scale are still imperative for better crop management.

5. Conclusions

In this study, we assessed the performance of three single-source energy balance models available on public domains that require minimum calibrations (i.e., SEBS, QWaterModel, and METRIC-EEFlux). They were tested on an irrigated maize rotational field, which was a registered EC site on ICOS and FLUXNET platform (the global network of micrometeorological tower sites). Results revealed that the installed micro-meteorological and soil sensors adequately sensed the field variations that occurred during crop cycles (i.e., irregularities in irrigation and resultant low biomass production), though water use efficiency varied significantly per crop, like fennel showed higher efficiency under rain-fed conditions than maize. Moreover, easy-to-use models can also estimate in-field ET_a variabilities in the case of tall crops like maize, which was not visible for short crops (i.e., fennel or rye-grass in ET_a maps).

Overall, SEBS and METRIC-EEFlux showed almost the same ET_a trend with varying magnitudes, except for the QWaterModel. It showed the lowest fitness r² = 0.09 and highest ET_a difference with RMSE = 2.62 mm day⁻¹, suggesting that this application needs the manual calibration of Landsat images before processing through the model. Though METRIC-EEFlux provided reasonable ET_a estimates with medium uncertainty (RMSE = 1.2 mm day⁻¹), the execution code seriously needs an L7 gap-filling protocol because its absence affects product quality over the site. On the other hand, as SEBS developed on PCRaster gave better results with the lowest RMSE = 0.71 mm day⁻¹, it would be an advantage if it gets a graphical user interface like METRIC-EEFlux. The study endorses the idea of automated model calibrations and user-friendly ET_a estimation tools, however, their product validation is also important for various agricultural sites to gain confidence for precise crop water estimations at the global or the regional scale.