Accounting for Almond Crop Water Use under Different Irrigation Regimes with a Two-Source Energy Balance Model and Copernicus-Based Inputs

: Accounting for water use in agricultural ﬁelds is of vital importance for the future prospects for enhancing water use efﬁciency. Remote sensing techniques, based on modelling surface energy ﬂuxes, such as the two-source energy balance (TSEB), were used to estimate actual evapotranspiration (ET a ) on the basis of shortwave and thermal data. The lack of high temporal and spatial resolution of satellite thermal infrared (TIR) missions has led to new approaches to obtain higher spatial resolution images with a high revisit time. These new approaches take advantage of the high spatial resolution of Sentinel-2 (10–20 m), and the high revisit time of Sentinel-3 (daily). The use of the TSEB model with sharpened temperature (TSEB S2+S3 ) has recently been applied and validated in several study sites. However, none of these studies has applied it in heterogeneous row crops under different water status conditions within the same orchard. This study assessed the TSEB S2+S3 modelling approach to account for almond crop water use under four different irrigation regimes and over four consecutive growing seasons (2017–2020). The energy ﬂuxes were validated with an eddy covariance system and also compared with a soil water balance model. The former reported errors of 90 W/m 2 and 87 W/m 2 for the sensible (H) and latent heat ﬂux (LE), respectively. The comparison of ET a with the soil water balance model showed a root-mean-square deviation (RMSD) ranging from 0.6 to 2.5 mm/day. Differences in cumulative ET a between the irrigation treatments were estimated, with maximum differences obtained in 2019 of 20% to 13% less in the most water-limited treatment compared to the most well-watered one. Therefore, this study demonstrates the feasibility of using the TSEB S2+S3 for monitoring ET a in almond trees under different water regimes.


Introduction
Almond is the most cultivated tree nut crop worldwide. The largest almond-producing countries are Spain and the US, followed by other Mediterranean or semi-arid countries. Since 2014, Spain almond production has increased by an average of 10.3% each year. At 421,610 metric tons in 2020, it is the third-highest almond-producing country, after the US and Australia [1]. In terms of cultivated area, almond is also the third crop in Spain, with 844,244 ha planted [2]. In Spain, almonds have typically been cultivated in marginal areas under rainfed conditions. Nevertheless, the surface area of irrigated almond orchards has increased considerably in recent years, converting it into a high-value crop, which currently represents about 16% of woody crops in Spain.
Some studies have shown that almond yield is positively related to the amount of water applied or consumed [3][4][5]. Several studies have reported that the annual irrigation requirements of intensive mature almond orchards located in California and Spain range Remote Sens. 2022, 14, 2106 3 of 25 sophisticated methods involve data fusion approaches to compute daily ET a at 30 m using the ALEXI/DisALEXI scheme [32,[38][39][40]. The performance of this approach has been successfully applied in different agricultural fields, for example, with Landsat, MODIS or even with the ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) mission [41,42]. On the other hand, the European Commission's Copernicus Program provides operational Earth Observation imagery in near real time with long-term continuity of its services and therefore it is interesting to continue exploring the feasibility of new similar methods.
The recent launch of Sentinel-2 (S2) and Sentinel-3 (S3) satellites opens up the possibility of producing operational agricultural products, since they contain most of the required spatial, temporal and spectral characteristics [43]. Although S3 acquires daily thermal infrared (TIR) data at a coarse spatial resolution of 1000 m at nadir, sharpening approaches have been proposed to downscale the LST to 20 m spatial resolution [44] and thus generate ET a maps at S2 scale. This approach relies on downscaling S3 thermal bands to S2 spatial resolution using a data mining sharpener (DMS; [45]). Consequently, in combination with meteorological data forcing from the Copernicus Climate Data Store (CDS) and S2-derived vegetation biophysical and structural properties, it is possible to obtain daily ET a maps at farm scale. Guzinski et al. [44] showed that the use of sharpened high-resolution LST data as input of the two-source energy balance (TSEB) model [30] was a sound approach, giving a relative root mean square error (RMSE) of instantaneous latent heat flux of around 30% in agricultural areas and, therefore, demonstrating the potential of this approach to estimate ET a time series. Although ET estimates using this approach were validated in different land cover types against in situ eddy covariance (EC) towers, validations corresponding to agricultural crops were mostly performed in fields with hardly any water stress conditions and thus the accuracy of ET a estimates under more severe water stress conditions remained uncertain. More recently, Guzinski et al. [46] extended the range of agricultural sites for model validation, including a wider range of herbaceous and woody crops under rainfed and irrigated conditions in semi-arid fields of Spain and Tunisia. Their results showed that the TSEB model combined with Copernicus S2+S3 data was able to accurately track ET a with minimal bias (~0 mm/day) and low RMSE (<1 mm/day). However, Bellvert et al. [47] pointed out that one of the limitations in accurately estimating the crop's water use using this approach in grapevines under water stress conditions was the inability of the DMS approach to capture the full range of surface temperatures, particularly when short-term water stress did not reduce the amount of biomass. Accurate estimates of ET a in heterogeneous row crops with limited water availability are also a major challenge with a still high degree of uncertainty. This is explained by the fact that differences in spacing distances, training systems or canopy management affect both the turbulent exchange of heat and water vapor and the accurate estimation of the vegetation parameters.
The aim of this study is to evaluate the performance of the TSEB model with sharpened LST images from S2 and S3 to assess the seasonal evolution of ET a in almond trees irrigated under different water regimes, and to do so with different crop water consumption and states, during four growing seasons (2017-2020). Modelled fluxes obtained with the TSEB model using S2+S3 were validated with an eddy covariance system. In addition, TSEB estimates from Landsat-8 TIR images were also produced as the benchmark for thermalbased ET a product, and the ET a of each irrigation treatment was also compared against the estimates obtained through a soil water balance approach.

