Evapotranspiration Acquired with Remote Sensing Thermal-Based Algorithms: A State-of-the-Art Review

Abstract

Almost fifty years have passed since the idea to retrieve a value for Evapotranspiration (ET) using remote sensing techniques was first considered. Numerous ET models have been proposed, validated and improved along these five decades, as the satellites and sensors onboard were enhanced. This study reviews most of the efforts in the progress towards providing a trustworthy value of ET by means of thermal remote sensing data. It starts with an in-depth reflection of the surface energy balance concept and of each of its terms, followed by the description of the approaches taken by remote sensing models to estimate ET from it in the last thirty years. This work also includes a chronological review of the modifications suggested by several researchers, as well as representative validations studies of such ET models. Present limitations of ET estimated with remote sensors onboard orbiting satellites, as well as at surface level, are raised. Current trends to face such limitations and a future perspective of the discipline are also exposed, for the reader’s inspiration.

1. Introduction

Evapotranspiration (ET) is a process that includes evaporation from the surface, such as open water bodies, soil and vegetation, and transpiration as the water released mostly by the plant leaves transported from the root system. ET is a key factor in the hydrological cycle, since it describes the mechanism and energy needed to transport the liquid water stored in the soil-watershed-canopy system to the atmosphere, converted into water vapour.

In the broader context, irrigated agriculture still accounts for 70% of freshwater withdrawals [1]. This is a serious problem noticed by the Horizon Europe program. In this framework, ET plays an important role as a controlling factor of the water cycle and energy transport among the biosphere, atmosphere and hydrosphere. The importance of ET has been evidenced from the very beginning of the Remote Sensing discipline, in application fields such as hydrology, for instance for the scheduling of field-scale irrigations and tillage [2] or in meteorology and climatology, for prediction of natural hazards such as floods and droughts [3] or in agriculture, where for vegetated surfaces, ET can be used as an indicator of plant water stress [4].

In situ ET field retrievals from ground instrumentation, based on conventional techniques (e.g., Eddy covariance, Bowen ratio, lysimetry, etc.), offer the advantage of a direct determination of the value of ET for a measuring site [5]. However, these instruments are often costly, their data require an adequate treatment which is quite time-consuming, and instrumental problems may happen [6]. A commonly-used alternative for the estimation of ET is the Penman-Monteith equation, which uses monitored variables (e.g., air temperature or humidity), very convenient with basic meteorological stations, especially in the form suggested by Food and Agriculture Organization (FAO) [7,8]. Furthermore, for an area which is spatially heterogeneous, measurements made in a single location may not be representative of the full area [9]. For a given area in a humid climate, annual ET can be as large as half of the precipitation, whereas for arid and semi-arid areas ET is almost equivalent to the total annual rainfall [10].

In summary, even if accurate ET estimates can be obtained for the instrumented sites, these values are not useful at larger scales due to the usual heterogeneity of the landscape and the complexity of the hydrological processes [2]. Conventional ground measurement techniques do not appreciate spatial distribution at large-scales, especially in regions composed of well-contrasted areas where advection may be important.

Remote sensing technology is a globally consistent and economically feasible means to estimate ET values at regional and meso-scales on the Earth’s surface [11], since the approach directly links surface radiances and the components of the surface energy balance [12,13,14]. Over the last 30 years, combined use of satellite remote sensing data from optical and thermal infrared sensors has provided substantial progress in the estimation of ET [15,16,17,18]. Based on the concept of surface energy balance and net radiation, most remote sensing models have estimated ET for application studies such as water consumption, water resources planning and management over watersheds [19] or modeling ecological processes and analyzing biophysical characteristics of landscape [20].

The present document is a review report of the most popular remote sensing Land surface temperature (LST)-based ET algorithms used during the last three decades. This article revisits and updates of the state of the art on this topic reported in previous review articles [17,21,22,23,24,25], incorporating recently published works showing modifications and improvements of the models. This review aims to combine the current state of the art on the subject under study, making it suitable for both experts and users of ET estimated from satellite. Nevertheless, reading the aforementioned review articles is highly recommended for a more complete and up-to-date idea of the subject under study.

A brief description of the surface energy budget is provided in Section 2, while some of the most popular remote-sensing approaches to estimate ET are described in Section 3, indicating their evolution with time and the representative validation studies in the last twenty years. Also advantages and disadvantages of each ET model are discussed. Section 4, it is shown an in-depth review of the limitations and current trends in remote sensing techniques to estimate ET, as well as the future prospects to be addressed by the remote sensing ET community. Finally, Section 5 resumes the main conclusions of this review study.

2. Surface Energy Balance (SEB) Theory

It is a common approach to consider that an energy balance equation can be applied to any point of the Earth Surface [26]. The Net Radiation is the flux term that warms the surface during the day and cools it at night.

Instead ET is a process related with the quantity of energy used by the surface for the water phase change. When condensation occurs at the surface (usually at night) energy is transferred to the surface. Ignoring at this point any other effects (e.g., advection, heat stored in the upper soil or the plant energetics), the surface energy balance equation can be written as:

R_n = G + H + LE

(1)

where the net radiation R_n (in W/m²), represents the total heat energy considered in the balance between incoming and outgoing short-wave (0.15 to 3 μm) and long-wave (3 to 100 μm) radiation under steady atmospheric condition, and it can be expressed as:

R_n = (1 − α_s) R_s + ε_s ε_aσT_a⁴ − ε_sσT_s⁴

(2)

