Towards a global actual evapotranspiration product for the Copernicus Land Monitoring Service

Copernicus Land Monitoring Service (CLMS) produces biogeophysical maps of the global land surface. The CLMS portfolio so far did not include actual evapotranspiration (ETa), despite it being a direct link between the energy, water and carbon cycles and its importance for global food security, efficient water resources management and weather forecasting. However, a global CLMS ETa product was recently developed and entered operational production at the end of 2025. It consists of water flux (ETa, evaporation and transpiration) and energy flux (latent and sensible heat) sub-products, has a spatial resolution of 300 m, dekadal (10 d mean) temporal resolution for water fluxes and daily temporal resolution for instantaneous energy fluxes, is produced in near-real-time, and (like all other CLMS products) is distributed under free and open data policy. It is based mainly on Copernicus input data with primary satellite imagery coming from the observations of OLCI and SLSTR sensors on board of Sentinel-3 satellites. Such product fills a gap in previously existing global and operational ETa products, thus satisfying a wide range of potential users' needs. In this paper, we describe, evaluate and justify the various design choices taken during the development of the ETa product, ranging from cloud masking and gap-filling, through derivation of biophysical traits, radiation components and weather forcings to spatial sharpening of the land surface temperature observations. Those data were then used to drive two ETa models: TSEB-PT and ETLook. A prototype implementation of the ETa processing chain was used to produce ETa data (using both models individually and their ensemble) across a globally representative range of climatic zones and plant functional types, which was validated against measurements from 206 Eddy Covariance flux tower sites. The ensemble resulted in the all-site root mean squared error (RMSE) of 0.87 mm d⁻¹ (relative RMSE of 47 %, site range: 0.31–1.96 mm d⁻¹), bias of 0.01 mm d⁻¹ (relative bias of 0.6 %, site range: −1.39–1.50 mm d⁻¹) and correlation coefficient of 0.81 (site mean of 0.77), which compares well with WaPOR global ETa dataset and therefore supports the operational production of global, near-real time, Copernicus-based 300 m ETa maps. Finally, we propose a number of potential evolutions of the CLMS ETa product, including enlarging and enhancing the model ensemble, producing a reanalysis version of the product with more robust meteorological forcing and improved gap-filling, and focusing on aspects of modelling related to influence of surface roughness and soil moisture on energy and water fluxes.

1 Introduction

Copernicus (https://www.copernicus.eu, last access: 1 September 2025) is the Earth observation (EO) program of the European Union that provides free and open access to data acquired by the Sentinel satellites (and other contributing missions) and to higher level and non-space products through the Copernicus services. One of the main objectives of the program is to bring operational mindset to EO by: (1) providing data in near-real-time (NRT) with regular updates and derived using robust and validated methods; (2) guaranteeing long-term future continuity of this data and predictable planning of future evolutions. The Copernicus Land Monitoring Service (CLMS – https://land.copernicus.eu, last access: 1 September 2025) produces a series of qualified global biogeophysical products on the status and evolution of the land surface. The products are used to monitor vegetation, water cycle, energy budget and terrestrial cryosphere. Production and delivery are carried out in a timely manner and are complemented by the constitution of long-term time series.

Actual evapotranspiration (ETa) can be defined as the phase change from liquid to gaseous water that is transferred from the land surface towards the atmosphere. It can be partitioned between evaporation and transpiration through plants, with the former component implicitly including also evaporation from intercepted water in the canopy. It has historically not been included in the CLMS portfolio. This changed at the end of 2025 when the CLMS near-real-time (NRT) ETa product entered operational production (https://land.copernicus.eu/en/products/evapotranspiration, last access: 28 January 2026). This product will be of importance for many applications and research questions. The ETa is one of the Essential Climate Variables (ECV) as defined by the Global Climate Observing System (Bojinski et al., 2014). It binds different processes occurring at the Earth's surface (related to energy, water and carbon cycles) and evolves coherently with other biophysical variables, including those which are part of the CLMS portfolio (e.g. leaf area index, gross primary productivity, soil moisture, surface temperature). Since actual evapotranspiration is a direct proxy of plant water use it can be utilized for consistent agricultural water use monitoring across natural and political boundaries, and is therefore essential for Sustainable Development Goal (SDG) reporting, e.g. of SDG indicators 6.4.1 and 6.4.2 (O'Connor et al., 2020). ETa can also be useful for forest fire risk/spread forecasting (Vidal et al., 1994), drought monitoring (Anderson et al., 2011), hydrological modelling (Zhang et al., 2020; Larsen et al., 2016), irrigation accounting (Zhang and Long, 2021) and yield modelling (Jurečka et al., 2021; Gómez-Candón et al., 2021). In addition, spatially distributed fields of ET at an adequate spatial resolution can help to improve weather forecasting (Boone et al., 2025), in particular in irrigated fields in semi-arid climates where the additional water supply in the soil affects the surrounding microclimate and thus the atmospheric boundary layer conditions (Udina et al., 2024; Lunel et al., 2024).

One of the initial primary users of the CLMS ETa product is the Food and Agriculture Organisation (FAO) of the United Nations. FAO has been producing and disseminating ETa datasets, including a global product, through its “WAter Productivity through Open-access of Remotely sensed derived data” (WaPOR – https://www.fao.org/in-action/remote-sensing-for-water-productivity/en, last access: 30 July 2025) project. The WaPOR data has been used in multiple applications such as improving water use productivity or monitoring agricultural practices and water consumption (Chukalla et al., 2022; Hajirad et al., 2023; Seijger et al., 2023). Since there is no guarantee of long-term continuity of the WaPOR project, FAO expressed strong interest to the European Commission for the introduction of a global ETa dataset in the CLMS portfolio (https://www.fao.org/in-action/remote-sensing-for-water-productivity/news-and-events/news/news-detail/copernicus-global-land-service-launches-new-global-evapotranspiration-product-based-on-the-wapor-etlook-model/en, last access: 1 May 2026).

Table 1 Specifications (requirements) of the Copernicus Land Monitoring Service actual evapotranspiration product based on call for tenders JRC/2023/OP/0273 (https://ec.europa.eu/info/funding-tenders/opportunities/portal/screen/opportunities/tender-details/13795, last access: 11 August 2026).

Developing a global, operational product is constrained by certain programmatic, financial and technical limitations. The CLMS products are designed for long-term continuity and are funded by the European Commission. Therefore, issues like compute efficiency, data storage and dissemination costs or long-term availability of sources of core input data need to be considered. Taking those constraints into account, and in order to satisfy a wide range of potential users needs while ensuring consistency with vegetation, primary productivity and other global CLMS products, the European Commission set the requirements for the CLMS ETa product of 300 m spatial resolution, dekadal temporal resolution for water fluxes sub-products and data availability in NRT (within 2 d after the end of each dekad). The product requirements are summarized in Table 1. Those requirements in turn impact the design choices presented in this study for the first version of the CLMS ETa product.

Those requirements also closely match the specifications of the WaPOR global ETa dataset and fill a gap in currently available satellite-based ETa products. Apart from WaPOR (the long-term continuity of which is not guaranteed), global operational ETa datasets produced in NRT with closest matching spatio-temporal specifications are the MODIS and VIIRS ETa products (Román et al., 2024). They have a higher temporal resolution (8 d) but lower spatial resolution (500 m) and use a modeling approach which does not make direct use of land surface temperature (LST) measurements (Mu et al., 2011). Another operational and global product which utilizes MODIS and VIIRS data is produced by United States Geological Survey for their Famine Early Warning Systems Network (FEWS NET) using SSEBop energy balance model (Senay et al., 2020) with dekadal temporal resolution and 1 km spatial resolution. Other ETa datasets have either much lower spatial resolutions (e.g. ETa product based on geostationary observations produced by the EUMETSAT Satellite Application Facility on Land Surface Analysis – Barrios et al., 2024, or the microwave-based GLEAM model – Miralles et al., 2011), do not have global coverage (e.g. OpenET – Melton et al., 2022) or are not produced in NRT (e.g. ETMonitor – Zheng et al., 2022).

Preparatory activities required to develop an operational CLMS ETa product recommended that ET models from two modelling frameworks should be further investigated. The first one is the Sen-ET framework (Guzinski et al., 2020, 2021) developed to model ETa with Copernicus data at various spatial scales and using the Two-Source Energy Balance Priestley-Taylor (TSEB-PT) ET model (Norman et al., 1995; Kustas and Norman, 1999; Anderson et al., 2024). The second is the WaPOR framework developed by FAO through the WaPOR project and using the WaPOR ETLook ETa model (hereafter referred to as ETLook) closely based on Bastiaanssen et al. (2012). Both models, although conceptually different, estimate evaporation and transpiration and use LST as one of core input forcings. This recommendation was mainly based on the availability of mature open-source implementations of the two ET models, on previous studies demonstrating the applicability of both models with Copernicus data sources (i.e. Sentinel-3 imagery and meteorological data from European Center for Medium Range Weather Forecasts) (Guzinski et al., 2021) and on FAO's familiarity with both approaches. This does not imply that those two modelling frameworks clearly outperform all other approaches and indeed it has been shown that the performance of an individual model depends on the landcover and climatic conditions (Reitz et al., 2025). However, due to constraints on developing a publicly funded global and operational dataset (outlined previously) there was a need to limit the design of the first version of the CLMS ETa product to those two frameworks.

This paper aims to describe, evaluate and justify the design choices made for the CLMS ETa product. In Sect. 2, we present the design of the CLMS ETa processing chain, starting with a brief outline of the two ETa models, followed by the description of input data sources and their pre-processing and finishing with the method used for gap-filling of the resulting maps. This is followed by Sect. 3 in which a prototype ETa dataset, produced with both TSEB-PT and ETLook and their ensemble, is validated against measurements from 206 Eddy Covariance (EC) flux tower sites. In Sect. 4, we compare the modelled outputs with each other and with other similar ET datasets, evaluate and justify the design choices described in Sect. 2 and outline suggestions for product improvement. Finally, conclusions are presented in Sect. 5.

2 Data and methods

2.1 ET modelling

2.1.1 TSEB-PT model

The Two-Source Energy Balance (TSEB) modelling scheme was proposed by Norman et al. (1995) and afterwards refined and applied in a multitude of applications and studies (Anderson et al., 2024) including the Sen-ET framework (Guzinski et al., 2020) and implemented as open source in pyTSEB Python package (Nieto et al., 2025a). In this modelling scheme the directional radiometric LST (T_R(θ)) is split into the temperatures of vegetation (T_C) and soil (T_S) based on Leaf Area Index (LAI) and LST observation geometry.

$$\begin{array}{}\text{(1)}& {\displaystyle }{T}_{\mathrm{R}}\left(\mathit{\theta }\right)\approx {\left[f\left(\mathit{\theta }\right)T_{\mathrm{C}}^{\mathrm{4}}+\left[\mathrm{1}\mathrm{-}f\right(\mathit{\theta }\left)\right]T_{\mathrm{S}}^{\mathrm{4}}\right]}^{\mathrm{0.25}}\text{(2)}& {\displaystyle }f\left(\mathit{\theta }\right)=\mathrm{1}\mathrm{-}\exp \left({\displaystyle \frac{\mathrm{-}\mathrm{0.5}\mathrm{\Omega }\left(\mathit{\theta }\right)\mathrm{LAI}}{\cos \mathit{\theta }}}\right)\end{array}$$

where θ is the view zenith angle of the thermal observation and Ω(θ) is the clumping factor of the vegetation at view angle θ (Kustas and Norman, 1999) and has a value of less than 1 for clumped vegetation.

Based on this split, the energy fluxes of vegetation and soil (net radiation – R_n, sensible heat flux – H, latent heat flux – λE, soil heat flux – G) are estimated separately, before being combined to obtain the bulk surface fluxes.

$$\begin{array}{}\text{(3)}& {\displaystyle }{R}_{\mathrm{n},\mathrm{C}}=H_{\mathrm{C}}+\mathit{\lambda }{E}_{\mathrm{C}}\text{(4)}& {\displaystyle }{R}_{\mathrm{n},\mathrm{S}}=H_{\mathrm{S}}+\mathit{\lambda }{E}_{\mathrm{S}}+G\text{(5)}& {\displaystyle }{R}_{\mathrm{n}}=(H_{\mathrm{C}}+H_{\mathrm{S}})+(\mathit{\lambda }{E}_{\mathrm{C}}+\mathit{\lambda }{E}_{\mathrm{S}})+G=H+\mathit{\lambda }E+G\end{array}$$

The soil (canopy) sensible heat flux is computed from the gradient between the soil (canopy) temperature (T_S and T_C respectively) and the air temperature at the sink-source height. This transfer of heat between the two components and the atmosphere is modulated by resistances to heat exchange (r_a,S (Kustas and Norman, 1999; Kondo and Ishida, 1997; Sauer and Norman, 1995) and r_a,C (McNaughton and Van Den Hurk, 1995), respectively) organized in a series resistance network (in analogy to electrical systems) which depend on canopy structure (roughness) and weather (wind and atmospheric stability) conditions.

Since there are multiple solutions to T_C and T_S satisfying Eq. (1), an iterative approach is employed. The initial assumption is that green canopy transpires at potential rate based on Priestley-Taylor formulation (Priestley and Taylor, 1972):

$$\begin{array}{}\text{(6)}& \mathit{\lambda }{E}_{\mathrm{C}}={\mathit{\alpha }}_{\mathrm{PT}}{f}_{\mathrm{g}}{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}{R}_{\mathrm{n},\mathrm{C}}\end{array}$$

where α_PT is the Priestley-Taylor coefficient, f_g is the fraction of vegetation that is green and thus transpiring, Δ is the slope of the vapour pressure to air temperature curve (hPa K⁻¹) and γ is the psychrometric constant (hPa K⁻¹). In all land-covers α_PT has an initial value of 1.26, except for coniferous forests where it is lowered to account for reduction in stem conductivity (and therefore stomatal conductance) with height (${\mathit{\alpha }}_{\mathrm{PT}}=max(-\mathrm{0.269}\ln (h_{\mathrm{c}})+\mathrm{1.31},\mathrm{0.4})$) and broadleaved and mixed forests where it is assigned a value of 0.82, both following Komatsu (2005). If unrealistic fluxes are obtained (λE_S<0 indicating condensation in the soil during daytime) during this first iteration, then the canopy transpiration (i.e., α_PT) is sequentially reduced and soil and canopy temperatures and fluxes are recalculated until realistic values are obtained. This implementation of the TSEB scheme is called the TSEB-PT model.

TSEB-PT outputs instantaneous fluxes at the time of thermal image acquisition. The modelled instantaneous latent heat flux (λE_inst), calculated during clear-sky conditions, can be extrapolated to daily ET values as $\mathit{\lambda }{E}_{\mathrm{daily}}=\mathit{\lambda }{E}_{\mathrm{inst}}\times \frac{S_{\mathrm{daily}}^{\downarrow }}{S_{\mathrm{inst}}^{\downarrow }}$ (Cammalleri et al., 2014), where $S_{\mathrm{daily}}^{\downarrow }$ and $S_{\mathrm{inst}}^{\downarrow }$ are the daily and instantaneous shortwave irradiances, respectively.

2.1.2 ETLook model

ETLook model (Bastiaanssen et al., 2012) is used in the WaPOR framework and is described in detail in the Methodology section of the WaPOR ETLook wiki (https://github.com/un-fao/wapor-et-look/wiki/Actual%20Evapotranspiration%20and%20Interception#methodology, last access: 22 January 2026) (FRAME Consortium, 2024). Similarly to TSEB-PT, ETLook is a two-source model, meaning that it derives soil evaporation (E) and vegetation transpiration (T) as two separate fluxes, and it ensures conservation of energy at the land-surface. The model assumes potential rates of daily E and T based on the Penman–Monteith equation (Monteith, 1965) that are throttled down to actual E and T using stress factors:

$$\begin{array}{}\text{(7a)}& {\displaystyle }E={\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n},\mathrm{S}}-G)+\mathit{\rho }{c}_p\frac{{\mathrm{\Delta }}_{\mathrm{e}}}{r_{\mathrm{a},\mathrm{S}}}}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\mathrm{S}}}{r_{\mathrm{a},\mathrm{S}}}\right)}}\text{(7b)}& {\displaystyle }T={\displaystyle \frac{\mathrm{\Delta }\left(R_{\mathrm{n},\mathrm{C}}\right)+\mathit{\rho }{c}_p\frac{{\mathrm{\Delta }}_{\mathrm{e}}}{r_{\mathrm{a},\mathrm{C}}}}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\mathrm{C}}}{r_{\mathrm{a},\mathrm{C}}}\right)}}\end{array}$$

