geeSSEBI: Evaluating Actual Evapotranspiration Estimated with a Google Earth Engine Implementation of S-SEBI

Abstract

Quantifying and mapping evapotranspiration (ET) from land surfaces is crucial in the context of climate change. For decades, remote sensing data have been utilized to estimate ET, leading to the development of numerous algorithms. However, their application is still non-trivial, mainly due to practical constraints. This paper introduces geeSSEBI, a Google Earth Engine implementation of the S-SEBI (Simplified Surface Energy Balance Index) model, for deriving ET from Landsat data and ERA5-land radiation. The source code and a graphical user interface implemented as a Google Earth Engine application are provided. The model ran on 871 images, and the estimates were evaluated against multiyear data of four eddy covariance towers belonging to the ICOS network, representative of both forests and agricultural landscapes. The model showed an RMSE of approximately 1 mm/day, and a significant correlation with the observed values, at all the sites. A procedure to upscale the data to monthly is proposed and tested as well, and its accuracy evaluated. Overall, the model showed acceptable accuracy, while performing better on forest ecosystems than on agricultural ones, especially at daily and monthly timescales. This implementation is particularly valuable for mapping evapotranspiration in data-scarce environments by utilizing Landsat archives and ERA5-land radiation estimates.

1. Introduction

Evapotranspiration (hereafter ET) is a key component of the terrestrial water cycle [1] and serves as a link to the surface energy balance. Evapotranspiration is the sum of the water that evaporates from soil surfaces and transpiration from vegetation tissues. Several techniques have been developed to directly measure the fluxes from land surfaces: eddy covariance [2,3], Bowen ratio measurement [4,5], and weighting lysimeters [6,7]. All of them provide estimates over restricted areas and are costly; due to these constraints, their application is still restricted. For several decades, remote sensing-based methods have been under continuous development and refinement to obtain spatial estimates of ET [8]. Currently, several medium- and coarse-resolution ET products are available, and they have been extensively evaluated [9,10,11,12]. Although these products are invaluable for assessing evapotranspiration at regional or global scales, their coarse spatial resolution significantly limits their applicability in field-scale investigations [13,14]. This is especially true in landscapes with heterogeneous land cover, where the phenomena under investigation may occur at smaller spatial scales [15]. This limitation can be addressed by using higher-resolution imagery, such as that produced by the Landsat, Ecostress missions, or the ASTER sensor, or even by performing airborne remote sensing campaigns over the area of interest [13]. In these use cases, the application of specific models or algorithms enables ET estimation [16]. The Surface Energy Balance (SEB) models partition the incoming radiation into different energy fluxes by using remotely sensed land surface temperature data, along with additional information on the land surface properties and meteorological conditions. Among this class of models, there is the TSEB [17] (Two Source Surface Energy Balance), S-SEBI [18] (Simplified Surface Energy Balance Index), SEBS (Surface Energy Balance System), SEBAL [19] (Surface Energy Balance Algorithm for Land), METRIC (Mapping Evapotranspiration with Internalized Calibration) [20], SSEBop (Simplified Surface Energy Balance operational) [21], ALEXI (Atmosphere Land Exchange Inverse model) [22], and DisALEXI (Disaggregated Atmosphere Land Exchange Inverse model) [23]. The basic assumption of the surface energy balance approach is represented the equation:

$$Rn=G_0+H+LE$$

(1)

where Rn is the net radiation at the earth surface, G₀ is the heat flux to the soil, H is sensible heat, and LE is latent heat of evaporation of water. The equation refers to the instantaneous energy balance at the earth surface. The energy balance can be solved using a two-source approach (considering separately the contribution of vegetation and soil) such as the one of TSEB, or a one source approach, without separating the two components.

TSEB partitions the radiation budget between the canopy and the soil, and for the canopy, assumes the Priestley–Taylor model of evapotranspiration. Then, through an iterative process, it forces the components of the energy budget to sum 0 and uses as constraint to determine the actual evaporation and transpiration rate, the temperatures of soil and vegetation, obtained by disaggregation from the land surface temperature. It considers separately the resistances of soil and of vegetation [17].