where α_s is the surface shortwave albedo, calculated with satellite data as a combination of narrow spectral reflectance estimated from bands of several remote sensing instruments [27], however, it can also be estimated in other ways such as using a 4-component net radiometer. R_s is commonly called insolation and it is a measure of the incident solar radiation, either photosynthetically active radiation in the visible spectrum (400–700 nm) or in the total shortwave (150–3000 nm), received on a given surface area and time. There also exist several remote sensing products to estimate R_s. For instance, the geostationary satellite Geostationary Operational Environmental Satellite (GOES) [28] or polar satellites like the Moderate Resolution Imaging Spectroradiometer (MODIS) [29] or Advanced Very High Resolution Radiometer (AVHRR) [30]. Atmospheric emissivity (ε_a) can be estimated from vapor pressure values [31]. There are several ways to retrieve air temperature (T_a). Some of them are the remote sensing MODIS products MOD07 [32], or from synthetic atmospheric profiles [33] included in radiative transfer codes, which has been proven to be more accurate than MOD07 compared with meteorological data from a near station [34]. Also T_a can be retrieved from other datasets like tower-based measurements or reanalysis data such as European Centre for Medium-Range Weather Forecasts (ECMWF, https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5 (accessed on 1 July 2022). Surface emissivity (ε_s) and LST can be retrieved from remote sensing techniques separately [35,36,37] as well as at the same time [38,39]. σ is the Stefan-Boltzmann constant. Two last terms in Equation (2) can be substituted by downward and upward longwave radiation terms. Both longwave radiation variables can be estimated by remote sensing techniques. For the downward term there exist operational products from MODIS and GOES [40], AVHRR and Spinning Enhanced Visible and InfraRed Imager (SEVIRI) [30], while for the upward radiation term there are products and algorithms for MODIS, Advanced Spaceborne Thermal Emission and Reflection (ASTER), Thematic Mapper (TM), SEVIRI, AVHRR, Along Track Scanning Radiometer (ATSR), Advanced Along-Track Scanning Radiometer (AATSR), GOES, Geostationary Meteorological Satellite (GMS), Visible Infrared Imaging Radiometer Suite (VIIRS), and Atmospheric infrared sounder (AIRS) [41]. Reader is referred to review study of [42] for further details about R_n estimated from ground measurements to remote sensing tools.

The soil heat flux G (in W/m²), is the heat energy used for warming or cooling substrate soil volume. As of today, it is impossible to directly measure G from remote sensing data. However, in many studies G is retrieved as a pondered estimate of Rn [43], or related to other variables easily retrievable from remote sensing data, such as Leaf Area Index (LAI) [44], Normalized Difference Vegetation Index and LST (NDVI) [45], and Fraction Vegetation Cover (FVC) [46]. Reader is referred to review study of Gao et al. [47] where this issue is investigated and a variety of methods are intercompared.

The sensible heat flux H (in W/m²), is the driving mechanism through which the air above the surface warms or cool as heat is transported from or to the surface. In the SEB Equation (1), it is expressed as:

$$H=\rho C_P\frac{\left(T_o-T_a\right)}{r_{aH}}$$

(3)

where ρ is the density of air (kg/m³), C_P is the specific heat of air (1004 J/kgK), T_a is the air temperature usually at a height z of 2 m above the canopy [48], T_o is the aerodynamic temperature near the surface and r_aH is the aerodynamic resistance (s/m) obtained from the Monin-Obukhov Stability Theory (MOST) [49]. T_o is a parameter not easily estimable, so usually in remote sensing it is considered as the radiometric T_s retrieved from TIR sensors on board orbiting satellites, after some reformulation of the r_aH parameter to account for such approximation [43]. This issue is addressed in the next Section 2.2.

The Latent Heat Flux LE (in W/m²), governs the loss of latent heat from surface due to evapotranspiration. It is usually retrieved as a residual from the SEB Equation (1) and it is linearly related with ET:

$$LE=\lambda \mathit{ET}$$

(4)

where λ is the latent heat of vaporization of water, which is dependent on the temperature, and usually takes a value of 2.45 × 10⁶ J/kg at 20 °C.

The remote sensing LST-based ET models produce instantaneous (ET_i) values, for which the algorithms will be described in the next section. This instantaneous estimation is of not much use for most hydrological and water resources management applications and it is customary to convert ET_i retrieved from overpass times of satellites to daily (ET_d) or longer time scales to make full usable application of the remote sensing data. Several authors have proposed approximation methods to derive ET_d values from ET_i estimated from satellite data. Some of these well-accepted methods are described below.

Jackson et al. [50] related ET_i to ET_d value according to the diurnal solar irradiance evolution, following the sine function equation:

$${\mathit{ET}}_d=\frac{2N}{\pi \sin \left(\frac{\pi t}{N}\right)}{\mathit{ET}}_i$$

(5)

$$N=0.945\left(n_1+n_2{\sin }^2\left(\frac{\pi \left(D+10\right)}{365}\right)\right)$$

(6)

$$n_1=12-5.69·{10}^{-2}lat-2.02·{10}^{-4}la{t}^2+8.25·{10}^{-6}la{t}^3-3.15·{10}^{-7}la{t}^4$$

(7)

$$n_2=0.123 lat-3.1·{10}^{-4}la{t}^2+8·{10}^{-7}la{t}^3+4.99·{10}^{-7}la{t}^4$$

(8)

where t is the time of day, D day of year and lat is the geographical latitude between 60°S and 60°N. For cloud free or a relatively constant cloud cover days, the sine function (5) has produced reliable estimates of ET_d [11,51].

Another method assumes that evaporative fraction (EF) remains constant during the daylight hours [52]. EF is the rate between ET_i and available instantaneous energy (Rn−G)_i. However, it has been found that during nighttime ET could reach as many as about 10 percent of the daily totals [53]. So, the night induces a percentage that is added to EF, and after that assumes EF constant during the 24 h of the day [52,53]. Then ET_d can therefore be written as:

$${\mathit{ET}}_d=1.1\frac{R_{nd}}{{\left(R_n-G\right)}_i}{\mathit{ET}}_i$$

(9)

where R_nd is the daytime available energy, actually it is (Rn−G)_d but it is generally assumed that daily G is negligible (see Section 2.2).

Allen et al. [54] considered that EF was not able to capture the impact of advection, changing wind and humidity conditions along the day. So, they defined the constant reference ET fraction (ET_rF) as the ratio between ET_i and reference ET over the standardized 0.5 m tall alfalfa (ET_r) [7]. For horizontal areas (no significant slope terrain effects) Allen et al. proposed the next ET_d approximation:

$${\mathit{ET}}_d=\frac{E{T}_{r}F}{E{T}_r}=\frac{{\mathit{ET}}_i}{E{T}_{r}E{T}_{r,d}}$$

(10)

where ET_r,d is the cumulative daily reference of ET [54].

2.1. Aerodynamic Resistance to Heat Transfer Based on MOST

Aerodynamic resistance to heat transfer (r_aH in Equation (3)), is an ecosystem parameter which is affected by the surface roughness (vegetation height and structure), wind speed (u) and atmospheric stability. Increasing the intensity of the wind or the roughness of the surface imply a decrease in the value of r_aH [48]. It must be adjusted according to surface characteristics and based on the calculation of momentum flux and the observed logarithmic profile of wind speed close to surfaces, it is given by:

$$r_{aH}=\frac{\left[ln\left(\frac{z_u-d}{z_{0M}}\right)-{\psi }_M\right]\left[ln\left(\frac{z_T-d}{z_{0H}}\right)-{\psi }_H\right]}{{\kappa }^{2}u}$$

(11)

where z_u and z_T are the heights at which u and T_a were respectively measured, d is the zero-displacement height, z_0M and z_0H are the roughness length of momentum and sensible heat exchange, respectively. κ = 0.41 is the dimensionless von Karman constant, Ψ_M and Ψ_H are the Monin-Obukhov stability functions [49] for momentum and sensible heat turbulent fluxes. Usually, two levels of temperature and wind are required to compute these functions, typically at 2 m above the canopy height for the atmospheric variables, while wind is set to zero at z_0M and the value of LST is taken at z_0H.

Atmospheric stability correction on r_aH is an important issue since not considering it leads to an overestimation of ET when the canopy-air temperature difference is greater than about ±2 °C [48]. To determine the stability functions the Obukhov length [49] is defined as:

$$\mathsf{\Lambda }=\frac{u^{*}{}^3\rho C_P{T}_a}{g\kappa H}$$

(12)

where u* is the friction velocity and g the gravity acceleration. It must be solved with H and u* [55]. When $\mathsf{\Lambda }$ < 0, means atmospheric instability, and for $\mathsf{\Lambda }$ > 0 stable conditions. It is worth to note that both T_a and u* are external parameters which retrieval is non-dependent on remote sensing techniques [56,57,58,59,60].

In case of atmospheric thermal instability and no predominant free convection conditions (usually at daytime), Ψ_M and Ψ_H can be expressed as [61]:

$${\mathsf{\psi }}_M=2ln\left(\frac{1+x}{2}\right)+ln\left(\frac{1+x^2}{2}\right)-2arctan\left(x\right)+\frac{\pi }{2}$$

(13)

$${\mathsf{\psi }}_H=2ln\left(\frac{1+x^2}{2}\right)$$

(14)

$$\mathrm{x}={\left(1-16\frac{z_u-d}{\mathsf{\Lambda }}\right)}^{0.25}$$

(15)

Under atmospheric stable conditions (prevailing at night) [62,63] proposed the same formulation for the Monin-Obukhov stability functions, as:

$${\mathsf{\psi }}_M={\mathsf{\psi }}_H=-5\frac{z_u-d}{\mathsf{\Lambda }}$$

(16)

The effective roughness length for momentum (z_0M) is a distance above the zero-plane displacement height where the wind speed is assumed to be zero when log-profile wind speed is extrapolated downward, rather than at true ground surface [64]. z_0M and d are complex functions of the surface roughness (vegetation height and structure), estimated from models that consider canopy height, width, and element spacing into account [65,66]. There are also alternative methods that integrate the LAI in the estimation of z_0M and d [67].

The roughness length for heat exchange (z_0H) can be expressed as a function of z_0M as [67]:

$$ln\left(\frac{z_{0M}}{z_{0H}}\right)=\frac{\kappa }{B}$$

(17)

where B the Stanton number (dimensionless). κ/B rate is a complex function, also called “excess resistance”, which is dependent on time of the day, weather, and the type of vegetation, but constant values for z_0H/z_0M of 0.1 and 0.2 are often used for dense and sparse crops, respectively [68].

In summary, surface roughness plays an important role in the estimation of H and it is worth noting it still remains a challenging issue for large scale retrieval of the turbulent heat fluxes. However, resistance Equation (11) is applied since early seventies over local and regional scales and for varied vegetation covers [48,69,70]. We will avoid further details on MOST at this point not to distract the reader from the main aim of this review. For more information about surface roughness formulation, reader is referred to reviews of [71,72] and an in-depth study by [73].

2.2. Satellite Estimation of T_o

Regarding T_o (Equation (3)), it represents the temperature of the apparent source/sink of sensible heat flux [74], and it is computed at the level of d + z_oH, averaged from the weighted contribution to the aerodynamic conductance of each element composing the full canopy scene. According to [75] T_o is neither measured nor easily estimated, while T_s is relatively easily measured with infrared thermometers and scanners. T_o has been usually substituted by T_s in Equation (3) since it is very difficult to measure [69,76]. Min et al. [77] indicated that for homogeneous and isothermal surfaces the definition of aerodynamic and thermodynamic (canopy or surface radiometric) temperatures are equivalent, however, assuming T_o equal to T_s at canopy level could be problematic since it is known that T_s is 2–3 °C higher than T_o for uniform canopy covers [78], and up to 10–15 °C higher for incomplete canopy covers [79]. T_s is also influenced by the viewing angle, so if it is measured at nadir T_s tends to be larger than T_o due to the direct incident sunlight. Nevertheless, at larger viewing angles between 50° and 70° from nadir, T_s coincides with T_o [80]. Difference between T_s and T_o leads to errors in the estimation of H which in turn leads to errors in the estimation of LE (ET). T_o account for those differences, researchers commonly have parameterized, with mixed results, the H equation in order to be able to use T_s instead of T_o. However, there is an effect of vegetation three-dimensional (3D) structure and its interaction with directional wind speed (wind direction, angle of attack), on the generation of T_o, roughness length (for heat and momentum transfer), friction velocity; which in turn affect H and LE.

More recently, [81] improved a T_o model for deficit irrigated corn by including a turbulence mixing-row resistance that attempts to incorporate the interactions between crop rows orientation and wind direction (the angle of attack). However, the approach is only applicable to the corn canopy structure and site conditions. In terms of vegetation canopy structure effects on energy balance models, [82] found out that, in the case of the unique canopy structure and configuration of vineyards, these SEB models were not able to adequately describe the physical processes driving evapotranspiration from vineyards.

In general, research results reported in the literature show that the surface aerodynamic temperature can be locally parameterized for different vegetation surfaces, surface homogeneity and environmental and climatological conditions, with relative success. However, researchers still use surface radiometric temperature images, from satellites, instead of the appropriate aerodynamic temperature because currently there is no single robust and comprehensive T_o model, for the application of the one-source energy balance of land surfaces model, that works for multiple surface cover types and environmental conditions. Thus, there is a need to study the mechanisms that generate the surface aerodynamic temperature, for different vegetation/crop surface types and environmental conditions.

2.3. Energy Balance Closure Flux Terms Frequently Not Considered

The hypothesis that the four terms described above close the surface energy budget is usually not fulfilled [26], with imbalances amounting typically between 10 and 40% depending on the terrain heterogeneity and complexity [83]. The most relevant motions seem to take place at the hectometre scale [84,85].

The main reasons for this lack of closure usually mentioned are the uncertainty of the measurements, the horizontal transport by advection, the unadequate computation of heat stored in the upper soil or the canopy and the lack of consideration of the energy used in the plant physiological processes for dense canopies. Several studies using drone or satellite-derived LST fields show that the advection of heat may be important [86,87,88] and well correlated with the imbalance.

Forcing the closure of the four-term surface budget equation to obtain LE, implies that the imbalance is considered entirely as latent heat flux, for which there is no clear evidence, although some studies point in that direction [89]. Instead most authors tend to distribute the observed imbalance between the latent and the sensible heat fluxes proportionally to the observed Bowen ratio (H/LE) [90]. In the rest of this work, we will consider that the surface energy budget is closed unless stated otherwise.

3. Remote Sensing LST-Based ET Algorithms

ET algorithms based on thermal remote sensing data can be divided into three categories: (i) single-source models, (ii) VI-TS Triangle/trapezoidal models, and (iii) dual (two)-source models.

3.1. Single-Source Energy Balance Models

3.1.1. Surface Energy Balance Index (SEBI) and Simplified-SEBI (S-SEBI) Models

Originally proposed by [91], the physically-based assumption is that ET varies with the T_s for a homogeneous surface. For a given surface albedo and boundary layer conditions, there exists a temperature gradient between the surface and the atmosphere (T_s−T_a) at some crop height, that ranges from a maximum value for which ET is assumed to be zero, to a minimum value for which the potential ET (ET_p) is taken, as given by the Penman-Monteith equation [7].

Figure 1 shows a complete scheme of the SEBI principles. There exists a correlation between T_s and the surface reflectance (albedo), where for low reflectance values T_s remains constant with the progressive increase of such reflectance. This phenomenon is attributed to presence of body water or well-watered surfaces, where all the available energy is used for the evaporation. From certain reflectance threshold T_s increases with the reflectance, inducing a decrease in the evaporation, in the so-called “evaporation controlled” range. From another higher reflectance threshold, T_s decreases with the increase of the reflectance, because the R_n decreases due to the high reflection of the energy received, a process known as “radiation controlled”.

Based on the previous assumptions, for a specific surface given Ts, [91] provided a value for ET by means of an interpolation expression:

$${\mathit{ET}}_i=\left(1-\frac{\frac{\left(T_s-T_{pbl}\right)}{r_{aH}}-\frac{\left(T_{s,min}-T_{pbl}\right)}{r_{aH,min}}}{\frac{\left(T_{s,max}-T_{pbl}\right)}{r_{aH,max}}-\frac{\left(T_{s,min}-T_{pbl}\right)}{r_{aH,min}}}\right){\mathit{ET}}_p=\left(1-SEBI\right){\mathit{ET}}_p$$

(18)

where T_pbl is the average air potential temperature at higher elevation or at the top of the Planetary Boundary Layer (pbl) and r_aH,min and r_aHmax are the minimum and maximum aerodynamic resistance in the sensible heat flux H [49].

A simplified version of the SEBI method (S-SEBI) was proposed by [92] with the advantage that no additional meteorological data are needed if the hydrological conditions of the surface are extreme enough to observe a significant range of gradient temperatures for the corresponding reflectance spectra.

Therefore if for a given α_s value the two extremes of T_s exist, that is for H_max (T_s,max, α_s) and LE_max (T_s,min, α_s), [92] expressed the ET value as:

$${\mathsf{\lambda }\mathit{ET}}_i=\frac{T_{s,max}-T_s}{T_{s,max}-T_{s,min}}\left(R_n-G\right)$$

(19)

The method is called simplified if the wet and dry pixel conditions are observed for a specific albedo value in a satellite scene. But [92] simplified even more Equation (19) by adjusting linearly the trends of maximum and minimum T_s with the albedo (see Figure 1), according the following equations:

$$T_{s,max}=a_{max}+b_{max}{\alpha }_s$$

(20)

$$T_{s,min}=a_{min}+b_{min}{\alpha }_s$$

(21)

where a_max and b_max are the slope and offset coefficients for the linear regression of the maximum temperature values with the albedo and a_min and b_min are the corresponding coefficients for the linear regression of the minimum temperatures, also with the albedo. These four parameters must be adjusted in each individual satellite scene. Including Equations (20) and (21) in Equation (19), ET expression from S-SEBI model can be defined then as:

$${\mathsf{\lambda }\mathit{ET}}_i=\frac{a_{max}+b_{max}{\alpha }_s-T_s}{a_{max}-a_{min}+\left(b_{max}-b_{min}\right){\alpha }_s}\left(R_n-G\right)$$

(22)

Very few modifications have been proposed in the literature for the S-SEBI model. One of them is made by [93], which in order to avoid using the LE and H fluxes, it is proposed using instead the evaporative fraction (EF), that can be extracted from the T_s-α scatterplot shape.

A recently published study [94] proposed an adaptation of EVapotranspiration Assessment from SPAce (EVASPA) tool [95] to S-SEBI model. This model, denoted as E3S (for EVASPA S-SEBI Sahel), proposes two relevant changes. The first one is related to the handling of the algorithms determining the wet and dry edges, since in some regions like the Sahel with strong seasonal contrast only one edge may be identified. This led to the introduction of twelve additional seasonally-dependent algorithms to estimate the dry and wet points, using an ensemble member weighting. Secondly a cloud edge filtering method was implemented, to avoid that pixel outliers can distort the T_s scatterplot shape. The cloud edge filtering was developed in two levels: (1) the first level filter was based on the quality flag provided by the MODIS LST products. Thus, all pixels bordering clouds and associated with an LST error >1 K were removed; (2) a second level filter consisted in extending the first level filter to cloud-bordering pixels with LST below the first LST quartile in the entire image. This quantile level selection could be site and climate-specific.

Table 1. Representative chronological compendium of validation studies of the S-SEBI method. Reference, spatial resolution (SR), number (N) of satellite scenes, Sensor Name, Uncertainty and Type of Surface are detailed (uncertainty statistics are averaged values of all type of surface studied). Statistical uncertainties on ET_d (mm/d) comprise Root Mean Square Error (RMSE) and bias.

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2005, [44]	20 m	19	Airborne	RMSE/bias = 1/0.2	corn, alfalfa, sunflower and wheat
2005, [96]	1100 m	520	AVHRR	RMSE = 1.2	forest
2011, [97]	90 m	11	ASTER	RMSE = 0.8	vineyard
2014, [45]	90 m	7	ASTER	RMSE/bias = 4/1	wheat, brocoli, frijoles chili, potatoes, chickpea, sunflower, orange and corn
2016, [98]	30 m	149	Landsat 5-TM Landsat 7-ETM	RMSE = 0.9	Citric Orchad, grazing, swamps, lakes
2017, [99]	30 m	19	Landsat 7-ETM Landsat 8-TIRS	RMSE = 0.9	sorghum
2017, [100]	1000 m	248	MODIS	RMSE/bias = 2/0.2	Pine trees, Green manure–weed– mustard (irrigated), Rice–rice (irrigated), Soybean–wheat (irrigated), Mixed crops (sugarcane, vegetables, turmeric, maize)
2019, [101]	30 m	52	Landsat 8-TIRS	RMSE = 1.12	Cropland (corn and muskmelon)

shows a review of the validation studies carried out from the S-SEBI proposal up to date. The validation studies are listed chronologically, indicating the sensor used and the spatial resolution. The uncertainties for each study are indicated in terms of Root Mean Square Error (RMSE), or bias. The corresponding type of surface (i.e., crop, forest, plateau, etc.) is also included in Table 1. The main advantage of the S-SEBI method is that no ground-based measurements are required. Instead, all the four input parameters (α_s, T_s, R_n and G) can be retrieved from satellite data (G is related with R_n, see Section 2.2). Nevertheless, a disadvantage is that extreme temperatures are mandatory and these are site-specific. Daily ET values for the S-SEBI ETi, are estimated in each study of Table 1, using Equations (5)–(10).

It is worth to note, that for Table 1,

Table 2. Representative chronological compendium of validation studies of the SEBS algorithm. Reference, spatial resolution (SR), number (N) of satellite scenes, Sensor Name, Uncertainty and Type of Surface are detailed (uncertainty statistics are averaged values of all type of surface studied). Statistical uncertainties are Root Mean Square Error (RMSE), bias, and Mean Absolute Difference (MAD), are shown for ET_d (mm/d). Different uncertainty rows inside the cell means uncertainty associated to the satellite sensor.

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2002, [46]	Tower footprint scale	620	Field Campaign Data	RMSE/MAD = 1 to 3/0.9 to 1.9	cotton shrub, brush grass
2005, [110]	30 m	1	Lansat 7-ETM	RMSE/bias = 1 to 3/0.3 to 0.6	corn soybean
2006, [14]	60 m	1	Landsat 7-ETM	bias = 0.2	corn, soybean
	90 m		ASTER	bias = 0.6
	1000 m		MODIS	bias = 2
2007, [111]	7–10 km	935	GSM 5-VISSR	RMSE = 5	Tibetan plateau
2007, [112]	1000 m	164	MODIS	RMSE/bias = 0.7 to 4/−1 to −0.17	Meadow, forest and corn
2010, [113]	1000 m	-	MODIS	RMSE/bias/MAD = 2/0.2 to 1.1/1.6	wheat, corn
2011, [114]	1000 m	33	MODIS	RMSE/bias = 3/1.7	wheat, corn
2014, [45]	90 m	7	ASTER	RMSE = 5	wheat, brocoli, frijoles chili, potatoes, chickpea sunflower, orange and corn
2015, [115]	1000 m	-	AATSR	RMSE = 1.2	varied non-especified surfaces (crop, forest, etc.)
2016, [98]	30 m	149	Landsat 5-TM Landsat 7-ETM	RMSE = 0.74	Citric Orchad, grazing, swamps, lakes
2017, [99]	30 m	19	Landsat 7-ETM Landsat 8-TIRS	RMSE = 1.1	sorghum
2018, [116]	30 m	11	Landsat 8-TIRS	RMSE/MAD = 0.22/0.21	double-cropped rice, peanut/sweet potato rotation, and orange groves
2019, [117]	30 m	22	Landsat 5-TM Landsat 7-ETM	RMSE = 1.8	wheat (dominant crop), barley or cotton in winter
2019, [101]	30 m	52	Landsat 8-TIRS	RMSE = 1.3	Cropland (corn and muskmelon)
2020, [118]	30 m	27	Landsat 8-TIRS	RMSE = 0.8	processing tomatoes and maize
2021, [119]	30 m	42	Landsat 7-ETM Landsat 8-TIRS	RMSE/bias/MAD = 1.3/−0.2/0.9	winter wheat

,

Table 3. Representative chronological compendium of validation studies of the SEBAL algorithm. Reference, spatial resolution (SR), number (N) of satellite scenes, Sensor Name, Uncertainty and Type of Surface are detailed (uncertainty statistics are averaged values of all type of surface studied). Statistical uncertainties Root Mean Square Error (RMSE), bias, and Mean Absolute Difference (MAD), are shown for ET_d (mm/d). Different uncertainty rows inside the cell means uncertainty associated to the satellite sensor.

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2005, [132]	90 m	1	ASTER	RMSE/bias = 1.3/0.04	corn, soybean
2007, [133]	6–12 m	5	airborne	RMSE/MAD = 2.4/2.0	shrubs, meadow wheat
2009, [134]	15 m	1	airborne	bias = 0.02	olive, vineyard and citric orchads
2009, [135]	30 m	10	Landsat 5-TM	RMSE = 1.2	Mango, vineyard, vegetation
			Landsat 7-ETM	RMSE = 0.4
2012, [122]	30 m	3	Landsat 5-TM	RMSE/bias = 1.9/−0.5	corn, soybean
			Landsat 7-ETM	RMSE/bias * = 1.4/−0.14
2012, [136]	1000 m	302	MODIS	RMSE = 0.5	wheat, corn, sunflower
2015, [137]	1000 m	7	MODIS	RMSE = 1.5	several non-especified crops
2016, [98]	30 m	149	Landsat 5-TM Landsat 7-ETM	RMSE = 0.8	Citric Orchad, grazing, swamps, lakes
2017, [99]	30 m	19	Landsat 7-ETM Landsat 8-TIRS	RMSE = 0.97	sorghum
2018, [116]	30 m	11	Landsat 8-TIRS	RMSE/MAD = 0.4/0.5	double-cropped rice, peanut/sweet potato rotation, and orange groves
2019, [101]	30 m	52	Landsat 8-TIRS	RMSE = 1.3	Cropland (corn and muskmelon)
2020, [118]	30 m	27	Landsat 8-TIRS	RMSE = 1.3	Almond, processing tomatoes and maize
2021, [119]	30 m	42	Landsat 7-ETM Landsat 8-TIRS	RMSE/bias/MAD = 1/−0.4/0.8	winter wheat

* Statistical uncertainties after applying the M-SEBAL method.

,

Table 4. Representative chronological compendium of validation studies of the METRIC algorithm. Reference, spatial resolution (SR), number (N) of satellite scenes, Sensor Name, Uncertainty and Type of Surface are detailed (uncertainty statistics are averaged values of all type of surface studied). Statistical uncertainties Root Mean Square Error (RMSE), and bias, are shown for ET_d (mm/d). Different uncertainty rows inside the cell means uncertainty associated to the satellite sensor.

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2009, [141]	60 m	2	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 2.4/0.2	corn, soybean
2009, [142]	120 m 60 m	2 1	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 0.6/−0.3	corn, soybean
2014, [143]	120 m 60 m	16	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 0.7/0.5	olive orchads
2015, [144]	30 m	12	Landsat 5-TM	RMSE/bias = 0.8/−0.1	cocoa, cotton, wheat, soybean
			MODIS	RMSE/bias * = 0.5/−0.3
2015, [145]	30 m	34	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 1.2/−0.3	cotton
2016, [98]	30 m	149	Landsat 5-TM Landsat 7-ETM	RMSE = 0.95	Citric Orchad, grazing, swamps, lakes
2017, [99]	30 m	19	Landsat 7-ETM Landsat 8-TIRS	RMSE = 1.5	sorghum
2019, [117]	30 m	22	Landsat 5-TM Landsat 7-ETM	RMSE = 1.6	wheat (dominant crop), barley or cotton in winter
2020, [118]	30 m	27	Landsat 8-TIRS	RMSE = 1.4	Almond, tomatoes and maize
2021, [119]	30 m	42	Landsat 7-ETM Landsat 8-TIRS	RMSE/bias/MAD = 1.2/0.4/0.9	winter wheat

* Statistical uncertainties derived from METRIC model after downscalling MODIS resolution to Landsat 30 m scale.

,

Table 5. Representative chronological compendium of validation studies of the T_s−VI Triangular/Trapezoidal methods. Reference, spatial resolution (SR), number (N) of satellite scenes, Sensor Name, Uncertainty and Type of Surface are detailed (uncertainty statistics are averaged values of all type of surface studied). Statistical uncertainties Root Mean Square Error (RMSE) and bias are shown for ET_d (mm/d).

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2001, [154]	1000 m	6	AVHRR	RMSE/bias = 3/0.3	crop, grass, forest
2003, [155]	1000 m	26	AVHRR	RMSE/bias = 1.6/025	forest, shrub, wheat, corn, soybean
2006, [156]	1000 m	15	AVHRR MODIS	RMSE/bias = 1.9/−0.5	wheat, cotton
2008, [149]	5000 m	123	SEVIRI	RMSE/bias = 1.4/−0.04	grazing
2009, [157]	1000 m	730	AVHRR	RMSE/bias = 1.2/−0.5	swamp, sugar cane water
2009, [141]	60 m	2	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 4/0.8	soybean, corn
2011, [114]	1000 m	33	MODIS	RMSE/bias = 3/−1.7	wheat, corn
2014, [158]	1000 m	52	MODIS	RMSE/MAD = 0.9/0.7	grazing, shrub
2017, [100]	1000 m	248	MODIS	RMSE/bias = 1.4/0.07	Pine trees, Green manure–weed– mustard (irrigated), Rice–rice (irrigated), Soybean–wheat (irrigated), Mixed crops (sugarcane, vegetables, turmeric, maize)
2019, [117]	1000 m	-	MODIS	RMSE/bias = 3/0.9 to 1.3	corn fields, sandy deserts, desert steppe, Gobi Desert, wetlands, and orchards
2021, [119])	30 m	42	Landsat 7-ETM Landsat 8-TIRS	RMSE/bias/MAD = 1.1/−0.4/0.8	winter wheat