where Δ_e (hPa) is vapor pressure deficit, ρ (kg m⁻³) is the air density, c_p (J kg⁻¹ K⁻¹) is specific heat of dry air, r_a,S and r_a,C are aerodynamic resistances for heat turbulent transport for soil and canopy respectively calculated following Allen (1998) and adjusted for buoyancy using Monin-Obukhov similarity theory in unstable conditions, and r_S and r_C are resistances of water transfer in respectively the soil and canopy. All resistances are in s m⁻¹.

The resistance of soil to water vapour transfer is calculated as

$$\begin{array}{}\text{(8)}& r_{\mathrm{S}}=r_{\mathrm{s},\min }{\left(S_{\mathrm{e}}^{\mathrm{top}}\right)}^c\end{array}$$

where r_s,min is minimum soil resistance with constant value of 800 (s m⁻¹), c is a power factor with a constant value of −2.1 and $S_{\mathrm{e}}^{\mathrm{top}}$ is relative top-soil soil moisture (Camillo and Gurney, 1986).

The canopy stomatal resistance is affected by air temperature stress (S_t), vapour pressure stress (S_v), radiation stress (S_r) and root-zone soil moisture stress (S_m)

$$\begin{array}{}\text{(9)}& r_{\mathrm{C}}=\left({\displaystyle \frac{r_{\mathrm{st},\min }}{{\mathrm{LAI}}_{\mathrm{eff}}}}\right)\left({\displaystyle \frac{\mathrm{1}}{S_{\mathrm{t}}{S}_{\mathrm{v}}{S}_{\mathrm{r}}{S}_{\mathrm{m}}}}\right)\end{array}$$

where r_st,min (s m⁻¹) is the minimum stomatal resistance and LAI_eff is effective leaf area index (Jarvis, 1976; Stewart, 1988).

In addition, ETLook estimates evaporation from rainfall intercepted by the canopy (I in mm d⁻¹):

$$\begin{array}{}\text{(10)}& I=\mathrm{0.2}\mathrm{LAI}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{\mathrm{FVC}\cdot P}{\mathrm{0.2}\mathrm{LAI}}}}\right)\end{array}$$

where P is precipitation, and FVC is fractional vegetation cover. The energy used in this process is subtracted from net radiation (R_n) before it is split into the soil and canopy components (R_n,S and R_n,C respectively) and used to estimate E and T in Eq. (7).

The main difference between WaPOR implementation of ETLook and Bastiaanssen et al. (2012) is the method of obtaining soil moisture required for E and T stress factors. In Bastiaanssen et al. (2012) it is based on microwave observations while in WaPOR ETLook it is derived using a trapezoid constructed in the LST – vegetation fractional cover (f_C) space (Yang et al., 2015). The trapezoid corner values are set based on theoretical calculations by inverting the Penman-Monteith equation for both dry and moist bare soil and vegetated conditions. Then the soil moisture of a pixel (representing both top-soil and root-zone) is estimated using the relative location of LST and f_C of that pixel within that trapezoid. The CLMS implementation of ETLook adopts the WaPOR methodology. The one major difference between the WaPOR and CLMS ETLook-derived fluxes is that the former includes I as a separate sub-product while the latter does not. To ensure consistency of ETLook model formulations, we pass rainfall estimates to the ETLook and let it calculate I internally, before adding it to E calculated in Eq. (7):

$$\begin{array}{}\text{(11)}& E_{\mathrm{ETLook},\mathrm{CLMS}}=E+I\end{array}$$

This aligns also with TSEB-PT model in which I is implicitly included in evaporation. E_ETLook,CLMS is referred to as evaporation from ETLook in the rest of this manuscript.

Apart from land-cover dependent parameters (see Sect. 2.2.6), ETLook uses other spatially-distributed but temporally constant parameters (e.g. theoretical albedo of bare soil and full vegetation cover used during soil moisture estimation) which are provided as maps distributed by FAO.

2.2 Input data sources

Modelling actual evapotranspiration with satellite observations and TSEB-PT or ETLook models is a complex task requiring diverse set of input forcing variables (Table 2). Those are derived from shortwave optical imagery needed to estimate biophysical properties of the surface (e.g. leaf area index – LAI and albedo), thermal infrared observations of land surface temperature which is the boundary condition for the land surface – air energy exchange and a proxy for root-zone soil moisture, weather forcing, which drive (e.g. solar irradiance) and modulate (e.g wind speed) the energy exchange between the land surface and the air, and ancillary data (e.g. digital elevation models and canopy height maps) that cannot be derived from the other data sources.

Table 2 Input forcing variables for Copernicus Land Monitoring Service actual evapotranspiration product. PAR indicates photosynthetically active radiation part of the spectrum (400–700 nm) while NIR indicates near infra-red part of the spectrum (700–2500 nm).

Copernicus products should, to the largest extent possible, be based on other Copernicus data. This is to ensure the free and open license conditions, long-term future continuity and consistency across the CLMS portfolio. In case of the CLMS ETa product, this means observations from optical (both shortwave and thermal infrared) sensors on board of Sentinel-3 satellites and Copernicus Digital Elevation Model (DEM – European Space Agency and Airbus, 2022), both available from the Copernicus Data Space Ecosystem (https://dataspace.copernicus.eu/, last access: 25 August 2025), and products provided by CLMS and Copernicus Atmosphere Monitoring Service (Peuch et al., 2022) (Table 3). The only exception is the canopy height map, which is derived from fusion of GEDI LiDAR measurements and Sentinel-2 imagery (Lang et al., 2023).

Table 3 Input data sources for the Copernicus Land Monitoring Service actual evapotranspiration near-real-time product.

The subsections below describe the pre-processing methods selected to convert the input data from Table 3 into ET model input forcing variables shown in Table 2.

2.2.1 Sentinel-3 cloud-masking and gap-filling

Biophysical characterization of the land-surface is based on imagery obtained by two sensors on-board Sentinel-3 satellites: shortwave Ocean and Land Colour Intrument (OLCI); and combined shortwave and thermal infrared Sea and Land Surface Temperature Radiometer (SLSTR). Both sensors operate in the optical spectral domain, which is blocked by clouds and therefore cloud detection and masking needs to be performed before further analysis of these data.

Atmospherically corrected and geolocated Top-Of-Canopy (TOC) reflectances derived from shortwave optical detectors on OLCI and SLSTR sensors are generated and distributed by the CLMS (Copernicus Land Monitoring Service, 2025b). The cloud mask of this product is based on quality flags from Level 1 SLSTR and OLCI products and IDEPIX approach (Wevers et al., 2022). No further cloud masking is performed during the ETa production and instead the recommendations on using the annotation flags from Sect. 5.1 of Product User Manual (Copernicus Land Monitoring Service, 2025d) are followed.

Land Surface Temperature (LST) product (SL_2_LST___) is a Level 2 product based on thermal observations by the SLSTR sensor. It comes with quality layers indicating the presence of clouds. As recommended in the Copernicus Sentinel-3 SLSTR Land User Handbook (https://sentiwiki.copernicus.eu/__attachments/a_0b084365 0153c5e0933699cf9259d22d5995022524c99c86c744555f7 ebe3af4/OMPC.ACR.HBK.002%20-%20Sentinel%203%20 SLSTR%20Land%20Handbook%202024%20-%201.4.pdf? cb=e7c5d241da616f170b23f3a61f552a4e, last access: 27 August 2026), the probabilistic cloud mask is used during ETa production. This cloud mask uses a semi-Bayesian approach by estimating a probability of a clear-sky using radiative transfer modelling and meteorological conditions at the time of satellite overpass and observational climatology (Bulgin et al., 2014).

The TOC reflectance is used to derive land surface biophysical traits (Sect. 2.2.2) and albedo (Sect. 2.2.3) and to sharpen the LST (Sect. 2.2.4). The reflectance values, as well as biophysical traits, usually show a clear seasonal cycle and spatial similarity, and have same or similar values in both cloudy and sunny conditions (e.g. leaf area index does not change day to day depending on cloudiness). Therefore, gaps in this data are highly suitable for filling using spatio-temporal gap-filling, e.g. StarFM (Gao et al., 2006), and smoothing methods, e.g. Whittaker-Eilers smoother used by WaPOR (Eilers, 2003). On the other hand, LST is highly variable in both time (the temperature can change significantly at short time intervals even though it also has a strong seasonal cycle) and space. In addition, LST under clouds is different to LST in clear-sky conditions (Abbasi et al., 2020). This makes gap-filling of LST particularly difficult and using inaccurately gap-filled values can lead to energy imbalance at the land surface. Therefore, LST is usually not gap-filled, especially if it is to be used as input into ETa models.

Figure 1 Widths and locations of swaths of Sentinel-3 OLCI and SLSTR instruments (courtesy of Donlon et al., 2012) (left) and maximum view zenith angles (VZA) of those swaths overlaid on an example OLCI (colour) and SLSTR (greyscale) image (right). Also indicated are potential limits of VZA on the western side of the swath (45 and 35°).

In the case of the CLMS ETa product, we are using TOC reflectance and LST observations acquired by sensors onboard the same satellite platform and thus observing the land surface (and the clouds) at the same time. Therefore, since LST is not gap-filled, there is also no need to gap-fill TOC reflectance acquired at the same time and location. However, the Sentinel-3 OLCI sensor has a swath width of 1270 km, while the nadir pointing scan of SLSTR sensor has a swath width of 1400 km (Donlon et al., 2012). Both sensors are tilted to reduce the sun glint effect and therefore larger part of the swath is located westwards of the satellite orbit path (Fig. 1, left).

The sensors are positioned such that their swaths align at the western edge, while on the eastern edge the coverage of OLCI ends at view zenith angle (VZA) of 24° while coverage of SLSTR ends at VZA of 35° (Fig. 1, right). When modelling ETa, LST observations acquired at high VZA (e.g. above 45° as in the Sen-ET approach; Guzinski et al., 2020) are usually omitted due to increased uncertainty caused by longer atmospheric path, LST anisotropy and larger pixel footprint. However, SLSTR observations from the eastern edge of the swath with VZA between 24 and 35° could still potentially be used if OLCI data for that part of the swath was gap-filled.

Figure 2 Mean number of monthly cloud-free SLSTR observations for period 2020–2024 (apart from March 2020 and December 2024) for March (top-left), June (top-right), September (bottom-left) and December (bottom-right). The titles of the sub-plots indicate view zenith angle (VZA) limits on the western and eastern edges of the swath (e.g. W45–E35 indicates a VZA limit of 45° on the western edge and 35° on the eastern edge). The last column of each panel shows the difference in cloud-free observations per month between limiting eastern swath to either 35° or 24°. Less than 3 observations per month means on average less than one observation per 10 d aggregation period (dekad), less than 6 per month means less than 2 per dekad, less than 12 per month means less than 4 per dekad.

To assess whether it is worth to gap-fill the missing part of OLCI's eastern swath, and whether maximum VZA should be limited to 45° or 35° (threshold mission requirement value for planned Copernicus Land Surface Temperature Monitoring mission; Koetz et al., 2021), we performed an analysis of number of monthly cloud-free Sentinel-3 LST acquisitions over Europe and Africa (Fig. 2). This analysis was performed using 5 years of data (2020–2024) for months of March, June, September and December (apart from March 2020 and December 2024) and while limiting VZA of the western edge of swath to either 45° or 35° and eastern edge to either 35° or 24°.

Figure 2 shows a clear seasonal trend in the number of cloud-free SLSTR observations. In June and September, most of Europe has sufficient number of cloud-free observations for each dekad (at least 2 but often more than 4), while in March and, especially in December, there are areas where less than two or even less than one cloud free observations are available in each dekad. In the equatorial region of Africa, there are predominantly cloudy conditions throughout the year but they shift north in June and September and south in December and March. The difference between number of monthly cloud-free observations with eastern VZA limited to 35° or 24° is in large majority of cases less than 2, except in the Sahara region and southern Africa where there is sufficient number of observations in either case. The difference between limiting western swath to either 45° or 35° is more pronounced and can push some regions to having less than 1 observation per dekad.

Based on this analysis, it was decided to limit the western VZA to 45° and eastern VZA to 24°. This means that the extra processing that would be required to gap-fill the OLCI swath is avoided. Even though gap-filling of input imagery will not be performed, the filling of cloud gaps in the produced ETa maps is required to fulfill the requirements of the ETa product. Therefore gap-filling of outputs will be performed as described in Sect. 2.3.

2.2.2 Biophysical traits

A sound and generalizable approach to obtain the biophysical traits required by the ETa models is the inversion of canopy radiative transfer models (RTM), such as the combined PROSPECT (Jacquemoud and Baret, 1990; Feret et al., 2008; Féret et al., 2017, 2021) and 4SAIL (Verhoef et al., 2007) RTMs. Together those models can simulate the canopy spectra from 400 to 2500 nm at every nm based on the inputs listed in Table 4 (Jacquemoud et al., 2009). One computationally efficient inversion method consists of training a regression model based on large number of PROSPECT and 4SAIL simulations (Weiss et al., 2000, 2002; Verrelst et al., 2012). Then this regression model is applied to the spectral imagery to obtain the biophysical traits maps.

Campbell (1990)

Table 4 List of biophysical traits and ancillary information required by PROSPECT-D+4SAIL radiative transfer models together with range of values which are simulated or set during the simulation. Ranges are based on information extracted from the LOPEX database (Hosgood et al., 1993) and the Sentinel-2 Biophysical ATBD (Weiss et al., 2020).

In case of CLMS ETa, we are using a vectorized version of Féret et al. (2017) PROSPECT-D and Verhoef et al. (2007) 4SAIL models (Nieto, 2025). This vectorized version allows us to run the RTMs very efficiently and thus enables us to generate specific regression models for each scene, given their actual observation and illumination conditions. This contrasts with other similar algorithms, such as the BiophysicalOp in the Sentinel-2 toolbox of Sentinel Application Platform (Djamai et al., 2019), which applies the same generic hybrid inversion approach to all scenes. The biophysical processing framework is summarized as:

1. Estimating the proportion of spectral direct and diffuse irradiance using the 6S model (Vermote et al., 1997) and mean sun zenith and azimuth angles together with mean aerosol optical thickness and total column water vapor of Sentinel-3 TOC scene.

2. Generating ca. 40 000 PROSPECT-D+4SAIL simulated spectra by:

   • a.

     Creating the same amount of Monte Carlo random samples (Saltelli et al., 1999) of biophysical traits and observation angles, based on prescribed plausible ranges listed in Table 4.

   • b.

     Running PROSPECT-D+4SAIL for each of these samples together with the generated proportion of spectral direct and diffuse irradiance and the mean solar angles of the scene.

   • c.

     Convolving the simulated narrowband spectra using the Sentinel-3 OLCI and SLSTR spectral response functions in order to obtain a set of 40 000 simulated Sentinel-3 TOC reflectances.

3. Training a random forest regression in which the dependent variables are the randomly generated biophysical traits and the explicative variables the correspondent simulated TOC reflectances and observation angles.

4. Applying the random forest regression to the actual TOC reflectance scene.

As outputs we obtain 8 products: leaf area index (LAI), Campbell (1990) mean leaf angle, leaf chlorophyll a+b concentration, leaf carotenoids concentration, leaf antocyanin concentration, brown pigments, leaf dry matter content, and leaf water content. The retrieved leaf pigments concentration gives us an idea about the canopy greenness, and we have observed an increase of Sentinel-3 estimated antocyanin concentration during vegetation curing in summer and the leaf senescence in fall. Indeed antocyanins are red pigments and thus may become dominant over other leaf pigments in those situations (Féret et al., 2017). For that reason, we can express the fraction of LAI that is green (f_g) using an empirical piecewise linear relation to the antocyanins (Ant):

$$\begin{array}{}\text{(12)}& f_{\mathrm{g}}=\left\{\begin{array}{ll}\mathrm{1},& \text{if Ant }\le \mathrm{5}\mathrm{\mu }{\text{g\hspace{0.17em}cm}}^{-\mathrm{2}}\\ \mathrm{1}-\mathrm{0.8}\frac{\mathrm{Ant}-\mathrm{5}}{\mathrm{20}},& \text{if Ant }>\mathrm{5}\mathrm{\mu }{\text{g\hspace{0.17em}cm}}^{-\mathrm{2}}\text{ and}\\ & \text{ Ant }\le \mathrm{25}\mathrm{\mu }{\text{g\hspace{0.17em}cm}}^{-\mathrm{2}}\\ \mathrm{0.2},& \text{if Ant }>\mathrm{25}\mathrm{\mu }{\text{g\hspace{0.17em}cm}}^{-\mathrm{2}}\end{array}\right.\end{array}$$

2.2.3 Albedo and net shortwave radiation

Net radiation (R_n) provides the energy that drives all other energy fluxes at the land-surface (including ET) and can be approximated as:

$$\begin{array}{}\text{(13)}& R_{\mathrm{n}}=S_{\mathrm{n}}+L_{\mathrm{n}}=S^{\downarrow }\left(\mathrm{1}-\mathit{\alpha }\right)+\mathit{ϵ}\left(L^{\downarrow }-\mathit{\sigma }{\mathrm{LST}}^{\mathrm{4}}\right)\end{array}$$

where S_n (S^↓) and L_n (L^↓) are the shortwave and longwave net radiation (incoming irradiances) respectively, α and ϵ are surface albedo and emissivity, respectively, LST is the Land Surface Temperature, and $\mathit{\sigma }\approx \mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ (W m⁻² K⁻⁴) is the Stefan-Boltzmann constant.

Considering the larger magnitude of shortwave irradiance (S^↓) compared to the longwave irradiance (L^↓), and the fact that L_n is usually computed internally by each ET model, we will focus on method for deriving S_n. In particular, this study concentrates on the canopy and leaf properties that influence albedo and radiation partitioning between soil and canopy. The albedo (α) is defined as the proportion of incident shortwave radiation that is reflected by the surface. The shortwave net radiation (S_n) is therefore the balance between the incident shortwave irradiance (S^↓) and the reflected shortwave radiance ($S^{\uparrow }=\mathit{\alpha }{S}^{\downarrow }$).

The spectral properties of the surface are key in determining the albedo. Leaves, due to their photosynthetic activity, absorb a large proportion of light due to the presence of chlorophylls, as well as other leaf pigments. On the other hand, soils can have a large range of albedo values, depending on their mineral composition, texture and topsoil moisture. Therefore, the albedo of a vegetated surface, and also radiation partitioning between soil and canopy, will depend not only on leaf chlorophyll concentration but also on canopy density and, in a lesser degree, on the soil albedo in situations of sparse vegetation or initial growth stages. Indeed, most of the Earth's surface show certain anisotropic behaviour when reflecting radiation i.e. it scatters different amounts of radiation depending on the scattering direction. Vegetation, as it is mainly composed by an array of leaves, is also affected by this anisotropic behaviour. Therefore, plants will reflect radiation differently depending on their structural characteristics as well as the illumination geometry (i.e. the solar position) and the scattering direction (i.e. the sensor position), changing their albedo with time.

To compute the net shortwave radiation and its partitioning between the canopy (Sn_C) and the soil (Sn_S) the model of Campbell and Norman (1998) is used. The key aspect of this model is the calculation of the transmitted shortwave radiation through the canopy (τ_C), which is wavelength dependent due to vegetation absorbing a greater portion of photosynthetically active radiation (PAR – 400–700 nm) than near infra-red (NIR – 700–2500 nm) wavelengths. τ_C is partitioned into two components (direct/diffuse) and in two spectral band (PAR/NIR).

$$\begin{array}{}\text{(14a)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{n},\mathrm{C}}=& \left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{C},\mathrm{DIR},\mathrm{PAR}}\right)\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{C},\mathrm{DIR},\mathrm{PAR}}\right){\mathrm{PAR}}_{\mathrm{DIR}}+\\ & \left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{C},\mathrm{DIR},\mathrm{NIR}}\right)\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{C},\mathrm{DIR},\mathrm{NIR}}\right){\mathrm{NIR}}_{\mathrm{DIR}}+\\ & \left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{C},\mathrm{DIF},\mathrm{PAR}}\right)\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{C},\mathrm{DIF},\mathrm{PAR}}\right){\mathrm{PAR}}_{\mathrm{DIF}}+\\ & \left(\mathrm{1}-{\mathit{\tau }}_{\mathrm{C},\mathrm{DIF},\mathrm{NIR}}\right)\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{C},\mathrm{DIF},\mathrm{NIR}}\right){\mathrm{NIR}}_{\mathrm{DIF}}\end{array}\text{(14b)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{n},\mathrm{S}}=& {\mathit{\tau }}_{\mathrm{C},\mathrm{DIR},\mathrm{PAR}}\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{S},\mathrm{PAR}}\right){\mathrm{PAR}}_{\mathrm{DIR}}+\\ & {\mathit{\tau }}_{\mathrm{C},\mathrm{DIR},\mathrm{NIR}}\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{S},\mathrm{NIR}}\right){\mathrm{NIR}}_{\mathrm{DIR}}+\\ & {\mathit{\tau }}_{\mathrm{C},\mathrm{DIF},\mathrm{PAR}}\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{S},\mathrm{PAR}}\right){\mathrm{PAR}}_{\mathrm{DIF}}+\\ & {\mathit{\tau }}_{\mathrm{C},\mathrm{DIF},\mathrm{NIR}}\left(\mathrm{1}-{\mathit{\rho }}_{\mathrm{S},\mathrm{NIR}}\right){\mathrm{NIR}}_{\mathrm{DIF}}\end{array}\end{array}$$

