Evapotranspiration Estimates at High Spatial and Temporal Resolutions from an Energy–Water Balance Model and Satellite Data in the Capitanata Irrigation Consortium

: The feasibility of combining remotely sensed land surface temperature data (LST) and an energy–water balance model for improving evapotranspiration estimates over time distributed in space in the Capitanata irrigation consortium is analysed. The energy–water balance FEST-EWB model (ﬂash ﬂood event-based spatially distributed rainfall–runo ﬀ transformation—energy–water balance model) computes continuously in time and is distributed in space soil moisture (SM) and evapotranspiration (ET) ﬂuxes solving for a land surface temperature that closes the energy–water balance equations. The comparison between modelled and observed LST was used to calibrate the model soil parametres with a newly developed pixel to pixel calibration procedure. The e ﬀ ects of the calibration procedure were analysed against ground measures of soil moisture and evapotranspiration. The FEST-EWB model was run at 30 m of spatial resolution for the period between 2013 and 2018. Absolute errors of 2.5 ◦ C were obtained for LST estimates against satellite data; while RMSE around 0.06 and 40 Wm − 2 are found for ground measured soil moisture and latent heat ﬂux, respectively.


Introduction
Distributed pixelwise knowledge of evapotranspiration (ET) and soil moisture (SM) is of extreme importance for a wide range of applications, especially for estimates of irrigation water needs in agricultural areas. However, their direct measurement at large scale is still not feasible; there exist several estimates from modelling and remote sensing data that are affected by intrinsic errors.
In the last 40 years, several satellites data from passive or active microwave sensors are available for the retrieval of superficial soil moisture (few centimetres) [1][2][3], which is nowadays used for a wide range of applications, even in agriculture water management [4]. However, problems can arise if the satellite images need to be used in conjunction with hydrological models for operational water management applications. In fact, SM estimates from passive sensors have a very low spatial resolution in between 25 and 50 km, which is too coarse for agricultural applications. Soil moisture pixel disaggregation [5,6] based on thermal infrared data is a way to improve the spatial resolution. The second concern is related to the retrieval of the soil water content of a few centimetres depending on the measurements' wavelength, which is not congruent with the hydrological active soil for flood generation or plant root zone uptake. Moreover, SM estimates from active radars can be affected by problems linked to the saturation of the retrieval algorithms [7] and their ability to detect soil moisture over vegetated surfaces [8].
Satellite land surface temperature (LST) data have been widely used as the main input variable of energy balance models to compute evapotranspiration as the residual term of the energy budget [9][10][11][12][13]. These models have limitations related usually to no consideration of the water mass balance component, to the estimates of ET at satellite overpass time, which may also be influenced by the presence of clouds preventing its computation, and daily ET estimates based on net radiation weights. Another main issue is related to the satellite LST spatial resolution, mainly in between 1000 m for MODIS and 100 m for LANDSAT data, which is lower than Visible Near-Infrared bands. Hence, disaggregation techniques are needed to recover high spatial resolution information, which are able to capture the crop fields dimension [14,15].
Instead, hydrological models based on water and energy balances can provide estimates continuously in time of evapotranspiration and also of soil moisture distributed in space, using the satellite land surface temperature as a control variable [16]. These models can then overcome the limitations related to cloud coverage, typical of visible and thermal infrared satellite images, and to the instantaneous estimates at the time of the satellite overpass. Moreover, the distribution of irrigation volumes into the hydrological model is one of the most uncertain input variables [17][18][19]. In fact, most of the time irrigation over large areas, certainly in Italy where the paper case study is located, is an unknown quantity and, if available, only an integral value is known. In the last years, from the integration of satellite soil moisture data and hydrological modelling, few estimates of distributed irrigation volumes maps are available but still with some uncertainties [4,20,21]. The need for several soil hydraulic parametres, which play a relevant role in water and energy fluxes estimates, can pose another limitation to the application of these models. These parametres are usually available at low spatial resolution for large scales and defined using soil texture maps [22], but problems of representativeness arise due to the pixels' heterogeneity. The calibration of these parametres may be achieved using satellite data, which allows a pixelwise tuning of these parametres, overcoming the traditional calibration based on a single multiplicative value retrieved from the comparison between observed and simulated ground discharge [23,24] or local ground measurements of soil moisture or evapotranspiration. This is feasible thanks to the distributed models' structure, which allows a pixel by pixel control of water and energy fluxes [25] (e.g., soil moisture (SM), land surface temperature (LST) and evapotranspiration fluxes (ET)). Even though little effort has been made in this direction, some examples are available [26][27][28]. In particular, references [29][30][31] used remotely sensed LST for calibrating the parametres of the FEST-EWB energy-water balance model showing its potentiality with respect to the traditional calibration with ground river discharge data.
In this study, we further investigate the integration of satellite data and hydrological modelling trying to address the following research questions: (1) To what extent does the calibration based on LST allow a better representation of the hydrological fluxes (e.g., SM, ET) in terms of spatial and temporal resolutions? (2) How much do these improvements depend on the ability of the model to reproduce the satellite LST? (3) What is the effect of a precise irrigation distribution on the LST and ET estimates in highly cultivated and irrigated areas?
Hence, the main paper objective is to improve the accuracy of pixelwise evapotranspiration and soil moisture estimates using remotely sensed LST to calibrate soil and vegetation parametres of an energy-water balance model. Moreover, the effect of the irrigation distribution will also bex evaluated on the pixelwise correctness of land surface temperature estimates. The procedure is then validated against ground measurements of soil moisture and evapotranspiration. The distributed energy-water balance model, flash-flood event-based spatially-distributed rainfall-runoff transformation-energy-water balance model [16,32], will be implemented over the Capitanata Irrigation Consortium (Italy) for the calibration period, from 2013 to 2016, and for the validation period, from 2017 to 2018.
With respect to previous studies [16,[29][30][31], three important innovations are introduced: (1) irrigation volumes as input to the FEST-EWB model are distributed over the consortium area considering only the irrigated crops areas retrieved from satellite data, (2) validation of FEST-EWB results against several eddy covariance data in different crops fields (e.g., tomatoes, cabbage, asparagus, wheat), (3) a different study basin highly irrigated and cultivated with very dry summer seasons, which helps to prove the ability of the FEST-EWB model to be able to correctly simulate ET and SM fluxes over any region of the World, (4) accuracy of satellite downscaled LST for model calibration.

Case Study
The case study is the Capitanata Irrigation Consortium, specifically the Sud Fortore district, located in Southern Italy in the Puglia region ( Figure 1) and delimited by the Apennines on the west and by the Gargano Promontory on the east side. It covers an area of about 65,000 hectares. Forty-five per cent of the total area is irrigated through the Consortium water distribution network (56,700 ha), while the remaining areas are irrigated with private wells [33]. The role of irrigation is crucial, in fact the mean irrigation volume in the irrigation season from April to October is of about 600 mm, while the seasonal rainfall amount is about 150 mm. Daily irrigation volumes measured in the main aqueduct are available from 2013 to 2018 with values between 2.1 and 6 10ˆ7 m 3 with a mean value of 4.6 10ˆ7 m 3 .  [16,[29][30][31], three important innovations are introduced: (1) irrigation volumes as input to the FEST-EWB model are distributed over the consortium area considering only the irrigated crops areas retrieved from satellite data, (2) validation of FEST-EWB results against several eddy covariance data in different crops fields (e.g., tomatoes, cabbage, asparagus, wheat), (3) a different study basin highly irrigated and cultivated with very dry summer seasons, which helps to prove the ability of the FEST-EWB model to be able to correctly simulate ET and SM fluxes over any region of the World, (4) accuracy of satellite downscaled LST for model calibration.

Case Study
The case study is the Capitanata Irrigation Consortium, specifically the Sud Fortore district, located in Southern Italy in the Puglia region ( Figure 1) and delimited by the Apennines on the west and by the Gargano Promontory on the east side. It covers an area of about 65,000 hectares. Fortyfive per cent of the total area is irrigated through the Consortium water distribution network (56,700 ha), while the remaining areas are irrigated with private wells [33]. The role of irrigation is crucial, in fact the mean irrigation volume in the irrigation season from April to October is of about 600 mm, while the seasonal rainfall amount is about 150 mm. Daily irrigation volumes measured in the main aqueduct are available from 2013 to 2018 with values between 2.1 and 6 10^7 m 3 with a mean value of 4.6 10^7 m 3 .
The Sud Fortore district is an intensive cultivation area, mainly devoted to durum wheat and tomatoes during the spring-summer season and fresh vegetables (sown in late summer and harvested October-February).  The Sud Fortore district is an intensive cultivation area, mainly devoted to durum wheat and tomatoes during the spring-summer season and fresh vegetables (sown in late summer and harvested October-February).

Meteorological Data
Meteorological data are available from 21 stations from ARPA-Regione Puglia and from the association Meteonetwork ( Figure 1). Then three other stations managed by Politecnico di Milano and farm Agricola Guzzetti s.r.l. are also available. All the measured variables at the stations (precipitation, air temperature and relative humidity, incoming solar radiation and wind speed) are used as input to the model at 1-hour temporal resolution from 2013 to 2018.

Soil Data
Soil type information are available from Regione Puglia as a shape file, allowing the identification of only two main soil classes for all the analysed area of 55,000 ha: silty clay and gravel sandy soils. This coarse spatial resolution of the soil type data is then reflected in a coarse spatial resolution of the soil hydraulic parametres used as input of the FEST-EWB model. In fact, soil physical and hydraulic properties can have an extremely high variability over space (even in an area of 1 m 2 ) and also time [34] and their correct quantification at high spatial resolution is essential for a correct estimate of the water balance components.
In particular, the soil is characterised by: soil hydraulic parametres, as the saturated hydraulic conductivity (ksat), the field capacity (fc), wilting point (wp), residual (θr) and saturated (θs) soil water content, Brooks-Corey index (BC), bubbling pressure (bp). These properties are defined for the two identified classes according to [22], assigning to each pixel of the analysed area its respective value. Another important parametre to define in the hydrological model is soil depth, which has been modelled as a single layer with an initial mean area value of 0.5 m, considering the predominant growing zone of fresh vegetable roots, relevant for the evapotranspiration process.
In Table 1, the mean and standard deviation (std dev) values of the main parametres are reported. The DEM shows a very flat area in the eastern part of the basin and a mountainous area in the west part ( Figure 1).

Eddy Covariance Data
Two eddy covariance stations are installed in the Sud Fortore area. Each year from 2014, the stations are moved in neighbourhood fields to follow different crops cycles [35]: These stations allow to measure primarily the latent (e.g., evapotranspiration) and sensible heat fluxes, but usually also the net radiation, ground heat flux, soil temperature and moisture.
In particular, soil moisture measurements are performed by TDR (Time Domain Reflectometre) soil moisture sensors (CS616, Campbell Sci, UT, USA) that were installed close to the eddy station at depths of 15 cm, considering the predominant growing zone of crops roots.
The turbulent fluxes (e.g., latent and sensible heat) are usually retrieved at high frequency and need to be post-processed applying a series of corrections well assessed in the scientific community [36]. The PEC software (Polimi Eddy Covariance) [37] is used in this study to post-process the turbulent raw data applying all the instrumental and physical corrections. Moreover, latent and sensible heat fluxes usually underestimate the available energy, leading to a non-closure of the energy balance. Hence, these fluxes are modified according to the procedure developed by [38], based on the preservation of the Bowen-ratio method to reach the energy budget closure.
Moreover, due to the intrinsic measurement technique of the eddy covariance station, the footprint (e.g., the source areas) of the measured vertical H and LE fluxes have been computed, following the two-dimensional footprint model of [39]: where f is the footprint dimension, L is the Obukhov length, D and P are similarity constants for unstable, neutral and stable atmospheric conditions, x is the footprint in the upwind direction and y in the orthogonal one, k is von Karman's constant, In particular, soil moisture measurements are performed by TDR (Time Domain Reflectometre) soil moisture sensors (CS616, Campbell Sci, UT, USA) that were installed close to the eddy station at depths of 15 cm, considering the predominant growing zone of crops roots.
The turbulent fluxes (e.g., latent and sensible heat) are usually retrieved at high frequency and need to be post-processed applying a series of corrections well assessed in the scientific community [36]. The PEC software (Polimi Eddy Covariance) [37] is used in this study to post-process the turbulent raw data applying all the instrumental and physical corrections. Moreover, latent and sensible heat fluxes usually underestimate the available energy, leading to a non-closure of the energy balance. Hence, these fluxes are modified according to the procedure developed by [38], based on the preservation of the Bowen-ratio method to reach the energy budget closure.
Moreover, due to the intrinsic measurement technique of the eddy covariance station, the footprint (e.g., the source areas) of the measured vertical H and LE fluxes have been computed, following the two-dimensional footprint model of [39]: where f is the footprint dimension, L is the Obukhov length, D and P are similarity constants for unstable, neutral and stable atmospheric conditions, x is the footprint in the upwind direction and y in the orthogonal one, k is von Karman's constant, y is the standard deviation of the cross-wind direction, zm is the height of measurement and zu is a length scale. Hence, FEST-EWB estimates of latent and sensible heat fluxes are integrated as a weighted sum over the footprint area in order to be comparable to eddy covariance measures [39]: where i is the position of a pixel in an image, F is the flux value in each pixel and ̅ is the average flux in the footprint area.

Methodology
The developed methodology improves FEST-EWB estimates of soil moisture and evapotranspiration and their spatial representativeness. The different followed steps are: (1) retrieval of high resolution satellite data of vegetation parametres and LST, (2) calibration of soil and vegetation parametres in each pixel of the area independently one from the other by comparing satellite observed and simulated land surface temperature from 2013 to 2016, (3) model validation at regional scale against LST from 2017 to 2018, and also at the local scale by comparison with measured fluxes from eddy covariance stations (e.g., soil moisture, evapotranspiration).

Land Surface Temperature
LST is computed with two different techniques: Single Channel (SC) and Split Window (SW) algorithms, which have been applied to Landsat-7 Enhanced Thematic Mapper Plus (ETM+) and Landsat-8 Thermal InfraRed Sensor (TIRS) respectively.
The SC algorithm structure has been extracted from [40]. LST for ETM+ is estimated using band 6, as: where Tsen is the at-sensor brightness temperature, ε is the Land Surface Emissivity (LSE), bγ = c2/λ is equal to 1277 K , and 1 , 2 and 3 are the so-called atmospheric functions, which are directly related to the upwelling and downwelling atmospheric radiances and the atmospheric is the standard deviation of the cross-wind direction, z m is the height of measurement and z u is a length scale. Hence, FEST-EWB estimates of latent and sensible heat fluxes are integrated as a weighted sum over the footprint area in order to be comparable to eddy covariance measures [39]: where i is the position of a pixel in an image, F is the flux value in each pixel and F is the average flux in the footprint area.

Methodology
The developed methodology improves FEST-EWB estimates of soil moisture and evapotranspiration and their spatial representativeness. The different followed steps are: (1) retrieval of high resolution satellite data of vegetation parametres and LST, (2) calibration of soil and vegetation parametres in each pixel of the area independently one from the other by comparing satellite observed and simulated land surface temperature from 2013 to 2016, (3) model validation at regional scale against LST from 2017 to 2018, and also at the local scale by comparison with measured fluxes from eddy covariance stations (e.g., soil moisture, evapotranspiration).

Land Surface Temperature
LST is computed with two different techniques: Single Channel (SC) and Split Window (SW) algorithms, which have been applied to Landsat-7 Enhanced Thematic Mapper Plus (ETM+) and Landsat-8 Thermal InfraRed Sensor (TIRS) respectively.
The SC algorithm structure has been extracted from [40]. LST for ETM+ is estimated using band 6, as: Remote Sens. 2020, 12, 4083 where T sen is the at-sensor brightness temperature, ε is the Land Surface Emissivity (LSE), b γ = c 2 /λ is equal to 1277 K, and Ψ 1 , Ψ 2 and Ψ 3 are the so-called atmospheric functions, which are directly related to the upwelling and downwelling atmospheric radiances and the atmospheric transmissivity. Fitting these parametres against the total atmospheric water vapour content (w), the SC water vapour approximation is computed as a direct relation between Ψ parametres against w: where matrix coefficients have been retrieved through direct relation between Ψ parametres against w. See more information in [15,41].
LST is retrieved from TIRS data following the algorithm developed by [41], which is based on the Split-Window (SW) technique: where T i and T j are the at-sensor brightness temperatures at the SW bands i and j (in K), ε is the mean LSE, ε = 0.5 (ε i + ε j ), ∆ε is the LSE difference, (∆ε = ε i − ε j ), w is the total atmospheric water vapor content (in g·cm −2 ) and a 0 to a 6 are the SW coefficients to be determined from simulated data.
For both algorithms, atmospheric water vapour is estimated by the computation of MOD07 product in MODTRAN radiative transfer code [44]. MOD07 is an atmospheric profile product measured by the MODIS sensor and available at Landsat pass time with a maximal time difference of 1:30 h.
The land surface emissivity (LSE) is computed using the NDVI Thresholds Method (NDVI-THM), originally developed by [45] and revised by [46] to estimate the LSE by the use of Fractional Vegetation Cover (FVC). The soil and vegetation emissivity values have been retrieved from [47,48].
A summary of formulas based on the NDVI-THM method for each Landsat thermal band is shown in Table 2. Table 2. Emissivity values for Landsat thermal bands using the Normalised Difference Vegetation Index Thresholds Method (NDVI-THM) method.

Landsat-8 B10
Landsat-8 B11 Landsat-7 ETM+ Downscaling Based on Nearest Neighbor Temperature Sharpening NNTS Methodology Land surface temperature data are usually available at lower spatial resolution than data in the visible or near-infrared bands. Especially for agricultural applications, downscaling procedures need to be implemented to increase LST spatial resolution; hence, coarse thermal data are disaggregated to a finer resolution data providing sub-pixel temperatures, using sharpening or unmixing algorithms.
In this paper, LST has been downscaled with the Nearest Neighbor Temperature Sharpening (NNTS) proposed in [15]. The method is based on the local variant model proposed by [49] called Temperature Sharpening (TsHARP), which takes into account similar pixel properties and its distance, both over a sliding N×N window over an LST image using relationships of Normalised Difference Water Index (NDWI) and Normalised Difference Vegetation Index (NDVI) indices versus the LST.
Scheme of downscaling procedure is shown in Figure 2. The method is based on the assignation of two weights: one, related to the direct relationship between LST and NDVI or NDWI (wfit); and the other one related to proximity of similar NDVI or NDWI pixels (wdist). When poor relation is found in the sliding window, the neighborhood gains more importance in the retrieval of LST thin pixels and when strong relation is obtained, the regression acquires more importance (see Figure 2). Once weights have been computed, temperatures associated to each weight must be estimated as: where Tdist is the LST value retrieved in the function of pixel inverse distance (di) of sliding window (Equation (6)) and Tfit is the LST retrieved with the relation between the best spectra index fit (NDVI or NDWI highest correlation) versus the LST. If the correlation is poor (r 2 < 0.6) Tfit is simply computed as the average value of LST coarse pixels included in the N×N window (Equation (7)). The averages performed on Equation (6) and 1 in Equation (7) have been computed only for the coarse NDVI/NDWI pixels with values ±20% of thin pixel value to be estimated (similar spectra index pixels). Finally, thin pixel LST is obtained as:

Vegetation Parametres
FVC, as the fraction of the total pixel that is covered by canopy, is estimated from NDVI, according to [50]: where NDVIs and NDVIv are representative values for bare soil and green vegetation, respectively. These values were estimated as 0.15 and 0.9, respectively. Values of bare soil (FVC = 0) and full vegetation cover (FVC = 1) have been obtained from the NDVI histograms retrieved during the period of maximal vegetation grow. As crop type does not vary a lot year over year, these soil and vegetation The method is based on the assignation of two weights: one, related to the direct relationship between LST and NDVI or NDWI (w fit ); and the other one related to proximity of similar NDVI or NDWI pixels (w dist ). When poor relation is found in the sliding window, the neighborhood gains more importance in the retrieval of LST thin pixels and when strong relation is obtained, the regression acquires more importance (see Figure 2). Once weights have been computed, temperatures associated to each weight must be estimated as: where T dist is the LST value retrieved in the function of pixel inverse distance (d i ) of sliding window (Equation (6)) and T fit is the LST retrieved with the relation between the best spectra index fit (NDVI or NDWI highest correlation) versus the LST. If the correlation is poor (r 2 < 0.6) T fit is simply computed as the average value of LST coarse pixels included in the N×N window (Equation (7)). The averages performed on Equation (6) and 1 in Equation (7) have been computed only for the coarse NDVI/NDWI pixels with values ±20% of thin pixel value to be estimated (similar spectra index pixels). Finally, thin pixel LST is obtained as:

Vegetation Parametres
FVC, as the fraction of the total pixel that is covered by canopy, is estimated from NDVI, according to [50]: where NDVIs and NDVIv are representative values for bare soil and green vegetation, respectively. These values were estimated as 0.15 and 0.9, respectively. Values of bare soil (FVC = 0) and full vegetation cover (FVC = 1) have been obtained from the NDVI histograms retrieved during the period of maximal vegetation grow. As crop type does not vary a lot year over year, these soil and vegetation NDVI limits tend to be stable every year. FVC estimates are then employed to compute LSE and the Leaf Area Index (LAI).
The LAI is calculated as [51]: where k(θ) is the light extinction coefficient, which is a measure of attenuation of radiation in the canopy for a given solar zenith angle (θ). This last depends on terrain geometry, solar declination and elevation angle, latitude and day of the year. In our case, it was estimated equal to 0.5.

Albedo
Albedo (α) is computed as the integration of at-surface reflectance across the shortwave spectrum in a different way for Landsat 7 (L7), Landsat 8 (L8) and Sentinel 2 (S2) because of the differences of the three sensors. As S2 and L8 bands have good agreement and high band value correlation [52], albedo was extracted in the same way for both sensors, following the relation proposed by [53]: where b1 to b12 represents the surface reflectance for each band and, the coefficients, are representative of the surface solar radiation fraction within the spectral range for each band.

FEST-EWB Hydrological Modelling
FEST-EWB is a distributed hydrological model based on the energy-water balance system [16], which has been developed from the simple water balance FEST-WB [32], allowing to continuously compute all the main fluxes of the hydrological cycle distributed in space (e.g., evapotranspiration, soil moisture, infiltration, deep percolation, snow dynamic and flow routing). FEST-EWB has been tested in several case studies proving accurate estimates of soil moisture and evapotranspiration against ground measurements from eddy covariance stations [16] and also in agricultural basin scale in comparison with satellite data of LST at low spatial resolution (1 km) [54]. The FEST-EWB model has also been calibrated and validated at the basin scale by comparison with satellite LST and ground discharge data for the Po and the Yangtze river basins [29][30][31]. In this paper, new issues are analysed: (1) providing irrigation volumes as input to the model, which is distributed over space considering only the irrigated crops areas retrieved from satellite data, (2) using high spatial resolution LST data (30 m) for model calibration and validation.
The energy and mass balance equations are solved looking for the representative equilibrium temperature (RET), the temperature that closes the energy budget equation and that regulates the water and energy fluxes. The system of equations, which are linked through evapotranspiration, is the model core: where is the sensible heat and LE (Wm −2 ) the latent heat flux, and ∆W/∆t (Wm −2 ) encloses the energy storage terms. The energy budget equation estimates all the energy fluxes as a function of RET. Then, the latent heat flux (e.g., evapotranspiration equal to less than two parametres) is used as input to the water balance equation besides the soil moisture at the previous time step, to compute the new soil moisture update. The main input parametres to the model are: (1) meteorological forcings: air temperature and relative humidity, incoming shortwave radiation, wind speed and rainfall; (2) soil hydraulic parametres, as the saturated hydraulic conductivity (ksat), the field capacity (fc), wilting point (wp), residual (θr) and saturated (θs) soil water content, Brooks-Corey index (BC), bubbling pressure (bp) and (3) vegetation parametres, as leaf area index (LAI), vegetation fraction and minimum stomatal resistance (rsmin); and (4) the digital elevation model (DEM) and land use/cover map.
The simulations are performed continuously in time for the period between 2013 and 2018 at 1-hour temporal resolution and at 30 m spatial resolution. The FEST-EWB model is run in a simplified version without considering the runoff and groundwater flow, which are considered to be not relevant in the case study area. In fact, no major rivers are present in the case study area and the evapotranspiration and SM dynamic processes are not directly related to flow routing; while the groundwater table is not emerging and interacting with the superficial soil layer [33].

Irrigation Distribution into FEST-EWB
The Capitanata irrigation consortium, as previously described, is a highly irrigated area; hence irrigation water distribution needs to be included in the FEST-EWB model to correctly simulate the water and energy fluxes. However, only daily volumes data measured in the distribution aqueduct of Sud Fortore district are available for each of the analysed years. During the different years the volumes ranges over the whole season between 2.1 and 6 ×10 7 m 3 with a mean value of 4.6 × 10 7 m 3 . These volumes are distributed over the consortium area as mm/day of irrigation considering only the possible irrigated crop areas. In fact, during spring and summer, the irrigated crops are mainly tomatoes and a few plurennial crops (asparagus, vineyard and olives trees); while during late summer and the beginning of autumn, besides the plurennial crops, the irrigated crops are all the vegetables. These possible irrigated crops areas are then identified using the satellite data of the vegetation fraction to discriminate between bare soil and vegetated areas (FVC > 0.05) and then the NDVI maps are employed to discriminate between tomatoes and wheat. In fact, between April and May NDVI values of a wheat field are much higher (0.7-0.9) than a tomato field with small plants (0.1-0.3) (Figure 3). In fact, wheat is not irrigated and it is harvested at the beginning of June. Hence, the cultivated areas that are irrigated are about 40 % of the whole area during summer. Similar results are obtained for all the four analysed years.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 26 W/t (Wm −2 ) encloses the energy storage terms. The energy budget equation estimates all the energy fluxes as a function of RET. Then, the latent heat flux (e.g., evapotranspiration equal to less than two parametres) is used as input to the water balance equation besides the soil moisture at the previous time step, to compute the new soil moisture update. The main input parametres to the model are: (1) meteorological forcings: air temperature and relative humidity, incoming shortwave radiation, wind speed and rainfall; (2) soil hydraulic parametres, as the saturated hydraulic conductivity (ksat), the field capacity (fc), wilting point (wp), residual (r) and saturated (s) soil water content, Brooks-Corey index (BC), bubbling pressure (bp) and (3) vegetation parametres, as leaf area index (LAI), vegetation fraction and minimum stomatal resistance (rsmin); and (4) the digital elevation model (DEM) and land use/cover map.
The simulations are performed continuously in time for the period between 2013 and 2018 at 1hour temporal resolution and at 30 m spatial resolution. The FEST-EWB model is run in a simplified version without considering the runoff and groundwater flow, which are considered to be not relevant in the case study area. In fact, no major rivers are present in the case study area and the evapotranspiration and SM dynamic processes are not directly related to flow routing; while the groundwater table is not emerging and interacting with the superficial soil layer [33].

Irrigation Distribution into FEST-EWB
The Capitanata irrigation consortium, as previously described, is a highly irrigated area; hence irrigation water distribution needs to be included in the FEST-EWB model to correctly simulate the water and energy fluxes. However, only daily volumes data measured in the distribution aqueduct of Sud Fortore district are available for each of the analysed years. During the different years the volumes ranges over the whole season between 2.1 and 6 10^7 m 3 with a mean value of 4.6 10^7 m 3 . These volumes are distributed over the consortium area as mm/day of irrigation considering only the possible irrigated crop areas. In fact, during spring and summer, the irrigated crops are mainly tomatoes and a few plurennial crops (asparagus, vineyard and olives trees); while during late summer and the beginning of autumn, besides the plurennial crops, the irrigated crops are all the vegetables. These possible irrigated crops areas are then identified using the satellite data of the vegetation fraction to discriminate between bare soil and vegetated areas (FVC > 0.05) and then the NDVI maps are employed to discriminate between tomatoes and wheat. In fact, between April and May NDVI values of a wheat field are much higher (0.7-0.9) than a tomato field with small plants (0.1-0.3) (Figure 3). In fact, wheat is not irrigated and it is harvested at the beginning of June. Hence, the cultivated areas that are irrigated are about 40 % of the whole area during summer. Similar results are obtained for all the four analysed years.

Calibration Methodology with Satellite Land Surface Temperature
The satellite-based calibration procedure is a physically based process where each single pixel is modified independently from the other; hence, in each pixel of the analysed domain the minimum of the pixel difference between modelled RET and satellite observed LST is searched over the period of calibration. In the case study analysed in this paper, the calibration period is in between 2013 and 2016, meaning that 22 satellite images are employed for calibrating more than 600,000 pixels (e.g., an area of about 55,000 ha at 30 m of spatial resolution). This procedure improves the traditional calibration based on local soil moisture or discharge data, which allows a global calibration of the parametres that are modified by a single factor over the entire basin [29]. The satellite-based technique applied to calibrate the FEST-EWB model guarantees a pixelwise improvement of soil moisture and evapotranspiration estimates.
The calibration is performed with a "trial and error" approach, where the soil hydraulic and vegetation parametres are tuned in each pixel independently from the other according to the mean pixel error for all the analysed images. A sensitivity analysis performed by [29,30] highlighted that the FEST-EWB parametres that have a relevant impact on the water fluxes and thus subjected to calibration are: soil hydraulic conductivity, soil depth, minimum stomatal resistance, soil evaporation resistance, aerodynamic resistance, Brooks-Corey index.
The tuning of the parametres is done step by step modifying each single parametre or combination of parametres until the difference between RET and LST is minimised. The initial run of the FEST-EWB model is performed using literature-suggested values (Section 2.2).
The parametres are then modified in each pixel independently from the other by a percentage that varies in each pixel according to the pixel by pixel difference values between RET and LST. This percentage is modified randomly in each simulation trial. The parametres are always modified keeping their values within their physical ranges of each soil class, as defined by [22]. The simulations are repeated with several combinations of parametre changes until a parametres' values combination reaches the minimum pixel by pixel difference between RET and LST [55].
The calibration effectiveness of the FEST-EWB model will be evaluated using the common statistical indices: the mean bias error (MBE), the absolute mean bias error (AMBE), the root mean square error (RMSE).
Moreover, the improvement in spatial variability representation of RET, evapotranspiration and the connected water fluxes, and hence of irrigation spatial distribution, is analysed considering the mutual relationship between its values in each pixel through the spatial autocorrelation function (AC). AC allows an analysis of this relationship between different pixel values at a fixed distance: where µ is the mean and σ 2 is the LST variance in the stationary hypothesis, so that a stochastic process, of which the joint probability distribution does not change in time or space, is considered. X 1 and X 2 are the generic positions at a fixed distance d. The autocorrelation function has been studied under an isotropy hypothesis so that d is a function only of the distance between two points and not of the direction. The goodness of spatial variability of FEST-EWB estimates is also evaluated studying their temporal relationships using the Pearson correlation analysis (r) (e.g., a measure of linear relationship between two variables [56]).

Satellite Estimates and Validation
Satellite data time series have been processed from year 2011 to 2018. Landsat 7 and Landsat 8 data are available every 16 days, while Sentinel 2 data every five days. Visible and Near-Infrared data have been obtained from Operational Land Imager (OLI) and Multispectral Instrument (MSI) sensors performing 16-day composite images for monitoring Albedo, FVC and LAI variables. All images have been downloaded from NASA and ESA servers and have been atmospherically corrected with 6S (Landsat-8 case) or Sen2Corr (Sentinel-2 case) software to estimate the bottom of atmosphere reflectance at a spatial resolution of 30 m. Thermal Infrared (TIR) data have also been extracted from TIRS and ETM+ to estimate the LST two times per 16 days. To meet the VNIR spatial resolution data, a sharpening process was performed to carry the resolution from 100 m (TIRS) and 60 m (ETM+) to 30m. Figure 4 shows the behaviour of those variables for the period 2011-2016 averaged over the whole Capitanata area (e.g., case study area). The evolution shows FVC and LAI peaks on winter and until the beginning of April, which is expected, as this period is the most intensive cultivation time. Regarding the albedo, the highest values are during the May-June period, which coincides with the wheat ripening. Then of course, the LST peak values are between June and August.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 26 have been obtained from Operational Land Imager (OLI) and Multispectral Instrument (MSI) sensors performing 16-day composite images for monitoring Albedo, FVC and LAI variables. All images have been downloaded from NASA and ESA servers and have been atmospherically corrected with 6S (Landsat-8 case) or Sen2Corr (Sentinel-2 case) software to estimate the bottom of atmosphere reflectance at a spatial resolution of 30 m. Thermal Infrared (TIR) data have also been extracted from TIRS and ETM+ to estimate the LST two times per 16 days. To meet the VNIR spatial resolution data, a sharpening process was performed to carry the resolution from 100 m (TIRS) and 60 m (ETM+) to 30m. Figure 4 shows the behaviour of those variables for the period 2011-2016 averaged over the whole Capitanata area (e.g., case study area). The evolution shows FVC and LAI peaks on winter and until the beginning of April, which is expected, as this period is the most intensive cultivation time. Regarding the albedo, the highest values are during the May-June period, which coincides with the wheat ripening. Then of course, the LST peak values are between June and August.

Satellite Product Validation
The reliability of remotely sensed land surface temperature and albedo estimates was analysed by comparison with ground measures performed at the location of the two eddy covariance stations. The main obstacle in this type of analysis is given by the resolution of satellite imagery, which is lower than the punctual in situ data. The field radiometres have a footprint of about 2-5m 2 according its instrument height. This implies the strong hypothesis that in situ data have the same homogenous values in the closer 30m, independently from the conditions of the closer area.
In Figure 5, the comparison between ground and satellite LST and albedo is reported. In the case of LST, MBE values are close to zero with an RMSE of 3.0 K and R 2 0.88 and an RMSE of 2.0 K and a R 2 of 0.92 for ETM+ and TIRS estimates, respectively. These values are in concordance with other validations performed in [15,52]. In the case of the accuracy of remote sensing albedo, the data were validated against a ground net radiometre. The assessment shows a good agreement with ground data with an MBE of −0.007 and RMSE of 0.033 with a R 2 of 0.1.

Satellite Product Validation
The reliability of remotely sensed land surface temperature and albedo estimates was analysed by comparison with ground measures performed at the location of the two eddy covariance stations. The main obstacle in this type of analysis is given by the resolution of satellite imagery, which is lower than the punctual in situ data. The field radiometres have a footprint of about 2-5m 2 according its instrument height. This implies the strong hypothesis that in situ data have the same homogenous values in the closer 30m, independently from the conditions of the closer area.
In Figure 5, the comparison between ground and satellite LST and albedo is reported. In the case of LST, MBE values are close to zero with an RMSE of 3.0 K and R 2 0.88 and an RMSE of 2.0 K and a R 2 of 0.92 for ETM+ and TIRS estimates, respectively. These values are in concordance with other validations performed in [15,52]. In the case of the accuracy of remote sensing albedo, the data were validated against a ground net radiometre. The assessment shows a good agreement with ground data with an MBE of −0.007 and RMSE of 0.033 with a R 2 of 0.1.

FEST-EWB Calibration and Validation
RET estimates from FEST-EWB are calibrated against the remotely sensed LST and the effect of the pixel by pixel calibration is also evaluated on the calibrated model parametres. The goodness of the spatial variability of the estimated RET, evapotranspiration and soil moisture are analysed with the AC function and the Pearson coefficient, also considering the effect of irrigation distribution in the model.
Finally, the FEST-EWB model is validated by comparing the simulated energy and water fluxes (LE, SM, H, Rn) and the land surface temperature with the observed ones at the eddy covariance stations location.

Calibration and Validation against Satellite LST
RET from FEST-EWB with the original parametres generally highly underestimates the remotely sensed LST, with a mean absolute difference pixel by pixel of 5 °C; while after the calibration procedure it is equal to 2.7 °C. In Figure 6a these differences pixel by pixel between LANDSAT LST and the calibrated and not calibrated RET are shown, highlighting the improvement in the LST estimates after the calibration. In Figure 6b, the mean aerial value over the whole Capitanata area of each LANDSAT image is reported along with the mean value of RET before and after the calibration procedure. The two LANDSAT sensors are also considered separately, and similar differences with FEST-EWB RET are obtained with L7 AMBE of 2.9 °C and L8 AMBE of 2.6 °C.

FEST-EWB Calibration and Validation
RET estimates from FEST-EWB are calibrated against the remotely sensed LST and the effect of the pixel by pixel calibration is also evaluated on the calibrated model parametres. The goodness of the spatial variability of the estimated RET, evapotranspiration and soil moisture are analysed with the AC function and the Pearson coefficient, also considering the effect of irrigation distribution in the model.
Finally, the FEST-EWB model is validated by comparing the simulated energy and water fluxes (LE, SM, H, Rn) and the land surface temperature with the observed ones at the eddy covariance stations location.

Calibration and Validation against Satellite LST
RET from FEST-EWB with the original parametres generally highly underestimates the remotely sensed LST, with a mean absolute difference pixel by pixel of 5 • C; while after the calibration procedure it is equal to 2.7 • C. In Figure 6a these differences pixel by pixel between LANDSAT LST and the calibrated and not calibrated RET are shown, highlighting the improvement in the LST estimates after the calibration. In Figure 6b, the mean aerial value over the whole Capitanata area of each LANDSAT image is reported along with the mean value of RET before and after the calibration procedure. The two LANDSAT sensors are also considered separately, and similar differences with FEST-EWB RET are obtained with L7 AMBE of 2.9 • C and L8 AMBE of 2.6 • C.

FEST-EWB Calibration and Validation
RET estimates from FEST-EWB are calibrated against the remotely sensed LST and the effect of the pixel by pixel calibration is also evaluated on the calibrated model parametres. The goodness of the spatial variability of the estimated RET, evapotranspiration and soil moisture are analysed with the AC function and the Pearson coefficient, also considering the effect of irrigation distribution in the model.
Finally, the FEST-EWB model is validated by comparing the simulated energy and water fluxes (LE, SM, H, Rn) and the land surface temperature with the observed ones at the eddy covariance stations location.

Calibration and Validation against Satellite LST
RET from FEST-EWB with the original parametres generally highly underestimates the remotely sensed LST, with a mean absolute difference pixel by pixel of 5 °C; while after the calibration procedure it is equal to 2.7 °C. In Figure 6a these differences pixel by pixel between LANDSAT LST and the calibrated and not calibrated RET are shown, highlighting the improvement in the LST estimates after the calibration. In Figure 6b, the mean aerial value over the whole Capitanata area of each LANDSAT image is reported along with the mean value of RET before and after the calibration procedure. The two LANDSAT sensors are also considered separately, and similar differences with FEST-EWB RET are obtained with L7 AMBE of 2.9 °C and L8 AMBE of 2.6 °C.  The model is then validated during 2017 and 2018 against L7 and L8 land surface temperature data. A mean absolute difference pixel by pixel of 3.5 °C is obtained, with similar values between the two LANDSAT sensors. In Figure 8 the mean spatial value of each LANDSAT LST image is reported along with the mean value of RET and also the mean absolute differences pixel by pixel confirming the goodness of the FEST-EWB model also in the validation phase. The model is then validated during 2017 and 2018 against L7 and L8 land surface temperature data. A mean absolute difference pixel by pixel of 3.5 • C is obtained, with similar values between the two LANDSAT sensors. In Figure 8 the mean spatial value of each LANDSAT LST image is reported along with the mean value of RET and also the mean absolute differences pixel by pixel confirming the goodness of the FEST-EWB model also in the validation phase. The model is then validated during 2017 and 2018 against L7 and L8 land surface temperature data. A mean absolute difference pixel by pixel of 3.5 °C is obtained, with similar values between the two LANDSAT sensors. In Figure 8 the mean spatial value of each LANDSAT LST image is reported along with the mean value of RET and also the mean absolute differences pixel by pixel confirming the goodness of the FEST-EWB model also in the validation phase.

Soil Hydraulic and Vegetation Parametres
Applying the pixel by pixel calibration procedure, the parametres change their spatial variability increasing their variability. In Table 1, the values of the main parametres are reported before and after the calibration process showing an increase of the standard deviation values after the calibration.
In Figure 9, an example of soil hydraulic and vegetation parametres before and after the calibration is shown, reporting a general decrease in the saturated hydraulic conductivity, an increase in soil depth and a decrease in the minimum stomatal resistance with a global increase in the spatial standard deviation for all parametres: from 2.8*10 8 to 4.5*10 7 ms −1 for the saturated hydraulic conductivity, from 0.5 to 0.78 m for the soil depth and from 40.8 to 75 sm −1 for the minimum stomatal resistance.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 26 Figure 8. During the validation phase, the mean spatial value of LANDSAT LST and RET after the calibration.

Soil Hydraulic and Vegetation Parametres
Applying the pixel by pixel calibration procedure, the parametres change their spatial variability increasing their variability. In Table 1, the values of the main parametres are reported before and after the calibration process showing an increase of the standard deviation values after the calibration.
In Figure 9, an example of soil hydraulic and vegetation parametres before and after the calibration is shown, reporting a general decrease in the saturated hydraulic conductivity, an increase in soil depth and a decrease in the minimum stomatal resistance with a global increase in the spatial standard deviation for all parametres: from 2.8*10 8 to 4.5*10 7 ms −1 for the saturated hydraulic conductivity, from 0.5 to 0.78 m for the soil depth and from 40.8 to 75 sm −1 for the minimum stomatal resistance.

Spatial Variability and Correlation
AC functions are computed for all images to understand improvements in fluxes' spatial variability representation. Figure 10 reports as an example the AC values for LST and NDVI from LANDSAT 8 images, and then for RET, SM and ET from FEST-EWB for a summer and a winter image: 08/07/2017 and 31/12/2017. The autocorrelation functions of satellite and modelled land surface temperature are almost the same till 600 m of distance, confirming the ability of FETS-EWB in reproducing satellite data. The AC values are equal to 1 at a 0 m distance and decreases till values near zero as the distance between the two pixels increases. The autocorrelation is computed till 600 m of distance, because at higher distances a scarce number of points couples are present. The obtained results confirm that the relationship between pixels at different LST is related to the fields dimensions, which are well delineated, and to the presence of bare soil or vegetation types at different growth stages and also to different soil moisture conditions. This is observable either for summer images and for winter images.

Spatial Variability and Correlation
AC functions are computed for all images to understand improvements in fluxes' spatial variability representation. Figure 10 reports as an example the AC values for LST and NDVI from LANDSAT 8 images, and then for RET, SM and ET from FEST-EWB for a summer and a winter image: 08/07/2017 and 31/12/2017. The autocorrelation functions of satellite and modelled land surface temperature are almost the same till 600 m of distance, confirming the ability of FETS-EWB in reproducing satellite data.
The AC values are equal to 1 at a 0 m distance and decreases till values near zero as the distance between the two pixels increases. The autocorrelation is computed till 600 m of distance, because at higher distances a scarce number of points couples are present. The obtained results confirm that the relationship between pixels at different LST is related to the fields dimensions, which are well delineated, and to the presence of bare soil or vegetation types at different growth stages and also to different soil moisture conditions. This is observable either for summer images and for winter images. Similar autocorrelation functions are obtained for the modelled evapotranspiration which have a spatial variability correlated to the RET one.
In Figure 10, as an example, the autocorrelation function is also reported for the saturated hydraulic conductivity ksat, which has been previously calibrated. AC of ksat was less correlated than LST with lower values of AC, as expected. In fact, the soil structure variability is in general independent from the fields shapes and the available information are at a lower scale than the 30 m of resolution of the satellite data. After the pixel by pixel calibration procedure of the FEST-EWB model, the autocorrelation of ksat is correctly increased. Pearson correlation analysis (r) is then computed among the satellite LST and NDVI with the modelled RET, soil moisture and evapotranspiration ( Table 3). As expected high correlation is found between modelled variables, which derives from the interconnected water-energy balances equations. Hence, high correlation (0.75-0.6) is obtained between RET and ET and SM, respectively, during the winter season; while during the irrigated season (e.g., from April to September) the variables are negatively correlated, as expected. Similar correlations are also found between LST from remote sensing the modelled ET and SM, confirming the goodness of the FEST-EWB model in reproducing the LST temporal variability. The high spatial correlation between the modelled SM and respectively the LST and RET is shown in Figure 11 highlighting the ability of FEST-EWB model in also reproducing the spatial variability of LST, with a high positive correlation in almost all vegetated pixels between either modelled and observed LST with soil moisture during the winter season and instead a lower correlation during the summer season. This is mainly due to the main cultivated crops in the Capitanata consortium: the irrigated tomatoes and the non-irrigated wheat. Table 3. Pearson correlation coefficient between modelled and observed variables during the irrigation and no irrigation seasons. Similar autocorrelation functions are obtained for the modelled evapotranspiration which have a spatial variability correlated to the RET one.

No Irrigation Season
In Figure 10, as an example, the autocorrelation function is also reported for the saturated hydraulic conductivity ksat, which has been previously calibrated. AC of ksat was less correlated than LST with lower values of AC, as expected. In fact, the soil structure variability is in general independent from the fields shapes and the available information are at a lower scale than the 30 m of resolution of the satellite data. After the pixel by pixel calibration procedure of the FEST-EWB model, the autocorrelation of ksat is correctly increased.
Pearson correlation analysis (r) is then computed among the satellite LST and NDVI with the modelled RET, soil moisture and evapotranspiration ( Table 3). As expected high correlation is found between modelled variables, which derives from the interconnected water-energy balances equations. Hence, high correlation (0.75-0.6) is obtained between RET and ET and SM, respectively, during the winter season; while during the irrigated season (e.g., from April to September) the variables are negatively correlated, as expected. Similar correlations are also found between LST from remote sensing the modelled ET and SM, confirming the goodness of the FEST-EWB model in reproducing the LST temporal variability. The high spatial correlation between the modelled SM and respectively the LST and RET is shown in Figure 11 highlighting the ability of FEST-EWB model in also reproducing the spatial variability of LST, with a high positive correlation in almost all vegetated pixels between either modelled and observed LST with soil moisture during the winter season and instead a lower correlation during the summer season. This is mainly due to the main cultivated crops in the Capitanata consortium: the irrigated tomatoes and the non-irrigated wheat. Table 3. Pearson correlation coefficient between modelled and observed variables during the irrigation and no irrigation seasons.

Effect of Irrigation Distribution on RET Estimates
The goodness of the irrigation distribution methodology into the FEST-EWB model is evaluated considering the differences between the measured and modelled LST patterns. Besides the simulation performed for the calibration and validation (irrigation distributed considering FVC and NDVI), two other simulations are performed with: (1) a uniform distribution of irrigation over the whole consortium area, (2) a distribution of irrigation considering only the irrigated areas with a vegetation fraction higher than 0.05. The averaged spatial value of the pixel by pixel difference between satellite LST and RET is computed over the irrigation season also for these two irrigation distribution methodologies and values of 4.7 and 4°C are obtained. The AMBE obtained for the original selected irrigation distribution based on both FVC and NDVI was of 2.7 °C. In Table 4, the AMBE values are reported considering vegetated and not vegetated areas separately: in not vegetated areas similar values are obtained for the three configurations (AMBE about 3 °C); while for the vegetated areas high differences among the irrigation configurations are obtained (AMBE from 7.9 to 2.5 °C).
As an example in Figure 12, RET estimates with the different irrigation distribution techniques into the FEST-EWB model are compared with satellite LST for the date of 9 July 2017. It is clearly visible that the uniform distribution of irrigation produces an RET estimate higher than the observed one of 7.5 °C, which is due to the less amount of water distributed in the crop pixels; while if irrigation water is distributed only over the vegetated fields and even more considering NDVI, in the RET maps the cooler fields are clearly visible as in the satellite map with a mean difference of 5.7 and 2.6 °C, respectively.

Effect of Irrigation Distribution on RET Estimates
The goodness of the irrigation distribution methodology into the FEST-EWB model is evaluated considering the differences between the measured and modelled LST patterns. Besides the simulation performed for the calibration and validation (irrigation distributed considering FVC and NDVI), two other simulations are performed with: (1) a uniform distribution of irrigation over the whole consortium area, (2) a distribution of irrigation considering only the irrigated areas with a vegetation fraction higher than 0.05. The averaged spatial value of the pixel by pixel difference between satellite LST and RET is computed over the irrigation season also for these two irrigation distribution methodologies and values of 4.7 and 4 • C are obtained. The AMBE obtained for the original selected irrigation distribution based on both FVC and NDVI was of 2.7 • C. In Table 4, the AMBE values are reported considering vegetated and not vegetated areas separately: in not vegetated areas similar values are obtained for the three configurations (AMBE about 3 • C); while for the vegetated areas high differences among the irrigation configurations are obtained (AMBE from 7.9 to 2.5 • C). As an example in Figure 12, RET estimates with the different irrigation distribution techniques into the FEST-EWB model are compared with satellite LST for the date of 9 July 2017. It is clearly visible that the uniform distribution of irrigation produces an RET estimate higher than the observed one of 7.5 • C, which is due to the less amount of water distributed in the crop pixels; while if irrigation water is distributed only over the vegetated fields and even more considering NDVI, in the RET maps the cooler fields are clearly visible as in the satellite map with a mean difference of 5.7 and 2.6 • C, respectively.

Evapotranspiration and Energy Fluxes Validation at Local Scale
The calibration of FEST-EWB parametres lead to changes not only in the RET estimates, but also in all the connected variables through the energy-water balance system, as SM and ET. Hence, FEST-EWB is validated by comparing the estimated energy and water fluxes with the measured fluxes at the eddy covariance stations. A general good agreement is observable confirming the ability of the model in producing high spatial resolution and continuous in time estimates of evapotranspiration and soil moisture. In Table 5, the results are shown for the six considered fields (e.g., crop seasons), reporting the RMSE and, the angular coefficient and the coefficient of determination of the regression line between observed and modelled variables.
For the plurennial crop asparagus field, a good agreement is obtained for LE with a coefficient of determination equal to 0.85 and an RMSE of 33.3 W m −2 , while for H a R 2 of 0.76 and an RMSE of 32.7 W m −2 are obtained. Soil moisture is also correctly reproduced with an RMSE of 0.04. For the cabbage field, an overall good agreement is also obtained not only for LST, but also for the other water and energy fluxes. In particular, soil moisture is correctly reproduced with an RMSE of 0.1, showing the SM increments in correspondence of irrigations and rainfall. The agreement is confirmed for LE by a coefficient of determination equal to 0.7 and an RMSE of 49 W m −2 , while for H with R2 of 0.8 and RMSE of 78 W m −2 . The final error between observed and simulated cumulated ET is equal to 2.7 %. In the tomato field with sandy soil for the 2016 season, the FEST-EWB model is able to correctly reproduce SM with an RMSE of 0.07. The good performances are also obtained for all the energy fluxes with an RMSE around 50 W m −2 except LE with an RMSE of 43 W m −2 . Similar results are also obtained for the tomato field with silty clay soil for the same 2016 season. In fact, an RMSE

Evapotranspiration and Energy Fluxes Validation at Local Scale
The calibration of FEST-EWB parametres lead to changes not only in the RET estimates, but also in all the connected variables through the energy-water balance system, as SM and ET. Hence, FEST-EWB is validated by comparing the estimated energy and water fluxes with the measured fluxes at the eddy covariance stations. A general good agreement is observable confirming the ability of the model in producing high spatial resolution and continuous in time estimates of evapotranspiration and soil moisture. In Table 5, the results are shown for the six considered fields (e.g., crop seasons), reporting the RMSE and, the angular coefficient and the coefficient of determination of the regression line between observed and modelled variables. For the plurennial crop asparagus field, a good agreement is obtained for LE with a coefficient of determination equal to 0.85 and an RMSE of 33.3 W m −2 , while for H a R 2 of 0.76 and an RMSE of 32.7 W m −2 are obtained. Soil moisture is also correctly reproduced with an RMSE of 0.04. For the cabbage field, an overall good agreement is also obtained not only for LST, but also for the other water and energy fluxes. In particular, soil moisture is correctly reproduced with an RMSE of 0.1, showing the SM increments in correspondence of irrigations and rainfall. The agreement is confirmed for LE by a coefficient of determination equal to 0.7 and an RMSE of 49 W m −2 , while for H with R2 of 0.8 and RMSE of 78 W m −2 . The final error between observed and simulated cumulated ET is equal to 2.7%. In the tomato field with sandy soil for the 2016 season, the FEST-EWB model is able to correctly reproduce SM with an RMSE of 0.07. The good performances are also obtained for all the energy fluxes with an RMSE around 50 W m −2 except LE with an RMSE of 43 W m −2 . Similar results are also obtained for the tomato field with silty clay soil for the same 2016 season. In fact, an RMSE of 0.06 is obtained for SM, while for LST 3.5 • C and for H and LE 75 and 40 W m −2 , respectively. In the same field, for the 2017 tomatoes season, an overall good agreement is obtained not only for LST with the angular coefficient of one and a coefficient of determination of one, but also on the other water and energy fluxes. In particular, SM has an RMSE of 0.04, showing the SM increments in correspondence of irrigations and rainfall. The agreement is confirmed for LE by an angular coefficient equal to 0.92 and an RMSE of 43 W m −2 , while for H with m of 0.8 and RMSE of 65 W m −2 . Finally, the results for the wheat field confirm the ability of the FEST-EWB model to correctly reproduce SM with an RMSE of 0.04 even in the winter season, LST has an RMSE of 4.1 • C, and LE and H an RMSE of 48 and 67 W m −2 , respectively.
In Figure 13, the scatterplots between ground measured LST and the calibrated RET are reported for all the stations. A good agreement is found for the asparagus field with an angular coefficient (m) of the linear regression forced through the origin between RET (x) and LST (y) equal to 0.94 and a R 2 of 0.96; m is equal to 1 for the cabbage field with R 2 of 0.9; an m of 0.98 is obtained for the tomato field with silty clay soil (R 2 = 0.83) while a lower m value of 0.90 is obtained in the tomato field with sandy soil (R 2 = 0.90). Both tomato fields in 2017 and the wheat field have an angular coefficient of 1, while R 2 is of 1 and 0.4, respectively. The lowest agreement obtained in the wheat field may be due to the different time period of the wheat crop season (from October to May) with respect to the other summer crops and irrigation volumes. The comparison between observed and simulated soil moisture before and after the calibration is shown in Figure 14 for all the validation stations. Similar results are reported for all the fields confirming the ability of FEST-EWB, after its calibration with LST data, to correctly reproduce the observed SM in different soil types, period of the year and irrigation volumes. In fact, cabbage, wheat and partially the asparagus fields are winter crops while tomatoes are summer ones. Similar errors are obtained with RMSE between 0.1 and 0.04. Similar results are obtained from the comparison between RET and LST from satellite data in the station pixel with m and R 2 around 0.9-1, except for the two tomato fields during 2016 with R2 around 0.4. This may be due to the very few numbers of available satellite images.
The comparison between observed and simulated soil moisture before and after the calibration is shown in Figure 14 for all the validation stations. Similar results are reported for all the fields confirming the ability of FEST-EWB, after its calibration with LST data, to correctly reproduce the observed SM in different soil types, period of the year and irrigation volumes. In fact, cabbage, wheat and partially the asparagus fields are winter crops while tomatoes are summer ones. Similar errors are obtained with RMSE between 0.1 and 0.04.
The model is also able to follow the soil type variability being capable of correctly reproducing observed SM in sandy soils and as well as silty clay soils for the same crop (tomatoes) with very low RMSE differences (between 0.04 and 0.07). In Figure 15, observed and modelled latent heat fluxes are shown for all the analysed fields, showing a good correspondence. Ground data have some missing values due to the correction procedures applied to eddy covariance data (paragraph 2.4). As observed for SM, FEST-EWB is able to correctly reproduce LE under different irrigation regimes and season of the year. In particular, low LE is correctly modelled during winter time for the cabbage and wheat fields, while very high LE is computed during the summer period for the tomato fields.
Moreover, if the cumulated fluxes are computed, the final error between observed and simulated cumulated ET ranges between 10% to 1.9%.  In Figure 15, observed and modelled latent heat fluxes are shown for all the analysed fields, showing a good correspondence. Ground data have some missing values due to the correction procedures applied to eddy covariance data (paragraph 2.4). As observed for SM, FEST-EWB is able to correctly reproduce LE under different irrigation regimes and season of the year. In particular, low LE is correctly modelled during winter time for the cabbage and wheat fields, while very high LE is computed during the summer period for the tomato fields. In Figure 15, observed and modelled latent heat fluxes are shown for all the analysed fields, showing a good correspondence. Ground data have some missing values due to the correction procedures applied to eddy covariance data (paragraph 2.4). As observed for SM, FEST-EWB is able to correctly reproduce LE under different irrigation regimes and season of the year. In particular, low LE is correctly modelled during winter time for the cabbage and wheat fields, while very high LE is computed during the summer period for the tomato fields.
Moreover, if the cumulated fluxes are computed, the final error between observed and simulated cumulated ET ranges between 10% to 1.9%.  Moreover, if the cumulated fluxes are computed, the final error between observed and simulated cumulated ET ranges between 10% to 1.9%.

Discussion
Evapotranspiration estimates computed in this paper with the FEST-EWB model are continuous in time (hourly time step) and at high spatial resolution (30 m). These estimates are at the same resolutions of the ET estimates provided by several remotely-sensed single and two source models which compute evapotranspiration as the residual term of the energy balance equation [9][10][11][12][13]. However, these last models are not able to provide high temporal resolution estimates, but only an instantaneous ET information at satellite overpass time, which may also be influenced by the presence of clouds preventing its computation, or a daily ET estimate based on net radiation weights.
The model has been calibrated against remote sensing data of land surface temperature, which have been preferred over soil moisture information. These seem to be more recommended to calibrate model soil parametres, thanks to their direct estimate of soil moisture and the fact that they do not suffer from cloud coverage problems, but remote sensing SM are usually limited by the few centimetres of soil depth representativeness, the complex reading geometry, the low sensitivity to soil water content above a significant threshold [57], and a low spatial resolution for passive microwave [3,58].
Furthermore, LST image uncertainty should also be taken into account, especially over non-homogeneous areas, where the spatial resolution, scan angle of view of the sensor and emissivity should be considered. For example, a mean bias with ground data equal to 2 K has been evaluated for MODIS land surface temperature data [59], which is comparable to modelled RET errors.
Regarding the FEST-EWB soil parametres, a possible uncertainty on an ET estimate may arise from the soil depth definition, which is a difficult variable to determine. The FEST-EWB model is a single-layer model and in the analysed case study the initial mean area value has been set to 0.5 m, considering the predominant growing zone of fresh vegetable roots, which is relevant for the evapotranspiration process. However, the soil depth is a fixed value over the year, which may lead to an overestimation of evapotranspiration during bare soil period.
Different calibration and optimisation techniques, based on different objective functions as RMSE or the Nash-Sutcliffe criterion and on different minimisation techniques with a local gradient or a global method, may modify the accuracy of the model calibration. The FEST-EWB model has been calibrated with a non-automatic "trial and error" procedure, which is based on several repetitions until the errors do not change anymore but only congruent parametres combinations are tested within their physical values. The topic of model calibration may highlight the unresolved issue of equifinality, which is embedded in all hydrological models [60].
Another remarkable point of attention is the distribution of irrigation into the hydrological model, which is one of the most uncertain input variables in hydrological models. In fact, irrigation plays a crucial role in the Capitanata case study with a relevant impact on soil moisture and in an indirect way on evapotranspiration and on land surface temperature. In this work, irrigation has been distributed according to the satellite FVC and NDVI, thanks to the knowledge of the crop types, showing a better agreement of LST estimates than not considering them, in terms of absolute values and also temporal and spatial variability.

Conclusions
The paper showed the feasibility of the combined advantage of satellite data and distributed hydrological model for improving water resources assessment, in particular evapotranspiration estimates at high temporal and spatial resolution. In particular, satellite LST are employed for a pixel by pixel control of water and energy fluxes accuracy into FEST-EWB model. The modelled RET has been calibrated against LST data from LANDSAT data providing a good agreement in the calibration phase (AMBE 2.6 • C) and a slight increase in the validation phase (AMBE 3.5 • C). No relevant differences are found in the RET calibration if L7 or L8 data are used, confirming also the goodness of the retrieval and downscaling algorithms, which are verified against ground data, returning an RMSE of 0.3 ±2 K.
The model has then been validated against soil moisture and evapotranspiration data measured by six ground eddy covariance stations. Similar results to LST are reported for all the fields confirming the ability of FEST-EWB after its calibration with LST data, to correctly reproduce the observed SM and ET in different soil and crop types, period of the year and irrigation importance. An RMSE between 0.1 and 0.04 is found for soil moisture, and an RMSE around 40 Wm-2 for latent heat flux.