,

Table 6. Representative chronological compendium of validation studies of the TSEB model. Reference, spatial resolution (SR), number (N) of satellite scenes, Sensor Name, Uncertainty and Type of Surface are detailed (uncertainty statistics are averaged values of all type of surface studied). Statistical uncertainties Root Mean Square Error (RMSE), bias and Mean Absolute Difference (MAD), are shown for ET_d (mm/d). Different uncertainty rows inside the cell means uncertainty associated to the satellite sensor.

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2005, [171]	60–120 m	3	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 1.7/1.4	corn, soybean
2005, [132]	90 m	1	ASTER	RMSE/bias = 3	corn, soybean
2007, [133]	6–12 m	5	airborne	RMSE/MAD = 2/1.9	shrub, grazing, wheat
2008, [172]	30–120 m	3	Landsat 5-TM	RMSE = 1.2	grass, shrub
2009, [141]	60 m	2	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 2.4/1.3	corn, soybean
2009, [142]	60–120 m	3	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 0.6/−0.06	corn, soybean
2009, [134]	15 m	1	airborne	bias = −0.04	olive, vineyard and citric orchads
2011, [114]	1000 m	33	MODIS	RMSE/bias = 1.7/0.4	wheat, corn
2012, [173]	30 m	2	Landsat 5-TM	RMSE/bias/MAD = 1.9/0.5/1.4	cotton
2012, [174]	30 m	3	Landsat series 5 and 7	RMSE/bias = 1.6/0.4	corn, soybean
	90 m	1	ASTER	RMSE/bias = 2.2/−0.17