where ρ_C,DIR is the canopy directional-hemispherical reflectance, ρ_C,DIF is the canopy bihemispherical reflectance, τ_C,DIR canopy directional-hemispherical transmittance, τ_C,DIF canopy bihemispherical transmittance, which depend on leaf spectral properties (absortance) and canopy structure (LAI and leaf inclination distribution parameter; Campbell, 1990). On the other hand, ρ_S is the soil bihemispherical reflectance. In all cases τ and ρ are separated between the PAR and NIR regions of the solar spectrum. The calculation of canopy transmittances and reflectances is shown in Appendix A.

Campbell RTM requires the leaf absorptance ($\mathit{\zeta }=\mathrm{1}-\mathit{\rho }-\mathit{\tau }$) as input for calculating canopy transmittance and reflectance. To estimate this parameter, we are using the leaf traits retrievals from the biophysical processor described in Sect. 2.2.2. With the information of these traits, we can run the PROSPECT-D leaf RTM in forward mode to get a spectral reflectance (ρ) and transmittance (τ) that could then be integrated to the required broadband regions. However, running PROSPECT-D for a large array of pixels is too computationally expensive. For that reason, a PROSPECT-D emulator (Guzinski et al., 2021) was developed (Rivera et al., 2015): we generated broadband leaf reflectances and transmittances in the PAR and NIR from a large range of PROSPECT-D simulations, covering all plausible range of leaf traits from Table 4, which was then used to train a random forest model that relates the leaf traits to the broadband PAR and NIR leaf reflectance and transmittance. Figure 3 depicts the importance that each leaf trait has on the broadband reflectance emulator, showing that the most important pigment is the chlorophyll a+b concentration, due to its strong absorption of PAR radiation, followed by the leaf dry matter and water contents, which are the main absorbers of the NIR/SWIR radiation.

Figure 3 Mean permutation importance for leaf biochemical components in the PROSPECT-D leaf radiative transfer emulator. Cab – chlorophyll a+b concentration, Cw – leaf water content, Cm – dry matter content, Ant – antocyanins concentration, Cbrown – brown pigments concentration, Car – carotenoids concentration.

Soil reflectance has a smaller role in most situations, since the radiation reaching the ground is smaller due to the interception of light by the canopy. However, in sparse and semi-arid conditions, where vegetation is scarce or even not present, soil albedo plays a significant role. Furthermore, these semi-arid areas are usually characterized by brighter soils since they are composed by sands, salts and with a small fraction of organic matter. For that reason, Guzinski et al. (2023) developed a method to unmix the broadband soil reflectance (ρ_soil,λ) from the TOC reflectance:

$$\begin{array}{}\text{(15)}& {\mathit{\rho }}_{\mathrm{soil},\mathit{\lambda }}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{surface},\mathit{\lambda }}-\left[\mathrm{1}-P_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{v}}\right)\right]{\mathit{\rho }}_{\mathrm{leaf},\mathit{\lambda }}}{P_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{v}}\right)}}\end{array}$$

where P₀(θ_v) is the canopy gap fraction at the satellite observation angle θ_v computed by the Beer-Lambert law ($e^{-\mathit{\kappa }\left({\mathit{\theta }}_{\mathrm{v}}\right)\mathrm{LAI}}$), ρ_surface,λ is the broadband surface reflectance, and ρ_leaf,λ is the leaf broadband reflectance computed previously.

Table 5 Narrowband to Broadaband conversion coefficients for Sentinel-3A and Sentinel-3B platforms.

In order to convert the narrowband reflectances of the TOC product into broadband reflectances we have applied the approach proposed by Liang (2001), who derived a linear regression model between the satellite spectral bands and the broad bands corresponding to the PAR region (400–700 nm), NIR region (700–2500 nm) and the solar spectrum (SW – 400–2500 nm). Similarly to Liang (2001), to overcome the limited amount of in situ measurements of albedo, we again made use of the PROSPECT-D+4SAIL canopy RTM coupled with the 6S atmospheric RTM (Vermote et al., 1997) to run a large number of simulations covering all plausible illumination, atmospheric and canopy conditions. Once the RTM simulations are performed (at 1 nm step), we convolved the results to both Sentinel-3A and Sentinel-3B specific spectral response functions, and integrated these simulated spectra to the broadband regions defined above. Then, one multivariate linear regression model per broadband region and per platform is fitted with the broadband reflectances as dependent variables and the Sentinel-3 simulated TOC reflectances as explanatory variables. The resulting coefficients are shown in Table 5, showing the weight of each Sentinel-3 band has on the broadband reflectances. It is worth noting that the coefficients are consistent with the expected behaviour, since the spectral bands located in the PAR region have no influence on the NIR broadband reflectances, nor the spectral bands in the NIR/SWIR influence on the PAR broadband model.

2.2.4 Land surface temperature sharpening

The CLMS ETa product specification states a spatial resolution of 300 m. However, the spatial resolution of the SLSTR LST product is 1 km. Therefore, a method to perform thermal sharpening of the LST is required. Most such methods are based on machine learning approaches, of various complexities, which capture the relationships between shortwave spectral reflectance (and possibly other ancillary data) and the LST. The shortwave reflectance is of higher spatial resolution and is resampled to the resolution of LST to be used as explanatory variables for model training. Once trained, the model is applied to the reflectance at its original resolution to obtain a representation of LST at that resolution. In most approaches, this is followed by a bias correction step to ensure conservation of energy between the LST maps at both the original and sharpened spatial resolutions (Agam et al., 2007; Gao et al., 2012; Sánchez et al., 2020).