Among the “one source” models, SEBS, SEBAL, and METRIC work by first determining the sensible heat flux, each following a specific approach, and then use it for estimating the latent heat flux, and finally evapotranspiration. Specifically, SEBS, after the G₀ computation, assumes for each pixel two opposite boundary conditions: maximum LE, corresponding to the maximum possible evapotranspiration rate, that is potential evapotranspiration, and LE equal to 0. Actual H is computed based on the Monin–Obukov similarity theory, and based on that, the evaporative fraction is computed. On the other hand, SEBAL and METRIC rely on Equation (2) for estimating the sensible heat flux:

$$H=\rho c_p{\displaystyle \frac{dT}{r_a}}$$

(2)

where $\rho$ is air density, c_p is the specific heat of air, dT is the temperature difference between the surface and air, and r_a is aerodynamic resistance. dT is then determined based on surface temperature, calibrating the relationship between surface temperature and dT on extreme boundary conditions.

S-SEBI and SSEBop use the land surface temperature to determine the evaporative fraction (ratio of LE to the sum of H plus LE in the case of S-SEBI, and ratio of ET to ET_reference, respectively). S-SEBI assumes two opposite possible conditions: radiation-limited evapotranspiration, where there is no moisture limitation and all radiation is allocated to ET, associated with the coldest pixels, and a dry limit, where no evapotranspiration occurs, associated with the hottest pixels. The position of each pixel between these two extremes on the albedo LST scatterplot is used for deriving its evaporative fraction.

While a more complex model with a better description of physical processes should lead to a higher accuracy, in practice this is not always true. The main reason is often the insufficient accuracy of the data regarding atmospheric conditions, which overrides the advantages of a detailed model [24,25]. Thus, in many data-scarce areas, semi-empirical approaches may still constitute the optimal choice [24,26,27]. Among these, in 2000, the S-SEBI model was proposed by Roerink et al. [18]. The main strength of this model is that it requires only radiation data, besides the image acquired in the shortwave and infrared wavelengths, thus making it suitable for areas where there is large uncertainty in the meteorological conditions, as the evaporative fraction is derived from the relationship between land surface temperature (LST) and albedo [28].

Usually, for practical applications, total daily evapotranspiration estimates are more meaningful than instantaneous ones [29]. The instantaneous values retrieved with SEB models can be scaled to daily by employing an appropriate strategy [30]. Specifically, it has been demonstrated that the evaporative fraction (the ratio between latent heat of evaporation, and the sum of latent and sensible heat) varies very little during a single day [31], and thus it can be coupled with radiation estimates of the whole day to scale up ET to daily values. Another possibility is to use the ratio between actual and potential evapotranspiration [32], or to use the ratio between daily and instantaneous radiation as a scaling factor [33].

The availability of the Google Earth Engine cloud computing infrastructure greatly increased the computational capabilities for the geoscientific community [34]. The reimplementation of some of the surface energy balance models enabled a thorough validation over a large number of sites, and even the creation of ensemble products, averaging the estimates of different models [35,36,37,38]. A large-scale evaluation of models can give us a more general picture of their performance under different conditions, and the causes underlying their inaccuracies, enabling the scientific community to improve and refine them, leading finally to a deeper understanding of the processes [39]. In recent years, the METRIC, SEBAL, SEBALI, SEBU, and SSEBop models were implemented in Google Earth Engine [40,41], which enhanced their usability among the scientific community, and led to the creation of evapotranspiration products, like the openET dataset for CONUS (the Conterminous United States) [35]. However, for other regions of Earth, a similar, high-resolution product based on LANDSAT thermal imagery does not exist. These elements led to the necessity to (i) develop a scalable and widely available infrastructure for the S-SEBI model as well, and (ii) to test its accuracy at European eddy covariance sites, on multiannual timeseries. This will deepen our understanding of the S-SEBI performance under different hydrological and meteorological conditions, and show us whether the model, applied to LANDSAT imagery and ERA5-land radiation data, can reliably reproduce evapotranspiration. This is a crucial step towards the creation of operational evapotranspiration products.