Study Site
The experiment was carried out in a 7.5 ha almond orchard located in Maials, Spain (41 • 22.9 172 N, 0 • 31.27 619"E) from 2017 to 2020 (see Figure 1). The climate in the region is typically Mediterranean with a mean annual temperature of 14.5 • C and average annual precipitation of 307.9 mm, mainly occurring during autumn (October-December) and spring (March-June). The accumulated average reference evapotranspiration (ET 0 ) during the growing season (March-October) was 960 mm. Soil type is loamy with a slightly calcaric influence due to the petrocalcic origins of the rock. Soil depth varied from 120 to 180 cm, depending on the area. Almond tree management for pruning; disease and pest control; and soil management and fertilisation was based on Spanish integrated production management practices [48].

Study Site
The experiment was carried out in a 7.5 ha almond orchard located in Maials, Spain (41°22.9′172″N, 0°31.27′619″E) from 2017 to 2020 (see Figure 1). The climate in the region is typically Mediterranean with a mean annual temperature of 14.5 °C and average annual precipitation of 307.9 mm, mainly occurring during autumn (October-December) and spring (March-June). The accumulated average reference evapotranspiration (ET0) during the growing season (March-October) was 960 mm. Soil type is loamy with a slightly calcaric influence due to the petrocalcic origins of the rock. Soil depth varied from 120 to 180 cm, depending on the area. Almond tree management for pruning; disease and pest control; and soil management and fertilisation was based on Spanish integrated production management practices [48].

Field Experimental Design
Almond trees (cv. Vairo) were grafted onto GF-677 rootstock and planted in 2015 in a 6 × 6 grid (277 trees/ha). Two drip irrigation laterals were placed 80 cm away from the tree rows in 2018. In 2019, the laterals of treatment T4 were moved to 70 cm from the tree rows. Drippers were separated and placed 75 cm apart and discharged 25.6 l/h through pressure-compensating emitters per tree. Irrigation was scheduled on a weekly basis and conducted on a daily basis. Crop water requirements (ETc) were calculated following the methodology proposed by Girona et al. [15]. The Penman-Monteith ET0 was obtained from the closest weather station from the Meteorological Service of Catalonia (SMC), which was located 3.5 km from the orchard (www.ruralcat.net/web/guest/agrometeo.estacions, accessed on 15 July 2021). The crop coefficients (Kc) were adapted from Gold-

Field Experimental Design
Almond trees (cv. Vairo) were grafted onto GF-677 rootstock and planted in 2015 in a 6 × 6 grid (277 trees/ha). Two drip irrigation laterals were placed 80 cm away from the tree rows in 2018. In 2019, the laterals of treatment T4 were moved to 70 cm from the tree rows. Drippers were separated and placed 75 cm apart and discharged 25.6 L/h through pressurecompensating emitters per tree. Irrigation was scheduled on a weekly basis and conducted on a daily basis. Crop water requirements (ETc) were calculated following the methodology proposed by Girona et al. [15]. The Penman-Monteith ET 0 was obtained from the closest weather station from the Meteorological Service of Catalonia (SMC), which was located 3.5 km from the orchard (www.ruralcat.net/web/guest/agrometeo.estacions, accessed on 15 July 2021). The crop coefficients (K c ) were adapted from Goldhamer and Girona [49] and the shading factor (K r ) was taken from Fereres et al. [50]. Therefore, irrigation was performed to replace crop evapotranspiration (ET c ) minus effective precipitation (P eff ).
A randomized complete block design with four block replicates was used. Each block housed four experimental plots with 80 trees in each. Four differential irrigation treatments were adopted: Treatment 1 (T1) consisted of simulating a scenario with seasonal water availability of 9000 m 3 /ha. During the four years of this study, these trees never reached this water demand, so T1 was always fully irrigated. Treatment 2 (T2) consisted of simulating a scenario with seasonal water availability of 6500 m 3 /ha. During the first two years of this study, the trees did not reach the pre-established water demand and so were fully irrigated. However, during 2019 and 2020, sustained deficit irrigation was adopted Remote Sens. 2022, 14, 2106 5 of 25 which consisted of gradually reducing the crop water requirements to 25% ETc over the course of the growing season. Treatment 3 (T3) consisted of simulating a scenario with seasonal water availability of 4500 m 3 /ha. In 2017, almond water demand was below this threshold, and so in this year the trees were fully irrigated. In 2018, a moderate sustained deficit was adopted which consisted of reducing water prescriptions by 10% ETc. In 2019 and 2020, irrigation was scheduled based on a reduction by 35% ETc. Finally, treatment 4 (T4) consisted of simulating a scenario with seasonal water availability of 2500 m 3 /ha. In 2017, almond water demand was below this value, and so the trees were fully irrigated. However, sustained deficit irrigation was applied in all the other years, which consisted of reducing the amount of water to below 2500 m 3 /ha. To achieve this, reductions by 40% ETc in 2018 and by 60% ETc in both 2019 and 2020 were applied. P eff was calculated as half of the rainfall for a single event day with more than 10 mm of precipitation, or was otherwise considered to be zero. The total amount of water applied was measured weekly using digital water meters (CZ2000-3M, Contazara, Zaragoza, Spain).

TSEB Scheme
Most energy balance approaches assume that the total available energy in the surface (net radiation (Rn)-soil heat flux (G)) is partitioned into sensible heat flux (H) and latent heat flux (LE), neglecting other sources of heat dissipation such as biochemical changes (i.e., photosynthesis, respiration) or heat stored in the canopy due to their low contribution in proportion to the overall energy budget.
The TSEB is a model formulated originally by Norman et al. [30] and was afterward modified and improved by Kustas and Norman [51]. In TSEB, each of the components of the energy balance equation is partitioned into canopy and soil, except for G, as it assumes negligible heat storage at the canopy layer.
Equation (2a-c) represent the partition of the latent heat of evaporation (LE), where the subscripts c and s refer to canopy and soil, respectively. G is estimated as a fraction of soil net radiation [52].
The fraction of net radiation which is stored in the soil c G is dependent on the type of soil, soil moisture and time of the day, but typically is set as 0.35 for near-noon conditions [30]. There are other approaches that set G as a function of daytime [53], but for this study it is assumed that using a constant fraction, as in Equation (3), is still valid as the satellite overpasses are within the 2-3 h around solar noon.
Sensible heat flux is partitioned into soil and canopy in which fluxes between soil and canopy are connected in series, as shown in Figure 2, following an analogy of Ohm's law for electric transport.
In Equation (4a-c) the sensible heat flux is partitioned into soil and canopy, where ρ is the density of air, C p is the specific heat of air, T s is the soil temperature, T a is the air Remote Sens. 2022, 14, 2106 6 of 25 temperature, T ac is the temperature in the canopy air space (equivalent to the aerodynamic temperature), r s and r x are the surface boundary layer resistances to heat transfer from the soil and canopy, respectively, to the air-canopy layer, and r a is the aerodynamic resistance to heat transport between the air-canopy layer and the overlying air layer.
These resistances represent the impediment or the enhancement of vertical heat and water transport related to the degree of air turbulence. Therefore, they mainly depend on wind speed, atmospheric stability (i.e., convective transport) and canopy/soil structure (i.e., surface roughness and wind attenuation through the canopy). Although there are several approaches to quantify these resistances [51,54,55], and several parameters that could be fitted [56], we used the original resistance formulation proposed in TSEB [51] with their standard coefficients in order to guarantee the operational aspect of the TSEB S2-S3 approach as proposed by Guzinski et al. [43,44,46].

=
(4b) In Equation (4a-c) the sensible heat flux is partitioned into soil and canopy, where is the density of air, is the specific heat of air, is the soil temperature, is the air temperature, is the temperature in the canopy air space (equivalent to the aerodynamic temperature), and are the surface boundary layer resistances to heat transfer from the soil and canopy, respectively, to the air-canopy layer, and is the aerodynamic resistance to heat transport between the air-canopy layer and the overlying air layer.
These resistances represent the impediment or the enhancement of vertical heat and water transport related to the degree of air turbulence. Therefore, they mainly depend on wind speed, atmospheric stability (i.e., convective transport) and canopy/soil structure (i.e., surface roughness and wind attenuation through the canopy). Although there are several approaches to quantify these resistances [51,54,55], and several parameters that could be fitted [56], we used the original resistance formulation proposed in TSEB [51] with their standard coefficients in order to guarantee the operational aspect of the TSEBS2-S3 approach as proposed by Guzinski et al. [43,44,46]. Soil and canopy temperatures cannot be directly retrieved from coarse-resolution satellite-derived images. However, the directional radiometric temperature observed by the sensor can be considered to be a mixture of warm soil and cool vegetation: Soil and canopy temperatures cannot be directly retrieved from coarse-resolution satellite-derived images. However, the directional radiometric temperature observed by the sensor can be considered to be a mixture of warm soil and cool vegetation: with f c (θ) as the fraction of vegetation observed by the thermal infrared sensor at a zenith angle of θ.
In order to derive T c and T s , the Priestley-Taylor (PT) approach is applied [30,58] and an initial maximum potential rate of transpiration is set.
where ∆ is the slope of the saturated vapor pressure to the air temperature, γ p is the psychometric constant, and f g is the fraction of leaf area index (LAI) that is green and is therefore transpiring. The α PT coefficient is set to 1.26. With this first guess of LE c known, together with Equations (2b), (4b) and (5), it is possible to obtain an initial estimate of T s and therefore of LE s , using Equations (4a) and (2c). LE s at this stage usually yields For this study, high-resolution shortwave observations acquired by the Multispectral Instrument (MSI) mounted on Sentinel-2A and 2B tandem satellites were used. MSI acquires shortwave radiance information in 13 spectral bands with a spatial resolution of 10 m, 20 m or 60 m, depending on the spectral band, and there is a global geometric revisit of at least five days when both satellites are considered, and higher as latitude increases. The MSI differs from previous open data medium-resolution multispectral sensors by including three spectral bands in the red edge part of the electromagnetic spectrum. The red edge responds highly to chlorophyll content as these bands are sensitive to leaf pigments. It also includes the traditional visible and near-infrared bands and two bands in the shortwave infrared, which are sensitive to water content. Sentinel-2 images at Level-2A (L2A; bottom of atmosphere reflectance) were directly retrieved from CREODIAS (https://creodias.eu/, accessed on 20 October 2021), downloading four tiles for this study (T31(TBF, TCG, TBG, TCF)). In order to have a larger dataset of an agricultural area for the sharpening algorithm and improve the accuracy of the machine learning model, the whole area of Lleida (northeast of Spain) (273,940, 4,573,440, 359,220, 4,653,320 m UTM 31 N) was processed and mosaicked to retrieve the biophysical parameters required by TSEB at 20 m spatial resolution. Since the 10 m resolution bands of S2 are not required in this study, band 8, which is a broader bandwidth version of band 8A, was excluded for this study.
The S3 SLSTR sensor was used to obtain the thermal data. This includes three thermal infrared (TIR) channels with 1 km nominal spatial resolution and less than two days temporal resolution with one satellite at the equator and half a day when Sentinel 3A and B are available (April 2018). The LST was obtained directly as an L2A product of S3 SLSTR at 1 km resolution. All the available images acquired in the 10 days ahead of an S2 acquisition in the first semester of 2017 and all SLSTR LST data with 5 days ahead from the second semester of 2017 to the end of 2020 were downloaded. A total of 175 S2 dates and 714 S3 cloud-free scenes were fetched from March 2017 to October 2020.
Landsat-8 is currently the satellite LST source at high spatial resolution, since it includes several bands in the shortwave spectrum at field-scale resolution plus a highresolution thermal band (100 m). The main issue with Landsat-8 is that it only overpasses every 16 days. Temporal resolution is essential when monitoring ET, mainly because the ET dynamics can change within short temporal windows. If Landsat-8 imagery is available in a turnaround time of 16 days, relevant ET information is likely to be lost. Moreover, there is the risk of cloud occlusion during the satellite overpass cycle, which may lead to not having a clear image in more than a month. Therefore, in the present study, Landsat-8 LST was used as the benchmark TSEB product for which results from TSEB S2+S3 could be compared against.

Sentinel-2 Biophysical Parameters of the Vegetation
From the S2 images, the vegetation biophysical parameters were computed through the Biophysical Processor [59] available in the SNAP software v8.0 (https://step.esa.int/ main/download/snap-download/-last accessed on 10 June 2021). This processor relies on building a randomized dataset of vegetation biophysical variables from which reflectance simulations by the radiative transfer models (RTM) PROSPECT [60] and SAIL [61] are produced. Then, the simulated reflectance values are trained against the biophysical parameters through a machine learning approach. The trained model is hence used to predict effective values of green LAI, fractional vegetation cover (FVC), fraction of absorbed photosynthetically active radiation (FAPAR), canopy chlorophyll content (CC) and canopy water content (CW) over the actual S2 image. The derived values of CC and CW are then used to derive the leaf optical properties, using CC for the transmittance in the visible spectrum and CW for the transmittance in the infrared spectrum. These transmittances are used to estimate net radiation and its partitioning into canopy and soil using the Campbell and Norman model [51,62]. As an additional outcome from these outputs, the fraction of green vegetation (fg) is derived by combining information of LAI and FAPAR [45]. In order to avoid inconsistencies in estimation and cloud-related issues, the multi-step Savitzky-Golay filtering algorithm [63] was used to smooth the S2 LAI data, and the upper envelope of LAI time series was generated.

Sharpening Land Surface Temperature (LST)
The S3 SLSTR LST images were disaggregated using the DMS approach [43,45]. The sharpening objective is to downscale the 1 km LST imagery from the S3 to the S2 spatial scale (20 m). This methodology relies on the empirical relation that exists between shortwave and thermal data. Similar approaches describe some relations between these two regions, such as TsHARP, which relates the normalized difference vegetation index (NDVI) with thermal data [64] or even with other vegetation indexes (VIs). The main limitation to these previous studies is that the full multispectral information was not represented in VIs but only a small part of it, since they only use two or three bands of the sensors [45]. Often they assume linear relations and, in most cases, these relations are non-linear. The DMS technique instead uses a non-parametric approach that considers a full set of spectral bands and ancillary information (e.g., elevation, exposure, etc.). In addition, it uses a machine learning algorithm that permits non-linear relations between shortwave and thermal data and ensures the conservation of emitted thermal energy between the sharpened and the original LST data by a bias correction approach. DMS aggregates the high-resolution image to the thermal image coarse resolution, and afterwards builds a machine learning model using the most homogeneous coarse-resolution pixels in the scene. In that way, urban areas, roads and highly heterogeneous landscapes can be discarded in the training dataset [65].
The machine learning algorithm used in this study for the non-parametric relation was a decision tree as it is fast and at the same time reliable in comparison with other algorithms such as neural networks [44]. The pyDMS Python module (https://github.com/radosuav/ pyDMS last accessed on 20 October 2021) was used to apply this algorithm.

Meteorological Data
Meteorological data from the ERA5-Land reanalysis dataset of the Copernicus Climate Change Service (C3S) (2019) were used. These data are extremely useful for running energy balance models at regional and basin scale as they consist of hourly gridded information at 0.1 • of spatial resolution and contain several layers of meteorological data. For the current study, the variables used were air temperature at 2 m, dew point temperature at 2 m, wind speed at 10 m, surface pressure and total column water vapor (TCWV).

Ancillary Data
Vegetation height and leaf inclination angle cannot be retrieved from remote sensing passive data directly. The leaf angle distribution parameter indicates how leaves are oriented in the vertical plane. This is important when accounting for determining the beam shortwave net radiation partitioning [66]. For almond, it is assumed that leaf angle follows a spherical distribution [67] in which the leaves point in all directions equally. Canopy height was measured periodically per plot, ranging from 2.3 m in the stressed plots (T3 and T4) to 5 m in the most irrigated ones during the last year of the study (T1 and T2). This Remote Sens. 2022, 14, 2106 9 of 25 variable, together with LAI, has an impact on the aerodynamic roughness [68], which in the TSEB scheme determines the atmospheric stability and aerodynamic resistance used to estimate the sensible heat flux.

Daily ET Upscaling
Instantaneous energy fluxes at the satellite overpass time were estimated using TSEB with the previously described inputs. Then, the latent heat flux was integrated into the daily evaporation rate. To achieve this, the instantaneous value, in W/m 2 , was upscaled to daily water fluxes, expressed in units of mm/day, by multiplying the instantaneous ratio of latent heat flux over solar irradiance by the average daily solar irradiance [69].

Gap Filling
It is a well-known issue that clouds can partially or totally occlude satellite optical images. Due to the S2+S3 synergy methodology, it is possible to have occlusion either on S2 high-resolution images or on S3 coarse-resolution images. If there is an occlusion with S2 data, it will extend for the following 5-day window or until an unoccluded image is available. If the occlusion corresponds to clouds in an S3 image, the period of occlusion will be shorter, but each masked area will correspond to 1 km pixel.
To obtain a continuous dataset with daily ET a values, a crop stress coefficient (K cs ) was retrieved from the valid ET estimates by computing the ratio between actual and reference ET, with the latter computed from the meteorological gridded data.
This ratio (K cs ), which represents the combined effect of the FAO-56 crop and stress coefficients (K cs = K c * K s ), is used to fill the cloud gaps present in the following days, under the assumption that K cs will remain constant for a short-term period [46].

Eddy Covariance (EC) Fluxes
To validate the modelled energy fluxes, an eddy covariance (EC) system was set up at the study site, equipped with an open path infrared gas analyzer (LICOR-7500), a sonic anemometer (CSAT3), a four-component net radiometer (NR01 Campbell) and four soil heat flux plates (Hukseflux) located at 0.5, 1, 1.5 and 2 m from the tree line. A CR3000 Campbell datalogger was used to store the 20 Hz measurements as well as the half-hourly aggregated fluxes. The fluxes were corrected by the Webb, Pearman and Leuning (WPL) method described in Foken et al. [70] to correct density fluctuations in open path gas analyzers. Then, the 30-min period fluxes in which the S3 overpassed (~10:30 local time) were used. For solving the non-closure of the energy fluxes, one closure forcing option is to apply the Bowen ratio closure method, which assumes that the ratio between H and LE (β = H/LE) is correctly measured by the EC system so that individual values of these fluxes can be adjusted to balance Equation (1) [71,72].
To estimate the relative contribution of each pixel to the measured fluxes, it is mandatory to calculate the tower footprint. This footprint will depend on the wind direction and speed and the surface roughness, which is mainly driven by the vegetation structure in the surroundings. The footprints were estimated using a two-dimensional parametrization of the flux footprint predictions ( [73]; www.footprint.kljun.net last accessed on 19 September 2021). Then, the footprints were resampled to 20 m to match S2 resolution ( Figure 3). Finally, the latent heat flux and the sensible heat flux were aggregated as a weighted sum of the pixel contribution to be compared against the tower fluxes. However, the R n footprint from the net radiometer covers only an area encompassed by the sensor viewing angle, which in this case falls in the site where the eddy covariance is set (i.e., T1).
in the surroundings. The footprints were estimated using a two-dimensional parametr zation of the flux footprint predictions ( [73]; www.footprint.kljun.net last accessed on 1 September 2021). Then, the footprints were resampled to 20 m to match S2 resolution (Fi ure 3). Finally, the latent heat flux and the sensible heat flux were aggregated as weighted sum of the pixel contribution to be compared against the tower fluxes. Howeve the footprint from the net radiometer covers only an area encompassed by the sens viewing angle, which in this case falls in the site where the eddy covariance is set (i.e., T1 Figure 3. Eddy covariance tower footprint for 21st July 2020, representing the fetch and the relati contribution for a typical day, at high resolution (left) and resampled at 20 m to match Sentinel resolution (right). The overall spot represents an 80% contribution for both plots.

Model ETa Intercomparison through a Soil Water Balance Approach
Average evapotranspiration for the period between two soil water content (SWC measurements was estimated from the following equation: where WB is water balance, IR is irrigation, Peff is effective precipitation and ΔSWC is th difference of soil water content between two measurements. Since the soil has a high i filtration rate and null slope, runoff was assumed to be zero and, therefore, Peff was calc lated as described above. Then, in order to compare values of ETa between irrigation trea ments and with the TSEBS2+S3 approach, ETa WB was normalized by dividing the total E of a specific period by the number of days between two consecutive measurements. A neutron probe (Campbell Pacific Nuclear Scientific, Model 503) was used to mea ure SWC down to 150 to 180 cm in one plot in each irrigation treatment. In each plot, thr tubes were installed, one in the emitter-wetted area, a second in the middle of the lan and a third in an intermediate location between rows. Readings were taken at 20 cm i tervals every 15 to 20 days throughout the growing season. In total, the number of mea urements ranged from eight to twelve, depending on the growing season. The neutro probe was calibrated for the experimental soil by taking soil samples for volumetric moi ture content (ɵ, cm 3 of water/cm 3 of soil) when the access tubes were installed. . Eddy covariance tower footprint for 21st July 2020, representing the fetch and the relative contribution for a typical day, at high resolution (left) and resampled at 20 m to match Sentinel-2 resolution (right). The overall spot represents an 80% contribution for both plots.

Model ET a Intercomparison through a Soil Water Balance Approach
Average evapotranspiration for the period between two soil water content (SWC) measurements was estimated from the following equation: where WB is water balance, IR is irrigation, P eff is effective precipitation and ∆SWC is the difference of soil water content between two measurements. Since the soil has a high infiltration rate and null slope, runoff was assumed to be zero and, therefore, P eff was calculated as described above. Then, in order to compare values of ET a between irrigation treatments and with the TSEB S2+S3 approach, ET a WB was normalized by dividing the total ET a of a specific period by the number of days between two consecutive measurements. A neutron probe (Campbell Pacific Nuclear Scientific, Model 503) was used to measure SWC down to 150 to 180 cm in one plot in each irrigation treatment. In each plot, three tubes were installed, one in the emitter-wetted area, a second in the middle of the lane, and a third in an intermediate location between rows. Readings were taken at 20 cm intervals every 15 to 20 days throughout the growing season. In total, the number of measurements ranged from eight to twelve, depending on the growing season. The neutron probe was calibrated for the experimental soil by taking soil samples for volumetric moisture content ( September 2021). Then, the footprints were resampled to 20 m to match S2 resolution (Figure 3). Finally, the latent heat flux and the sensible heat flux were aggregated as a weighted sum of the pixel contribution to be compared against the tower fluxes. However, the footprint from the net radiometer covers only an area encompassed by the sensor viewing angle, which in this case falls in the site where the eddy covariance is set (i.e., T1). . Eddy covariance tower footprint for 21st July 2020, representing the fetch and the relative contribution for a typical day, at high resolution (left) and resampled at 20 m to match Sentinel-2 resolution (right). The overall spot represents an 80% contribution for both plots.

Model ETa Intercomparison through a Soil Water Balance Approach
Average evapotranspiration for the period between two soil water content (SWC) measurements was estimated from the following equation: where WB is water balance, IR is irrigation, Peff is effective precipitation and ΔSWC is the difference of soil water content between two measurements. Since the soil has a high infiltration rate and null slope, runoff was assumed to be zero and, therefore, Peff was calculated as described above. Then, in order to compare values of ETa between irrigation treatments and with the TSEBS2+S3 approach, ETa WB was normalized by dividing the total ETa of a specific period by the number of days between two consecutive measurements. A neutron probe (Campbell Pacific Nuclear Scientific, Model 503) was used to measure SWC down to 150 to 180 cm in one plot in each irrigation treatment. In each plot, three tubes were installed, one in the emitter-wetted area, a second in the middle of the lane, and a third in an intermediate location between rows. Readings were taken at 20 cm intervals every 15 to 20 days throughout the growing season. In total, the number of measurements ranged from eight to twelve, depending on the growing season. The neutron probe was calibrated for the experimental soil by taking soil samples for volumetric moisture content (ɵ, cm 3 of water/cm 3 of soil) when the access tubes were installed.
, cm 3 of water/cm 3 of soil) when the access tubes were installed. Both TSEB S2-S3 and the SWB approach were intercompared by computing their correlation as well as their root mean squared deviation (RMSD = ). RMSD is analogous to RMSE but it denotes deviations between two independent estimates (x and y), rather than errors from one estimate to a reference measurement.

Validation with the Eddy Covariance Flux Tower
The eddy covariance closure was analyzed in order to observe the reliability of measured fluxes. Figure 4 shows the closure at daytime hours (8:00-18:00) out of this time range. The value of the slope of the closure was 0.76, excluding night, dusk and dawn observations, which falls within the range of typical values found, for instance, in the FLUXNET database [74], of between 0.56 and 0.97 [75]. Usually, EC-derived fluxes bring inconsistencies, especially at dawn and dusk. There is better closure when the available energy is lower than 300 Wm −2 associated with morning hours, while at high irradiance, estimates of H+LE were up to 25% lower than Rn-G. The average Bowen ratio at the satellite overpasses was 0.74, being higher in July-August, when it rose above unity, due to the arid conditions of these months when temperatures can reach up to 40 • C or more.
observations, which falls within the range of typical values found, for instance, in the FLUXNET database [74], of between 0.56 and 0.97 [75]. Usually, EC-derived fluxes bring inconsistencies, especially at dawn and dusk. There is better closure when the available energy is lower than 300 Wm −2 associated with morning hours, while at high irradiance, estimates of H+LE were up to 25% lower than Rn-G. The average Bowen ratio at the satellite overpasses was 0.74, being higher in July-August, when it rose above unity, due to the arid conditions of these months when temperatures can reach up to 40 °C or more. The energy balance fluxes from TSEB using the S2+S3 method (TSEBS2+S3) were evaluated against eddy covariance measurements from 2018 to 2020 during the vegetative growing period (March-October) ( Figure 5, Table 1). Net radiation appears to be slightly underestimated by TSEBS2+S3 with an RMSE of 63 W/m 2 . Sensible heat flux has an RMSE of 90 W/m 2 , LE of 87 W/m 2 , and G an error of 37 W/m 2 . Turbulent fluxes are the most scattered ones. Rn shows a bias of −46 W/m 2 that corresponded to the −43 W/m 2 underestimation of sensible heat flux from TSEBS2+S3 over the EC. High values in LE measured by the EC tower are related to dates after a precipitation event where the tower measured high LE in contrast with TSEBS2+S3; this sharpening approach may not identify the sharp surface temperature decrease due to the larger evaporation from the soil. Soil heat flux appears to be overestimated by TSEBS2+S3 at lower fluxes and overestimated for higher net radiation days.
The energy fluxes estimated from TSEB using Landsat LST (TSEBLandsat) were also compared against the eddy covariance data. The error fluxes derived from TSEBLandsat (RMSE; Rn: 103 W/m 2 , H: 80 W/m 2 , LE: 108 W/m 2 , G: 43 W/m 2 ) did not differ significantly from TSEBS2+S3, showing the same pattern for net radiation. Even a higher error can be The energy balance fluxes from TSEB using the S2+S3 method (TSEB S2+S3 ) were evaluated against eddy covariance measurements from 2018 to 2020 during the vegetative growing period (March-October) ( Figure 5, Table 1). Net radiation appears to be slightly underestimated by TSEB S2+S3 with an RMSE of 63 W/m 2 . Sensible heat flux has an RMSE of 90 W/m 2 , LE of 87 W/m 2 , and G an error of 37 W/m 2 . Turbulent fluxes are the most scattered ones. Rn shows a bias of −46 W/m 2 that corresponded to the −43 W/m 2 underestimation of sensible heat flux from TSEB S2+S3 over the EC. High values in LE measured by the EC tower are related to dates after a precipitation event where the tower measured high LE in contrast with TSEB S2+S3 ; this sharpening approach may not identify the sharp surface temperature decrease due to the larger evaporation from the soil. Soil heat flux appears to be overestimated by TSEB S2+S3 at lower fluxes and overestimated for higher net radiation days.
Remote Sens. 2022, 14, 2106 12 observed, which is hard to determine given the fewer observations compared to TSE Observing the distribution pattern, it could be assumed that there is no significant provement in using a higher-resolution thermal sensor at this scale, where spatial he geneity is under the 100 m resolution, or the sharpening technique can identify these ferences by determining the actual relationships between stressed and non-stressed etation for the training area in the corresponding time window.   The energy fluxes estimated from TSEB using Landsat LST (TSEB Landsat ) were also compared against the eddy covariance data. The error fluxes derived from TSEB Landsat (RMSE; Rn: 103 W/m 2 , H: 80 W/m 2 , LE: 108 W/m 2 , G: 43 W/m 2 ) did not differ significantly from TSEB S2+S3 , showing the same pattern for net radiation. Even a higher error can be observed, which is hard to determine given the fewer observations compared to TSE S2+S3 . Observing the distribution pattern, it could be assumed that there is no significant improvement in using a higher-resolution thermal sensor at this scale, where spatial heterogeneity is under the 100 m resolution, or the sharpening technique can identify these differences by determining the actual relationships between stressed and non-stressed vegetation for the training area in the corresponding time window.

Water Applied and Biophysical Parameters of the Vegetation
Reference ET (ET 0 ) is a proxy for atmospheric demand, as it is usually modelled over well-irrigated grass, depending only on meteorological conditions. In the orchard location, ET 0 was higher in 2017 and 2019 by 100 and 90 mm, respectively, in comparison to 2018 and 2020. The two latter years corresponded to years with higher precipitation ( Table 2). The irrigation period started the second week of March for all treatments and years. The amount of water applied in 2017 was similar for all treatments, except for T4 which received 28% less water in comparison to other treatments. In 2018, differences between T4 and other treatments were also significant. Irrigation differences for T3 in comparison to other treatments can be observed in 2018, probably due to the high rainfall, but were not significant. In 2019 and 2020, the differences were significant between all treatments. The amount of water applied in T2, T3 and T4 in comparison to T1 was, respectively, 26%, 43% and 63% less in 2019, while in 2020 it corresponded to 10%, 25% and 55% less.
With respect to the biophysical parameters of the vegetation, the LAI is one of the most important when accounting for canopy size and vigor, and is used as input in TSEB. Figure 6 shows the time series of S2 LAI values after applying the Savitzky-Golay interpolation filter for each irrigation treatment during the four-year study. There is a clear increase in vegetative growth from the first year of study to 2020, increasing from averaged maximum LAI values of 0.75 in 2017 to 1.80 in 2020 (Figure 6a). A seasonal trend of LAI can also be observed representing crop development at the beginning of each growing season and leaf senescence at its end. A peak can be seen by the end of spring, when vegetative growth of the almond trees coincided with maximum vegetative development of the cover crop in the interrow. Differences between treatments were also noticeable in some years. In 2017, no significant differences were observed between treatments. In 2018, T4 had significantly smaller LAI values in comparison to other treatments (p < 0.01). In 2019, differences in LAI between treatments were more evident and significant differences were observed between all treatments (p < 0.01). In comparison to T1, the seasonal average LAI decreased 18%, 25% and 40% for T2, T3 and T4, respectively. In 2020, significant differences were observed between T1 and T3 and T4, and also between T4 and T2. In comparison to T1, the seasonal average LAI decreased 18% and 33% in T3 and T4, respectively. Although FVC is not used by TSEB, as vegetation gap fraction is modelled in TSEB from LAI and the clumping index [51], several studies have demonstrated high correlation with crop evapotranspiration [4,76]. Therefore, in this study, we also wanted to assess the seasonal pattern of FVC for different irrigation treatments. As was expected, results showed a similar pattern to that observed for LAI (Figure 6b). Sharpened LST at 20 m resolution was also compared between treatments (Figure 6c). Non-significant differences in LST were observed in the first two years of the experiment. However, LST of T1 and T4 showed significant differences, respectively showing the lowest and highest values in 2019 and 2020 in comparison to other treatments. These differences mostly occurred during the summer time, when there is a higher atmospheric demand and differences between well-watered and stressed treatments were more evident.   Figure 7 shows the seasonal evolution of the modelled ETa with TSEBS2+S3 for each irrigation treatment, using the gap-filling procedure. ETa increased sharply from early March to mid-June, when maximum values of ~4 to 5 mm/day were attained, depending on the year. After mid-July, ETa progressively dropped until the end of the season as ET0 and LAI also decreased. It can also be seen how ETa values drop for those days with precipitation events, due to lower solar irradiance under cloudy conditions (lower ET0), after which ETa rapidly ramps up in sunny days as soil moisture is replenished. As expected, Letters show Tukey's HSD analysis for years with significant differences (p ≤ 0.05). Lines show the mean and the shading region shows the error bounds for the four plots within each irrigation treatment. * means significant differences at p ≤ 0.05, while n.s. means not statistically significant differences at p ≤ 0.05. Figure 7 shows the seasonal evolution of the modelled ET a with TSEB S2+S3 for each irrigation treatment, using the gap-filling procedure. ET a increased sharply from early March to mid-June, when maximum values of~4 to 5 mm/day were attained, depending on the year. After mid-July, ET a progressively dropped until the end of the season as ET 0 and LAI also decreased. It can also be seen how ET a values drop for those days with precipitation events, due to lower solar irradiance under cloudy conditions (lower ET 0 ), after which ET a rapidly ramps up in sunny days as soil moisture is replenished. As expected, ET a increased year after year as vegetative growth increased, but also depending on seasonal ET 0 . On average, the highest ET a values were observed in 2020 due to the higher LAI values. In 2017, no significant differences between treatments were observed, except at the end of the growing season when ET a of T4 was significantly lower. This coincided with a decrease in the irrigation amount of water from mid-July until the end of the growing season. ET a values were all below 4 mm/day. In 2018, significant differences between T4 and other treatments were also observed from 15 July. A similar pattern was observed in 2019 and 2020. and other treatments were also observed from 15th July. A similar pattern was observed in 2019 and 2020. With respect to total cumulative ETa throughout the growing season, differences between treatments were significant from 2018 (Figure 8a). In 2018, a total of 682 mm was evaporated in T4, which was 9-14% less in comparison to other treatments. In that year, the sum of irrigation and rainfall water applied in T4 was also ~20% less than in the other treatments (Figure 8b). In 2019, differences in cumulative ETa between treatments were also significant. Despite the significant differences in the amount of water applied between all treatments, a significant decrease in ETa was only observed in T4, ranging between 20% and 13% less evaporated water. Although it was expected that the highest differences between treatments would be observed in 2020, these were not as evident as in 2019, probably due to a lower ET0 and higher number of precipitation events. Cumulative ETa for T4 was ~8% less in comparison to other treatments. Results also showed that T4 almond trees evaporated 8%, 10%, 22% and 36% more than the amount of water applied respectively for 2017 to 2020, while for T1 and T2 the total amount of water applied was higher than that consumed by the crop. Therefore, differences in ETa between treatments and years demonstrate the ability of this approach to identify periods of water stress and peaks over the growing season.  With respect to total cumulative ET a throughout the growing season, differences between treatments were significant from 2018 (Figure 8a). In 2018, a total of 682 mm was evaporated in T4, which was 9-14% less in comparison to other treatments. In that year, the sum of irrigation and rainfall water applied in T4 was also~20% less than in the other treatments (Figure 8b). In 2019, differences in cumulative ET a between treatments were also significant. Despite the significant differences in the amount of water applied between all treatments, a significant decrease in ET a was only observed in T4, ranging between 20% and 13% less evaporated water. Although it was expected that the highest differences between treatments would be observed in 2020, these were not as evident as in 2019, probably due to a lower ET 0 and higher number of precipitation events. Cumulative ET a for T4 was~8% less in comparison to other treatments. Results also showed that T4 almond trees evaporated 8%, 10%, 22% and 36% more than the amount of water applied respectively for 2017 to 2020, while for T1 and T2 the total amount of water applied was higher than that consumed by the crop. Therefore, differences in ET a between treatments and years demonstrate the ability of this approach to identify periods of water stress and peaks over the growing season.

Comparison of ETa Estimates with TSEBS2+S3 and a Soil Water Balance Approach
While the eddy covariance footprint represents an integrated measurement of the energy fluxes from all four treatments, a way to assess the ETa estimates with TSEBS2+S3 for each individual irrigation treatment, ranging from stressed to well-watered conditions, is through an intercomparison with a soil water balance model. This was obtained for weekly or bi-weekly periods, on the basis of soil water content reads with a neutron probe.
A scatter plot encompassing all ETa data obtained through the water balance and estimated with the TSEBS2+S3 approaches showed an average seasonal root mean square deviation (RMSD) and bias of 1.54 and −0.02 mm/day, respectively (Figure 9a). Results did not show significant differences between irrigation treatments. Although the joint data analysis does not show good results, when this was analyzed individually for each irrigation treatment and year, the time series of both methodologies followed the same trend throughout the growing seasons ( Figure 10). The largest sources of error came from longer periods between neutron probe measurements and due to precipitation events. The Letters show Tukey's HSD analysis for treatments with significant differences (p ≤ 0.05). Lines show the mean and the shading region shows the error bounds for the four plots within each irrigation treatment. * means significant differences at p ≤ 0.05, while n.s. means not statistically significant differences at p ≤ 0.05.

Comparison of ET a Estimates with TSEB S2+S3 and a Soil Water Balance Approach
While the eddy covariance footprint represents an integrated measurement of the energy fluxes from all four treatments, a way to assess the ET a estimates with TSEB S2+S3 for each individual irrigation treatment, ranging from stressed to well-watered conditions, is through an intercomparison with a soil water balance model. This was obtained for weekly or bi-weekly periods, on the basis of soil water content reads with a neutron probe.
A scatter plot encompassing all ET a data obtained through the water balance and estimated with the TSEB S2+S3 approaches showed an average seasonal root mean square deviation (RMSD) and bias of 1.54 and −0.02 mm/day, respectively (Figure 9a). Results did not show significant differences between irrigation treatments. Although the joint data analysis does not show good results, when this was analyzed individually for each irrigation treatment and year, the time series of both methodologies followed the same trend throughout the growing seasons ( Figure 10). The largest sources of error came from longer periods between neutron probe measurements and due to precipitation events. The RMSD ranged from 0.6 to 2.5 mm/day, but when the periods with precipitation events corresponding to 15 May 2019 and 18 May 2020 were not considered, maximum RMSD dropped to 2.0 mm/day. Similarly, average seasonal RMSD decreased to 1.19 mm/day when removing these dates (Figure 9b). In both approaches, a clear increase can be observed in ET a from 2017 to 2020. RMSD and bias were higher in 2020, when seasonal cumulative ET a reached its maximum values. For that year, ET a TSEB S2+S3 values for T3 and T4 were significantly higher than those obtained with ET a WB, suggesting a slight overestimation in periods of severe water stress. However, this cannot be confirmed in this study since T4 trees were also stressed in 2018 and 2019 and this overestimation was not observed.

Validation of Energy Fluxes
The imbalance between Rn-G and LE+H represents the non-closure of the surface energy balance which might be attributable to differences between the sampling frequencies or the footprint measurements. Moreover, other energy sources or sinks are not considered in these energy fluxes, such as canopy storage [77]. In another context, flux uncertainties may indicate an underperformance of TSEB since typical low errors for sensible heat flux (H) usually fluctuate around 50 W/m 2 and the ones obtained in this study are slightly higher. However, the footprint of the present study had a high degree of heterogeneity due to differences in almond water status and vegetative growth. This makes these results harder to compare with those from other studies in which fluxes are usually measured in well-irrigated homogeneous fields. Despite this, Guzinsky et al. [44] reported an error in H of 94 W/m 2 in woody crops using the TSEB-PT approach. On the other hand, He et al. [34] evaluated the METRIC modelling approach using Landsat in an almond orchard and obtained an RMSE of approximately 75 W/m 2 for LE. These results were not far from the RMSE of 90 and 80 W/m 2 for H and of 87 and 108 W/m 2 for LE observed in our study, respectively, for TSEB S2+S3 and TSEB Landsat .
The eddy covariance footprint may also vary from day to day. Although the main wind direction in the orchard was NW, it can vary throughout the day. In part, this could affect the measured fluxes according to the relative contribution of each treatment. On the other hand, the bias observed in Rn (−46 W/m 2 ) may be related to the fact that the net radiometer is placed over a fully irrigated tree of T1, which means that it is above the highest and densest canopy in the orchard. Therefore, this bias is actually related to the canopy-intercepted radiation and therefore can lead to higher observed Rn values than those modelled with TSEB but as well in larger EC closure errors. It should also be noted that the inputs for the net radiation calculation are obtained from a 20 × 20 m pixel which integrates both bare soil and canopy. Observed overestimations in H may be related to horizontal advection of hot, dry air. During these days, temperatures were over 35 • C, and the surrounding bare soil can heat up this air, provoking such conditions. These advective conditions can significantly increase the localized evaporative demand over an irrigated field, deriving larger latent heat fluxes. In such cases, the additional energy required for evaporation is extracted from the surface itself, resulting in negative sensible heat, increased closure errors [78] and more difficulties of TSEB-PT in accurately estimating sensible heat flux [79].
Orchards often present highly clumped canopies and therefore a significant area of the ground has no or very little leaf biomass. This causes distinct wind attenuation towards the soil surface, with larger attenuation within the canopy and lower in-ground exposed areas. This issue alters the momentum, heat and vapor exchange in the soil-plant-atmosphere continuum [80]. It has been proven that using wind profile models adapted to vertically heterogeneous canopies can improve the accuracy of surface energy flux estimations [81,82] compared to other profiles designed for homogeneous canopies [83]. However, these wind attenuation profile models require additional parameters related to the canopy structure that may not always be available, in particular for operational and large-scale applications such as regional water accounting with satellite imagery [46].
Landsat-8 was also used to analyze how the sharpening technique applied to Sentinel imagery can influence the flux estimations under water stress conditions. However, results from Figure 5 showed that flux from the sharpened S2+S3 LST were not significantly different from those of the lone Landsat 8 LST. An even higher-resolution sensor may be required for this study site, mainly because the 100 m resolution TIR on board Landsat-8 might not be enough to acquire representative LST pixels per plot. In fact, data from the Ecosystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS), using the International Space Station (ISS) [41], were also evaluated for the present study, but the lack of valid observations made this comparison unviable and it was therefore not presented in this article.

Land Surface Temperature (LST) and Biophysical Parameters of the Vegetation
The DMS approach used in this study to sharpen low-resolution LST with highresolution shortwave data is based on the assumption that there exists a statistical relation between the two spectral regions. The decision tree machine learning algorithm used allows a high level of complexity in non-linear relations between the shortwave reflectance bands and the LST [45]. However, there is some controversy about the accuracy of high-resolution LST estimates using this approach in highly variable landscapes, where there is an interaction between different land uses (e.g., rainfed and irrigated fields), an agricultural-to-urban landscape continuum, or in cases where crops are stressed but there is no reduction in biomass [47]. Until the launch of higher-resolution LST missions, the use of data sharpening methods is necessary. Some other machine learning ensemble algorithms exist which have also demonstrated high accuracy in downscaling LST and, at the same time, less overfitting. For instance, a decision-treebased ensemble machine learning algorithm that uses a gradient boosting framework (XGBoost) has recently become popular mainly due to its high accuracy and ease of implementation [84,85]. However, the computation time can be high when scaling it to a bigger dataset or a high number of variables. Other approaches based on use of a light gradient boosting machine (LightGBM) [86] have reported similar accuracy to that of XGBoost but with lower computation time. Adoption of such approaches could probably improve the accuracy of LST estimates. Although there were no in situ measurements to validate actual sharpened LST, a spatiotemporal pattern in LST was clearly observed between treatments ( Figures 6 and 11), with T4 showing significantly higher LST values during 2019 and 2020. Since sharpened LST was obtained from machine learning algorithms which relate S3 LST with S2 shortwave bands at coarse resolutions, it is probable that differences in sharpened LST between irrigation treatments can mainly be attributed to differences in vegetative growth variables. This has previously been observed in grapevines by Bellvert et al. [47], who reported that sharpened LST for severely stressed vines tends to be underestimated, especially if there is no evidence of a decrease in biomass.
Radiative transfer models for retrieving biophysical parameters are already being operationally produced and integrated into models for decision making [60]. Among them, LAI is one of the main drivers needed in TSEB. Burchard-Levine [87] demonstrated that a change of 0.1 in LAI is associated with a 2.9% change in simulated H with TSEB. However, the retrieval of accurate LAI and canopy structure estimates, especially in heterogeneous row crops, continues to be a major challenge. Nonetheless, the current study was able to detect spatiotemporal differences in LAI and FVC between irrigation treatments (Figures 6 and 11). model and clearly fits already existing datasets. As most of the current available st were conducted in mature almond trees with an FVC > 60%, the combination of the able datasets with the results obtained in the current study in young almond trees all completion of the regression. It is also interesting to see that despite the fact that ceived less water and transpired significantly less during the last three years of s points followed the same trend as other treatments. Therefore, the lower ETa observ T4 was basically explained by trees being smaller. Values with symbols corresponded to well-watered treatments (T1 and T2), and in red circ T4. x corresponds to the following references: [4,5,[7][8][9]20,34].
Finally, a spatial pattern was also observed for the biophysical parameters LS ETa, with differences observable between irrigation treatments. On average, maxi differences in FVC and LAI between treatments were observed during the last two of the experiment (Figure 12a). The highest spatial variability in ETa was also observ 2019 and 2020. Figure 12b shows in more detail the spatial variability of LST, LA daily ETa for three representative dates along the growing season. This figure also r sents the spatial resolution with respect to the plots. It can be observed that only o two pixels fall entirely into each treatment. LST at the beginning of the season was ar 295K in contrast with mid-July LST where the soil can be up to 320K or even high bare areas. LAI differences were also lower at the beginning of the season and be higher at mid-summer and at the end of the season. ETa peak was at mid-summer reaching around 5 mm/day in T1. Values with symbols corresponded to well-watered treatments (T1 and T2), and in red circle was T4.
It is worth noting that several assumptions are taken into account when predicting the biophysical parameters of vegetation from Sentinel-2, most importantly that the radiative transfer model that is used assumes horizontally and vertically homogeneous vegetation. Therefore, in orchards and row crops or in transition zones between different vegetation patches, inversion of this RTM may yield larger uncertainty. Moreover, any retrieval algorithm is susceptible to satellite signal uncertainty due to geometric misregistration, anisotropic reflectance effects, electronic errors, artifacts due to data resampling, atmosphere and clouds [88]. In those cases, noisy retrievals can be observed, such as unexpected peaks inherent to the vegetation phenology or other external factors such as pruning or other agronomic practices that may affect vegetation vigor. For this reason, in the present study, the biophysical parameters were filtered using a Savitzky-Golay filter, which smoothed the temporal series without losing the actual information [89]. Another assumption to be considered is that biophysical retrievals may be affected by cover crop in the interrow, especially at the beginning of each growing season. In this study, an increase in LAI was observed at the beginning of each growing season (until late June), probably attributable to grass and weeds growing in the interrow. This effect cannot be separated from a single-view satellite observation (for the S2 tile T31TBF, the VZA is 2.28 degrees in the study site). Three-dimensional RTM models, such as FLIGHT [90] or DART-3D [91], where the simulated scenes can include heterogeneous canopy structures based on plant architecture and other aspects such as weeds in the interrow, may be required for retrieving biophysical parameters in sparse woody crops. It would be interesting to explore the use of these models in future works since more efficient solutions have recently been implemented. For example, DART-3D applied an unbiased bidirectional path tracing method called DART-Lux, which can be a hundred times more efficient than the previous version [92]. However, the greater the complexity of the model, the lower the suitability of its inversion, as it could lead to many likely solutions even when different inputs are used.

Seasonal Actual Evapotranspiration (ET a )
For many years, different studies have shown the feasibility of using different SEB models to estimate energy fluxes in different crops, using data from different satellites. In general, satellite-derived energy fluxes have been validated with eddy covariance or surface renewal instrumentation [8,46,78,93]. Common to most of these studies is that fluxes are assessed temporarily during one or several growing seasons, and the spatial variability in terms of the crop's vegetative growth or water status of energy fluxes within the footprint is low. Moreover, from our knowledge, no study has assessed ET a estimates through surface energy models such as the TSEB S2+S3 in heterogeneous row crops irrigated at different levels, thus forcing eventual differences in water status.
This study demonstrates the usefulness of the TSEB S2+S3 approach to assess ET a in an almond orchard with a high degree of spatial variability in terms of water status. The seasonal average RMSD ranged from 1.4 to 1.6 mm/day and the bias from −0.53 to 0.53 mm/day in comparison to the water balance approach, showing similar values for all irrigation treatments.
The maximum cumulative ET a values obtained in this study for each specific year ranged from 760 to 915 mm, corresponding with a maximum FVC of 40% and 61%, respectively. These seasonal values were comparable with other studies also conducted in a semi-arid Mediterranean climate. For instance, Lopez-Lopez et al. [9] reported an annual ET a ranging from 923 to 1220 mm in fully irrigated almond trees with an FVC of 58 and 75%, respectively. The same study reported an average seasonal ET a of 699 mm for stressed almonds, which on average corresponded to 36% less in comparison to the fully irrigated ones. In our study, maximum differences in ET a between stressed (T4) and well-watered almonds (T1) were observed in 2019 and were as high as 20%. Similarly, Stevens et al. [20] observed a cumulative ET a of 1257 mm in almond trees with an FVC of 65%. In this case, however, the seasonal ET 0 (ET 0 = 1365 mm) and the total sums of irrigation and precipitation (I + P = 1409 mm) were higher than in our study. Goldhamer and Fereres [7] also estimated similar seasonal ET a values in almond trees with an average FVC of 80% and under ten irrigation levels, which ranged from 1100 to 1350 mm (ET 0 = 1170 mm). In these studies, calculations of ET a were obtained through a water balance approach, where soil water content was either monitored with neutron probe measurements or time domain reflectometry (TDR) sensors. Since neutron probe measurements are scarce in space and time, mainly because measurements are performed manually and a certified person is needed, measurements cannot be carried out periodically. Likewise, installing TDR sensors may not be an appropriate methodology due to the large number of sensors required to characterize total soil water content within tree spacing distance. For both reasons, commonly used soil water balance approaches may have some error in accounting for ET a . Studies such as that of Bellvert et al. [8], which assessed seasonal ET a on the basis of eddy covariance or surface renewal instrumentation, also reported annual water consumption of 1194.8 mm for a dense (80% FVC) high canopy of mature almonds. In Figure 11, we plotted average FVC with cumulated ET a and crop coefficients (Kc) only for those treatments that were transpiring at potential rates (T1 and T2) and were compared with those of T4 as well as with data from other available studies. Since the almond trees of our study were planted in 2016, the maximum average FVC values were only 50% in 2020, which corresponded with a seasonal ET a of 900 mm. The regression follows an exponential model and clearly fits already existing datasets. As most of the current available studies were conducted in mature almond trees with an FVC > 60%, the combination of the available datasets with the results obtained in the current study in young almond trees allowed completion of the regression. It is also interesting to see that despite the fact that T4 received less water and transpired significantly less during the last three years of study, points followed the same trend as other treatments. Therefore, the lower ET a observed in T4 was basically explained by trees being smaller.
Finally, a spatial pattern was also observed for the biophysical parameters LST and ET a , with differences observable between irrigation treatments. On average, maximum differences in FVC and LAI between treatments were observed during the last two years of the experiment (Figure 12a). The highest spatial variability in ET a was also observed in 2019 and 2020. Figure 12b shows in more detail the spatial variability of LST, LAI and daily ET a for three representative dates along the growing season. This figure also represents the spatial resolution with respect to the plots. It can be observed that only one or two pixels fall entirely into each treatment. LST at the beginning of the season was around 295K in contrast with mid-July LST where the soil can be up to 320K or even higher in bare areas. LAI differences were also lower at the beginning of the season and became higher at mid-summer and at the end of the season. ET a peak was at mid-summer 2020, reaching around 5 mm/day in T1.

Conclusions
The Copernicus program provides a wide range of services and datasets with capabilities that can be used for several operational agricultural applications. For irrigation management, however, a high-resolution TIR sensor is still a missing feature for operational estimates of crop evapotranspiration. Consequently, while waiting for the launch of higher-resolution LST missions, the use of modern data fusion techniques is required to compensate for the current lack of high-resolution TIR observations. This study demonstrates that use of the TSEB model with sharpened Sentinel-2 and Sentinel-3 imagery is suitable to monitor water use in almond trees irrigated under different water regimes. Sensible heat flux showed an RMSD of 90 W/m 2 , LE of 87 W/m 2 , and G an error of 37 W/m 2 . The assessed Sentinel-2 biophysical parameters such as LAI and FVC were able to distinguish differences in vegetative growth between irrigation treatments. Moreover, sharpened LST showed significant differences between the most extreme irrigation treatments (T1 and T4) during the last two years of the study period. The comparison of ETa with the soil water balance model showed an RMSD ranging from 0.6 to 2.5 mm/day. The maximum differences in cumulative ETa between irrigation treatments were obtained in 2019, amounting to between 20% and 13% less when comparing T1 and T4.

Conclusions
The Copernicus program provides a wide range of services and datasets with capabilities that can be used for several operational agricultural applications. For irrigation management, however, a high-resolution TIR sensor is still a missing feature for operational estimates of crop evapotranspiration. Consequently, while waiting for the launch of higher-resolution LST missions, the use of modern data fusion techniques is required to compensate for the current lack of high-resolution TIR observations. This study demonstrates that use of the TSEB model with sharpened Sentinel-2 and Sentinel-3 imagery is suitable to monitor water use in almond trees irrigated under different water regimes. Sensible heat flux showed an RMSD of 90 W/m 2 , LE of 87 W/m 2 , and G an error of 37 W/m 2 . The assessed Sentinel-2 biophysical parameters such as LAI and FVC were able to distinguish differences in vegetative growth between irrigation treatments. Moreover, sharpened LST showed significant differences between the most extreme irrigation treatments (T1 and T4) during the last two years of the study period. The comparison of ET a with the soil water balance model showed an RMSD ranging from 0.6 to 2.5 mm/day. The maximum differences in cumulative ET a between irrigation treatments were obtained in 2019, amounting to between 20% and 13% less when comparing T1 and T4.