Both the Sen-ET and WaPOR (version 3) modelling frameworks use a Data Mining Sharpener (DMS) thermal sharpening approach (Gao et al., 2012) as implemented in the pyDMS Python package (https://github.com/radosuav/pyDMS, last access: 30 July 2025). It is a quite complex method that however works well even for sharpening by a ratio of 50 i.e. sharpening Sentinel-3 LST from 1 km to 20 m using Sentinel-2 reflectance (Guzinski and Nieto, 2019; Guzinski et al., 2023; Sánchez et al., 2024). Based on this work, we are applying it to sharpen Sentinel-3 LST using CLMS TOC product with 300 m spatial resolution. The DMS regression models are trained on the whole Sentinel-3 TOC tile (10° by 10°) as well as on subsets of 30 by 30 LST pixels in a moving window fashion.

In Sen-ET framework all relevant reflectance bands from Sentinel-2 or Sentinel-3 (depending on the target spatial resolution) covering visible, near-infrared and shortwave-infrared parts of the spectrum were used as explanatory variables. For Sentinel-3 this means the following 17 bands from CLMS TOC product: Oa02, Oa03, Oa04, Oa05, Oa06, Oa07, Oa08, Oa09, Oa10, Oa11, Oa12, Oa16, Oa17, Oa18, Oa21, S5N, S6N. This configuration of variables, in addition to DEM and cosine of solar inclination angle, is called “DMS – Reflectance” in Sect. 4.3.2. During the WaPOR project, feature engineering was performed and spectral bands were converted to indices before evaluating their usefulness in DMS (FRAME Consortium, 2024). This resulted in following 11 spectral bands, reflectance indices and DEM-related products being used as explanatory variables: Oa04 (blue), Oa17 (NIR), Modified Normalized Difference Water Index (MNDWI), Plant Senescence Reflectance Index (PSRI), Normalized Difference Moisture Index (NMDI), Visible Atmospherically Resistant Index Red Edge (VARI_RED_EDGE), Bare Soil Index (BI), elevation, cosine solar inclination angle, aspect and slope. More details and the list of evaluated indices are available in the WaPOR wiki (https://github.com/un-fao/wapor-et-look/wiki/Land%20Surface%20Temperature#processing-approach, last access: 22 January 2026). This combination of explanatory variables is called “DMS – WaPOR” in Sect. 4.3.2. Finally, since the main influence of aspect and slope on LST is already captured in the solar inclination angle variable and their additional impact on sharpened LST accuracy is negligible (see Sect. 4.3.2), those two variables were removed from the WaPOR list. The resulting combination of 9 variables (called “DMS – WaPOR selected” in Sect. 4.3.2) is used in the ETa processing chain to sharpen the 1 km Sentinel-3 LST to the required 300 m spatial resolution.

2.2.5 Weather forcing

Weather forcing is critical for accurate estimation of ET. Due to the 2 d timeliness requirements for the CLMS ETa product, the weather data source is Copernicus Atmosphere Monitoring Service (CAMS) forecasts (Peuch et al., 2022), produced by the European Center for Medium Range Weather Forecasts and distributed freely and openly through the CAMS Data Store. CAMS data contain surface meteorological parameters with forecast runs every 12 h covering the whole Earth on a 0.4° grid and hourly temporal resolution. In this study we used CAMS forecasts with 12 h lead time also for the production of historical data to ensure fair comparison with the NRT dataset.

Table 6 List of required CAMS global atmospheric composition forecast fields and their topographic correction status.

Instantaneous weather forcing at the satellite overpass are used to drive both ET models and include air temperature, vapor pressure, wind speed, surface pressure, and clear-sky solar irradiance (Table 6). All instantaneous data were obtained by linear interpolation between two CAMS hourly forecasts to the time of Sentinel-3 SLSTR acquisition over the area of interest. Daily weather forcing is used to drive the ETLook ET model and to extrapolate and interpolate the instantaneous estimates of ET and include solar irradiance as well as air and dew temperatures, wind speed, and pressure, which are then used to calculate the FAO-56 reference ET (Allen, 1998) required for the gap filling (see Sect. 2.3). They are being integrated over a 24 h period starting at midnight local time.

The pre-processing of CAMS weather forcing, including topographic correction, was done using the open source Python software meteo_utils (Nieto et al., 2025b) and as described in Guzinski et al. (2021). The only differences from Guzinski et al. (2021) being the use of Copernicus DEM (COP-DEM_GLO-90-DTED – European Space Agency and Airbus, 2022) resampled to a resolution of 300 m for topographic correction and the use of REST2 model (Gueymard, 2008) to estimate instantaneous clear sky solar irradiance and its PAR and NIR, beam and diffuse components. In addition daily total precipitation is calculated using 24 h integration of CAMS hourly total precipitation forecasts.

2.2.6 Structural and ancillary parameters

Resistance energy balance models need additional ancillary inputs, such as canopy height or roughness. The latter influences the efficiency of the turbulent transport of heat and water between the land surface and the overlying air (Raupach, 1994; Alfieri et al., 2019). Vegetation structure and density are thus important for estimating turbulent transport of momentum, heat and water vapour in the canopy air space (Garratt and Hicks, 1973; Thom, 1972; Raupach, 1994; Shaw and Pereira, 1982).

Table 7 Land Cover (LC) Look-Up-Table for ancillary and structural parameters required by ET models, adapted from Guzinski et al. (2020). h_min (m) is the minimum canopy height; h_max (m) is the maximum canopy height occurring when leaf area index (LAI) reaches its optimal maximum LAI_max (for annual plant functional types only); f_c is the at-nadir fraction of the ground occupied by a clumped canopy (f_c=1 for a homogeneous canopy); $w_{\mathrm{c}}/h_{\mathrm{c}}$ is the canopy shape parameter, representing the canopy width to canopy height ratio; l_w (m) is the average leaf size; and r_st (s m⁻¹) is the minimum stomatal resistance. LC classes come from the CLMS Dynamic Land Cover map. Note that snow/ice and water surfaces are masked in the current implementation.

For that reason, the CLMS Global Dynamic Land Cover map (Copernicus Land Monitoring Service, 2015) is used to assign vegetation parameters, which are difficult to estimate directly from other Earth Observation data (Guzinski et al., 2021). Those parameters, and values assigned to different plant functional types, are listed in Table 7. The TSEB-PT model requires all of the parameters, apart from stomatal resistance, while ETLook requires only vegetation height and stomatal resistance. Values of all parameters, except for vegetation height of forested land covers, were adapted from Guzinski et al. (2020).

Vegetation height is one of the most important of the ancillary parameters, especially for the TSEB model (Burchard-Levine et al., 2020), as it influences the aerodynamic resistance to heat transport. In order to better estimate the obstacle (canopy) height in different land covers we use a framework that differs depending on predominant plant functional type of each land cover. For annual plant functional types, canopy height is dynamically computed considering its growth as:

$$\begin{array}{}\text{(16)}& h_{\mathrm{c}}=h_{min}+(h_{max}-h_{min})\times min\left({\displaystyle \frac{\mathrm{LAI}}{{\mathrm{LAI}}_{max}}},\mathrm{\hspace{0.17em}}\mathrm{1}\right)\end{array}$$

where the symbols are described in Table 7.

For forest plant functional types we set a static canopy height based on a 10 m spatial resolution (averaged to 300 m resolution) global forest canopy height map developed by combining the Global Ecosystem Dynamics Investigation (GEDI) LiDAR and Sentinel-2 observations using a probabilistic deep learning model (Lang et al., 2023). Whenever the land cover map indicated a forest while the GEDI-based canopy height map was below the h_min parameter, the minimum value was enforced.

2.3 Output gap-filling and temporal compositing

The CLMS ETa product consists of 5 sub-products: instantaneous sensible and latent heat fluxes (in W m⁻²), and dekadal (10 d mean) evapotranspiration and its components evaporation and transpiration (in mm d⁻¹). The instantaneous heat fluxes represent values modeled at the time of Sentinel-3 satellite overpass and, therefore, do not undergo any gap-filling. On the other hand, the daily water fluxes (used for the production of the dekadal sub-products) need to undergo gap-filling. Otherwise, the dekadal mean would only take into account fluxes modeled during clear-sky conditions within the aggregation period. This would firstly result in frequent gaps and secondly in a systematic overestimation of the dekadal aggregate. The data availability is strongly time and region dependent, but the majority of the dekades have fewer than 6 daily ET estimates per pixel. To illustrate this, Fig. 4 shows the number of valid observations in the first dekad of July 2021 across the globe.

Figure 4 Number of land pixels (y-axis) having 1–10 observations on the first dekad of July 2021 per region of the world.

The gap-filling of ET product, i.e. estimation of ET during cloudy conditions, is usually performed using a reference quantity that can be derived for any date regardless of cloud conditions. This implies that this reference quantity is mostly dependent on weather forcing. In the Sen-ET approach the choice was made to use reference evapotranspiration calculated using the FAO-56 method (Allen, 1998). A ratio of modeled daily ETa to reference ET (called crop-stress coefficient K_s,c) is calculated on dates for which daily ETa is available and is used to recreate ETa on the target date which needs to be gap-filled, as described in Guzinski et al. (2021). This method was further developed in Guzinski et al. (2023) to better account for soil drying by performing linear interpolation of K_s,c (applicable only in non-NRT modelling when ETa after the target date is available) and for soil wetting through rainfall by setting K_s,c to the maximum value observed during the gap-filling period if a simple water-balance approach indicates sufficient moisture in the soil. Taking rainfall into account is important especially for longer gap-filling periods (Delogu et al., 2021) and in climates in which rainfall initiates the growing season (e.g. rainy season in the Sahel and other semi-arid areas). For the CLMS ETa processing chain we further improved the robustness of this method, especially for longer gap-filling windows, by replacing the maximum K_s,c observed within the gap-filling window with the 80th percentile of ratios observed within the gap-filling period. If the latest know K_s,c before the target date was larger than 80th percentile then the value of that last know K_s,c was preserved. This was done to ensure that taking rainfall into account would only increase, and not reduce, the gap-filled ET values. To accommodate longer periods with few satellite observations (see Fig. 2) a 60 d gap-filling window is used during ETa production.

Once the gap-filled ET is estimated, the split into evaporation and transpiration is performed using the ratio of evaporation or transpiration to ET from the closest non-gap-filled preceding date. The all-sky gap-filled daily estimates of evapotranspiration, evaporation and transpiration are then aggregated to dekadal timesteps by taking the mean of all valid values within the aggregation period. Finally the Ensemble dataset is computed by averaging the gap-filled evapotranspiration, evaporation and transpiration estimates from the TSEB-PT and ETLook models. The degree of gap-filling for a given dekad is indicated in the quality layers as the number of valid modelled clear-sky (i.e. not gap-filled) ET values within each dekad and the average temporal distance (in days) between the gap-filled ETs within a given dekad and the closest previous valid clear-sky modelled ET.

3 Prototype product validation

The main aim of the validation is to evaluate the performance of the prototype product through comparison of dekadal ETa values, as modelled by the TSEB-PT and ETLook algorithms, with ETa data derived from in-situ measurements in the period 2020–2022, from stations covering main worldwide climatic zones and plant functional types. In addition to assessing the performance of the two ETa models driven by input forcing as described in Sect. 2, the statistical metrics were also computed for an ETa dataset composed by the average of the TSEB-PT and ETLook estimates. This additional dataset is referred to as Ensemble ETa (Volk et al., 2024). The statistical metrics used in the analysis were: bias, root mean squared error (RMSE) and Pearson correlation coefficient (r). This comparison focused only on the variables that could be extracted from measurements at eddy covariance stations (ETa, λE, H). Evaporation and transpiration were generated by the models but could not be routinely contrasted against in situ observations and therefore accuracy of partition of ET into E and T is uncertain.

3.1 Validation data collection and preparation

The search for eddy covariance sites was intended to collect data representing the largest possible diversity of climate regions and plant functional types. Different eddy covariance networks and datasets were explored in search for validation data for recent years. The years 2020–2022 were selected for this analysis considering the overall in situ data availability.

When building the reference database we used the energy balance closure (EBC) (Foken et al., 2011) correction provided by the dataset, in order to ensure transparency and minimize scientific bias. Therefore, priority was given to datasets in which λE and H data had been corrected for EBC problem. A commonly applied procedure delivering corrected values of λE and H is the ONEflux processing which uses the Bowen Ratio correction method (Pastorello et al., 2020). Data generated by this procedure could be obtained from the AmeriFlux (FLUXNET dataset) (Novick et al., 2018), ICOS (Integrated Carbon Observations System – Heiskanen et al., 2022) and OzFlux (Isaac, 2014) networks. In search for more diversity in the reference datasets, other networks/datasets were queried even if the processing chain was not accounting for the EBC. Thus, data from the European Fluxes Database Cluster (EFDC – https://www.europe-fluxdata.eu/, last access: 1 September 2025), l'observatoire AMMA-CATCH (Analyse Multidisciplinaire de la Mousson Africaine – Couplage de l'Atmosphère Tropicale et du Cycle Hydrologique – Galle et al., 2018; AMMA-CATCH, 1990), AsiaFlux (Mizoguchi et al., 2009), JapanFlux (Ueyama et al., 2025), South African Environmental Observation Network (SAEON, https://observationsmonitor.saeon.ac.za, last access: 15 January 2026) and AmeriFlux (BASE dataset) (Chu et al., 2023) networks were also included in the analysis.

Figure 5 Eddy covariance sites and networks used for ETa validation.

Figure 6 Number of eddy covariance sites grouped per PFT (left) and per Köppen climate region (right). CRO: Cropland; CSH: Closed Shrubland; CVM: Crop Vegetation Mosaic; DBF: Deciduous Broadleaved Forest; EBF: Evergreen Broadleaved Forest; ENF: Evergreen Needleleaved Forest; GRA: Grassland; OSH: Open Shrubland; SAV: Savanna; URB: Urban; WET: Wetland; WSA: Woody Savanna.

Figure 5 shows the location of the eddy covariance sites considered in this study and the network/dataset the data were taken from, while Fig. 6 provides an overview of the representation of climate regions and plant functional types (PFT) in this set of eddy covariance sites. Further details of the location, climate, PFT and network of each site can be found in Appendix B.

Figures 5 and 6 show that the available eddy covariance networks are not equally distributed across the globe which results in imbalanced representativeness of climate regions and plant functional types. However, there are stations located in all major climate zones and PFTs and across all continents (except Antarctica). Despite lack on in-situ measurements in some larger geographical areas (mainly in Asia, Africa and South America), there are other sites which are located in similar climatic zones, e.g. tropics (northern Australia) as well as arid and semi-arid (Spain, southern Australia and western US). Therefore, we believe that the validation should be representative also of the geographical areas for which in-situ ET data is currently missing.

The data acquired from the AmeriFlux, ICOS and OzFlux networks was available at half-hourly and daily time steps and format aspects like variable naming, units, quality flags, etc. were uniform. Daily λE data were discarded if the value of the quality flag was lower than 0.6 (i.e. less than 60 % of sub-daily data was measured or had a good quality gap filling). A dekadal ETa value was computed if the dekad was composed of at least 7 valid daily ETa values. The computation of ETa from λE and air temperature was conducted as follows:

$$\begin{array}{}\text{(17)}& {\displaystyle }lv=\left(\mathrm{2.501}\mathrm{-}\mathrm{0.00237}{T}_{\mathrm{a}}\right)\times {\mathrm{10}}^{\mathrm{6}}\text{(18)}& {\displaystyle }\mathrm{ET}=\mathrm{3600}\times \mathrm{24}\times {\displaystyle \frac{\mathit{\lambda }E}{lv}}\end{array}$$

where lv is the latent heat of vaporization, T_a is the average daily air temperature, λE is the daily average latent heat flux and ET is the daily evapotranspiration in mm d⁻¹.

The data from other networks was delivered at half-hourly time step only. Therefore, an additional aggregation step needed to be considered to obtain daily ETa values. In doing this, the half-hourly data were filtered using quality flag (provided with the data) to remove invalid and poor quality values. For each day, the number of valid timeslots between sunrise and sunset was computed (the sunrise and sunset times change as function of geographic location and time of the year). Missing half-hourly data during the day were computed by linear interpolation if the number of valid timeslots during daytime was at least 50 % of the total number of timeslots in that period. Otherwise, the day was discarded. The criterion for aggregating the daily ETa values to dekadal values was the same as indicated in the previous paragraph.

The last preparatory step was matching the modelled and reference values on the basis of timeslot and location. The temporal matching is a straightforward step for the dekadal ETa values. The data at satellite overpass time (λE and H) were matched to the nearest timeslot of the eddy covariance half-hourly datasets. The satellite overpass time, in local time, varied along the year for the different sites. For instance, the AU-Cum site registered values between 09:05 and 10:52 local time (LT); US-Ton, between 10:03 and 11:52 LT; ES-LM1, between 11:27 and 13:13 LT; etc. The analysis for λE and H at satellite overpass time was conducted for TSEB-PT only as the ETLook algorithm does not generate those variables. Spatial matching was performed on the basis of selecting the 300 m pixel within which the flux tower is located. The spatial scale mismatch between the flux tower footprint and the pixel size could affect the reported accuracy of the products, especially at heterogeneous sites and during stable atmospheric conditions.

3.2 Validation results

The procedure described in the previous section resulted in a dataset of over 13 000 records (206 sites) of dekadal ETa values obtained from eddy covariance measurements. A first appraisal on the performance of the ETa models can be obtained from the statistical scores shown in Table 8. A graphical representation of the performance of the models under consideration is presented in the Taylor plot of Fig. 7. The overall statistics (i.e. not at site level) were generated from the full dataset of observed and modelled ETa pairs across all the validation sites.

Table 8 Bias (modelled − measured: mm d⁻¹), Root Mean Squared Error (RMSE: mm d⁻¹), and Pearson correlation coefficient (r) scores for dekadal ETa per model (all sites combined), and summary statistics at site level. rBias and rRMSE are relative metrics (i.e., divided by mean measured ET). SD stands for standard deviation.

Figure 7 Taylor-plot overview of ETa accuracies of TSEB-PT, ETLook and Ensemble dataset. Azimuthal angle represents the Pearson r, radial distance from the origin is the standard deviation (mm d⁻¹). Circle marker at the horizontal axis represents perfect model performance, markers closer to this point indicate better model performance. Radial distance from this point is equivalent to centered RMSE (contour values shown in mm d⁻¹). The statistics shown in the Taylor diagram were generated from the full dataset of observed and modelled ETa pairs across all the validation sites.

The points to be highlighted in these statistical scores are:

• All ETa models exhibit high correlation with the in situ observations (r above 0.75).

• Pairwise t-tests of the site-level validation results indicated that the Ensemble model had a significantly (p<0.05) higher r, lower RMSE and smaller bias, compared to TSEB-PT and ETLook.

• The bias values of all three datasets are small (relative bias below 4 %) with Ensemble showing negligible bias.

Around 82 % of the data records used in calculating the scores of Table 8 and Fig. 7 had been corrected for the EBC problem. The same analysis was conducted for the subset of data records composed only of EBC-corrected data. The overall results were very similar to those presented above (Fig. C1 and Table C1), although not addressing the EBC issue could have large impact when looking at individual sites.

Figure 8 Boxplots of r (left), RMSE (center) and bias (right) between in-situ observations and values of dekadal ET as computed by ETLook and TSEB-PT and the Ensemble dataset for years 2020–2022, for sites grouped per climates (top) or per plan functional types (bottom). The horizontal line indicates the median value, the box limits show the first and third quartiles, while the whiskers give the minimum and maximum value.

The plots of Fig. 8 show a more detailed view of these results, split into climate regions (based on Köppen climate classification system) and plant functional type classes. Complementary Taylor plots are available in Appendix C (Fig. C2 for climate regions and Fig. C3 for PFTs).

While ETLook and TSEB-PT have a different modelling approach, outputs of both models achieved comparable results in many cases. The CRO, GRA and WET and shrublands (CSH, OSH) Taylor plots reveal that both models exhibited similar performance in the generation of dekadal ETa values in those groups. The five classes together account for 49 % of the evaluated sites spread across different climate regimes and represent a wide range of conditions. The Ensemble ETa series exhibits slightly better correlation scores than the individual models.

Differentiating between climates, the r-score of the Ensemble was significantly better in continental and temperate climate regions, compared to dry and tropical regions (t-test p value < 0.05). In terms of bias, Ensemble ET estimates in dry climate regions had a significantly larger overestimation, compared to continental and temperate sites. This was also the case for tropical sites, though not significantly.

The sites in tropical and dry climates demonstrated the largest variability in model output performance and inter-model disagreement. This can also be observed in the classification according to PFT, where EBF, SAF, OSH and WSA exhibit a low and highly variable r. In tropical sites (and EBF sites), this is partly explained by the lack of pronounced seasonal cycle. In dry climates, this underlines the challenge to estimate ET under water-restricted conditions and the impact of different modelling approaches. Overall, the best scores were obtained with the Ensemble in tropical and dry areas, though a notable exception is the significantly smaller bias of TSEB-PT in the dry climate (0.17 mm d⁻¹) compared to ETLook (0.31 mm d⁻¹) and the ensemble (0.23 mm d⁻¹). TSEB-PT also has the smallest bias in the tropics, though not significantly better than the other models.

Looking at PFTs, all models achieved the highest correlation in DBF (Ensemble median r: 0.92). However, this is not fully confirmed in the other validation metrics. Notably, the RMSE of TSEB-PT for this class is high (0.93 mm d⁻¹), and significantly higher compared to ETLook and the Ensemble (0.72 and 0.74 mm d⁻¹, respectively). Similar behavior was found in MF, where ETLook had a significantly lower RMSE, compared to TSEB-PT. This aligns with the underestimated variability of the fluxes (standard deviation in Fig. C3) in DBF and MF by TSEB-PT.

Conversely, ETLook had a significantly higher RMSE in SAV, compared to TSEB-PT. This corresponds to a large (positive) bias in ETLook for this class, whereas TSEB-PT has a small (negative) bias.

Some other classes also exhibit interesting behavior (e.g. contrasting RMSE and bias in DNF and WSA between models) but might be misleading due to small sample size (see Fig. 6) and are therefore not discussed here.

Table 9 Bias (modelled − measured: W m⁻²), Root Mean Squared Error (RMSE: W m⁻²), and Pearson correlation coefficient (r) scores for instantaneous TSEB-PT λE and H (all sites combined), and summary statistics at site level. rBias and rRMSE are relative metrics (i.e., divided by mean measured λE or H).

Figure 9 Boxplots of r (left panels) and bias (right panels) values of λE (top panels) and H (bottom panels) at overpass time as modelled by the TSEB-PT model. The horizontal line indicates the median value, the box limits show the first and third quartiles, while the whiskers give the minimum and maximum value.

Table 9 shows the overall accuracy statistics and summary scores at site levels for λE and H, while Fig. 9 shows the range of r and bias values obtained for each climate and PFT class. A number of aspects in this plot can be connected to the analysis of dekadal ETa presented above. For instance, the large positive bias and low correlation of λE in the evergreen broadleaved forest (EBF) classes that is translated to the dekadal ET performance of Fig. 8. However, looking at the performance of H in EBF we see that TSEB-PT is explaining fairly well the variability in H, with mean correlation values of 0.6 and bias near 0 W m⁻². An opposite behavior occurs for ENF, with a systematic underestimation of H but no bias in λE. It is also notable that the median bias in λE was positive and small (below 50 W m⁻²) in the majority of the sites while H estimates show a dominant negative bias with a larger magnitude, as can also be seen in timeseries modelled at the individual sites (Fig. C4). The length of the boxes in the boxplots of Fig. 9 suggests an important degree of variability in the correlation and error of the modelled values among the sites of each PFT class. Apart from modelling uncertainty, the differences in the time window at which satellite overpass takes place, the uncertainty associated with the in situ energy balance closure correction, and the large heterogeneity in ecosystem properties (e.g. heat storage) within each class can partly explain this variability and the differences between biases seen in H and λE.

Validation was not performed in urban areas since neither of the models is designed to estimate fluxes in those conditions. Urban areas present specific challenges when applying energy balance models due to, among other things, strong shadowing effects of buildings surrounding vegetated patches, increased complexity in modelling turbulent transport between surface and atmospheric boundary layer caused by the interaction between buildings of different heights and vegetated surfaces, and additional heat and water vapour sources interacting between themselves and with the atmosphere: buildings, paved roads, parks with different canopies coexisting (grass/trees). Therefore, a specific modelling approach for urban areas would be required to accurately account for those processes.

4 Discussion

4.1 Comparison between the TSEB-PT, ETLook and Ensemble datasets

The validation assessment presented in Sect. 3 was based on the comparison of TSEB-PT, ETLook and Ensemble ETa estimates with measurements at eddy covariance towers and presented as one set of statistics for the whole validation period. However, the view of the spatial and temporal patterns in the ETa calculation by the ETLook and TSEB-PT models can give interesting insights towards the design of an operational ETa product.

Figure 10 Differences (ETLook − TSEB-PT: mm d⁻¹) between estimates of ET, E and T of the two models on the first dekad of January and July 2021.

The global distribution of differences between the two models during two selected dekads (northern hemisphere winter and northern hemisphere summer) can be seen in Fig. 10 with absolute bias and values of ET on the same dekads shown in Appendix C (Figs. C5 and C6 respectively). A complementary view of the differences and similarities between the two models is also presented in Appendix C (Fig. C7), which shows histograms of the dekadal ETa values (for the same two dekads) as modelled by ETLook and TSEB-PT for a number of 10° lat/lon tiles capturing different climatic zones.

Looking at the figures it is clear that the differences between ETLook and TSEB-PT are more pronounced in July. This is expected as northern hemisphere summer is the main global growing season and thus the magnitude of the fluxes is the largest. During January the main differences appear in savannas, shrublands and grasslands in the southern parts of the globe, e.g. Sahel and southern Africa, parts of Brazil and Argentina and Australia. In those areas, ETLook yields higher ET values due to a higher estimation of transpiration by that model. Observable differences, but of smaller magnitude, are also present in a band spanning latitudes of around 20 to 40° N and the Andes, where TSEB-PT produces higher ET values due to higher estimation of E. In July the picture is mixed with larger ET estimates by TSEB-PT at both higher latitudes (both northern and southern) and higher altitudes (e.g. Tibetan plateau and Andes), and ETLook estimating higher rates in semi-arid areas. The higher ET estimated by TSEB-PT can be mainly attributed to higher evaporation estimates by that model, and conversely higher ET estimates by ETLook can be attributed to higher transpiration estimates. In tropical rainforests (Amazon, Congo, Borneo) and agricultural areas (e.g. central and eastern Europe, parts of India and China) the differences between the ET (and E and T) from the two models are minimal in both dekads.

The above observations are supported by the tower-based validation results in Sect. 3.2. It was demonstrated that ETLook and TSEB-PT generally achieve similar performance, and that the main differences emerge in dry climates and forests in temperate and continental climates.

Dry climate regions and savannas sites (SAV, WSA) were generally characterized by a stronger positive bias in ETLook, compared to TSEB-PT. Furthermore, TSEB-PT displayed higher correlation and lower RMSE than ETLook. Still, the Ensemble was generally outperforming TSEB-PT. Those regions are water, not radiation, limited and therefore in ETLook the soil moisture stress factor plays an important role. Previous studies have shown that the WaPOR ETa dataset, based on the same ETLook model implementation as used in this study, tends to overestimate ET in such regions due to likely overestimation of the internally derived soil moisture (Blatchford et al., 2020; Cogill et al., 2025). Those studies also highlighted the difficulty of ETLook in capturing the ET dynamics in tropical (humid) regions due to poorer performance of Penman–Monteith based models in areas with low vapour pressure deficit. This can be observed in the current study too with much lower correlation of ETLook in tropical climate compared to TSEB-PT, despite having comparable RMSE and (to smaller extent) bias in this climate.

In forest sites (DBF, ENF, MF), ETLook in general showed higher variability (standard deviation) of fluxes and closer to that of flux-towers, resulting in a higher accuracy. In most cases (DNF, EBF, ENF) the Ensemble performed better than the two individual models but in DBF and MF ETLook showed the best overall model performance. The lower performance of TSEB-PT in those cases might be caused by two reasons: (i) its relative higher sensitivity to canopy height (Burchard-Levine et al., 2020), and (ii) its dependence of a first-guess potential transpiration based on the Priestley-Taylor parameter, which in case of broadleaved and needleleaved forests might deviate from (respectively) the fixed value of 0.82 and the empirical relationship with canopy height described by Komatsu (2005). In deciduous needleleaved forest (DNF) in particular the bias and RMSE of TSEB-PT is much higher than that of ETLook and this could point to poor applicability of Komatsu (2005) relationship in that type of forest caused by only two out of 45 needeleaved sites used to derive it being located in this PFT.

Figure 11 Dekadal ET, E and T as modelled by TSEB-PT and ETLook and Ensemble product and dekadal ET from the eddy covariance towers in five validation sites covering different PFTs. The Tower line represents EBC corrected value.

Looking at the timeseries of selected sites (Fig. 11), covering different PFTs, again supports the above conclusions. At the two of the sites where TSEB-PT and ETLook show similar performance against in-situ measurements (MF BE-Vie and CRO FR-EM2) both models also show very good agreement in the T estimates, although at the FR-EM2 site the E estimate of TSEB-PT is significantly higher than that of ETLook. Conversely, at the three sites where models show divergent behavior against EC measurements (EBF Au-Whr, SAV ES-LM1 and GRA US-Var), there is a large disagreement in T estimates with ETLook producing significantly higher values. In addition, at ES-LM1 the E estimate of TSEB-PT is significantly higher and at both ES-LM1 and US-Var the temporal patterns of E produced by the two models are inverted (TSEB-PT shows higher values in winter and spring and ETLook shows higher values in summer). In case of ES-LM1 and US-Var, TSEB-PT is able to capture a sudden drop of ET, which could be a result of agricultural practices or dry conditions, while ETLook continues to estimate high T and ET for the duration of the growing season. On the other hand, at AU-Whr TSEB-PT models a similar drop in ET and T, which results in negative correlation with tower measurements, while ETLook overestimates ET and T but captures well the relative temporal changes. Finally, Au-Whr and ES-LM1 are good examples in which the errors of the two models cancel each other out in the Ensemble dataset that results in good agreement with EC measurements, while in US-Var the Ensemble shows worse performance compared to TSEB-PT outputs.

Detailed attribution of the behavior described above to the different model assumptions and formulations is beyond the scope of this study. However, some general observations can already be made. Firstly, the partition of ET into E and T is to a large extent guided by the partition on net radiation into the soil and canopy components respectively. Here, TSEB-PT employs a more complex radiation transfer model (Kustas and Norman, 1999; Campbell and Norman, 1998) that accounts for differential spectral behavior in the PAR, NIR spectral regions, while ETLook relies on a simpler partition based on LAI and Beer-Lambert law (FRAME Consortium, 2024). In addition the longwave component of net radiation (L_n) is also computed differently. In the TSEB-PT model, L_n is computed according to Eq. (2a) from Kustas and Norman (1999) using CAMS “surface thermal radiation downwards” field together with internally computed leaf and soil temperatures, constant emissivity values for the leaf and soil components (0.98 and 0.95 respectively) and effective leaf area index. ETLook estimates it following the FAO-56 approach and daily values of air temperature, vapour pressure and transmissivity (FRAME Consortium, 2024).

Another aspect is the partition of the net radiation of soil and canopy into latent and sensible heat fluxes (and ground heat flux in case of soil). In TSEB-PT the transpiration is initially assumed to be at the potential rate and this is later throttled down if invalid values of other fluxes are obtained (see Sect. 2.1.1). In forests the potential ET is set to a lower values, depending on forest height and type, based on the empirical model developed by Komatsu (2005). This model might not be applicable to all the different forest types and during all phenological stages which might lead to low variability (standard deviation) of the observed fluxes in most forest classes (see Fig. C3). Taking AU-Whr as an example, it is an open eucalyptus forest and classified as evergreen broadleaved forest and thus assigned a constant α_PT parameter of 0.82. However, this value might not capture the different transpiration rates of the understory or specific adaptations of tree-grass ecosystems in dry climates. In the TSEP-PT model, evaporation is estimated as a residual of the other fluxes. For example in FR-EM2 TSEB-PT E has high values in spring and autumn, leading to overestimation of TSEB-PT ET during those seasons, and this could be caused by increased uncertainty in modelling the other energy fluxes of the soil at the time when it is exposed.

In ETLook the partition between latent and sensible heat fluxes of soil and canopy depends on the environmental stress factors and resistances which are based on meteorological and structural conditions as well as soil moisture (SM). The soil moisture is derived using LST (see Sect. 2.1.2) and therefore in vegetated conditions (when satellite mainly observes the canopy) it is more closely connected to root-zone moisture while in bare soil conditions it is connected mainly to top-soil moisture. However, the same SM estimate is used for both E and T calculations and therefore could lead to increased uncertainty of E or T depending on the canopy cover. In Fig. 11 it can be observed that T is overestimated in the three sites located in arid conditions (AU-Whr, ES-LM1 and US-Var) which could point to problems with the root-zone SM estimation in such conditions. The evaporation estimated by ETLook tends to display both low variability and low values. This again could be linked to uncertainties in the estimation of top-soil moisture and to the use of constant parameterization (regardless of soil type) when calculating soil resistance (Eq. 8).

Since the two ET models take different approaches, with different assumptions and approximations, they perform differently in various PFTs, climates and even phenological stages. The Ensemble (average) of the two provides more accurate ET than either of the models in most cases due to cancellation of random (and sometimes systematic) errors. However, there are also sites and PFTs where individual models might obtain better accuracy statistics. Therefore, as part of the CLMS ETa product users have access to individual TSEB-PT and ETLook model outputs as well as the Ensemble dataset. In addition, per-pixel standard deviation between the dekadal ET estimated by the two models is also included in the product as an ancillary quality layer since it can provide information on the uncertainty of the ET estimates due to model formulations. Finally, the TSEB-PT model also produces instantaneous latent and sensible heat fluxes which are available as part of the CLMS ETa product. Those heat fluxes do not undergo gap-filling and are distributed with a daily time-step and therefore might be preferred for applications such as fire risk/spread monitoring which are time-critical and in which gap-filling might produce increased uncertainty due to potential of significant changes at the land surface.

4.2 Comparison with WaPOR and OpenET ETa products

WaPOR global ETa dataset is the only existing operational and global ETa product with specifications very similar to the CLMS ETa. Furthermore, given the particular interest of FAO in the quality of the upcoming CLMS ETa product, WaPOR dekadal ETa values (WaPOR Level 1 version 3 – https://data.apps.fao.org/catalog//iso/7f4a7339-d56e-4393-8712-a8ffeffe2731, last access: 4 August 2025) were extracted for the study sites and timeslots and included in the analysis presented in Table 10. Based on those figures, the accuracy of the CLMS Ensemble ETa is significantly improved compared to WaPOR ETa (especially bias), while the individual model runs of TSEB-PT or ETLook have smaller bias but slightly lower correlation than WaPOR. It is also worth noting that the bias values of TSEB-PT and the Ensemble ETa are positive whereas ETLook and WaPOR exhibited a negative bias.

Table 10 Bias (modelled − measured: mm d⁻¹), Root Mean Squared Error (RMSE: mm d⁻¹), and Pearson correlation coefficient (r) scores for dekadal ET for the CLMS products and WaPOR. rRMSE and rBias are the relative RMSE and relative Bias (i.e. RMSE and Bias divided by mean measured ET).

When comparing global spatial patterns (Fig. C8), it can be observed that compared to the CLMS datasets, WaPOR generally produces higher T estimates and lower E estimates. As expected the differences are smallest when comparing ETLook and WaPOR since the model setup should be fairly similar (the operational production setup of WaPOR is not public although FAO has created equivalent open-source package which produces similar results when run with the same input data). Therefore, the differences can be attributed mainly to the input data. Here, four possible reasons appear likely:

• While CLMS relies on DMS to improve the spatial resolution of 1 km Sentinel-3 LST, WaPOR (version 3) uses the LST observations acquired by the VIIRS sensor on board of Suomi-NPP satellite with 375 m spatial resolution.

• The CLMS ETa dataset was produced in NRT mode, meaning that forecast meteorological forcing were used and gap-filling was performed using only preceding dates. On the other hand, historical WaPOR ETa used in this study is a reanalysis product which means that reanalysis meteorological forcing was used and gap-filling was performed using both preceding and succeeding dates.

• WaPOR ETLook implementation relies of some temporarily-static – spatially-distributed layers which parameterize the model. Those layers were used in CLMS ETLook implementation whenever available but there are still some layers which are not publicly disclosed in which case default constant values were used.

• WaPOR dataset explicitly separates evaporation from canopy interception into a sub-product while in the CLMS dataset (for both ETLook and TSEB-PT) it is bundled together with E.

OpenET (Melton et al., 2022) is another dataset with which CLMS ETa can be compared. Although OpenET is produced with much higher spatial (30 m) and temporal (daily) resolutions and only in the United States, similarly to CLMS dataset it contains a product which is an average (ensemble) of individual ET models (6 in case of OpenET). In a recent validation study, the RMSE of monthly ensemble product was between 12 % and 30 % lower than that of individual models while r² increased from 0.83–0.87 for individual models to 0.9 for ensemble (Volk et al., 2024). In our case the RMSE of dekadal Ensemble product was 10 % lower than that of TSEB-PT and ETLook models, while r increased from 0.76/0.77 for individual models to 0.81 for Ensemble product. In another study, the daily mean Ensemble OpenET product had a RMSE of 0.96 mm d⁻¹, bias of −0.2 mm d⁻¹ and r² of 0.84 (Melton et al., 2022) which is very similar to the results presented in Table 8 for the Ensemble dekadal ETa despite the increased uncertainty when validating 300 m product due to spatial-scale mismatch between flux tower footprint and pixels size.

Further comparison of CLMS ETa product, and its E and T components, and other global ETa datasets can be found in Sect. 4.5 of the product Validation Report (Barrios et al., 2025).

4.3 Evaluation of CLMS ETa modelling framework components

Section 3 focused on the validation of the final CLMS ETa model outputs. During the development of the ETa modelling framework various design choices were made and intermediate products evaluated to justify those choices. In this section, we briefly present this evaluation.

4.3.1 Biophysical traits, albedo and net shortwave radiation

The biophysical traits were obtained by random-forest inversion of PROSPECT-D+4SAIL RTM (Sect. 2.2.2). In order to test the sensitivity and robustness of such inversion method, we ran an independent set of PROSPECT-D+4SAIL simulations to compare how the regression model predicts the biophysical traits. Random white noise was added to this test dataset, considering that retrieved TOC have relative uncertainty of 10 % (${\mathit{\rho }}_{\mathrm{test}}={\mathit{\rho }}_{\mathrm{ProSAIL}}[\mathrm{1}+\mathcal{N}(\mathrm{0},\mathrm{0.1}\left)\right]$). This is shown as an example in Table 11.

Table 11 Evaluation performance for the Random Forest hybrid inversion of the PROSPECT-D+4SAIL radiative transfer model. The “observed” dataset corresponds to 40 000 independent simulation of PROSPECT-D+4SAIL for a VZA = 0°, SZA = 37.5°, a standard atmosphere, and a relative uncertainty of 10 % in the TOC reflectance retrievals. Cab, Car, Ant, Cbrown, are respectively the leaf concentrations of chlorophyll a+b, carotenoids, anthocyanins and brown pigments; Cm and Cw are respectively the leaf dry matter and water contents; LAI is the leaf area index and Leaf Angle is the Campbell (1990) mean leaf inclination angle. Physical units in the error metrics are consistent with those on Table 4.

A more computationally efficient alternative for the retrieval of LAI and f_g (but not pigments used for albedo estimation) would be to use the daily estimates of LAI and fraction of absorbed photosynthetically active radiation (fAPAR), which are internal CLMS datasets used to produce the 300 m dekadal biophysical CLMS product. In this case we could use a simple relationship between f_g, LAI and fAPAR (Fisher et al., 2008) in order to derive the fraction of LAI that is green (f_g).

$$\begin{array}{}\text{(19)}& {\displaystyle }{f}_{\mathrm{g}}={\displaystyle \frac{\mathrm{fAPAR}}{\mathrm{fIPAR}}}\text{(20)}& {\displaystyle }\mathrm{fIPAR}=\mathrm{1}-e^{K\left({\mathit{\theta }}_{\mathrm{s}}\right)\mathrm{LAI}}\end{array}$$

where LAI is the total Leaf Area Index $\left(\mathrm{LAI}=\frac{\mathrm{gLAI}}{f_{\mathrm{g}}}\right)$, fIPAR is fraction of intercepted photosynthetically active radiation and K(θ_s) is the shortwave beam coefficient of extinction at solar zenith angle θ_s.

Since CLMS LAI product actually represents the green LAI (Copernicus Land Monitoring Service, 2022), and with f_g and total LAI unknown, an iterative procedure proposed by Guzinski et al. (2020) is performed to find the optimal value of f_g based on gLAI and fAPAR. An initial LAI is assumed equal to gLAI (i.e. f_g=1), from which fIPAR is computed from Eq. (20) and then f_g recalculated with Eq. (19). This process is repeated until the f_g value converges between iterations.

Figure 12 Density scatterplot intercomparison between the biophysical processor LAI, green LAI (gLAI) and f_g (x-axis) and those derived from the alternative use of CLMS LAI/FAPAR 300 m global products, version 1 (y-axis).

In order to evaluate whether both approaches provide consistent data, we ran the biophysical processor over sites included in the ICOS WarmWinter2020 database (Warm Winter 2020 Team et al., 2022). These selected sites are listed in Table B5. We thus compared both LAI and f_g retrieved between 2019 and 2021 using the method described in Sect. 2.2.2 against the Fisher et al. (2008) and Guzinski et al. (2020) LAI/fAPAR approach using the CLMS FAPAR (Copernicus Land Monitoring Service, 2017a) and LAI (Copernicus Land Monitoring Service, 2017b) global products at 300 m, version 1. The density plots of Fig. 12 shows that overall the LAI products agree well, in particular at values lower than 2–3, with best agreement for green LAI (gLAI). Nonetheless, the larger scatter at higher LAI (i.e. denser vegetation) is not of great concern, as over these very dense canopies the interception (transmission) of radiation is already close to the maximum (minimum), i.e. near the light saturation, and thus these uncertainties have a minimal effect on ETa modelling. It is worth noting these results do not evaluate whether any of the two approaches is better than the other. Indeed, the CLMS gLAI has been intensively validated with a wide dataset of in situ LAI measurements (Copernicus Land Monitoring Service, 2025a) and thus it can be trusted with great confidence.

Figure 13 Timeseries intercomparison between the LAI (in black), gLAI (in green in LAI sub-panel) and f_g (in green in FG sub-panel) products retrieved from the biophysical processor (plain line) and those derived from the alternative use of CLMS LAI/fAPAR products (hollow circles).

However, the fraction of LAI that is green (f_g) shows larger scatter and with the point cloud in Fig. 12 far from the 1 : 1 line. This deviation has indeed an effect on the comparison of total LAI with larger scatter and a slight positive bias towards the total LAI derived from the CLMS LAI/fAPAR. In order to better understand the different behavior of f_g, Fig. 13 shows the timeseries for LAI, gLAI and f_g over selected representative ICOS sites. These sites include both temperate and semi-arid conditions and a wide range of biomes: broadleaved and conifer forests, croplands, grasslands, savannas, and orchards/vineyards.

These timeseries confirm that LAI trends are mostly in agreement, with only some bias present in very dense vegetation: croplands during their peak development, broadleaved forests in summer and conifer forests. However, the behavior of f_g shows discrepancies in certain cases. For instance, we would expect that herbaceous croplands (Fig. 13a) remain mostly green during the growing phase (spring) with f_g values very close to 1 which then decrease during crop senescence in summer. However, thef_g data derived from the CLMS LAI/fAPAR products shows values significantly lower than one during spring, and sometimes values close to 1 when LAI decreases during senescence. Another inconsistent behaviour was found in the savanna site (Fig. 13e), in which the CLMS LAI/fAPAR f_g even shows an opposite trend as one would expect i.e. decrease of f_g in late spring reaching a minimum in summer where most of the grass layer in this site is dead and only the evergreen oak canopy remains green and then a re-greening with the first rains of autumn. In addition, the f_g derived with the CLMS LAI/fAPAR product seems to be also underestimated in the temperate grassland (Fig. 13b) and evergreen forest (Fig. 13d), as over these two temperate biomes the canopy should remain mostly green all year round.

PROSPECT-D model, through an emulator, was also used to derive leaf bihemispherical reflectances and transmittance (Sect. 2.2.3). The evaluation of the performance of this emulator is shown in Fig. 14. There is a larger scatter for the NIR reflectances and transmittances, likely due to the fact that we are intentionally excluding the PROSPECT-D leaf structural parameter, that basically controls the multiple scattering within the leaf tissues, which is of a larger magnitude in the NIR region than in the PAR. However, the uncertainties when deriving the reflectances and transmittances seem to be cancelled out when computing the leaf absorptance and therefore the conversion from pigments to leaf spectra seems sufficiently robust.

Figure 14 Emulator of PROSPECT-D leaf radiative transfer model for the retrieval of broadband leaf reflectance, transmittance and absorptance factors.

To derive soil bihemispherical reflectance we converted narrowband reflectance values measured by Sentinel-3 satellites into broadband reflectance in PAR and NIR spectral regions using a linear regression approach (Sect. 2.2.3). We evaluated the goodness of this approach by running an independent set of simulations and compared the retrieved broadbands using the coefficients of Table 5. For both platforms, the evaluation confirms the robustness of conversion from the Sentinel-3 TOC reflectances to the broadband PAR, NIR and SW spectral regions, with negligible mean bias (≈ 0), very low RMSE and very good correlation (≈ 1) (Fig. 15).

Figure 15 Evaluation of the narrowband to broadband conversion for the estimated Sentinel-3A (top) and Sentinel-3B (bottom) coefficients using the Liang (2001) method. These results are obtained after applying the coefficents in Table 5 to 40 000 PROSPECT-D+4SAIL independent simulations.

4.3.2 Land surface temperature sharpening

In the case of CLMS ETa product, the LST sharpening needs to be performed by a ratio of around 3 (i.e. from 1 km to 300 m). Because of this relatively small ratio, the sharpening could potentially be achieved with methods which are more computationally efficient than DMS. Therefore, we also evaluated the performance of a classic and simple (and therefore faster) sharpening method called TsHARP (Agam et al., 2007). It relies on finding a linear regression between the Normalized Difference Vegetation Index (NDVI) and LST. The regression is derived on the whole LST scene to be sharpened with NDVI resampled to the LST spatial resolution. Afterwards, the linear regression is applied to NDVI at its native resolution to estimate the representation of LST at this resolution. Finally, a bias correction step is applied to ensure the consistency of LST between the original and sharpened maps.

Figure 16 Locations of areas of interest (left) and outline of the framework (right) used to evaluate DMS and TsHARP approaches for sharpening SLSTR LST images.

The two LST sharpening methods were tested in a number of geographically distributed areas of interest (AOI) covering different climatic zones and PFTs (Fig. 16 left panel) and across different seasons (images from at least four dates were sharpened at each AOI) to ensure a robust comparison. SY_2_SYN___ (SYN) Sentinel-3 product was used as a proxy for CLMS TOC reflectance product since the latter was still in production at the time this analysis was performed. The two products share the same spectral bands and spatial resolutions. The DMS regression models are trained on the entire Sentinel-3 SYN 3 min product data unit (PDU) (around 1400 km by 1200 km) as well as on subsets of 30 by 30 LST pixels in a moving window fashion. The evaluation of the sharpening methods was performed using the Sentinel-3 LST (as shown in Fig. 16 right panel) due to the lack of in situ LST data which could be used to validate satellite LST with 300 m spatial resolution. The SLSTR LST product was first resampled from 1 to 3 km and the SYN reflectance product was resampled from 300 m to 1 km, as was the ancillary data. The two sharpening methods were then used to recreate the LST at 1 km resolution, and the resulting map was compared with the original LST product. This evaluation framework assumes that there are no significant differences in relations between the explanatory variables and the LST when sharpening from 3 to 1 km and when sharpening from 1 to 300 m.

The results of the comparison are summarized in Table 12. At all sites the DMS method produces more accurate sharpened LST compared to TsHARP (see supplementary Table C2). The largest difference is at the Central Europe AOI where RMSE of DMS is up to 0.3 K lower than that of TsHARP and Mean Absolute Error (MAE) is up to 0.24 K lower. Looking at all sites, the RMSE of DMS is around 0.1 K lower compared to TsHARP (12 % difference) and MAE is around 0.06 K lower (10 % difference). Bias is minimal for both methods because bias correction is incorporated in both of them. Both methods also have similar and very high r², with DMS being slightly better, and the slope of the linear regression between sharpened and original LST very close to 1.

Table 12 Accuracy statistics for thermal sharpening: coefficient of determination (r²), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), bias (modelled minus observed) and slope of the linear regression line between modelled and observed values. Model configurations are described in text.

Regarding the three sets of DMS explanatory variables (see Sect. 2.2.4 for details), the differences between them are negligible. Looking at the details, there is no site in which “DMS – WaPOR” performs better than the other two configurations, while “DMS – Reflectance” has slightly better performance at some sites and “DMS – WaPOR selected” at others. However, the reduction in explanatory variables from 19 (“DMS – Reflectance”) to 9 (“DMS – WaPOR selected”) results in speed-up of the execution by a factor of around 2.3.

Qualitative assessment of the sharpened LST maps reveals that, in some cases, when variability of NDVI is low, TsHARP fails to find a meaningful relation between NDVI and LST (i.e. the linear regression has a slope close to 0). In such cases, the resulting sharpened LST is mainly an output of the bias correction step and has a blurry appearance consistent with simple resampling of lower resolution LST to higher resolution (see Fig. C9). In those cases, the quantitative analysis might still result in good accuracy statistics due to a small difference in spatial resolution between 3 and 1 km LST. DMS is less sensitive to, but not fully unaffected by, this issue because it relies on a range of spectral bands, indices, and DEM-based datasets. It is also noticeable that, in some cases, when there is insufficient information in the spectral data, DMS relies too heavily on DEM-based datasets, which can result in artifacts in the sharpened LST. One of the root causes of this issue could be the insufficient or incorrect atmospheric correction applied to SYN spectral bands. For operational production of CLMS ETa (and for production of data validated in Sect. 3) the SYN product is replaced by CLMS TOC reflectance product. The latter has shown improved agreement with in-situ measurements (Copernicus Land Monitoring Service, 2025c) and could therefore lead to an improvement in LST sharpening in such cases.

4.3.3 Weather forcing and ETa gap-filling

The suitability of the topographically corrected CAMS forecast data (Sect. 2.2.5) for modelling of actual evapotranspiration was assessed by comparing modelled daily solar irradiance and reference ET against measurements from 43 EC flux towers located in western Europe (France, Belgium, western Germany, Switzerland Spain, and northern Italy), United States and Australia (Fig. C10) for the year 2020. Those towers represent various topographical conditions from flat and low-lying to mountainous terrain as well as different climates, from arid to temperate. The statistical results of this comparison are shown in Fig. 17 and confirm the applicability of CAMS forecast for ETa modelling and suitability of the correction methods.

Figure 17 Density scatter plots between daily solar irradiance (Daily SW Irradiance panel in W m⁻²) and reference ET (FAO56 ET_r panel in mm d⁻¹) modelled with topographically corrected CAMS forecast data (x-axis) and measured in 43 representative Eddy Covariance (EC) stations located in western Europe, United States and Australia (y-axis).

Apart from forcing the ET models, CAMS data was used for gap-filling of modelled ETa maps through the use of reference ET and precipitation (PR) (Sect. 2.3). We evaluated three gap-filling approaches: ignoring rainfall (called GF_ETr – Guzinski et al., 2021); taking rainfall into account as described in Guzinski et al. (2023) (called GF_ETr+PR); conservative modification of GF_ETr+PR described in Sect. 2.3 (called GF_ETr+PR80). We also evaluated near real-time (NRT) gap-filling periods with different maximum durations: 15 d (as used by default in Sen-ET approach), 30 and 60 d (as used in WaPOR during reprocessing). Finally, a non-time-critical (NTC) gap-filling with a maximum duration of 60 d (i.e. up to 60 d before and after target date) was also evaluated.

To evaluate the behaviour of the methods in diverse climates, we performed an analysis using in-situ ET measured at geographically distributed ICOS stations (Table B5). The daily ET at ICOS stations was calculated by summing up good quality instantaneous ET at hourly or half-hourly timesteps. Dates on which more than two measurements were of poor quality were ignored and no correction for lack of energy balance closure was performed. The determination of quality of ET data was based on ICOS quality flags associated with each measurement. Reference ET and rainfall sums were calculated from the stations' meteorological measurements. Cloudy days were identified as those for which measured solar irradiance was less than 80 % of theoretical surface solar irradiance in clear-sky conditions.

Table 13 presents the results on this evaluation. Statistics are calculated on daily basis and only for data points which were gap-filled and for which the rainfall adjustment is relevant. Neglecting rainfall can introduce an overall negative bias (underestimation) of up to −0.18 mm d⁻¹ on wet, gap-filled dates. When rainfall is taken into account the introduced bias can become negligible and this reduction in bias leads to a less significant reduction in introduced uncertainty. Method GF_ETr+PR performs well with shorter window size but performance degrades significantly as window size increases and for largest windows the uncertainty is higher than for method GF_ETr. Method GF_ETr+PR80 results in lowest introduced bias and uncertainty with 60 d window size, chosen for the CLMS ETa processing chain. Finally, the NTC gap-filling performs better than NRT gap-filling and with method GF_ETr+PR80 performing the best in this case. An example timeseries of ET gap-filled with methods GF_ETr and GF_ETr+PR80 is shown in Fig. 18.

Table 13 Accuracy statistics for gap-filling of evapotranspiration on cloudy days using ICOS in-situ data and different methods (described in text) and maximum window sizes (in NRT mode apart from last rows). Numbers highlighted in bold represent the best obtained statistics for each window size.

Figure 18 An example of gap-filled ET timeseries at Voulundgaard ICOS station in 2015. Maximum gap-filling window of 60 d was used with, at left, gap-filling with rainfall not taken into account (method GF_ETr) and, at right, gap-filling with rainfall taken into account using the best performing method (method GF_ETr+PR80). Most notable differences are in spring (March–April) and late summer (July–August). Shaded area represents the uncertainty (energy balance closure error) of the data.

4.4 Limitations and potential improvements to the CLMS ETa product

The CLMS ETa product entered operational production at the end of 2025. Therefore, given the very short timeframe, it was not feasible to introduce changes for the initial version of the product. However, like most CLMS products, ETa might undergo evolution and reprocessing in the coming years and, therefore, it is worth to highlight some potential areas of improvement which might be the focus of this evolution.

One of the main limitations in producing medium- and high spatial resolution ET products (below 1 km) is the lack of LST observations with high spatio-temporal resolution, especially within the Copernicus Sentinel satellites. This situation should be resolved by the end of the decade when Land Surface Temperature Monitoring (LSTM) mission, with a primary objective of frequent monitoring of field-scale ETa, will join the Copernicus constellation (Koetz et al., 2018). In the meantime, it might be advantageous to make use of well established non-Copernicus sensors, such as VIIRS on board of Suomi-NPP satellite which can provide daily observations of LST with spatial resolution of 375 m.

Another potential improvement would be to compute a temporal running mean of biophysical parameters and TOC reflectance using the latest week of data. This should result in two benefits. Firstly, inverting the PROSPECT-D+4SAIL RTM (see Sect. 2.2.2) is an ill-posed problem, meaning that multiple combinations of biophysical parameters can result in the same reflectance. By taking a mean of multiple inversions, under the assumption that the biophysical conditions of vegetation do not change significantly over short time periods, the robustness of this inversion can be improved. Secondly, this will allow the full use of SLSTR swath with VZA below 35° by gap-filling the OLCI image on the eastern edge of the swath (see Sect. 2.2.1).

Producing a reanalysis dataset, e.g. two months after the completion of NRT production, might also lead to improvements in the accuracy of the CLMS ETa product. Firstly, because it would allow the use of ERA5 reanalysis meteorological data (Hersbach et al., 2020) instead of CAMS forecast as ET model forcing. The former has improved accuracy and increased spatial resolution compared to the latter and is therefore expected to result in more robust ET estimates. Secondly, gap-filling could be performed using 60 d both preceding and succeeding the target date which will result in less cloud gaps in the final ETa maps but also in improved accuracy (compare two bottom rows of Table 13).

Currently the Ensemble is created by simple averaging of the ET (and E and T) modelled by TSEB-PT and ETLook. However, the two models might not perform equally well, or be equally suitable, for all climates and PFTs (see Sect. 4.1). Therefore, following the approach being suggested for other ET model ensembles (Reitz et al., 2025) it might be advantageous to use the assessment of individual model performance and uncertainty in different conditions when creating the Ensemble dataset. Those two ET models were chosen based on preliminary analysis funded by European Commission which led in their explicit inclusion in the call for tenders for the CLMS ETa product (JRC/2023/OP/0273 – https://ec.europa.eu/info/funding-tenders/opportunities/portal/screen/opportunities/tender-details/13795). This however does not imply that they are the best performing models in all circumstances. Therefore, enlarging the ensemble member size with additional ET models, e.g. those included in the OpenET ensemble (Melton et al., 2022) or other models which were shown to be applicable globally such as ETMonitor (Zheng et al., 2022), can further improve the robustness of the CLMS ETa product in future evolutions.

Adding an explicit sub-product on evaporation from canopy interception (currently implicitly included in the evaporation sub-product) could be another development area. Modeling this component robustly and operationally is not a straightforward task as the currently existing semi-empirical methods might not account for all the complexities of this process. Therefore, a novel multi-source energy balance model might need to be developed.

The modelling of evaporation of inland water bodies can be important for many applications, such as water resources management or SDG reporting. Therefore, a post-processing step could be introduced in which this flux could be estimated, e.g. using Penman model (Monteith and Unsworth, 2013) parameterized for water.

Finally, the implemented ET modelling workflows present some limitations compared to the state-of-the-art in ET modelling research, particularly in some climates and PFTs as discussed in Sect. 4.1. There are multiple studies looking into different aspects of ET modelling including more detailed radiation partitioning (Ponce De León et al., 2026), accounting for multiple sources of heat fluxes (Mwangi et al., 2022), improved estimation of aerodynamic resistances (Li et al., 2019) and surface roughness (Alfieri et al., 2019) and integration of direct soil moisture estimates into ET models (Huang et al., 2025). In an operational and global product there is an inherent trade-off between keeping up to date with latest research and maintaining consistency and accuracy in a global context within reasonable computational costs. The current version of the CLMS ETa product is the result of such trade-off analysis, assessing the potential pros and cons of research findings for a global product, executed within tight budget and time constraints. Upcoming evolutions will provide opportunities to re-assess this, guided by user feedback and new scientific developments.

In addition to the potential input data and model improvements presented above, the validation methodology could also be extended. The focus of this study was on point validation using EC sites spread around the globe. Despite our best efforts, there are still large areas (especially in Asia, Africa and South America) from which no EC measurements are available. In such areas, other quality assessment methods could be performed (Tran et al., 2023), including catchment scale water balance (Blatchford et al., 2020; Weerasinghe et al., 2020; Comini De Andrade et al., 2024), and statistical analysis of consistency and uncertainty using other global or regional products (He et al., 2025).

5 Conclusions

A global and operational actual evapotranspiration (ETa) product was introduced to the portfolio of the Copernicus Land Monitoring Service (CLMS) at the end of 2025. This addition is motivated by a request from the Food and Agriculture Organization of the United Nations but the new product will have a multitude of applications in water resources management, SDG reporting, food security, forest management and other fields. The product is designed to have a 300 m spatial resolution, dekadal (10 d average) temporal resolution for water fluxes (evapotranspiration and its components of evaporation and transpiration) and daily temporal resolution for instantaneous heat fluxes (latent and sensible heat), global extent and to be produced in near-real-time with 2 d delay (Table 1). As all other CLMS products, it is distributed with a free and open license and will have guaranteed long-term continuity.

The product is based almost exclusively on Copernicus data, ranging from imagery acquired by OLCI and SLSTR sensors on board of Sentinel-3 satellites through meteorological forecast data provided by Copernicus Atmosphere Monitoring Service (CAMS) to higher-level products such as land cover maps produced by Copernicus Land Monitoring Service (Table 3). Those data undergo significant pre-processing to turn them into input forcing (Table 2) for two ET models: TSEB-PT and ETLook (Sect. 2). The cloud masked top-of-canopy (TOC) reflectance from OLCI and SLSTR shortwave optical bands were obtained from CLMS TOC product, while LST from SLSTR thermal bands was based on the Sentinel-3 LST product. A decision was made not to gap-fill the cloud gaps in those input products since the TOC reflectances were acquired at the same time and location as LST and it was preferred not to gap-fill the latter in order to obtain most accurate energy balance at the surface. In addition, a maximum view zenith angle of 45° was enforced to limit uncertainty due to different directionality effects. The derivation of biophysical parameters (e.g. LAI, leaf pigments, fraction of LAI that is green, canopy albedo) was based on inversion of radiative transfer models driven by the TOC product spectral bands. This allowed us to obtain LAI comparable with the CLMS LAI product while at the same time to derive fraction of vegetation that is green which better aligns with the expected phenological patterns, and leaf pigments required for canopy albedo estimation. Applying efficient vectorization and inversion techniques allows even such complex approaches to be applied in NRT at global scale and daily timestep. Meteorological forcing was based on the CAMS forecast dataset and after undergoing rather simple topographical correction showed to be highly suitable for modelling of ET at 300 m resolution. The 1 km SLSTR LST product was sharpened to 300 m using Data Mining Sharpener, a complex machine learning method, and 300 m TOC product based spectral indices, reflectance and a DEM with a specific model trained and applied to each pair of scenes. The resulting LST recreates well the finer spatial details of the higher resolution although can also produce artifacts or blurry appearance when there is insufficiently strong relationship between the shortwave and thermal observations. Finally, the outputs of the two ET models were gap-filled using reference ET and rainfall to obtain all-weather ET estimates, with the latter adding important information on soil wetting and leading to improved accuracy of gap-filled ET. The selected input data are available globally and the pre-processing and ET gap-filling methods do not undergo any localized fine-tuning or specific parameterization and therefore are applicable globally (Sect. 4.3) and suitable as input also for other ET models.

A demonstration dataset of three years was produced using the prototype end-to-end processing chain for the CLMS ETa product and validated using three years of measurements from 206 globally distributed eddy covariance stations (Sect. 3). Despite our best efforts for EC data collection, there are areas of the globe which lack in-situ measurements and in which product performance is therefore more uncertain.

The two ET models individually proved to be capable of producing global ET estimates, both achieving RMSE below 1 mm d⁻¹, correlation coefficient above 0.75, and absolute bias below 0.1 mm d⁻¹ when compared to the in-situ measurements. Those results compare favourably with the WaPOR global ETa dataset produced by FAO, which has the same spatio-temporal characteristics as the Copernicus ETa products, and other higher-resolution ETa datasets (Sect. 4.2). However, when looking at different climates and plant functional types (PFTs), the differences between the model become apparent (Sect. 4.1). In particular ETLook model shows higher positive bias in semi-arid and arid areas, most probably due to the limitations of the LST – vegetation fractional cover trapezoid soil moisture estimation method. It also achieves lower correlation in tropical climate, possibly due to the underlying constraints of Penman Monteith model in humid environments. On the other hand, the TSEB-PT model struggles more in forested PFTs due to its higher sensitivity to vegetation height and unsuitability of default parameterization of Priestley-Taylor potential ET in forests. We tried to address those issues by using a LiDAR based vegetation height map and empirical correction factors but especially the latter could benefit from a more thorough physics-based approach.

Taking an ensemble (average) of the two models has shown to be an effective method to improve the ETa product performance. The Ensemble product achieved RMSE of 0.87 mm d⁻¹, bias of 0.01 mm d⁻¹ and correlation of 0.81 and outperformed the individual models in all climates and large majority of PFTs (Sect. 3.2). This was mainly due to cancellation of random and systematic errors when averaging the model outputs. However, in some circumstances either of the individual models outperformed the ensemble and therefore both individual model outputs and the ensemble values are distributed as part of the CLMS ETa product. A more intelligent method of creating the ensemble, based on specific model performances in different climate-PFT zones, as well as addition of new model members to the ensemble is expected to further improve product performance (Sect. 4.4).

The results of this study are highly encouraging for the production of global ETa dataset based on Copernicus data sources. The addition of operational and NRT ETa product in the CLMS portfolio should significantly enlarge the CLMS user community and encourage further development of applications leveraging ET maps from regional to global scales.

Appendix A: Estimation of canopy transmittance and reflectance

The direct-hemispherical spectral transmittance (${\mathit{\tau }}_{\mathrm{C},\mathrm{DIR},\mathit{\lambda }}$) at a given solar zenith angle (θ_S) is calculated following the equations of Campbell and Norman (1998) for a single layer canopy:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\mathit{\tau }}_{\mathrm{C},\mathrm{DIR},\mathit{\lambda }}\left({\mathit{\theta }}_{\mathrm{S}}\right)=\\ & {\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }{\left({\mathit{\theta }}_{\mathrm{S}}\right)}^{\mathrm{2}}-\mathrm{1}\right)\exp \left(-\sqrt{{\mathit{\zeta }}_{\mathit{\lambda }}}{\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right)\mathrm{LAI}\right)}{\begin{array}{c}\left({\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }{\mathit{\rho }}_{\mathrm{S},\mathit{\lambda }}-\mathrm{1}\right)+{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left(\mathit{\psi }\right)\left({\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right)-{\mathit{\rho }}_{\mathrm{S},\mathit{\lambda }}\right)\\ \exp \left(-\mathrm{2}\sqrt{{\mathit{\zeta }}_{\mathit{\lambda }}{\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right)\mathrm{LAI}}\right)\end{array}}}\end{array}\end{array}$$

with λ being either the PAR or NIR. ${\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left(\mathit{\psi }\right)$ is the beam spectral reflection coefficient for a deep canopy with non-horizontal leaves (see Eq. A2), ζ_λ is the leaf absorptivity, κ_b is the extinction coefficient for direct-beam radiation (per LAI unit), and ρ_S,λ is the soil spectral reflectance. The multiple scattering between the soil and the canopy is accounted for in the ${\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }$ and ρ_S,λ terms.

$$\begin{array}{}\text{(A2)}& {\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right)={\displaystyle \frac{\mathrm{2}{\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right){\mathit{\rho }}_{\mathit{\lambda }}^H}{{\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right)+\mathrm{1}}}\end{array}$$

${\mathit{\rho }}_{\mathit{\lambda }}^H=\frac{\mathrm{1}-\sqrt{{\mathit{\zeta }}_{\mathit{\lambda }}}}{\mathrm{1}+\sqrt{{\mathit{\zeta }}_{\mathit{\lambda }}}}$ is the reflectance factor for a canopy with horizontal leaves.

Finally, the canopy beam extinction κ_b(ψ) is calculated based on the ellipsoidal LADF of Campbell (1990):

$$\begin{array}{}\text{(A3)}& {\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right)={\displaystyle \frac{\sqrt{{\mathit{\chi }}^{\mathrm{2}}+{\tan }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{S}}}}{\mathit{\chi }+\mathrm{1.774}{\left(\mathit{\chi }+\mathrm{1.182}\right)}^{-\mathrm{0.733}}}}\end{array}$$

Diffuse spectral transmittance (${\mathit{\tau }}_{\mathrm{C},\mathrm{DIF},\mathit{\lambda }}$) is calculated by numerically integrating κ_b over the hemisphere:

$$\begin{array}{}\text{(A4)}& {\mathit{\kappa }}_{\mathrm{d}}=\mathrm{2}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}{\mathit{\kappa }}_{\mathrm{b}}\left(\mathit{\psi }\right)\sin \mathit{\psi }\cos \mathit{\psi }\mathrm{d}\mathit{\psi }\end{array}$$