and

Table 7. Same as Table VI but for studies dated from 2013 to 2021.

Year, Reference	SR	N	Sensor	Uncertainty (mm/d)	Type of Surface
2013, [175]	30 m	5	Landsat 5-TM Landsat 7-ETM	RMSE/bias/MAD = 1.4/−0.2/1.1	corn, soybean
2014, [176]	30 m	11	Landsat 5-TM Landsat 7-ETM	RMSE/bias/MAD = 1.5/−0.7/1.4	corn, cotton, soybean
2014, [45]	90 m	7	ASTER	RMSE = 4	wheat, brocoli, frijoles chili, potatoes, chickpea, sunflower, orange and corn
2015, [145]	30 m	34	Landsat 5-TM Landsat 7-ETM	RMSE/bias = 1/−0.9	cotton
2015, [177]	90 m	6	ASTER	RMSE/bias = 1.8/0.17	crops surrounded by desert
2016, [178]	0.06 m	10	Drone	RMSE/MAD = 1.8/1.5	Olives orchads
2016, [179]	30 m	22	Landsat 8-TIRS	RMSE/bias/MAD = 0.7/0.2/0.5	vineyard
2016, [180]	90 m	9	ASTER	RMSE/bias = 2.4/0.5	cron, vegetables and trees
2017, [181]	30 m	8	Landsat 7-ETM Landsat 8-TIRS	RMSE/bias/MAD = 0.9/−0.1/0.7	forest
2018 [160]	30 m	5	Landsat 5-TM	RMSE/bias = 0.2/−0.6	Olives orchads
2020, [182]	90 m	9	ASTER	RMSE/bias = 4/3	vegetables, maize and orchard
2021, [119]	30 m	42	Landsat 7-ETM Landsat 8-TIRS	RMSE/bias/MAD = 1.5/0.3/1.1	winter wheat

, sometimes the spatial resolutions for the Landsat series is 30 m. This is because authors are referring to the original scale of the VNIR sensor, not the original resolution of the thermal InfraRed (TIR) sensor of Landsat, which is 60 m (L7), 120 m (L5) and 100 m (L8−L9). The Landsat team offers to the user the TIR products of its sensors at 30 m resolution, by complementing the data from its Visible and Near InfraRed (VNIR) sensors.

3.1.2. Surface Energy Balance System (SEBS) Model