The aim of this paper is to introduce geeSSEBI, a Google Earth Engine code, with a Google Earth Engine application as graphical user interface, for increasing its usability, and integrating the long-term Landsat imagery archives with the ERA5-land reanalysis data. This study, besides introducing an accessible implementation of the S-SEBI model, constitutes an unprecedented effort of validation and evaluation of the model accuracy, as it is the first to implement and validate S-SEBI using multiannual datasets from European ICOS sites. The results of the application of the model to 871 images, acquired between 2005 and 2023 at four ICOS sites, are presented, using as source of radiation data the ERA5-land reanalysis product. The evaluation of the model at many times, representing different meteorological and hydrological conditions, is crucial for understanding in depth its points of strength and potential shortcomings, thus helping model improvement [42]. Integration of the available big datasets is a key issue in remote sensing, hydrology, and related disciplines, and it is a key step towards advancing our knowledge and understanding of the biogeochemical cycles and ecosystem processes, especially in data-scarce areas, where ground data are not available.

2. Materials and Methods

2.1. Model Structure

The S-SEBI model was first proposed by Roerink et al. in 2000 [18]. In the implementation presented hereafter, the incoming shortwave and longwave radiation are retrieved from ERA5-land as hourly data, and instantaneous values are computed by dividing them by 3600. Then, they are used for computing the net radiation, following Equation (3).

$$R_n=\left(1-\alpha \right)R_{sw\downarrow }+\epsilon R_{lw\downarrow }-\epsilon \sigma T_s^4$$

(3)