and replacing κ_b by κ_d in Eq. (A1).

Similarly the canopy direct spectral albedo is computed as:

$$\begin{array}{}\text{(A5)}& {\mathit{\rho }}_{\mathrm{C},\mathrm{DIR},\mathit{\lambda }}\left({\mathit{\theta }}_{\mathrm{S}}\right)={\displaystyle \frac{\begin{array}{c}{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right)+\left[\frac{{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right)-{\mathit{\rho }}_{\mathrm{S},\mathit{\lambda }}}{{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right){\mathit{\rho }}_{\mathrm{s},\mathit{\lambda }}-\mathrm{1}}\right]\\ \exp \left(-\mathrm{2}\sqrt{{\mathit{\zeta }}_{\mathit{\lambda }}}{\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right)\mathrm{LAI}\right)\end{array}}{\begin{array}{c}\mathrm{1}+{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right)+\left[\frac{{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right)-{\mathit{\rho }}_{\mathrm{s},\mathit{\lambda }}}{{\mathit{\rho }}_{\mathrm{C},\mathit{\lambda }}^{\ast }\left({\mathit{\theta }}_{\mathrm{S}}\right){\mathit{\rho }}_{\mathrm{S},\mathit{\lambda }}-\mathrm{1}}\right]\\ \exp \left(-\mathrm{2}\sqrt{{\mathit{\zeta }}_{\mathit{\lambda }}}{\mathit{\kappa }}_{\mathrm{b}}\left({\mathit{\theta }}_{\mathrm{S}}\right)\mathrm{LAI}\right)\end{array}}}\end{array}$$