Proposed by [46] the SEBS method is based, as it is S-SEBI, SEBAL and METRIC (to be described below), on the same main concept as SEBI. The difference of SEBS from the rest is that it deals with estimating the ET fluxes at large scales, therefore considering heterogeneous surfaces.

Therefore, the assumption of constant atmospheric conditions for the whole scene is not applicable. Moreover, spatial variation of roughness height for different surfaces is allowed.

The SEBS is an extended model that defines the surface roughness length in the heat fluxes, reformulating the EF at the limiting cases. To this end in the SEBS theory the dry and wet case limits are:

Dry limit:

$${\mathrm{LE}}_{dry}=0$$

(23)

$${\mathrm{H}}_{dry}=R_n-G$$

(24)

Wet limit:

$${\mathrm{LE}}_{wet}=R_n-G-H_{wet}$$

(25)

$${\mathrm{H}}_{wet}=\frac{R_n-G-\frac{\rho C_{P}VPD}{r_{aH}\gamma }}{\left(1-\frac{\mathsf{\Delta }}{\gamma }\right)}$$

(26)

where VPD is the vapor pressure deficit (difference between saturated and non-saturated vapor pressure), Δ is the slope of the curve relating temperature difference (T_s-T_a) with VPD, and γ is the is the psychrometric constant.

As an intermediate step SEBS method defines a relative EF (EF_r) as:

$$E{F}_r=1-\frac{H-H_{wet}}{H_{dry}-H_{wet}}$$

(27)

Once all the required parameters have been estimated, LE is expressed in SEBS as:

$$LE=E{F}_r·{\mathit{LE}}_{wet}$$

(28)

Some studies presented a modification of the SEBS inner structure in terms of: (i) momentum roughness length and zero plane displacement; (ii) an empirical relationship between land surface temperature and air temperature; and (iii) designing a new input interface of meteorological field at reference, with the objective of representing land surface and climatic changes of a specific region. These studies adapted SEBS to be applicable in different countries like Taiwan (SEBS-Taiwan, [102]), China (SEBS-China, [103]) or Iran (SEBS-Iran, [104]). SEBS was also modified at local scale (SEBS-Urban, [105]) after including the anthropogenic heat flux in the SEB equation to estimate the ET of Beijing city. This study hypothesized that ET was equal to the water consumption in the study area.

An improvement of SEBS method was suggested by [106] by correcting the excess resistance (Equation (15)) with a scale factor dependent on soil moisture (SM), following:

$${\left(\frac{\kappa }{B}\right)}_{corr}=\left[a+\frac{1}{1+exp\left(b-cS{M}_{rel}\right)}\right]\frac{\kappa }{B}$$

(29)

where a, b and c are coefficients of the sigmoid function and SM_rel is a relative soil moisture defined as:

$$S{M}_{rel}=1-\frac{SM-S{M}_{min}}{S{M}_{max}-S{M}_{min}}$$

(30)

SM_max and SM_min values can be obtained using a time series analysis of the soil moisture data on annual or long-term basis [106]. Wu et al. [107] suggested the use of the Normal Differential Water Index (NDWI) instead of SM_rel in Equation (30). Yi et al. [108] also recommended replacing SM_rel by the inverse of the modified perpendicular drought index (MPDI).

It is also worth mentioning the study of Chen et al. [109], in which the excess resistance is internally modified to be robust for heat flux calculation for high and low canopies, since it has been proved its depends on environmental conditions, surface types, canopy structure and vegetation surface geometry. kB⁻¹ structure was adjusted from a one-foliage layer to a multi-foliage layer, which provides the possibility of including the impact of vertical variations in foliage leaf area density, foliage shelter factor, and foliage heat transfer coefficient.

Table 2 shows a review of the validation studies carried out for the SEBS model up to date. The main advantages of the SEBS are that new limit cases in the SEB reduce uncertainty from T_s and meteorological variables, SEBS introduces a new formulation for the roughness length for heat transfer instead of using fixed values, and a priori knowledge of the actual turbulent heat fluxes is not required. However, the disadvantages are that SEBS requires too many parameters and offers relatively complex solution for the turbulent heat fluxes, which can be a source of inconveniences when data are not available.

3.1.3. Surface Energy Balance Algorithm for Land (SEBAL) Model

Proposed by [120], SEBAL is an iterative and feedback-based numerical procedure that deduces the radiation, heat and evaporation fluxes. SEBAL estimates LE with the residual of Equation (1) in a pixel by pixel satellite scene, just considering two anchor points (dry and wet points) for the full scene. At the dry point LE = 0 and H = R_n-G, and at the wet point H = 0 and LE can be considered either R_n-G or a reference ET [7].

The main consideration of SEBAL is to establish a linear regression between T_s and difference of T_s with the air temperature at the crop height z (δT_a), as follows:

$$\delta T_a=c_1{T}_s-c_0$$

(31)

where c₁ and c₀ are the empirical coefficients derived from anchored conditions of wet (δT_a,wet = 0, T_s,wet) and dry points (T_s,dry), for which the temperature difference is expressed as:

$$\delta T_{a,dry}=\frac{(R_n-G)r_{aH}}{\rho C_P}$$

(32)

These coefficients must be estimated for each individual satellite scene, independently of the surface type composing the image. Once the regressed equation, for each T_s value assigned to each pixel of the scene is retrieved, the H is estimated with the expression:

$$H(T_s)=\frac{\rho C_P}{r_{aH}}\left(c_1{T}_s-c_0\right)$$

(33)

A first modification of SEBAL was proposed by [121] after introducing Digital elevation model (DEM) data to account for impacts of slope aspect on incident solar radiation. But the main subject of modification of SEBAL model is related with the choice of endmembers or anchor points (i.e., “hot” and “wet” points in the satellite scene). A Modified version of SEBAL (M-SEBAL) was proposed by [122], that based the election of the anchor (dry/wet) points on the type of surface, defined by the FVC. They adjusted Equation (31) for every specific surface defined by FVC and subsequently estimated H and LE. Feng et al. [123] derived geo-mathematically four reference dry and wet limits, instead of the classical option of two anchor points.

Another modification is addressed in [124], where authors deal with the endmember pixel selection process, proposing an automated procedure by introducing an exhaustive search algorithm (ESA), in which all possible candidates are systematically considered and checked to determine if the solution is satisfied, based on two steps: (1) create a binary map of candidate pixels that are identified using a simple rule-based classifier, and (2) apply ESA to identify hot and cold pixels that meet the defined criteria. In this sense, [125] created the automated SEBAL (ASEBAL) model that allows automatic selection of the endmember pixels according to criteria based on the variation in biophysical parameter values retrieved from satellite products.

It is also worth mentioning the study of [126], where authors found more appropriate to combine the capabilities of Synthetic Aperture Radar (SAR) with optical remote sensing technique in the applications where soil moisture conditions to locate anchor pixels, since SAR is more suitable to detect soil moisture conditions.

Two interesting studies that suggested modifications of SEBAL ET model were related with changes in the SEB equation, by introducing new fluxes terms that must be accounted in certain situations. Mkhwanazi et al. [127] proposed a modified version of SEBAL model (SEBAL-A), that is capable of accounting for advection flux term, estimated with a wind function which included a daily mean or extreme temperatures parametrization, which is incorporated it into the original SEBAL daily ET sub-model. [128] also suggested enhancing the SEBAL model when it is applied in a heterogeneous urban environment (uSEBAL), after considering urban land surface parameters and including the urban anthropogenic heat flux in the SEB budget.

Finally, it is worth to note that recent open access software code versions of SEBAL are accessible. One of them is a new python version of SEBAL (pySEBAL), which incorporates an automated pixel selection procedure and is currently under development and testing at the IHE-Delft Institute [129]. In addition, a new tool based on the SEBAL algorithm, modified with a simplified version of the Calibration using Inverse Modeling at Extreme Conditions process [130] for the endmembers selection, was developed within the Google Earth Engine (geeSEBAL) environment [131].

Table 3 shows a review of the validation studies carried out by the SEBAL model up to date. The main advantages of SEBAL are the minimum ground data required, the unnecessary atmospheric correction thanks to the internal calibration. However, election of dry/wet anchor points is subjective, SEBAL applications over non-flat areas and uncertainties on T_s estimates have a great impact on H retrievals.

3.1.4. Mapping Evapotranspiration at High Resolution and with Internalized Calibration (METRIC) Model

METRIC was proposed by [54], based almost completely on SEBAL method. The main differences with SEBAL are that R_n and T_s are corrected topographically using a Digital Elevation Model. Both methods differ in the calculation of H as in METRIC the temperature gradient of the wet point (Equation (32)) is expressed as:

$$\delta T_{a,dry}=\frac{(R_n-G-1.05E{T}_r)r_{aH}}{\rho C_P}$$

(34)

The assumption that ET rate (ET_r) is about 5% above that of the standard ET_r in the METRIC model is assumed because, for a large population of fields, some fields will have a wet soil surface beneath a full vegetation canopy that will tend to increase the total ET rate [54].

Since METRIC is directly based on SEBAL model, all the modifications previously discussed in Section 3.1.3 can also apply to METRIC. However, it is interesting to note a recent study [138], where calibration variability in METRIC is examined through the anchor pixels selection. Following METRIC original procedure to detect “cold” and “hot” pixel in a scene, the authors detected a range of different “cold” and “hot” pixels that conformed a total of five sets of calibration pixel pairs (based on combination of maximum and minimum “cold” and “hot”, as well as values closest to the flux station). So, different slopes and intercept values were retrieved from the plot δT vs. Ts (Equation (31)). In summary, authors observed significant spatial differences in ET across the study area, which suggested that candidate calibration pixels cannot be completely quantified using only NDVI and T_s. Therefore, it can be presumed that the large differences in δT among the candidate pixels were not random noise, and that the inclusion of available energy into the calibration pixel selection process could potentially improve the calibration process and provide better estimates of ET.

