Evaluating the Two-Source Energy Balance Model Using MODIS Data for Estimating Evapotranspiration Time Series on a Regional Scale

Citation: Bozorgi, M.; Cristóbal, J.; Pàmies-Sans, M. Evaluating the Two-Source Energy Balance Model Using MODIS Data for Estimating Evapotranspiration Time Series on a Regional Scale. Remote Sens. 2024,16, 4587. https://doi.org/10.3390/ rs16234587 Academic Editor: Yaoming Ma Received: 22 October 2024 Revised: 30 November 2024 Accepted: 4 December 2024 Published: 6 December 2024 Copyright: © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Article Evaluating the Two-Source Energy Balance Model Using MODIS Data for Estimating Evapotranspiration Time Series on a Regional Scale Mahsa Bozorgi 1,2,*, Jordi Cristóbal 1and MagíPàmies-Sans 2 1Department of Geography, Universitat Autònoma de Barcelona, Campus de Bellaterra, Edifici B, Carrer de la Fortuna, s/n, 08193 Bellaterra, Spain; jor[email protected] 2Efficient Use of Water in Agriculture Program, Institute of Agrifood Research and Technology, Fruitcentre, Parc Científic i Tecnològic Agroalimentari de Lleida 23, 25003 Lleida, Spain; [email protected] *Correspondence: mahsa.bozor[email protected] Abstract: Estimating daily continuous evapotranspiration (ET) can significantly enhance the monitoring of crop stress and drought on regional scales, as well as benefit the design of agricultural drought early warning systems. However, there is a need to verify the models’ performance in estimating the spatiotemporal continuity of long-term daily evapotranspiration (ET d ) on regional scales due to uncertainties in satellite measurements. In this study, a thermal-based two-surface energy balance (TSEB) model was used concurrently with Terra/Aqua MODIS data and the ERA5 atmospheric reanalysis dataset to calculate the surface energy balance of the soil–canopy–atmosphere continuum and estimate ET at a 1 km spatial resolution from 2000 to 2022. The performance of the model was evaluated using 11 eddy covariance flux towers in various land cover types (i.e., savannas, woody savannas, croplands, evergreen broadleaf forests, and open shrublands), correcting for the energy balance closure (EBC). The Bowen ratio (BR) and residual (RES) methods were used for enforcing the EBC in the EC observations. The modeled ET was evaluated against unclosed ET and closed ET (ET BR and ET RES ) under clear-sky and all-sky observations as well as gap-filled data. The results showed that the modeled ET presented a better agreement with closed ET compared to unclosed ET in both Terra and Aqua datasets. Additionally, although the model overestimated ET d across all different land cover types, it successfully captured the spatiotemporal variability in ET. After the gap-filling, the total number of days compared with flux measurements increased substantially, from 13,761 to 19,265 for Terra and from 13,329 to 19,265 for Aqua. The overall mean results including clear-sky and all-sky observations as well as gap-filled data with the Aqua dataset showed the lowest errors with ET RES , by a mean bias error (MBE) of 0.96 mm.day −1 , an average mean root square (RMSE) of 1.47 mm.day −1 , and a correlation (r) value of 0.51. The equivalent figures for Terra were about 1.06 mm.day −1 , 1.60 mm.day −1 , and 0.52. Additionally, the result from the gap-filling model indicated small changes compared with the all-sky observations, which demonstrated that the modeling framework remained robust, even with the expanded days. Hence, the presented modeling framework can serve as a pathway for estimating daily remote sensing-based ET on regional scales. Furthermore, in terms of temporal trends, the intra-annual and inter-annual variability in ET can be used as indicators for monitoring crop stress and drought. Keywords: remote sensing-based evapotranspiration; two-source energy balance; MODIS; drought 1. Introduction According to the Intergovernmental Panel on Climate Change Assessment Report (IPCC), European Mediterranean areas have experienced an increase in the frequency of concurrent droughts and heatwaves since the 1950s [ 1 ], associated with rising extreme temperatures and decreasing precipitation [ 2 ]. The Mediterranean climate is characterized by water-limited crop production; therefore, mapping the long-term spatiotemporal variability Remote Sens. 2024,16, 4587. https://doi.org/10.3390/rs16234587 https://www.mdpi.com/journal/remotesensing Remote Sens. 2024,16, 4587 2 of 29 in evapotranspiration (ET) is critical for comprehensively understanding how ecosystems respond to climate changes. ET plays a critical role in the climate system by exchanging energy and mass between the land surface and atmosphere, which sequentially affects plant growth and yield [ 3 – 5 ]. The estimation of daily evapotranspiration (ET d ) is critical for achieving efficient water use in irrigated agricultural areas [ 6 ] by determining crop water requirements [ 7 ], scheduling irrigation, and optimizing water use, especially in water-limited regions [8]. Additionally, the spatiotemporal variability in ET d can reflect ecosystem responses to hazards such as droughts [9] and climate changes [10], aiding in the design of early warning systems [11]. ET can be estimated with ground-based and remote sensing-based observations. Conventional ground observation methods, including lysimeters [ 12 ], the Bowen ratio [ 13 ], and eddy covariance systems [ 14 ], rely on direct ground measurements. Although these pointscale measurements have shown relatively accurate ET measurements for homogeneous areas, they have a limited ability to quantify the spatial variability in ET on regional scales due to spatial heterogeneity. Remote sensing data have revolutionized water resource management by providing continuous spatial and temporal information, even from remote and inaccessible regions, at various resolutions, which can address the limitations of ground observations [15–18]. There are various methods for estimating ET using remotely sensed land surface temperature (LST) [ 19 ]. Thermal-based surface energy balance (SEB) models are widely used to map ET by estimating surface energy fluxes [ 20 – 24 ]. SEB models consider the latent heat flux (LE), which represents the energy used for soil evaporation and plant canopy transpiration, as the differences between available energy (i.e., net radiation (R n ) and soil heat flux (G)) and sensible heat flux (H) [ 25 ]. However, there are challenges associated with satellite retrievals of LST [ 26 ]. Remote sensing-based LST integrates thermal data from both soil and vegetation components. Dual-source models have been introduced to overcome this limitation by petitioning radiometric temperature between soil and vegetation. These models differentiate soil and vegetation energy fluxes when estimating H and LE, showing better performance on vegetation–soil-mixed surfaces compared to signal-source models that assume the land surface as a signal system [ 27 ]. An example of a dual-source model is the TSEB model developed by [ 28 , 29 ], which has been successfully applied in various environments and land covers worldwide. The TSEB model minimizes uncertainties related to user experience, similar to the Surface Energy Balance Algorithm for Land (SEBAL) [ 30 ] and Mapping Evapotranspiration at high Resolution with Internalized Calibration (METRIC) [ 31 ]. When comparing SEBAL with TSEB, [ 32 ] found that TSEB outperformed SEBAL in areas with bare soil and sparse vegetation. Modifications to SEBAL inputs improved accuracy for specific land cover types but exacerbated discrepancies for others, further highlighting the model’s limitations [ 32 ]. Additionally, [ 33 ] found that TSEB performed better compared to METRIC in estimating instantaneous ET. However, the accuracy of the TSEB model is determined by the quality of the input data. Moreover, the parameterization of vegetation in flux estimation is one of the challenges in the application of the TSEB model, which can significantly influence the outputs [ 34 , 35 ]. The Atmosphere–Land Exchange Inverse (ALEXI) model introduced by [ 36 ] applied TSEB at two times to improve soil and vegetation components. The ALEXI disaggregation (DisALEXI) approach [ 33 ] aimed to improve the spatial resolution of the coarse temporal scale of the ALEXI model. The SEN-ET modeling framework “https://www.esa-sen4et.org (accessed on 19 March 2023)” [ 37 ] is another recent development, with significant improvements in the implementation of the TSEB model using the European Union’s Copernicus Earth Observation data (Sentinel-3 and -2) and meteorological reanalysis from ERA5 to estimate land surface fluxes, including ETd, at a 20 m resolution. Numerous ET products have been developed based on remote sensing observations, SEB models, and machine learning, each offering different spatial and temporal resolutions [ 38 ]. Daily ET products such as GLEAM [ 39 ] and SSEBop [ 40 ] are available, at coarse spatial resolutions. However, finer-scale ET mapping is required to enhance water re- Remote Sens. 2024,16, 4587 3 of 29 source management and drought monitoring [ 7 ]. Conversely, products like MOD16 [ 41 ], GLASS [ 42 ], and BESS provide high spatial resolution (1 km) but only at 8-day intervals, limiting their ability to monitor the sub-seasonal (weekly) variability in ET. Consequently, although numerous ET products are available on multiple spatial and temporal scales, there is still a lack of spatiotemporal continuous estimation of long-term ET d at a 1 km spatial resolution across southern Europe. Furthermore, many published studies only evaluated the performance of ET models over the crop growth period of a single year or short periods [ 24 , 43 ]. However, ET estimation approaches should be capable of capturing ET spatiotemporal year-round variability. Hence, this study aims to address these gaps by (1) estimating the ET d time series at two satellite overpasses using Terra and Aqua MODIS images at a 1 km spatial resolution from 2000 to 2022, and (2) comparing the performance of the TSEB model with long-term ground-based measurements under various land covers and environmental conditions to evaluate year-round model performance. 2. Materials and Methods 2.1. Study Area The study area covers the entire of the Iberina Peninsula (IP) and the southern part of France, including regions with both Mediterranean and coastal climates. This area is among the most vulnerable regions in southern Europe due to its heterogeneous atmospheric influences [ 44 ], Mediterranean climate, and agricultural predominance [ 45 ]. The diverse atmospheric patterns in this area lead to significant variability in precipitation, contributing to recurrent droughts and a notable trend toward increasing dryness over the last few decades [ 46 – 48 ]. Figure 1shows the study area, and the 11 selected flux towers used to evaluate the model. Figure 1. Location of the study area and the 11 selected flux towers. Projection system in UTM-30N WGS-84. 2.2. Data 2.2.1. Remote Sensing and Reanalysis Dataset TSEB Inputs The remote sensing and climate reanalysis data (ERA5) used to implement the TSEB model in this study are described in Table 1. Remote Sens. 2024,16, 4587 4 of 29 Table 1. The input data used in the TSEB model. Data Products Variables Resolution Temporal Spatial MODIS/Terra MOD09G Surface reflectance, NDVI Daily 500 m MOD11A1 LST, emissivity Daily 1000 m MODIS/Aqua MYOD09 Surface reflectance, NDVI Daily 500 m MYD11A1 LST, emissivity Daily 1000 m MODIS/Combined MCD15A3H LAI and FPAR 4-Day 500 m MCD19A2 Aerosol optical depth (AOD) Daily 1000 m MCD12Q1 Land cover Yearly 500 m ERA5 ECMWF Air temperature (Ta) Incoming surface radiation Total columnar water vapor (TCWV) Atmospheric pressure (P) Wind speed at 10 m height (u) Daily 27,830 m SRTM Elevation Yearly 90 m 2.2.2. Flux Towers Eddycovariance(EC)fluxmeasurementsfromtheFLUXNET2015“http://fluxnet.fluxdata. org/data/fluxnet2015-dataset (accessed on 27 January 2022)” dataset and ICOS (Warm Winter 2020 ecosystem eddy covariance flux data in FLUXNET-Archive format—release 2022-1 (Version 1.0), “https://doi.org/10.18160/2G60-ZHAK (accessed on 1 January 2022)” were used to evaluate the model. From those datasets, 11 flux stations covering a wide altitudinal range were selected. The flux tower data were selected based on having required variables, including daily flux components (Rn, LE, H, and G) and precipitation. Daily LE was converted to units of mass flux (mm·day−1) for the ETdmodel comparison. The location of selected flux towers is shown in Figure 1, and additional details about each flux tower are presented in Table 2. Table 2. Detailed information of selected flux towers. The height of the temperature measurement (z_T). The height of the wind measurement (z_U). The International Geosphere–Biosphere Programme (IGBP) defines landcover type. Savanna (SAV), Woody savanna (WSA), cropland (CRO), evergreen broadleaf forest (EBF), open shrubland (OSH). Flux ID Latitude Longitude Altitude z_T z_U IGBF Study Period ES_Abr 38.70 −6.79 279 2/12 12 SAV 2015–2020 ES_LM1 39.94 −5.78 266 2/15 15 SAV 2014–2020 ES_LM2 39.93 −5.78 270 2/15 15 SAV 2014–2020 ES_Cnd 37.91 −3.23 366 6 9.3 WSA 2014–2020 FR_Aur 43.55 1.11 250 2 2 CRO 2005–2020 FR_Lam 43.50 1.24 181 2 2 CRO 2005–2020 FR_Pue 43.74 3.60 270 10 10 EBF 2000–2014 ES_LgS 37.01 −2.97 2267 1.2 2.8 OSH 2007–2009 ES_Agu 36.94 −2.03 195 3.28 3.28 OSH 2006–2020 ES_LJu 36.93 −2.75 1600 2.5 2.5 OSH 2004–2020 ES_Amo 36.83 −2.25 58 2.9 2.9 OSH 2007–2012 To evaluate the modeled ET d , flux tower data were screened according to a daily energy closure >0.65 and with no precipitation present. The energy balance closure (EBC) was used to evaluate the reliability of the EC measurements at flux towers [ 46 , 49 ].The EBC formulation is based on the first law of thermodynamics, representing that the sum of turbulent fluxes (LE, H) is equal to the available energy (Rn, G, and heat storage (S)), as follows: EBC =∑(H+LE) ∑(Rn−G−S)(1) Remote Sens. 2024,16, 4587 5 of 29 Therefore, the ideal EBC is close to 1 and the lack of closure increases uncertainty in ET measurements [50,51]. Numerous studies have shown imbalances in the ratio of turbulent fluxes to available energy in the EC technique [ 49 , 52 ]. The Bowen ratio method [ 53 ] and residual method are two standard methods used for enforcing the EBC in the EC observations by attributing available energy to turbulent fluxes. In the current study, a distribution of the residual according to the Bowen ratio method (LE BR ) was applied as suggested by [ 53 , 54 ]. In addition, LE was recalculated as the residual (LE RES ) of the surface energy budget used in previous TSEB evaluations [ 55 ]; and both closure methods were then used to evaluate ET d . On clear-sky days, the maximum shortwave radiation peak occurs at local solar noon and reduces to zero at night. Hence, clear-sky observations were identified by calculating the ratio of incoming shortwave radiation to potential incoming shortwave radiation. Days with a ratio exceeding 70% were visually inspected, and those meeting this criterion were selected as clear-sky observations [56]. 2.3. Overall Modeling Framework A modeling framework for estimating ET d based on the SEN-ET model (detailed description can be found at “https://www.esa-sen4et.org (accessed on 19 March 2023)”, calculating the SEB through the TSEB model [ 28 ] (see Section 2.4 “TSEB model overview”), was refined for Terra/Aqua MODIS data (Figure 2). Figure 2. Flowchart of the modeling framework estimating ET d time series. The orange rectangle denotes the pre-processing of MODIS vegetation indices through TIMESAT. The SEN-ET model was designed to use Sentinel-2 and Sentinel-3 data as well as ERA5 reanalysis data as TSEB model inputs. To refine this methodology for Terra/Aqua MODIS data, several MODIS products (see Table 1) were used to obtain biophysical variables, vegetation properties, and surface characteristics. TIMESAT 3.3 [ 53 ] was employed to temporally smooth and spatially complete the MODIS normalized difference vegetation index (NDVI), leaf area index (LAI), and fraction of photosynthetically active radiation (FPAR). This smoothing process transforms noisy signals (due to clouds and other atmospheric and image acquisition artifacts) into smooth seasonal curves through a weighting mechanism and quality control (QC) layers of MODIS products. A double logistic smoothing method was selected, as it provided a better fit to the data, as suggested by [57]. As in the original SEN-ET model, meteorological data including air temperature, atmospheric vapor pressure, wind speed above the canopy, atmospheric pressure, and Remote Sens. 2024,16, 4587 6 of 29 total column water vapor (see Table 1) were downloaded from the ERA5 reanalysis dataset. Meteorological data were interpolated at the satellite overpass through the view time layers of Terra/Aqua LST products. Land cover data were used to parameterize TSEB surface characteristics (see Table 3). Instantaneous and daily incoming shortwave radiation were estimated using DEM, water vapor, and aerosols from MERRA-2 [58]. Table 3. LUT utilized to run the TSEB model for extracting the structural parameters of vegetation according to land cover type. HC represents the maximum canopy height considered for this study (m) (see Appendix A); Ωdenotes the clustering level of the canopy. Category HCΩ Evergreen Needleleaf Forest 20 1 Evergreen Broadleaf Forest 10 1 Deciduous Broadleaf Forest 10 1 Mixed Forests 15 1 Closed Shrublands 1.5 1 Open Shrublands 1.5 0.15 Woody Savannas 8 0.75 Savannas 8 0.25 Grasslands 0.5 1 Permanent Wetlands 10 1 Croplands 1.2 1 Urban and Built-Up Land 20 0 Cropland/Natural Vegetation Mosaic 1.2 0.5 Snow and Ice 0 0 Barren or Sparsely Vegetated Land 0.05 1 Water 0 0 After resampling all input data to a 1 km spatial resolution, the instantaneous ET was estimated using the TSEB model. The instantaneous ET was then upscaled to ET d . Finally, gap-filling methods were applied to estimate ET d on clouds or low-quality data based on MODIS QC (2.5 “Temporal upscaling and gap-filling methods”). 2.4. TSEB Model Overview To estimate surface energy fluxes, the Two-Source Energy Balance (TSEB) model in its series version, based on the Priestley–Taylor approximation (TSEB-PT) introduced by [ 28 ] and further developed by [ 59 , 60 ], was used (see Figure 3). The model considers the surface radiometric temperature (T RAD ), either measured from field-based radiometers or from satellite thermal sensors (LST), to be a composite of both soil (T S ) and canopy (T C ) temperatures, weighted by the fraction of vegetation (f C ) observed at a certain radiometer or thermal sensor viewing angle (θ): TRAD(θ)≈fC(θ)T4 C+1−fC(θ))T4 Si0.25 (2) where fC(θ) is the fraction of vegetation observed at a zenith angle θ defined by the following: fc(q)=1−exp−0.5WLAI cosq (3) where the factor Ω indicates the clustering level of the canopy [ 27 ] (see Table 3for value definition used in this study). T represents temperature (K). The subscripts C and S denote the canopy and soil components, respectively. The surface structural characteristics were parameterized into the Look-Up Table (LUT) presented in Table 3to accurately model the different fractions of the surface covered by the canopy fc according to the land cover type. Remote Sens. 2024,16, 4587 7 of 29 Figure 3. TSEB modelling scheme (adapted from [61]). The model contributes to partitioning the surface energy fluxes and radiometric land surface temperature (T RAD ) between vegetation and soil sources to estimate instantaneous R n , H, and G. Consequently, instantaneous LE can be calculated as the residual of the SEB (Equations (4)–(6)): LE =Rn−H−G(4) Rn,C≈HC+LEC(5) Rn,S≈HS+LES+G(6) R n is estimated as the sum of net shortwave and longwave radiation above the canopy. Net shortwave radiation is determined using incoming solar radiation and surface albedo, while net longwave radiation is calculated using air temperature and LST via the Stefan– Boltzman equation based on atmospheric emissivity [ 62 ]. Additionally, G is estimated as 35% of the net soil radiation (Rn,S) [28,63]. His calculated from the temperature gradients between the aerodynamic temperature (T A ), the temperature of the soil ( TS ), and the temperature of the canopy ( TC ), respectively (Equations (7) and (8)): HS=ρCP TS−TAC RS (7) Hc=ρCpTC−TAC RX (8) And the total sensible heat flux is calculated by Equation (9): HS+HC=ρCP TAC −TA RA (9) where ρ is the air density (kg m −3 ), CP is the specific heat of air (assumed constant at 1013 J kg−1K−1) ,T S and T C are soil and canopy radiometric temperature (K), respectively. T A is the air temperature in the surface layer measured at some height above the canopy (K), T AC is the temperature in the canopy air space (K), R S is the resistance to heat flow in the boundary layer immediately above the soil surface (s m −1 ), R X is the total boundary layer resistance of the complete canopy leaves (s m −1 ), and R A is the aerodynamic resistance (s m−1). Further information on the resistance formulation is given in Appendix A. Remote Sens. 2024,16, 4587 8 of 29 The canopy latent heat (LE C ) is estimated following the Priestley–Taylor formula [ 33 , 64 ], initially assuming a potential rate for LE C , and soil latent heat (LEs) is solved as a residual as follows: LEC=αPT fg ∆ ∆+γp Rn,C(10) LEs=RN−H−G−LEc(11) where αPT is the Priestley–Taylor coefficient, initially set to 1.26 [ 64 , 65 ], fg is the fraction of LAI that is green or actively transpiring, ∆ is the slope of the saturation vapor pressure versus air temperature curve, and γpis the psychometric constant (≈66 Pa K−1). The green vegetation fraction ( fg ) is calculated according to [ 66 ] using LAI and FPAR biophysical variables: fg=f APAR f IPAR (12) where the fraction of intercepted photosynthetically active radiation ( f IPAR ) is estimated depending on the Sun zenith angle (SZA): f IPAR =1−exp−0.5 LAI cos(SZA)(13) Under stress conditions, LEC is iteratively recalculated to reduce the αPT value [ 28 , 29 ]. To estimate the LE flux from the canopy component, the TSEB internally changes its initial αPT values of 1.26 to allow an acceptable partitioning between LE C and LE S under stressed vegetation conditions. The initial value of αPTC is down-adjusted when the TSEB results in negative values for LE S , given that condensation on the soil is unlikely to occur during the day [ 28 ]. This condition typically occurs if the potential transpiration flux is too high to be consistent with the observed surface temperature. The canopy temperature estimation is too low, resulting in a high soil temperature derived from Equation (1). This leads to high soil sensible heat (Equation (5)) exceeding the available energy at the soil surface; hence, LE S < 0 result from Equation (13). When this non-physical condition is encountered midday, it is assumed that the canopy is stressed and αPT is iteratively reduced until the LE S values are higher than 0. This iterative scheme works well in ecosystems where canopy values of αPTC are relatively conservative under unstressed conditions [ 66 , 67 ], while for stressed canopies where the soil surface is usually dry, a LE S value close to zero is a reasonable assumption [40,68]. 2.5. Temporal Upscaling and Gap-Filling Methods The instantaneous ET was upscaled by means of the ratio between instantaneous LEi estimated by the model and the ratio between instantaneous ( SIi ) and daily ( SId ) solar irradiance: ETd=SId LEi SIi (14) ETdis converted to units of mass flux (mm·day−1) through Equation (15): ETd=LEi Pwλ10386,400 (15) where Pw (kg.m −3 ) is the density of air-free water at a pressure of 101.325 kPa, λ (J.kg −1 ) is the latent heat of vaporization, λ Eis in W · m −2 , and the conversion factors 10 3 , 86,400 are in mm·m−1and s·h−1, respectively. The gap in the ET d time series due to cloud presence or low-quality satellite observations according to MODIS QC flags was filled using the ratio of actual ET and potential ET (ET a /ET 0 ) at the objective date (i.e., the date of the pixel that needs to be gap-filled) and closest valid actual ET output within 10 days of the objective date [69]. Remote Sens. 2024,16, 4587 9 of 29 2.6. Model Evaluation The modeled ET d estimated through the modeling framework adapted for MODIS data was evaluated against data collected from 11 flux towers. Certain dates were omitted from the analysis due to the following reasons: (a) absence of MODIS LST data over the towers because of cloud cover or sensor malfunctions; (b) presence of rainfall; (c) gaps in the flux tower dataset; and (d) instances where energy closure fell below 0.65. The TSEB model performance was evaluated by the following metrics: Pearson correlation coefficient (r, Equation (16)), the root mean square error (RMSE, Equation (17)), and the mean bias error (MBE, Equation (18)). r=∑n i=1Ei−EOi−O h∑N I=1Oi−O2i0.5h∑N I=1Ei−E2i0.5 (16) RMSE =s∑n i=1(Ei−Oi)2 n(17) BIAS =∑n i=1(Ei−Oi) n(18) where n is the number of observations, Oi is the observed value (in situ measurements acquired by the flux stations), Ei is the estimated value with the model, O is the average of measured values, and Eis the average of estimated values. 2.7. LST Evaluation The LST is a critical variable for estimating ET using the thermal-based TSEB model [ 7 , 24 ]. Thus, to evaluate its impact on the TSEB model, the accuracy of LST products was evaluated by comparing the LST derived from MODIS with those calculated using halfhourly longwave incoming and outgoing radiance measured from the net radiometer sensor at each flux tower at the times of the Terra and Aqua satellite overpasses, approximately 11:00 a.m. and 1:00 p.m. local solar time, respectively. The LST at the flux tower sites was calculated using Equation (19): LSTp=Lup −(1−ε)Ldown εσ 1/4 (19) where L up and L down are the longwave outgoing and incoming radiance, respectively. σ is the Stefan–Boltzmann’s constant (5.67 × 10 −8 W · m −2· K −4 ) and pthe subscript refers to the pyrgeometer. ε is the broadband emissivity that was estimated using spectral emissivity LST in bands 31 and 32 from the LST and emissivity MODIS product (MOD11A1), according to [70] following: ε=0.4857ε31 +0.5414ε32 (20) 3. Results 3.1. Flux Tower EBC Screening Results Flux tower EBC and R 2 values were lower when considering the entire dataset, with the EBC ranging from 0.68 to 1 and R 2 from 0.41 to 0.95 (top left box in each panel of Figure 4). Additionally, at all sites, the EC method tended to underestimate the energy flux components (H+LE), with the largest underestimation—up to 0.43—occurring at ES_LgS (IGBP OSH). However, after screening data were used to impose an EBC > 0.65 and with no precipitation present, the mean EBC and R 2 values improved significantly to 0.99 and 0.94, respectively. This indicated that the screened EC flux tower data were robust for evaluating the proposed modeling framework. After screening the dataset, a total number of 19265 EC observations for the whole flux tower dataset were available for model evaluation. Remote Sens. 2024,16, 4587 16 of 29 Figure 10. Mean monthly variability in estimated ET using Terra dataset. Figure 11. Mean monthly variability in estimated ET using Aqua dataset. 3.6. Seasonal Variability in Estimated ET from 2000 to 2022 in the Study Area Figures 12 and 13 reveal that seasonal ET shows more variation during the growing seasons (spring and summer) across all sites in both datasets. Additionally, the lowest ET value during the growing seasons was observed at the ES_Amo site, while the highest values were recorded at the ES_LgS and ES_Lam sites. Remote Sens. 2024,16, 4587 17 of 29 Figure 12. Seasonal variation in estimated ET from 2000 to 2022 in the study area using Terra dataset. Figure 13. Seasonal variation in estimated ET from 2002 to 2022 in the study area using Aqua dataset. 3.7. Annual Cumulative of Estimated ET from 2000 to 2022 in the Study Area The highest and lowest mean annual cumulative ET were observed at the CRO and SAV sites, respectively (Figure 14). At the CRO site, the mean ET values were around 1045 mm · year −1 and 999 mm · year −1 for Terra and Aqua datasets, respectively. Conversely, at the SAV sites, the mean ET was 818 mm/year in the Terra dataset and 767 mm · year −1 Remote Sens. 2024,16, 4587 18 of 29 in the Aqua dataset. The cumulative ET values exhibited the greatest variability at the OSH sites, ranging from 665 to 1189 mm · year −1 for Terra and from 582 to 1221 mm · year −1 for Aqua. The mean annual cumulative value of the estimated ET in the WSA site was estimated to be approximately 945 mm · year −1 and 868 mm · year −1 for Terra and Aqua datasets, respectively, which were higher than the corresponding values at the SAV sites. Figure 14. Temporal variation in annual cumulative of estimated ET from 2000 to 2022 in the study area using Aqua (left) and Terra (right) dataset. 3.8. Spatiotemporal Variability in Estimated ET from 2000 to 2022 in the Study Area Figures 15 and 16 highlight significant inter-annual variability in ET from 2000 to 2022, reflecting the complex dynamics of water exchange across the study area. Notably, in 2005, 2009, 2012, and 2022, there were pronounced changes in ET, corresponding to years marked by severe heatwaves and droughts [ 48 ]. These fluctuations in ET can largely be attributed to plant water stress, driven by water scarcity, a critical factor in the hydrological cycle. Figure 15. Spatiotemporal variability in annual cumulative of estimated ET from 2000 to 2022 in the study area using Terra dataset. Remote Sens. 2024,16, 4587 19 of 29 Figure 16. Spatiotemporal variability in annual cumulative of estimated ET from 2000 to 2022 in the study area using Aqua dataset. Additionally, the figures reveal the extended duration of drought periods, particularly from 2005 to 2006 and 2009 to 2012, during which the region experienced prolonged dry conditions. 4. Discussion 4.1. Model Evaluation and Comparison with Existing Studies In this study, the proposed modeling framework, using Terra/Aqua MODIS data, LST measurements at different viewing times, and ERA5 reanalysis data, was evaluated with 11 EC flux towers located in various land cover types (i.e., savannas, woody savannas, croplands, evergreen needleleaf forest, evergreen broadleaf forest, open shrublands). The modeled ET d evaluated under clear-sky and all-sky observations as well as the gap-filled dataset. The total number of clear-sky observation days (with available LST) were 5704 and 5630 for Terra and Aqua, respectively (Table 4). The model evaluation with both closed and unclosed datasets showed better performance in comparison with all-sky observations and gap-filled data (see Tables 5and 6). These findings highlight the significant impacts of cloudy conditions on the model’s performance. Cloudy conditions can alter various input variables critical to ET estimation, such as solar radiation, LST, vapor pressure, and surface emissivity [3,73], all contributing to increased errors in the model’s predictions. Additionally, the model performance was evaluated using all-sky observations, including 13,761 days for Terra and 13,329 days for Aqua, based on EC data collected from 2004 to 2020. The model showed the strongest agreement with ET RES , yielding a mean RMSE of 1.53 mm · day −1 and 1.40 mm · day 1 for Terra and Aqua datasets, respectively. The Aqua dataset yielded the lower mean RMSE, demonstrating higher performance compared to the Terra dataset, mainly due to differences in satellite overpass times. The mean Aqua overpass time is around 13:00 solar local time, when plants are under more pronounced water stress, likely contributing to the higher performance of Aqua data in estimating ETd(Tables 5and 6). After applying the gap-filling method, the total number of days compared with flux measurements increased significantly, from 13,761 to 19,265 for Terra and 13,329 to 19,265 Remote Sens. 2024,16, 4587 20 of 29 for Aqua. However, the statistical metrics gained from comparing the gap-filled modeled ET d with ET RES showed only slight changes, with a mean RMSE of 1.85 mm · day −1 and 1.76 mm·day−1for Terra and Aqua datasets, respectively. The overall mean results including clear-sky and all-sky observations as well as gapfilled data with Aqua showed the lowest errors with ETRES, by an MBE of 0.96 mm.day−1 and an average RMSE of 1.47 mm.day −1 , while the equivalent figures for Terra were about 1.06 mm.day−1and 1.60 mm.day−1(Table 7). The gap-filling procedure showed limited effectiveness in complex systems such as FR_Pue with the IGBP classification of EBF. In contrast, the gap-filling model showed strong performance in simpler systems, like in savanna areas at the ES_LM1, ES_LM2 and ES_Abr sites. Additionally, the gap-filled data showed lower errors with the Terra dataset, indicating that the gap-filling method has limitations in estimating ET under water-stressed conditions. Overall, the mean RMSE and MBE values revealed no significant changes compared to all-sky observations, highlighting the consistency and robustness of the modeling framework. This consistency demonstrates that the gap-filling method effectively overcame missing LST data due to cloud cover yielding a spatiotemporally continuous and robust ET d time series. Although both Terra and Aqua models overestimated ET, the intra-annual and inter-annual variability in ET in different land cover types was captured. The modelling framework with clear-sky and all-sky observations as well as gap-filled data exhibited lower performance over forested land cover (FR_Pue). This outcome is consistent with previous studies, highlighting the limited effectiveness of the TSEB model in forested regions. Given that the TSEB model relies on turbulent transportation over croplands, the aerodynamic resistance parameterization shows a low ability in complex ecosystems such as forests [74]. Conversely, the model demonstrated higher accuracy in woody savanna and savanna areas using both clear-sky and all-sky observations, as well as gap-filled data. This result aligns with previous research, which showed smaller discrepancies between the TSEB model and ground observations in bare soil and sparsely vegetated areas [32]. The authors of [ 75 , 76 ] estimated ET in almond cropland using the TSEB model with Sentinel-2 and Sentinel-3 at a high spatial resolution. ET was evaluated against the soil water balance model within 148 observations, resulting in a root mean square deviation of 1.54 mm.day −1 and r value of 0.33. Compared to current study, in croplands (FR_Aur and FR_Lam sites), RMSE values ranged from 1.58 mm.day −1 to 1.87 mm.day −1 and r values were greater than 0.50. The comparison of Terra/Aqua MODIS LST estimates with flux tower T RAD observations (Figures 8and 9) yielded a RMSE > 2 K and an MBE < − 2 K at most of the flux stations, thus underestimating field estimates of LST, and it was particularly high in the FR_Lam site. Although these measurements did not precisely coincide with the satellite overpass times, as they were 30 min LST averages, a similar bias was found in boreal forest settings, in which LST field measurements were collected at 1 min timesteps [ 22 ]. The increase in error and bias compared to local estimates of LST using flux tower data can be attributed to several factors. These include biases in the LST due to atmospheric and emissivity corrections, sensor biases, or the unrepresentative nature of the LST at a 1 km resolution relative to the flux tower footprint. When applied regionally, the TSEB technique is typically performed in a time-differential mode as applied in this study, with the LST as one of the main input variables. Thus, the LST disagreement with field observations might have resulted in an increased bias and error in ET d estimation, lowering the final agreement. To reduce errors caused by uncertainties in atmospheric correction and emissivity when determining the absolute LST value, an alternative time-differential TSEB technique using a dual-time difference (DTD) model in the observed LST will be implemented, adding the model refinements found in this paper and combining Terra and Aqua LST observations [34,60]. Overall, the statistical metrics in the present study indicated that the differences observed were consistent with those found in previous research on estimating ET at Remote Sens. 2024,16, 4587 21 of 29 regional scales. Although many of these studies only evaluated the ET accuracy over a partial annual cycle, the present study focused on estimating ET d and evaluating the modeled ET d with ground observations to capture the temporal year-round variability in ET in different land cover types from 2000 to 2022. For instance, [ 18 ] employed the SEB model to estimate ET d using MODIS imagery and ground-based meteorological data. The model’s results were compared with data from three EC flux towers located in native and managed grasslands over two years. The results showed that the model overestimated ET d , with MBE and RMSE values of less than 30 W · m −2 and 50 W · m −2 (1.05 mm · day −1 and 1.75 mm · day −1 , respectively). The authors of [ 77 ] utilized the evaporative fraction method with NOAA–AVHRR dataset, which has a medium spatial resolution of approximately 1 km, and the Digital Airborne Imaging Spectrometer (DAIS) dataset, which has a high spatial resolution of approximately 5 m. The ET estimations from these two datasets were compared over 13 selected dates across different land cover types. The results showed a mean RMSE of 1.4 mm · day −1 . The authors of [ 7 ] presented a method for estimating continuous ET d at a 1 km spatial resolution by integrating MODIS products and reanalysis data. To fill the gaps in LST products under cloudy conditions, a data fusion approach was applied to estimate ET using the TSEB model over croplands from 2008 to 2017. The model was validated against a total of 5707 days from five EC flux towers all located in croplands. The results showed that RMSE values range from 0.87 mm · day −1 to 1.37 mm · day −1 . The authors of [ 78 ] used the Penman–Monteith model, combined with MODIS data, to estimate ET d at 1 km spatial resolution. The model was evaluated using data from 501 days from three EC flux towers in agricultural regions (one in wheat/summer maize and the other one in orchard and maize areas), achieving a mean RMSE of 0.7 mm · day −1 . The authors of [ 79 ] evaluated the accuracy of the MOD16 product, which estimates ET using the Penman– Monteith model. The ET product was validated using data from 51 days from four EC flux towers in various land cover types (vegetation, village, maize, and orchard areas). The result showed a mean RMSE of 1.2 mm · day −1 . After model improvements, the RMSE dropped to 0.9 mm · day −1 . The authors of [ 80 ] proposed a gap-filling method to estimate continuous ET d using the Penman–Monteith model and remote sensing data at a 1 km spatial resolution for 2016. The model was evaluated using data from 730 days from two EC flux towers located in croplands. The results showed a mean RMSE of 0.79 mm · day −1 and 0.62 mm · day −1 . Hence, unlike many studies in the literature, the present study has thoroughly evaluated the performance of the ET model across different land cover types and over long-term daily temporal scales. Accurate ET estimation approaches should also be capable of capturing the temporal variability in ET across diverse land cover types and over year-round periods to ensure robust and reliable assessments are also evaluated in this study. Although the RMSE values from most of the flux sites align with previous studies, some flux stations exhibit lower agreement. While most correlations are statistically significant at α = 0.05, certain flux stations showed a particularly low agreement. The authors of [ 7 ] reported similar agreement values when modeling daily ET over a long time series of MODIS data in crop areas, aligning well with the performance of certain flux stations in this study (e.g., ES_LM2, ES_LM1, ES_Abr). However, a markedly lower agreement is observed at most open shrubland flux sites, such as ES_Amo, ES_LJu, and ES_Agu. This behavior might be explained by the low ET d variability in these sites. Averaged ET d values of 0.43 ±0.33, 0.44 ±0.43 , and 0.86 ± 0.58 were found for ES_Amo, ES_LJu, and ES_Agu, respectively. Although the RMSE values fall within the range reported in previous studies, these sites exhibited MBE values of 0.59, 0.17, and 1.37 mm · day −1 , respectively (see Aqua results in Table 4, Model vs. Residual column), leading to greater disagreement, particularly at ES_Amo (see Figure 5). At these sites, the modeling framework failed to capture temporal dynamics. However, the high number of processed and analyzed images, combined with the heterogeneity of these surfaces and the period studies, provide consistency to the error estimates. Thus, further research will be needed at local scales using field data to enhance the model performance in these environments. Remote Sens. 2024,16, 4587 22 of 29 Additionally, although many flux sites selected were representative of at least one MODIS pixel, some flux sites in open shrublands, such as ES_Agu, and croplands, such as FR_Aur and FR_Lam, were relatively smaller (around 500 × 500 m). Typically, the TSEB model performs adequately in croplands. However, results yielded higher RMSE values and intermediate agreement compared to other flux sites. This may be explained by the presence of forest patches near both crop sites, which might overestimate or underestimate the LAI when the crop was active or harvested, respectively, decreasing the model agreement (see Figure 5). Finally, it is noteworthy that the ES_Cnd site showed better agreement and a similar RMSE compared to previous studies using Sentinel imagery and the Sen-ET modeling framework at a higher spatial resolution. 4.2. Application of Long Time Series of ETd The modeling framework applied in the study area successfully captured the spatiotemporal variation in ET within various types of land cover from 2000 to 2022. Some past studies reported considerable drought periods within the study area over the last two decades. In connection with these studies, changes observed in the intraannual and inter-annual timescale of the estimated ET d align with the drought events documented [48,76]. The annual cumulative ET underwent a significant reduction in 2005 in all sites except for the FR_Aur and FR_Lam sites with both Terra and Aqua datasets. Ref. [ 47 ] utilized the Standardized Precipitation Index (SPI) and the Standardized Precipitation– Evapotranspiration Index (SPEI) to identify drought events in the IP from 1961 to 2011, noting 2005 as one of the major drought events in this study area. Additionally, [ 48 ] characterized drought events across the IP using the SPI from 1916 to 2020, identifying 2005 as having notable drought periods. In the northern (FR_Aur and FR_Lam sites) and southeastern (ES_Amo and ES_Agu sites) parts of the study area, there were remarkable reductions in the growing seasons (Figures 12 and 13) and annual cumulative ET (Figure 14) in 2003. In addition, the outliers beyond the lower fence (Figures 10 and 11) observed in the FR_Lam and FR_Aur sites in May and August 2003 are associated with the significant ET changes in this area. [ 81 ] analyzed soil moisture during the 2003 drought event in Europe and indicated that 2003 was a significant drought event in the Euro-Mediterranean region. Notably, there is a declining trend in the ET time series across all sites in 2022 (Figure 14). Ref. [ 82 ] demonstrated that the Euro-Mediterranean region was under severe drought conditions in 2022. Figures 10 and 11 show outliers beyond the lower fence at the ES_Pue site in July 2022, illustrating abnormal changes in the ET time series that were also captured by the modeling framework. In the southeastern part of the study area (OSH sites), the minimum annual cumulative ET values were recorded in 2005, indicating significant water stress on vegetation (Figure 14). This finding aligns with [ 83 ], which used the SPI to monitor drought in this region. The result demonstrated that the southeastern part of the study area experienced a significant drought in 2005. For both Terra and Aqua datasets, outliers beyond the lower fence at the ES_LgS site (Figures 10 and 11) were observed in May, September, and March in 2005, 2012, and 2017, respectively, also aligning with [48], which reported drought events. The annual cumulative ET values at sites located in savanna and woody savanna areas (ES_LM1, ES_LM2, ES_Abr, ES_Cnd) showed less fluctuation over the time series, with notable reductions in 2005 and 2012 (Figure 14). This result agrees with [ 84 , 85 ], which reported 2005 as the driest year in the Guandiana River basin. Additionally, the annual cumulative ET experienced significant reductions at the ES_LM1, ES_LM2 and ES_Abr sites in 2009, and at the ES_Abr site in 2019, aligning with the findings of [ 48 ]. In Figures 10 and 11, outliers beyond the lower fence observed at savanna sites during the growing seasons in 2009 and 2012 indicate unexpected changes in the ET time series, corresponding with drought years in this area reported by [ 46 , 48 ]. Moreover, outliers beyond the lower fence were observed in 2005 and 2012 at the ES_Cnd site from February to July. Remote Sens. 2024,16, 4587 23 of 29 This research highlights that the study areas experienced significant water stress from 2000 to 2022, varying across the regions (Figures 15 and 16). This underscores the critical importance of monitoring ET to recognize vegetation water stress and drought, as evidenced by significant changes in the ET time series. In addition, the ET trends suggest that future droughts could exacerbate water scarcity in these vulnerable areas. This pattern could have implications for food security and ecosystem health in the region. 4.3. Limitations of the Proposed Modelling Framework and Future Implementation Without any site-specific adjustments and relying only on global datasets, the modeling framework produced ET d estimates comparable to those reported in other studies in terms of error, demonstrating robustness across various land cover types. However, agreement at certain flux sites could be increased by improving some of the input variables. Similarly to other studies, the model accuracy and agreement could be improved by using more accurate meteorological data (especially wind speed measurements). As [ 24 , 37 ] suggested, atmospheric forces from numerical weather prediction models may introduce some uncertainty to the TSEB model compared to using local meteorological data. Several studies [ 7 , 24 ] have relied on ERA5 reanalysis data to estimate ET. Despite the low spatial resolution of ERA5 (approximately 35 km), comparisons with meteorological instruments at several flux sites have shown strong agreement for most key meteorological parameters, except for wind speed and solar radiation. This discrepancy may directly affect the accuracy of the TSEB model in estimating Rn and, consequently, ET, since ET is calculated as a residual of the energy balance. Therefore, errors in ET d could be significantly reduced if more accurate solar radiation inputs were used, especially in temperate areas (i.e., radiationlimited) which are more sensitive to uncertainty in available energy[ 37 ]. Additionally, errors in both air temperature and wind speed primarily affect the retrievals of the sensible heat flux, which is more critical in the estimation of LE in semi-arid (i.e., water-limited) areas [37]. To address these issues, an alternative climate forcing dataset such as ERA-Land, with an improved spatial resolution of approximately 9 km, could be used as TSEB meteorological input. However, the primary drawback is that the ERA-Land product applies a land mask, potentially excluding a 9 km area in regions near the seashore, rendering the modeling framework inapplicable in those areas. Additionally, there is no available layer for wind speed. For solar radiation, other methodologies or products could be utilized. The authors of [ 86 ] proposed a DEM-based hybrid methodology, which has been successfully used to estimate Rn at daily scales [ 87 ] and can also be useful for ET daily upscaling. Other current solar radiation products, such as the MSG-SAF solar radiation product, have demonstrated accurate evaluation results with in situ measurements, proving to be an effective input for calculating evaporative fluxes [88]. Finally, values for H c and Ω in Table 3were meant to fit a global-based land cover legend and can significantly deviate from the actual values at specific flux sites. Thus, utilizing ancillary remote sensing-based products derived from LiDAR data [ 89 ] would increase the accuracy of surface characteristic parameterizations. 5. Conclusions Estimating ET at regional scales with daily temporal resolution can benefit drought monitoring in water-limited regions. Global ET datasets often lack a daily temporal resolution with a finer spatial scale. This research addresses these gaps by providing spatiotemporally continuous ET data, offering a significant advancement over existing ET products in this region. A modeling framework, based on the SEN-ET scheme, was used with Terra/Aqua MODIS data and the ERA5 atmospheric reanalysis dataset to estimate ET across various land cover types (i.e., savannas, woody savannas, croplands, evergreen broadleaf forest, open shrublands) from 2000 to 2022. The main conclusions are as follows: 1. The modeling framework was evaluated with clear-sky and all-sky observations as well as gap-filled data. Under clear-sky observations, the model demonstrated lower Remote Sens. 2024,16, 4587 24 of 29 errors and higher correlations with both closed and unclosed datasets compared to all-sky observations and gap-filled data. These findings highlighted the negative impact of cloudy conditions on the model’s performance. Moreover, the overall mean results, including clear-sky and all-sky observations as well as gap-filled data with the Aqua dataset, showed the lowest errors against ET RES , by an MBE of 0.96 mm.day −1 , an average RMSE of 1.47 mm.day −1 , and r value of 0.51. The equivalent figures for the Terra dataset were about 1.06 mm.day −1 , 1.60 mm.day −1 , and 0.52. The Aqua dataset showed better performance compared to the Terra dataset, primarily due to differences in satellite overpass times and a smaller MBE between Aqua LST products and in situ LST measurements. The overestimation errors obtained in this study were within the range of errors observed in previous studies that estimated ET at regional scales. Importantly, the model was able to capture both intra-annual and inter-annual variability in ET in different land cover types. By successfully gap-filling cloudy or low-quality satellite observations, the model achieved more comprehensive spatiotemporal coverage, solving a major limitation in consistent ET monitoring. Furthermore, the performance metrics (RMSE and MBE) remained consistent after applying gap-filling methods, further underscoring the robustness and reliability of the model’s ET estimations. However, for some flux stations, a low agreement but statistically significant correlation was found. Thus, further modelling framework and TSEB adjustments will be needed to improve the modelling results, such as using more accurate solar radiation products and climate forcing data. Moreover, some local TSEB mode evaluations, especially in semi-arid open shrublands, will be also needed. 2. The alignment between the spatiotemporal trends in the estimated ET time series and documented drought events underscores the reliability of ET as a robust indicator for drought monitoring. The modeling framework’s ability to reflect changes in ET in response to water stress and drought conditions highlights its potential for improving water resource management at regional scales. The findings suggest that the ET data generated from this modeling framework can serve as a valuable tool for drought early warning systems. Author Contributions: Conceptualization, M.B. and J.C.; methodology, M.B. and J.C.; software, M.B. and J.C.; data curation, M.B., J.C., and M.P.-S.; writing, review, and editing, M.B., J.C. and M.P.-S.; graphical deployment, M.B.; supervision, J.C. All authors have read and agreed to the published version of the manuscript. Funding: This study was funded by the project ET4DROUGHT (No. PID2021-127345OR-C31) and DigiSPAC [TED2021-131237B-C21] both funded by the Ministry of Science and Innovation (MICIN-AEI). Data Availability Statement: Restrictions apply to the datasets. The datasets presented in this article are not readily available because the data are part of an ongoing study. Requests to access the datasets should be directed to Jordi Cristóbal (jordi.cristó[email protected]). Acknowledgments: Thanks are given to the FLUXNET Network and the Integrated Carbon Observation System (ICOS) for the open access “Warm Winter 2020 ecosystem eddy covariance flux product for 73 stations in FLUXNET-Archive format—release 2022-1, “https://doi.org/10.18160/2 G60-ZHAK (accessed on 1 January 2022)”. Thanks are also given to Hector Nieto for developing SENT_TSEB, “https://github.com/hectornieto/pyTSEB (accessed on 4 March 2018)”. Conflicts of Interest: The authors declare no conflicts of interest. Appendix A. Summary of Equations Used to Estimate Aerodynamic Resistances in TSEB In this appendix, the resistances involved in Equations (5)–(7) are summarized [ 28 ]. R A is the aerodynamic resistance in the surface layer computed from the stability-corrected temperature profile equations (s·m−1), defined as follows: RA=([In((zU−d0)/z0M)−ΨM][In((zT−d0)/z0H]−ΨH)] k2u(A1) Remote Sens. 2024,16, 4587 25 of 29 where d O is the displacement height, U is the wind speed measured at height z U , k is von Karman’s constant ( ≈ 0.4), z T is the height of the T A measurement, ΨM and ΨH are the Monin–Obukhov stability functions for momentum and heat, respectively, and z OM is the aerodynamic roughness length. Both d O and z OM were estimated using the LAI and canopy height (HC, see Table 3) as follows [90,91]: dO=0.7 − 1 (5n)1−e(−3.3n)!HC(A2) where n is the within-canopy wind speed profile extinction coefficient, parameterized as follows: n=CdLAI 2u∗ u(h)2(A3) where Cd is the drag coefficient of the foliage elements with a value of 0.2, and the ratio u∗/u(h) is parameterized as follows: u∗ u(h)=0.360 −0.264 e(−15.1∗Cd∗LAI)(A4) Finally, zOM is defined as follows: zOM =(1−dO)e−0.4 / u∗ u(h)Hc (A5) R S is the resistance to the heat flow in the boundary layer above the soil (s · m −1 ), defined as follows: RS=1 a′+b′US (A6) where a’ ≈ 0.004 m · s −1 , a’ ≈ 0.012 m · s −1 , U S is the wind speed (m · s −1 ) at the height above the soil surface where the effect of the soil surface roughness is minimal (~5 cm) and estimated from the [ 82 ] wind attenuation model through the canopy layer (see Equation (A8)). Finally, R X is the total boundary layer resistance of the vegetation leaf layer (s · m −1 ), defined as follows: RX=C′ LAIS Udo+Zom 1 2(A7) where C’ is derived from weighting a coefficient in the equation for the leaf boundary layer resistance (assumed a constant ~90 s 1/2 m −1 ) and S is the leaf size (m). U d+Zm is given by the following: Udo+Zom =UCexp−a1−do+Zom HC (A8) where U C is the wind speed at the canopy height (Hc) and ais the wind attenuation coefficient of the [ 82 ] wind attenuation model, which is a function of the LAI, Hc, and leaf size. For a full description of the TSEB resistance formulations, see [27–29]. References 1. IPCC. IPCC Synthesis Report—Climate Change 2023. In An Assessment of the Intergovernmental Panel on Climate Change; IPPC: Geneva, Switzerland, 2022; Volume 335. 2. Noto, L.V.; Cipolla, G.; Pumo, D.; Francipane, A. Climate Change in the Mediterranean Basin (Part II): A Review of Challenges and Uncertainties in Climate Change Modeling and Impact Analyses. Water Resour. Manag. 2023,37, 2307–2323. [CrossRef] 3. Wang, K.; Dickinson, R.E. A Review of Global Terrestrial Evapotranspiration: Observation, Modeling, Climatology, and Climatic Variability. Rev. Geophys. 2012,50, RG2005. [CrossRef] 4. Jung, M.; Reichstein, M.; Ciais, P.; Seneviratne, S.I.; Sheffield, J.; Goulden, M.L.; Bonan, G.; Cescatti, A.; Chen, J.; de Jeu, R.; et al. Recent Decline in the Global Land Evapotranspiration Trend Due to Limited Moisture Supply. Nature 2010,467, 951–954. [CrossRef]