where α is the surface albedo, ε is the emissivity, R_n is net radiation (as in Equation (1)), and T_s is surface temperature. Albedo and emissivity, in the case of Landsat images, can be estimated from surface reflectance in the different spectral bands and fraction of vegetation cover (f_c), respectively (Equations (4)–(6); for Landsat 8 and 9 [43], the albedo equation for Landsat 5 is reported in Appendix A Equation (A1). The fraction of vegetation cover can be derived from the NDVI values (Equation (7) for Landsat 8 and 9, with the Landsat 5 equation in Appendix A Equation (A2) [43,44]. The coefficient values for albedo computation used in this implementation are derived from [45].

$$\alpha =0.13{\rho }_{B1}+0.115{\rho }_{B2}+0.143{\rho }_{B3}+0.18{\rho }_{B4}+0.281{\rho }_{B5}$$

(4)

$$\epsilon =0.971\left(1-f_c\right)+0.982f_c$$

(5)

$$f_c=\left\{\begin{array}{c}\begin{array}{cc}NDVI<0.2, & 0\end{array} \\ \begin{array}{cc}0.2<NDVI<0.8, & {\left({\displaystyle \frac{NDVI-0.2}{0.6}}\right)}^2\end{array}\\ \begin{array}{cc}NDVI>0.8,& 1\end{array}\end{array}\right.$$

(6)

$$NDVI={\displaystyle \frac{{\rho }_{B5}-{\rho }_{B4}}{{\rho }_{B5}+{\rho }_{B4}}}$$

(7)

Soil flux (G) is estimated as a fraction of net radiation, dependent on the fraction of vegetation cover in the pixel (Equations (8) and (9)). A value of 0.05 is assumed for the soil flux coefficient (G₀) under full vegetation cover, while a 0.315 ratio is assumed for bare soil surfaces [46].

$$G=G_0{R}_n$$

(8)

$$G_0=0.05f_c+0.315\left(1-f_c\right)$$

(9)

The evaporative fraction is computed from the surface temperature of each pixel, and the potential cold and hot temperatures, given its albedo. For estimating maximum and minimum potential temperatures for each albedo class, temperature values are aggregated for albedo intervals of width 0.01. For each of these albedo classes, maximum and minimum temperature are obtained. Then, linear regression between LST_max (corresponding to the temperature of the dry limit T_dry) and albedo is used for estimating the dry limit, and linear regression between LST_min (corresponding to the temperature of the wet limit T_wet) and albedo for the wet limit (yellow and red regression lines in Figure 1, respectively). Albedo is then used for obtaining the actual dry and wet temperature limit for each pixel of the area. Finally, evaporative fraction is computed following Equation (10).

$${EF}_i={\displaystyle \frac{T_{dry}-T_i}{T_{dry}-T_{wet}}}$$

(10)

Negative EF values are set to 0, while EF values exceeding 1 are set to 1.

Latent heat flux is computed from the net radiation by using this evaporative fraction value, and the remaining part of energy is allocated to the sensible heat flux.

$$LE={EF(R}_n-G)$$

(11)

$$H=(1-EF)(R_n-G)$$

(12)

Instantaneous evapotranspiration is computed by multiplying the latent heat flux by the latent heat of vaporization of water (2.46 × 10⁶ J/Kg), thus yielding ET in mm (1 Kg/m² = 1 mm, assuming water density to be 1 Kg/dm³).

Figure 1. Example of the S-SEBI diagram, showing the relationship between LST and albedo, and the wet and dry limits, employed in the evaporative fraction computation, represented by the red and yellow regression lines, respectively. Blue dots represent several sample pixels of the area of interest.

Figure 1. Example of the S-SEBI diagram, showing the relationship between LST and albedo, and the wet and dry limits, employed in the evaporative fraction computation, represented by the red and yellow regression lines, respectively. Blue dots represent several sample pixels of the area of interest.

2.2. Hourly, Daily, and Monthly Scaling

Daily ET is scaled assuming a constant evaporative fraction during daytime [31]. A conversion factor C_di is computed (Equation (13)) [47], which is the ratio between the instantaneous and daily radiation downwelling. Instantaneous ET is converted to daily by multiplying it by this conversion factor (Equation (14)). The net radiation of the 24 h is obtained from the ERA5-land daily aggregates and is computed as the sum of the net radiation in both the shortwave and the longwave.

$$C_{di}={\displaystyle \frac{R_{\downarrow day}}{R_{\downarrow inst}}}$$

(13)

$${ET}_{day}=C_{di}\ast {ET}_{inst}$$

(14)

For scaling the daily ET to monthly, the Allen’s theoretical framework is followed [20], adapting it to the context of the S-SEBI model: the ratio of daily ET to net radiation, for each Landsat overpass date is computed, and the timeseries is linearly interpolated in time. Then, each value of this timeseries is multiplied by the net radiation of the corresponding day. Finally, the monthly aggregated evapotranspiration is the sum of the daily values of each month. The steps of the procedure are graphically represented in Figure 2.

2.3. Code Structure and Availability

A JavaScript code implementing the aforementioned model, written for the Google Earth Engine interface, is made available, both as a public library and as a Google Earth Engine application. The workflow of the code is summarized in Figure 3, and the functions available in the library are enumerated, along with their parameters, in Appendix B. The code is divided into two main functions: preprocessInputs and runSSEBI. The preprocessInputs functions exists in three variants, specific for each satellite platform (Table S1), each of them accepted as parameters: startDate, endDate, and geometry. It gets the first Landsat image (Collection 2 Level 2) acquired during the time interval specified, that intersects the geometry identifying the study area. The image is masked, removing all of the pixels marked as cloudy according to the quality assessment band included in the product, and all of the adjacent pixels. Pixels with land surface temperature values lower than 0 °C or higher than 70 °C are masked out as well. Albedo, NDVI, fraction cover, emissivity, and NDWI (Normalized Difference Water Index [48], Equation (15)) are computed following Equations (3)–(6) and (14). All of the pixels with NDWI higher than 0.0 are masked out, considered as water bodies [49]. Instantaneous and daily radiation values for the area of interest, for the hour corresponding to the acquisition time of the Landsat image, are retrieved from the ERA5-land reanalysis product; net radiation and the scaling coefficient C_di are then computed, following the equations described in the Section 2.1.

$$NDWI={\displaystyle \frac{{\rho }_{B3}-{\rho }_{B6}}{{\rho }_{B3}+{\rho }_{B6}}}$$

(15)

Then, in the runSSEBI function, pixels are grouped to albedo classes spanning 0.01. For each albedo class, the maximum and minimum surface temperature value is retrieved from the image, and then, linear regression is performed for determining the regression lines representing the moisture-limited and radiation-limited conditions (Figure 1). Then, the regression line equations are used for determining T_dry and T_wet for each pixel, given its albedo, for computing the evaporative fraction according to Equation (10). The soil flux coefficient is computed according to Equations (8) and (9), and the soil flux is determined from net radiation. Then, the evaporative fraction is used for computing LE and H. Finally, the instantaneous LE value is converted to daily ET as described in Section 2.2.

The code can be used through the Google Earth Engine application, available at https://ee-jpk1996jpk-2024.projects.earthengine.app/view/geessebigui. The source code is available at the GitHub repository https://github.com/jpkabala96/geeSSEBI (accessed on 24 October 2024). For use in Google Earth Engine, the code is available as a public library as well: “users/jpk1996jpk/geeSSEBI”. It can be imported in a Google Earth engine script by executing the following line of code:

var geeSSEBI = require(’users/jpk1996jpk/geeSSEBI:geeSSEBI_library’);

which enables the users to use the functions later in their script.

A script, showing some examples of the application of the code, is stored in the library and in the GitHub repository (file named: “examples.js”), along with the source code.

2.4. Evaluation Metrics

The estimates of the model are evaluated at each site at instantaneous and daily timescales. For the Bosco Fontana and San Rossore2 sites, the ICOS Level 2 data product is used. Daily and hourly values corresponding to the satellite overpass are extracted. For computing ET, the corrected LE value is multiplied by the latent heat of vaporization of water (2.46 MJ/Kg). For each site, the following model performance metrics were computed: the root mean square error (RMSE), Pearson’s correlation coefficient (R), coefficient of determination (R²), KGE (Kling–Gupta Efficiency), NSE (Nash–Sutcliffe Efficiency), and percentage bias using the hydroGOF R package [50]; the equations used for computing them are shown in Appendix C. ET values were extracted from the evapotranspiration maps produced by the Google Earth Engine code by extracting the nearest neighbour value to the coordinates of the eddy covariance towers. Additionally, the accuracy of the estimates produced with geeSSEBI were compared against the accuracies of the 8-day global MODIS evapotranspiration product. To do this, the MODIS ET values were extracted from the product for the same time period (indicated in

Table 1. Overview of the ICOS sites where model evaluation was performed: site name, latitude (°), longitude (°), elevation (m a.s.l.), ecosystem, mean annual temperature (MAT; °C), mean annual precipitation (MAP; mm), date of start and end of the studied period, and number of images selected for running the model.

Site	Lat	Lon	Elev	Ecosystem	MAT	MAP	Date Start	Date End	Number of Images
San Rossore2	43.73°	10.29°	4 m	Pine forest	15.3 °C	950 mm	1 January 2013	31 December 2023	296
Bosco Fontana	45.19°	10.74°	37 m	Mixed broadleaved forest	14.5 °C	697 mm	1 January 2016	31 December 2023	161
Lison	45.74°	12.75°	1 m	Vineyard	13.0 °C	1100 mm	1 January 2016	31 December 2023	223
Lamasquere	43.49°	1.24°	121 m	Crop (maize and winter wheat)	13.4 °C	677 mm	1 January 2005	31 December 2023	175

) by the nearest neighbour algorithm. geeSSEBI daily estimates of ET and eddy covariance daily ET values were aggregated for the same 8-day time intervals as the MODIS data. Then, the accuracy of both geeSSEBI and MODIS were compared.

Finally, the correlation between the error and the number of images per month, actual ET, and radiation was evaluated with Kendall’s correlation test (significance level alpha = 0.05). This additional result is reported in Supplementary Table S2.

2.5. Validation Sites and Data

The validation of the model outputs was performed on four different study sites, representing different ecosystems (Figure 4). Two are forested areas: San Rossore2 ICOS station [51,52], where measurements are performed in a stone pine forest; and Bosco Fontana ICOS site [53], where a mixed broadleaved forest is monitored. The Lison ICOS site [54] is representative of a vineyard, while the Lamasquere ICOS study site [55] is representative of croplands. For each of the study sites, the Level 2 data products were downloaded.

Instantaneous latent heat flux estimates by S-SEBI were extracted from the rasters produced by using the nearest neighbour algorithm. Then, they were compared with the corresponding half-hourly LE values provided in the product, corrected for the surface energy balance closure. Daily evapotranspiration estimates by S-SEBI were compared with daily average LE, corrected for the surface energy balance closure, which was converted from radiation units to millimeters of water per day, following Equation (15).

$$ET={\displaystyle \frac{LE\ast 3600\ast 24}{2.46\ast {10}^6}}$$

(16)

where LE is the average latent heat of evapotranspiration in W/m², 3600 is the number of seconds in an hour, 24 is the number of hours in a day, and 2.46 × 10⁶ is the latent heat of vaporization of water, in J/kg.

The model was fitted to all of the Landsat images available for each study site, selecting the ones with less than 20% of cloud cover in the area of interest. The satellite images were cloud masked, and an additional 100 m buffer was applied to the cloud mask in order to ensure the exclusion of pixels possibly affected by clouds. Moreover, to exclude water bodies from the computation, NDWI was computed, and pixels having a value greater than 0.0 were masked out. On the other hand, from the ICOS Level 2 databases, rows with invalid LE_CORR values were removed. The IDs of the Landsat images used in the comparison are reported in Supplementary Table S1. Additionally, validation on the non-overpass dates, which were used for producing the monthly timescale data, was conducted, and the same accuracy metrics as for the overpass dates were computed.

3. Results

3.1. Instantaneous Timescale

S-SEBI estimates of latent heat flux (LE) exhibited a significant positive correlation with eddy covariance measurements across all evaluated sites (Figure 5A–D). In terms of accuracy, San Rossore2 exhibited an RMSE of 137 W/m², Bosco Fontana had 130 W/m², Lison recorded 124 W/m², and Lamasquere showed 92 W/m².

3.2. Daily Timescale

At all of the sites, the model showed an accuracy around 1 mm/day in terms of RMSE, similar with what was reported by previous studies, and a good correlation between modelled and observed values (Figure 6A–D). At the daily scale, the sites representative of forest ecosystem showed lower RMSE than the ones in land used for agriculture, mainly due to the higher bias for these two sites (

Table 2. Accuracy of the S-SEBI estimates of daily evapotranspiration for the study sites.

Site	RMSE	R	R2	KGE	NSE	% Bias	N° Images
San Rossore2	1.05	0.73	0.54	0.70	0.37	−7.2	296
Bosco Fontana	1.03	0.81	0.66	0.80	0.61	−4.5	161
Lison	1.37	0.86	0.74	0.58	0.53	−32.7	223
Lamasquere	1.16	0.80	0.64	0.69	0.47	−24.5	175

), while showing higher R and R squared.

For non-overpass dates, the relationship between modelled and measured ET was very similar to that for overpass dates (Figure 7), and also, the accuracies were very similar (

Table 3. Accuracy metrics for the non-overpass dates of the daily timeseries of evapotranspiration reconstructed with the S-SEBI model.

Accuracy Metric	Bosco Fontana	San Rossore2	Lison	Lamasquere
RMSE	0.79	1.03	1.32	1.3
% bias	3.9	10.6	−43.5	−45.2
NSE	0.75	0.43	0.49	0.2
r	0.89	0.8	0.88	0.74
R2	0.79	0.64	0.78	0.55
KGE	0.85	0.68	0.43	0.44

). Specifically, for the two forest sites, the slope of the regression line between modelled and observed values was very close to one, while for the two agricultural sites, it was lower, suggesting a negative bias, as confirmed by the accuracy metrics in Table 3. Moreover, at all of the sites, the S-SEBI modelled values showed a positive NSE, which was higher at the forest sites Bosco Fontana and San Rossore2, and lower at the agricultural ones. RMSE, r, R squared, and KGE were comparable to those computed for the overpass dates, with slight differences.

3.3. Monthly Scaling

The monthly scaling for the San Rossore2 site resulted in a 24.5 mm/month RMSE, with 0.92 correlation (p < 0.05, Figure 8A). In the Bosco Fontana, the RMSE was 16.9 mm/month, with 0.95 correlation (p < 0.05, Figure 8B). At the Lison site, RMSE was 35.5 mm/month, and correlation was 0.94 (p < 0.05). S-SEBI underestimated the eddy covariance evapotranspiration, similar to what occurred at the daily timescale. At Lamasquere, RMSE was 37 mm/month, with 0.82 correlation (p < 0.05). For both of the sites representative of agriculture, the negative bias found at the daily scale was still present at the monthly scale, while the high correlation values confirm that the model properly represents the ET trend (Figure 8C,D).

3.4. Landsat 5 Evaluation

For the Lamasquere study site, where the eddy covariance timeseries started in 2005, the analysis of the model performance with Landsat 5 images was possible. The equations used for determining albedo and NDVI from Landsat 5 imagery [56,57] are shown in Appendix A. A total of 58 Landsat 5 images, with more than 80% of cloud-free area in the area selected for the analysis, were selected, spanning a wide range of dates between 2005 and 2011. The model showed a 90 W/m² RMSE for the instantaneous LE estimate, and a 1.29 mm/day RMSE for the daily scaled evapotranspiration. In both cases, the estimates produced by the model had high correlation with the measured values (0.82 for instantaneous and 0.78 for daily, Figure 9A,B). Due to the low amount of good-quality images available for model application on the Landsat 5 dataset (average 7 images/year, less than 1/month), the monthly scaling was not evaluated on this dataset.

3.5. Evaluation of S-SEBI Estimates Against MODIS Evapotranspiration

The comparison of the accuracy of S-SEBI and the MOD16A2 8-day global evapotranspiration product showed that in general, the two have a comparable accuracy (Figure 10). Both showed the worst performance at the Lison site, where they accurately represented the trend (good correlation and coefficient of determination against the eddy covariance data), but strongly underestimated the absolute ET values (strong negative % bias, and so high RMSE and poor KGE and NSE compared to the other sites). At the two forest sites (Bosco Fontana and San Rossore2), the two datasets performed similarly, with S-SEBI outperforming MODIS in some cases (e.g., KGE for the San Rossore2 site). At the Lamasquere site, MOD16A2 outperformed S-SEBI, not showing the strong negative bias, and thus reproducing better the measured ET timeseries (better KGE and NSE).

4. Discussion

This paper outlines a methodology for estimating actual evapotranspiration by applying the S-SEBI model to high-resolution remotely sensed data and utilizing ERA5-land reanalysis as the source for radiation data. Simple models like S-SEBI, which require minimal meteorological input, can be particularly valuable for estimating evapotranspiration in regions with limited data availability. To facilitate computation and enable large-scale application of the model, a Google Earth Engine script was developed.

The model’s accuracy was extensively evaluated at four ICOS sites with long-term eddy covariance timeseries of evapotranspiration. This represents an unprecedented effort to assess the model’s performance, which in previous studies was often limited to short timeframes [58,59]. Evaluating the model under a wide range of conditions is essential to ensure its reliability and to account for potential errors and limitations in its estimates. The model’s accuracy at a daily timescale was approximately 1 mm/day RMSE across all sites, aligning with previous research findings on this topic [58,59]. Notably, the accuracy was slightly higher for forest sites compared to agricultural ones, potentially due to the greater homogeneity of forest canopies relative to croplands [60]. The model effectively captured trends in evapotranspiration over the years. These results align with findings from other applications of the S-SEBI model, further confirming its ability to reproduce seasonal trends in evapotranspiration over time.

The higher error and negative bias in the estimates over agricultural landscapes found in this study might be due to the inability of the model to represent phenomena such as advection, which can lead ET to exceed incoming radiation, in irrigated or wet ecosystems under dry or semiarid climate [61]. In fact, the high correlation coefficient and the significant negative bias indicate that while the trend (in terms of relative variation over time) was successfully detected by the model, the absolute ET values were underestimated. The issue might be mitigated by the application of different models such as METRIC, which employ reference evapotranspiration computed from in situ meteorological data for its calibration [62]. The model absolute error is positively correlated with incoming radiation and the actual ET value (Supplementary Table S2), while the percentage error is negatively correlated with the number of images, radiation, and actual ET. This suggests that the availability of cloudless images plays an important role in determining the quality of the output.

A monthly scaling procedure was tested to upscale the non-continuous daily estimates of ET. The results indicate that using ERA5-land data and the ET-to-net radiation ratio enables successful upscaling of S-SEBI estimates to monthly values. The availability of user-friendly implementations of surface energy balance models can enhance our understanding of evapotranspiration in data-scarce regions. This information could be valuable in the future for analyzing the impacts of forest or crop management practices [63], quantifying the effects of events such as wildfires on the water cycle [63,64,65], characterizing irrigation requirements [66,67], assessing the water budget of specific areas [68], or examining the effects of drought on vegetation [69]. The tests on Landsat 5 imagery showed that even the archives of data collected by old sensors might be used for evapotranspiration estimation, with an accuracy comparable to the one that can be obtained with Landsat 8 data, when using the S-SEBI model. This might be very useful for analysing the effects of past events on evapotranspiration and for performing retrospective analyses [70].

The model showed, at the evaluation sites, a comparable accuracy to that of the MOD16A2 evapotranspiration product. Both reproduced the trends in the ET timeseries well; however, in some specific sites, they showed a significant bias in the absolute value of ET, likely due to local-scale phenomena that could not be accounted for.

Recently, Senay et al. [38] applied SSEBop in the CONUS area, reporting a variability in accuracy over different ecosystems: specifically, SSEBop performed best on croplands, while it showed a significant positive bias in forest sites. The overall accuracy of the model, according to this study, was 30 mm/month, which is similar to the results of this study, while geeSSEBI performs better on forest sites than on agricultural ecosystems. Future comparative experiments running different models in the same areas would be advisable, in order to make a wise and thorough comparison of the available tools. geeSEBAL, when applied at different study sites in Brasil, coupled to ERA5-land data as meteorological inputs, showed an accuracy around 0.7 mm/day [37]. Finally, eeMETRIC and SEBALIgee were applied on other areas that do not overlap with the areas examined in the present study, and on very different land cover types, reasons for which we do not report the results of those studies [71,72].

Moreover, a massive application of the models and the analysis of their errors will improve them in order to get more accurate estimates (e.g., [73,74]). Finally, easily usable implementations of the models might allow their intercomparison, in order to better understand their advantages and disadvantages under specific conditions [60].

5. Conclusions

The model’s accuracy was thoroughly evaluated at four ICOS sites with long-term eddy covariance time series of evapotranspiration. This effort represents an unprecedented attempt to assess the model’s performance, as previous studies often focused on shorter timeframes [58,59]. Validation against eddy covariance tower data demonstrated that the model’s accuracy is acceptable and comparable to the MOD16A2 product at the evaluated sites.

In recent years, several surface energy balance (SEB) models have been implemented within the modern Google Earth Engine cloud computing platform, significantly increasing their accessibility to a broad user base and enabling large-scale evaluations of model accuracy under diverse conditions. The development of geeSSEBI aligns with this trend, providing the scientific community with a user-friendly tool for applying the S-SEBI model. By leveraging the computational power and data catalog of Google Earth Engine, geeSSEBI significantly reduces the technical barriers to applying the S-SEBI model, making it accessible to a broader user base, including practitioners and researchers in data-scarce regions.

Additionally, a novel procedure for upscaling daily ET estimates to monthly timescales was evaluated, enabling its use for long-term water resource monitoring and ecosystem studies. The accuracy of the model was rigorously evaluated against multiyear data from eddy covariance towers, providing insights into its performance across different ecosystem types.

Despite the extensive dataset analyzed in this study, it would be valuable to expand the comparison to additional ICOS sites and to evaluate the model against other widely used evapotranspiration models. Such efforts could further enhance the understanding of turbulent fluxes and contribute to advancing the modelling discipline. Finally, the widespread availability of the model code may empower future users to refine and improve the models further.