Is also worth to note that METRIC results can be directly obtained at any region of the world with the new Earth Engine Evapotranspiration Flux (EEFlux) portal [139,140]. Table 4 shows a review of the validation studies carried out for the METRIC model up to date. The main advantages and drawbacks of METRIC are the same as for SEBAL, Highlighting that the main advantage over SEBAL is its application over mountainous areas.

3.2. Land Surface Temperatures–Vegetation Index Triangle/Trapezoidal ET Models

Vegetation Indices (VIs) offer important information from T_s−VI triangle/trapezoidal relationship. These diagrams can be useful for estimating regional ET values under the condition that full ranges of SM content and vegetation are present in the remote sensing scene.

The Triangle method was described by [146], and provides a biophysical justification for the combination of T and vegetation indices (VI) time series for land-cover mapping and land-cover change analysis. Figure 2a shows the simplified T_s−VI triangle diagram where the so called “dry edge” is defined by a decrease of the T_s with the increase of the VI, this edge is defined under the assumption that highest T_s for different types of surfaces (different VIs) is associated to points with unavailability of SM content. The lower envelope represents the wet edge, which is nearly constant and represents the lowest T_s for the different types of surfaces, assuming here that the ET is potential. Theoretically the two envelopes (dry and wet edges) intersect in a point of full cover vegetation and maximum ET.

Based on the Priestley-Taylor formulation [147] and fully remotely sensed data, an expression proposed by [148] allows estimating the LE according to the T_s−VI triangular relationship as:

$$LE=\phi \left[\left(R_n-G\right)\frac{\mathsf{\Delta }}{\mathsf{\Delta }+\gamma }\right]$$

(35)

where φ varies from 0 to 1.26, and it is estimated with a two-step linear interpolation. This interpolation consists in: (1) once set a φ_min = 0 (driest bare soil) and φ_max = 1.26 (largest VI and lowest T_s), for each intermediate VI a φ_min,i is estimated by interpolation between φ_min and φ_max and a φ_max,i is assigned to the minimum T_s value of such VI (usually 1.26), (2) Different φ_i can be estimated from the T_s range at each VI. VI can be either the NDVI, FVC, the Soil Adjusted Vegetation Index (SAVI), Enhanced Vegetation Index (EVI), etc.

The T_s−VI Trapezoidal method based on the proposed Water Deficit Index (WDI) by [70], an index related to the ratio of actual and potential evapotranspiration, which evaluates the evapotranspiration rates of both full-cover and partially vegetated sites. WDI based T_s−VI Trapezoidal method considers that full covered vegetation can both water-stressed and well-watered in a satellite scene, so there is no intersect in a point of full cover vegetation and the T_s−VI diagram become a trapezoid (Figure 2b). Moran et al. [70], extended the ET estimation from full to partially vegetated surface areas, considering that temperature gradient (T_s−T_a) varies linearly with the vegetation status. So, WDI is linearly related at each VI to the maximum (dry bare soil) and minimum (probably well-watered full-cover vegetation) gradient temperature, taking values between 0 and 1 and LE is calculated similar to Equation (19) as the ratio of potential ET calculated from Penman-Monteith equation [7], LE for the trapezoidal method is expressed as:

$$\mathit{LE}=\left(1-WDI\right)L{E}_P=\left(1-\frac{\left[{\left(T_s-T_a\right)}_{min}-{\left(T_s-T_a\right)}_i\right]}{\left[{\left(T_s-T_a\right)}_{min}-{\left(T_s-T_a\right)}_{max}\right]}\right)L{E}_P$$

(36)

Stinsen et al. [149], combined the triangle method with thermal inertia information obtained from MSG-SEVIRI sensor, by changing the surface temperature over time. Authors also applied a slightly modified version of the two-step interpolation scheme to obtain φ, based on a interpolation where the decomposition is regarded as non-linear with the argument that non-linear intersections are observed in the dry edge and in the iso-lines of equal moisture availability of the T_s−VI tringle feature. Adopting equations from this study, [150] proposed a Time-Domain Triangle Method (TDTM), based on a pixel by pixel T_s–VI feature space, obtained by exploring the temporal domains over each single image pixel, since the triangle contains the whole climatic variability of the T_s–VI. Also, in [151], an algorithm to determine quantitatively the dry and wet edges for the T_s–VI triangular space in arid and semi-arid areas, where wet pixels are not generally easily identified, was developed with an iterative process that automatically retrieve the two edges in the T_s–FVC triangle space in ten steps. The process starts dividing the range of FVC in the T_s–FVC triangle space into M intervals evenly and ends with a linear regression to obtain the dry edge.

Another interesting modification of the T_s−VI triangle method is proposed by [152] that estimate daily ET directly using the top of atmosphere (TOA) radiances without performing atmospheric correction or other processes, to avoid the uncertainty associated with the estimation of ET, attributed to satellite data products. Furthermore, a recent study [153] estimates the soil heat flux (G) in the trapezoidal space of Ts-FVC, prior to apply the original T_s−VI trapezoidal method afterwards.

Table 5 shows a review of the validation studies carried out by the T_s−VI Triangular/Trapezoidal methods up to date. The main advantage of Triangular method is that no atmospheric or ground data are required, but the drawback is that large numbers of pixels with wide range of SM and FVC are required to assure dry and wet limits, a subjective procedure. For the Trapezoidal method one advantage over triangular version is that a large number of pixels is not required since intermediate values are determined by four limiting vertices. But a drawback is that ground measured data (vapor pressure, T_a, u, stomatal resistances, etc.) are required. A major disadvantage for both Triangular/Trapezoidal methods is that they do not detect immediately the water stress, since it does not have instantaneous impact on vegetation cover.

3.3. Dual (Two)-Source Energy Balance (TSEB) Model

A two-layer model of turbulent exchange that includes the view geometry associated with directional radiometric surface temperature was developed and evaluated by [43]. The main assumption of TSEB model is that T_s estimated from remote sensing measurements is the composite of soil and canopy surfaces, and it can be expressed as:

$$T_s=\left[FVC{T}_c^4+\left(1-FVC\right)T_{soil}^4\right]$$

(37)

where T_c is the temperature of the canopy and T_soil is the temperature of the non-vegetated soil. Note the term FVC in this Equation (37) is a function of LAI and sensor viewing angle [44]. The original TSEB does not consider the weighting by the emissivities of each component, but some authors have claimed for its importance [132,159,160].

The key of the TSEB model is the partition of sensible heat flux into the canopy and soil layers. If we assume that there is an interaction between the fluxes of both components, due to a heating of the in-canopy air by heat transport from the soil, the resistance framework in TSEB can be considered to be in series [43,161,162,163]. A parallel configuration was also introduced [43], based on the assumption that the resistance of each source interacts independently with respect to the surface-atmosphere, although the canopy still has an effect on the radiation transmitted towards the soil and the wind attenuation through the canopy to the soil surface. Schemes in Figure 3 shows both the series and parallel configurations [164]. Each component has its own SEB equation expressed as:

$$R_{n,soil}=H_{soil}+L{E}_{soil}+G$$

(38)

$$R_{n,c}=H_c+L{E}_c$$

(39)

A variety of radiative transfer models have been implemented within TSEB to estimate R_n components. Traditional versions of TSEB estimate R_n using the physically-based radiative transfer model described in [164], first implemented within TSEB by [165]:

In this revised two-source model, the separation of net radiation between canopy and soil can be expressed as:

$$R_{nc}=L_{nc}+\left(1-{\tau }_s\right)\left(1-{\alpha }_c\right)S$$

(40)

$$R_{ns}=L_{ns}+{\tau }_s\left(1-{\alpha }_s\right)S$$

(41)

where S is solar radiation, ${\tau }_s$ is solar transmittance in the canopy, ${\alpha }_s$ is soil albedo, ${\alpha }_c$ is canopy albedo, L_nS and L_nC are longwave radiation for soil and canopy, respectively, and are estimated via the following expression:

$$L_{ns}=\exp (-{\kappa }_L\mathsf{\Omega }LAI)L_{sky}+\left(1-\exp (-{\kappa }_L\mathsf{\Omega }LAI\right))L_C-L_s$$

(42)

$$L_{nc}=\left(1-\exp (-{\kappa }_L\mathsf{\Omega }LAI\right))\left(L_{sky}+L_s-2L_c\right)$$

(43)

where L_Sky, L_S, and L_C are longwave radiation from the sky, soil, and canopy, and can be calculated from air, soil, and canopy temperature, Ω is the clumping factor, and κ_L is an extinction coefficient. Campbell and Norman [164] described how to estimate the transmittance, albedo, and extinction coefficients for soil and canopy.

Soil and canopy sensible heat fluxes H of Equations (38) and (39), can be deduced easily from Figure 3 as:

$$H_{soil}=\rho C_P\frac{T_{soil}-T_a}{r_{aH}+r_s}$$

(44)

$$H_c=\rho C_P\frac{T_c-T_a}{r_{aH}}$$

(45)

where the parameter r_s is the resistance to the heat flow in the boundary layer immediately above the soil surface, and it can be expressed as:

$$r_s=\frac{1}{0.0025{\left(T_s-T_c\right)}^{1/3}+0.012u_s}$$

(46)

where u_s is the wind speed at a height above the soil surface where the effect of soil surface roughness is minimal (i.e., 0.05–0.1 m).

Latent heat fluxes for soil (LE_s) and canopy (LE_c) are calculated as residuals of Equations (44) and (45).

Another type of two-source model formulation is the so-called patch model where it is assumed that all the fluxes act vertically and that there is no interaction between soil and canopy components (i.e., a complete energy balance between the atmosphere and each element [165,166]).

There exist modifications on the TSEB model, like the Two-Source Time-Integrated Model (TSTIM) proposed by [53], lately called Atmosphere-Land Exchange Inverse (ALEXI) in [167], which relates time-differential T_s observed (1.5–5.5 h after sunrise) from geostationary satellites to the time-integrated energy balance within the surface-atmospheric boundary layer system.