and the diffuse canopy albedo (${\mathit{\rho }}_{\mathrm{C},\mathrm{DIF},\mathit{\lambda }}$) by replacing κ_b(θ_S) by κ_d.

Appendix B: List of flux towers

Tables B1, B2, B3 and B4 list all the EC flux tower sites used to validate CLMS ETa datasets. The ICOS datasets used in the analysis were the ETC L2 FLUXNET (ICOS RI et al., 2025) and the Warm Winter 2020 (ICOSww – Warm Winter 2020 Team et al., 2022). AsiaFlux site is described in Meijide et al. (2017). The DOIs of AmeriFlux sites are listed below:

• CA-DBB: https://doi.org/10.17190/AMF/1881565 (Christen and Knox, 2025)

• US-ARM: https://doi.org/10.17190/AMF/1854366 (Biraus et al., 2026)

• US-Bi1: https://doi.org/10.17190/AMF/1871134 (Rey-Sanchez et al., 2025)

• US-Bi2: https://doi.org/10.17190/AMF/1871135 (Rey-Sanchez et al., 2026)

• US-BZF: https://doi.org/10.17190/AMF/1881570 (Euskirchen, 2026)

• US-HB3: https://doi.org/10.17190/AMF/2229378 (Forsythe et al., 2026)

• US-Jo1: https://doi.org/10.17190/AMF/1902833 (Tweedie, 2026)