ALEXI has minimal reliance on absolute (instantaneous) air or surface temperature input data, and therefore provides a relatively robust flux determination at the coarse geostationary pixel scale (5–10 km). For finer scale ET applications, ALEXI flux fields can be spatially disaggregated using higher resolution T_s information from polar orbiting systems (e.g., Landsat or MODIS) in an algorithm referred to Disaggregated ALEXI (DisALEXI) model [168].

Another modified TSEB model is the Simplified TSEB (STSEB) model proposed by [169], which suggests dividing the R_n, H and LE fluxes in soil and canopy weighted values, according to the FVC. The scheme proposed by [170], which is based on canopy conductance by means of a formulation that has sound parameterizations of transpiration and evaporation, that links soil surface temperature and soil evaporation eliminating the need of using soil moisture as input, and allowing the model predicting surface temperature and energy fluxes simultaneously. Furthermore, it can be calibrated using remotely sensed LST data.

Table 6 and Table 7 show a list of the validation studies carried out using TSEB model up to date. TSEB is obviously a very well-accepted model, based on the amount of validation studies. The main advantages of TSEB model are: angular effects are into consideration, no “resistance excess” correction is required, atmospheric and emissivity induced errors on T_s are reduced (case of ALEXI/DisALEXI model).

In general terms, the validation studies of all the ET models shown in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6 and Table 7 establish the estimation error of the daily ET between 1–5 mm/d, with an average error of around 1 mm/d regardless of the sensor spatial resolution or the type of surface under study. However, it can be seen for ET models such as S-SEBI, SEBS, LST-VI Triangle/Trapezoidal or Two-Source, that at high spatial resolutions (60–100 m) the error is smaller (around 1 mm/d), and around 2 mm/d for kilometric resolutions. Regarding the one-source ET SEBAL and METRIC models, the validation studies have shown the smallest errors, with values ranging between 0.4–2.4 mm/d, although most of these studies have been carried out with high spatial resolution sensors. Finally, it cannot be stated that the error for all models has decreased over the years.

This review has shown, in line with conclusions from the recent Carlson’s essay [183], that mathematical models of natural land-living systems applied to satellite-based measurement systems, will never be able of full-describing and predicting all relevant processes related to the functioning of such systems, insomuch as these models also require knowledge of plants and soils. However, such models can be useful, in describing key driving variables under limited environment conditions for plant growth and development, and surface energy balance at spatial resolutions ranging from tens of meters, up to hundreds of kilometers.

4. Limitations, Current Trends and Future Prospects

4.1. Toward High Spatiotemporal Resolution ET Retrievals

Water resource management applications require spatial information about daily and seasonal ET at pixel scales resolving individual agricultural fields or even finer management zones within fields. Such applications include irrigation scheduling, monitoring compliance with water rights, negotiating water rights transfers, sales or leases, predicting changes in regional water use due to land-use and/or climate changes, and studying the impacts of agriculture on surrounding ecosystems and stream flow [173]. T_s is the most influential input parameter in deriving the energy fluxes, because it provides an indication of the net effect of land-atmosphere interactions, which results in the process of evapotranspiration. However, coarser-scale TIR imaging devices (e.g., MODIS) can provide T_s data for ET mapping on a near-daily basis for relatively large fields, these data being too coarse to resolve water use at the scale of individual users in most irrigation districts. On the other hand, some polar orbiting satellites (e.g., Landsat series, ASTER) are the only platforms that provide routine, global thermal imagery at scales (60–120 m) that resolve water use patterns over heterogeneous agricultural areas, but overpass frequency (16 days) is not optimal for ET monitoring. Hence, a desirable alternative technique is to derive high resolution T_s data using the fine resolution VIS/NIR/TIR bands that are available from either the same or different satellite systems, at daily frequencies of temporal resolution.

A strong inverse relationship exists between T_s and VIs (see Section 3.2.), and VIs can be derived from VIS/NIR data at high spatial resolution (60–120 m). Such a relationship can be exploited to improve the spatial resolution of TIR band imagery (at a coarser resolution of 1–10 km) than the one associated to VIS/NIR band instruments. This relationship reflects the fact that denser vegetation cover tends to be correlated with lower surface temperatures due to cooling by transpiration.

Generally, most methods first imply a coincident overpass of two (high and coarse) spatial resolution satellite imager (e.g., MODIS-Landsat/ASTER), both operating in the VIS/NIR/TIR regions. Second, the methods are based on a statistical regression between VI (e.g., NDVI, SAVI, EVI, Modified SAVI (MSAVI), Normalized Difference Moisture Index (NDMI), NDWI), etc.), which is aggregated to the original resolution of the T_s (e.g., Landsat NDVI at 30 m is aggregated to 960 m of the original MODIS T_s) and the actual T_s. The resulting regression is used to disaggregate the coarse resolution T_s into fine resolution T_s with the VI retrieved for other dates. The fundamental assumption in this method of disaggregation is that the T_s–VI relationship remains continual at multiple spatial resolutions.

Techniques and tools to disaggregate either T_s, VI or ET are varied. So, the pioneering work of [184] used the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [185], which uses comparisons of one or more pairs of observed Landsat/MODIS maps, collected on the same day, to predict maps at Landsat-scale on other MODIS observation dates.

To estimate from ALEXI/DisALEXI application, daily ET values at MODIS (1 km) and Landsat (30 m) scales, with continuous data retrieved from geostationary GOES sensor. Such STARFM- ALEXI/DisALEXI combined algorithm has been proven successful in successive studies [173,175,176]. Other techniques like the hybrid use of the hot edge model and artificial neural network [186], machine learning approaches [187] or data mining [188] are also valuable disaggregation alternatives. However, temperature variations due to subpixel moisture variability cannot be recovered unless they are well-correlated with vegetation features. Therefore, dominant moisture variations must be well-resolved at the native TIR resolution. These results suggest that while thermal sharpening is a valuable tool for enhancing spatial information content in TIR imagery, it does not replace the need for TIR data collection at the sub-field scale. Other differences between both sensors exist due to several effects such as acquisitions at different time, spectral resolution, atmospheric correction, viewing angle, and pixel footprint that have to be considered and resolved, usually by normalization procedures.

An actual prospect which has been recently addressed in [189,190,191,192,193], is the synergistic application of disaggregation techniques to retrieve ET values at high spatiotemporal resolution, using two different satellites, one of which is carrying a sensor without thermal bands at coarser (~ 1 km) spatial resolution, and the other carrying a VIS-NIR sensor (i.e., MODIS- Spot-5), mapping Earth surfaces at finer spatial resolution (~10–30 m) and with a revisit of few days (e.g., 5 days for Sentiel-2, combining platforms A and B).

In this context, the ESA has funded the Sentinels for Evapotranspiration (Sen-ET) project, which main objective was to develop an open-source software application for accurate, pre-operational modelling of instantaneous evapotranspiration, by means of an optimal methodology for estimating ET at tens of meters spatial scale, based on synergistic use of Sentinel 2 and Sentinel 3 satellites observations. Also, the continuing Increasing Crop Water Use Efficiency at Multiple Scales Using Sentinel Evapotranspiration (ET4FAO) project, will implement further refinements to the Sen-ET method. For further information readers are referred to Guzinski et al. [194].

4.2. Accurate and Spatial Representative Field Instrumentation

Common ET measuring systems include water balance, lysimeters, Bowen ratio, eddy covariance, or sap flow [5]. Inherent to all these methods is the reality that the instrumental uncertainties and the difficulties to set a fully consistent experimental setup to determine the energy and water budgets can lead to water use estimates far from reality.

Instrumentations like water balance, lysimeters, Bowen ratio or sap flow are point specific and could not represent the ET of medium (100 m) to coarse (1–10 km) scales estimated from satellite data. The Eddy-Covariance systems are typically representative of an area surrounding them of a few hundred meters of radius [195], but in the case of terrain heterogeneities of smaller size their representativity is questionable, as local circulations at the hectometer may advect significant amounts of water vapour [196]. New Large Aperture Scintillometers (LAS) devices provide an approach to validate ET at much larger scales. LAS estimate the sensible flux based on measurements of small fluctuations in the refractive index of air caused by temperature, humidity, and pressure induced variations in density. However, according to [6] all the above-mentioned instruments showed average typical errors of: 10% (lysimetry), 15% (Bowen Ratio), 20% (soil water balance), 25% (Eddy covariance and LAS) and 35% (sap flow).

Eddy covariance and LAS are well-accepted techniques to validate ET estimates at medium to coarse (100–1000 m) and large (up to 10 km) scales, because of a good ratio error/advantage. But accurate instrumentation is required for ET validations at different scales. In this sense Laser LAS devices [26,197] are lately being used to ensure residual ET estimations (through accurate H estimations) at high and medium scales (10–100 m), reducing the problem of the SEB non-closure induced by the Eddy correlation [198].

The use of unmanned aerial vehicles (UAVs) has given a boost to precision agriculture, offering an opportunity to provide multi-spectral imagery at the highest spatial resolution (few cm per pixel, [87]). UAVs are portable and cost-effective, and flight procedures over crop fields helps increasing the production. New VIS/NIR/TIR sensors onboard UAVs have the ability of observing the agricultural fields with great detail. UAV images can produce ortho-mosaics and 3D information products to be used in mapping structure and volume of plant canopy, which is directly correlated to plant biomass and LAI, [199]. A significant number of studies recently published [200,201,202,203,204,205,206] have taken advantage of the good performance of UAVs, applying remote sensing techniques to estimates the crop ET. It must be emphasized the care to be taken, when applying ET models to thermal sensors on board UAVs, since it is well-known that uncooled commercial thermal cameras induce very significant errors in LST, which is key in the application of energy balance models [87].

4.3. Operational ET Products