• US-Me6: https://doi.org/10.17190/AMF/2204871 (Hanson et al., 2026)

• US-Mo1: https://doi.org/10.17190/AMF/2229382 (Schreiner-McGraw, 2026a)

• US-Mo2: https://doi.org/10.17190/AMF/2229383 (Schreiner-McGraw, 2026b)

• US-Mo3: https://doi.org/10.17190/AMF/2229384 (Schreiner-McGraw, 2026c)

• US-ONA: https://doi.org/10.17190/AMF/1832163 (Silveira and Bracho, 2025)

• US-Rls: https://doi.org/10.17190/AMF/2229387 (Flerchinger, 2026)

• US-Ro4: https://doi.org/10.17190/AMF/1881589 (Baker and Griffis, 2026)

• US-Ro5: https://doi.org/10.17190/AMF/1818371 (Baker and Griffis, 2025a)

• US-Ro6: https://doi.org/10.17190/AMF/1881590 (Baker and Griffis, 2025b)

• US-Sne: https://doi.org/10.17190/AMF/1871144 (Shortt et al., 2026)

• US-SRG: https://doi.org/10.17190/AMF/2204877 (Scott, 2026a)

• US-Ton: https://doi.org/10.17190/AMF/2204880 (Ma et al., 2026a)

• US-Tw4: https://doi.org/10.17190/AMF/2204881 (Eichelmann et al., 2026)