Historically, the MODIS evapotranspiration products (MO/YD16) are the first 0.5 km operational land surface ET database for the global vegetated land surface at 8-day, monthly, and annual intervals [207]. This product is based on the Penman-Monteith approach [7], not on thermal data. This non-TIR based product has been widely evaluated against ground measurements, showing errors around 0.31 mm/d, or 24%, using data from 46 eddy covariance flux towers [208]. Other important non-thermal-based datasets products are the FLUXCOM [209], WaPOR [210], or the European Centre for Medium-Range Weather Forecasts (ECMWF) (GLDAS) ET images [211]. However, some initiatives assimilate the land surface temperature in their ET products. This is the case of the ECOSTRESS ALEXI/JPL ET [212,213], or the EUropean Organisation for the Exploitation of METeorological SATellites (EUMETSAT) Satellite Application Facility on Land Surface Analysis (LSA-SAF), ET product from Meteosat Second Generation (MSG) [214]. The comparison of these operational ET products with 15 eddy covariance flux towers located in Europe showed the good performance of the EUMETSAT over the MODIS ET [208,215,216,217] products [194].

The Global Land Evaporation Amsterdam Model (GLEAM) [218] is the only global ET operational product designed to be driven by remote sensing observations derived from microwave sensors, including soil moisture and vegetation optical depth. To assimilate microwave soil moisture observations—typically sensitive to the first few cm of the soil—GLEAM uses a simple Newtonian nudging algorithm, which minimizes the computational demands [219]. All output datasets have a daily resolution and a common 0.25° global grid.

As previously commented on Section 3.1.4., METRIC ET products can also be directly downloaded, for any region of the world, with the new Earth Engine Evapotranspiration Flux (EEFlux) portal [140,220], that includes a tool that processes ET maps from input images of Landsat series 7 and 8, i.e., ET maps at the spatial resolution of 30 m every 16 days. Initiatives like online platform (OpenET, [221]) for mapping daily, monthly and yearly ET values at the scale of individual fields (i.e., 30 m). This webtool offers easily accessible satellite-based ET across the western US. Also missions like GRAPEX [222] tries to apply a multi-scale remote sensing ET toolkit for mapping crop water use and crop stress for improved irrigation scheduling and water management in vineyards in the Central Valley of California, a region of endemic periodic drought.

Despite the enhancement of operational ET products, there is still a demand of higher resolution and more accurate operational ET products. As pointed out in Section 4.1., water resource management applications require spatial information about daily and seasonal ET at pixel scales of the hundreds of meters or lower. Operational high spatiotemporal ET products are being demanded specially from the agricultural sector and this is the goal to be achieved in the near future.

4.4. New Advent of Orbiting Sensors with Improved Spatial and Temporal Features

Since its beginning nearly four decades ago, the field of Remote Sensing has experienced an exponential increase in terms of satellites, products and users. However, the advance of the satellite engineering is always a compromise between high spatial/low temporal resolution (or vice versa) of the scenes acquired by the imager onboard. Fortunately, these advances tend to converge to instrument features that attend the actual requirements of remote sensing users. So, future Hyperspectral Infrared Imager (HyspIRI) [223], will combine near and shortwave infrared imaging spectrometer with a multi-channel TIR radiometer, designed to have a spatial resolution of 60 m and a revisit time of 5 days. Another future micro-satellite is the MIcro Satellite for Thermal Infrared GRound Surface Imaging (MISTIGRI) [224], which intends to estimate daily the TIR radiance from surface at resolutions of 50 m, also with a revisit period of 5 days. Currently is orbiting and acquiring data the ninth Landsat satellite [225], which is carrying, as well as its immediate predecessor Landsat 8, two sensors on-board, Operational Land Imager (OLI) detecting the VIS-NIR spectral reflectance at 30 m resolution and TIRS-2, measuring the spectral thermal emittance at 100 m. Landsat 9 revisit time is 16 days, like its predecessors.

In parallel, work is being carried out on the development of increasingly higher resolution multispectral thermal sensors, such as the ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) mission [226], the Land Surface Temperature Monitoring (LSTM) mission [227] or the Thermal infraRed Imaging Satellite for High-resolution Natural resource Assessment (THRISNA) mission of Centre national d’études spatiales- Indian Space Research Organisation (CNES-ISRO) [228], which will allow further progress in the applications of this parameter in the estimation of crop water requirements, among other applications.

Recently, the U.S. Geological Survey- National Aeronautics and Space Administration (USGS-NASA) consortium has launched an operational LST product for the entire Landsat time series globally, which will also apply to the new TIRS-2 sensor on board Landsat 9, thus continuing the entire Landsat legacy beyond 2030.

4.5. SEB Closure and Impact of Advection Term

All the remote sensing-based ET models described in here, presume that the Equation (1) is undoubtedly true and the SEB closure takes place as R_n is used only as the source of the energy fluxes H, G and LE. The analysis of the surface energy budget with eddy-covariance measurements has made evident that this presumed closure seldom takes place, instrumentation error aside, and there are other terms that at least at some time of the day, it must be taken seriously into consideration. Cuxart et al. [83] indicated that the four mains traditionally terms would be only balanced if no other thermal sources or sinks existed, implying no thermal advection and surface homogeneity. Therefore, once surface heterogeneity comes into the picture, thermal advection will take place and the usual 4-term energy budget will not close anymore, and the corresponding imbalance is the sum of several other terms like the heat storage in the volume, the biological thermal exchanges, the temperature tendency, the thermal advection and any other unaccounted factor [88].

In a recent study [87] it was observed from high resolution (2 m) T_s maps obtained from a TIR camera assembled in a drone, that after up-scaling the original resolution to scales of 70–200 m, the advection term can explain an important amount of the imbalance occurred at midday in a SEB station. Non-closure of SEB equation is an important issue historically left aside, but new tools and remote sensing sensors commit to ET community to study such non-closure causes in more detail.

4.6. Machine Learning ET Retrievals

There is no doubt that the T_s is the most influential indicator of the land-atmosphere evapotranspiration. However, in the last years there have been a noticeable increase of published papers [229,230,231,232,233,234,235,236,237] that fused remote sensing techniques with a high variety of machine learning (ML) methods to estimate ET values. Artificial intelligence approaches using the ML techniques like: neural network (NN), support vector machine (SVM) or random forest (RF) have tried to consider the non-linearity between ET and the variables for meteorological and land surface conditions. However, these methods can be time-consuming since in order to use properly a machine learning algorithm, it is necessary to train it with a different dataset, and once optimized, the ML algorithm is checked against a test dataset, different from the training and validation datasets. A new review study about this topic [238] is recommended to the reader for more detail.

4.7. Merging ET Models

In the last five years the number of studies related with the ET estimations with multi-model techniques have raised. The reason of such increase is because hybrid ET models may address landscape heterogeneity occurring in most ecosystems. These areas are partially covered with soil, trees, and grasses patches, so taking advantage of better performance of single source ET models in homogeneous areas and Two-source ET models in heterogeneous ones is of key importance for best ET retrievals. So, a new growing trend to work with three-source SEB models, to account for complexity of various contributions of different vegetation elements in the scene has been stated [239], results showed that three-source SEB models can be of great utility for irrigation management of patchy surface as vineyard. Some of these hybrid ET models are also based on machine learning techniques such SVM [240], also with the Bayesian model averaging (BMA) method, which is based on the basis of the posterior probability distribution of multiple candidate models [101,241]. BMA can not only predict more accurate values, but do so with smaller biases and variances, therefore this method usually outperforms individual candidate models. Finally, a recent study [119] converted the single-source SEBS into a dual-source model (called TS-SEBS) to partition ET into E and T for detailed studies of flux exchanges between soils and canopies.

5. Conclusions

Monitoring the evapotranspiration of a crop or vegetated surface, through remote sensing techniques, has been one of the goals of disciplines such as agriculture, hydrology or climatology and meteorology, for a variety of purposes such as optimizing the irrigation scheduling and the agronomic management, or for the prediction of natural hazards such as floods and droughts. Direct estimation is simply not possible, and ET is generally obtained as a residue from the energy balance equation, once the terms R_n, G and H have been estimated, assuming that other effects, such as advection, are negligible, which is often not the case. However, the original concept of using an aerodynamic temperature as a reference was substituted by the information provided by the surface radiometric temperature, more easily derived from satellite data. Although results in the literature show that the aerodynamic surface temperature can be parameterized locally for different vegetation surfaces, surface homogeneity, and environmental and climatological conditions, with relative success, the lack of a robust and comprehensive model for the application of this approach keeps most efforts focused on the use of satellite thermal images. This review article brings together and updates the main ET models, based on LST data obtained from satellite. These models fall into 3 categories: (i) one-source models, (ii) triangular/trapezoidal LST-VI relationship models, and (iii) two-source models. A comprehensive selection of validation studies for each particular model is included.

After updating the state of the art of the subject under study, no firm conclusions can be stated regarding the better performance of one model versus others. One can only extract partial findings such as the better suitability of the one-source approaches for full cover crops, whereas the two-source models can better represent ET in woody or row crops. In general terms, the validation studies of all the ET models shown in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6 and Table 7 establish the estimation error of the daily ET ranging 1–5 mm/d, with a mean error of around 1 mm/d regardless of the sensor, the spatial resolution or the type of surface under study. However, it is noticed that for ET models such as S-SEBI, SEBS, LST-VI Triangle/Trapezoidal or Two-Source, the error is smaller (around 1 mm/d) at high spatial resolutions (60–100 m), and degrades up to 2 mm/d for kilometer resolutions. As for the one-source ET SEBAL and METRIC models, validation studies report the smallest errors, with values ranging from 0.4–2.4 mm/d, although most of these studies focus and are based on high spatial resolution sensors. One of the main findings in this review is that new formulations, approaches or model updates do not succeed in reducing the uncertainty in ET modeling over the years. For these reasons, decisions on the models to be implemented to derive operational ET products based on coming missions such as ESA′s LSTM or NASA′s SBG, should rather stand on reliability and robustness under a variety of environmental and surface conditions. Joint efforts initiatives such as the recent OpenET should prevail in this search. Some important challenges that the field of ET monitoring from satellite data have been identified, with special attention on the use of machine learning methodologies to improve ET estimates.