• US-Var: https://doi.org/10.17190/AMF/1993904 (Ma et al., 2026b)

• US-Whs: https://doi.org/10.17190/AMF/1984574 (Scott, 2026b)

• US-Wkg: https://doi.org/10.17190/AMF/1984575 (Scott, 2026c)

• US-xAE: https://doi.org/10.17190/AMF/1985434 (NEON, 2025a)

• US-xDC: https://doi.org/10.17190/AMF/1985437 (NEON, 2025b)

• US-xDS: https://doi.org/10.17190/AMF/1985439 (NEON, 2025c)

• US-xJE: https://doi.org/10.17190/AMF/1985443 (NEON, 2025d)

• US-xSB: https://doi.org/10.17190/AMF/1985451 (NEON, 2025e)

• US-xSE: https://doi.org/10.17190/AMF/1985452 (NEON, 2026a)

• US-xST: https://doi.org/10.17190/AMF/1985454 (NEON, 2026b)

• US-xTA: https://doi.org/10.17190/AMF/1985455 (NEON, 2026c)

• US-xWD: https://doi.org/10.17190/AMF/2229412 (NEON, 2026d)

Table B5 lists selected sites from ICOS WarmWinter2020 database which were used to evaluate biophysical modelling and ET gap-filling approaches.

Table C1 and Fig. C1 show validation of modelled dekadal ETa (similarly to Table 8 and Fig. 7) but only against the flux towers for which ECB correction was applied.

Table B1 Geographical location, climate region (Köppen classification), plant functional type (PFT) and network/dataset of origin for EC flux towers used for CLMS ETa model validation – part 1/4.

Table B2 Geographical location, climate region (Köppen classification), plant functional type (PFT) and network/dataset of origin for EC flux towers used for CLMS ETa model validation – part 2/4.

Table B3 Geographical location, climate region (Köppen classification), plant functional type (PFT) and network/dataset of origin for EC flux towers used for CLMS ETa model validation – part 3/4.

Table B4 Geographical location, climate region (Köppen classification), plant functional type (PFT) and network/dataset of origin for EC flux towers used for CLMS ETa model validation – part 4/4.

Table B5 Selected sites from ICOS WarmWinter2020 database (Warm Winter 2020 Team et al., 2022) used to evaluate biophysical modelling and ET gap-filling approaches. In addition to geographical location the table shows climate region (Köppen classification) and plant functional type (PFT) of each site.

Appendix C: Supplementary plots, maps and tables

In this appendix we present supplementary figures and tables supporting the analysis conducted in the Sect. 3.2 (Validation results) and Sect. 4 (Discussion). Reader is referred to those section for more details.

Figure C1 Taylor-plot overview of dataset accuracies using only those sites in which energy balance closure (EBC) correction was applied (AmeriFlux, ICOS, OzFlux networks). Azimuthal angle represents the Pearson r, radial distance from the origin is the standard deviation (mm d⁻¹). Circle marker at the horizontal axis represents perfect model performance, markers closer to this point indicate better model performance. Radial distance from this point is equivalent to centered RMSE (contour values shown in mm d⁻¹). The statistics shown in the Taylor diagram were generated from the full dataset of observed and modelled ETa pairs across all EBC corrected validation sites.

Table C1 Bias (modelled − measured: mm d⁻¹), Root Mean Squared Error (RMSE: mm d⁻¹), and Pearson correlation coefficient (r) scores for dekadal ETa per model (all sites combined), and summary statistics at site level for sites in which energy balance closure correction was applied (AmeriFlux, ICOS, OzFlux networks). rBias and rRMSE are relative metrics (i.e., divided by mean measured ET).

Figure C2 Taylor plots on the performance of the dekadal ET from TSEB-PT, ETLook and Ensemble datasets. Sites grouped per climate region (Köppen climate classification system). For number of sites in each group see Fig. 6. Azimuthal angle represents the Pearson r, radial distance from the origin is the standard deviation (mm d⁻¹). Circle marker at the horizontal axis represents perfect model performance, markers closer to this point indicate better model performance. Radial distance from this point is equivalent to centered RMSE (contour values shown in mm d⁻¹).

Figure C3 Taylor plots on the performance of the dekadal ET from TSEB-PT, ETLook and Ensemble datasets. Sites grouped per Plant Functional Type. For number of sites in each group see Fig. 6. Azimuthal angle represents the Pearson r, radial distance from the origin is the standard deviation (mm d⁻¹). Circle marker at the horizontal axis represents perfect model performance, markers closer to this point indicate better model performance. Radial distance from this point is equivalent to centered RMSE (contour values shown in mm d⁻¹).

Figure C4 Timeseries plots of LE (top panels) and H (bottom panels) (W m⁻²) as measured at the eddy covariance stations (blue) and as modelled by TSEB-PT (orange) at satellite overpass time in AU-Cum (EBF), CL-SDF (EBF), US-Ton (WSA), US-Var (GRA), Ragola (GRA) and IT-Cp2 (EBF) sites.

Figure C5 Absolute differences (mm d⁻¹) between estimates of ET, E and T of the ETLook and TSEB-pT models on the first dekad of January and July 2021.

Figure C6 Global maps of the Ensemble, ETLook and TSEB-PT ET estimates for norther hemisphere winter (left panels – first dekad of January 2021) and northern hemisphere summer (right panels – first dekad of July 2021).

Figure C7 Histograms of TSEB-PT (blue) and ETLook (green) dekadal ET in first dekads of January (left panels) and of July (right panels) 2021 across 10° lat/lon tiles located in different climate regions according to the Köppen classification system. Climate types per region: A – Dfc; B – Dfa; C – mosaic of Csa, BSk, Csb, Cfb; D – BWh; E – Cfb; F – Af; G – mosaic of BWh, BSk, Cwb, Cfb; H – Af.

Figure C8 Differences (WaPOR − CLMS: mm d⁻¹) between estimates of ET, E and T of the two products on the first dekad of January and July 2021.

Figure C9 Example from Amazon rainforest (forested and cleared areas) of original 1 km resolution Sentinel-3 LST and Sentinel-1 LST resampled to 3 km and then sharpened back to 1 km using TsHARP and DMS (see Sect. 4.3.2 for details).

Figure C10 Location of sites used to evaluate CAMS meteorological data (Sect. 4.3.3). Background map shows Köppen climate classification.

Table C2 Accuracy statistics for thermal sharpening: coefficient of determination (r²), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), bias (modelled minus observed) and slope of the linear regression line between modelled and observed values. Model configurations are described in Sect. 2.2.4.