P-LSHv2: a multi-decadal global daily evapotranspiration dataset enhanced with explicit soil moisture constraints

Accurately quantifying the impact of soil water availability on evapotranspiration (ET) is crucial for improving ET retrieval accuracy. However, most global satellite-derived ET datasets do not explicitly incorporate soil moisture constraints, leading to significant uncertainties, particularly in water-limited regions. In this study, we propose an enhanced soil moisture constraint scheme that effectively captures soil moisture's influence on vegetation transpiration and soil evaporation using a quantile-based approach. Unlike previous methods, this scheme relies solely on soil moisture data, reducing uncertainties associated with heterogeneous soil hydraulic properties. We integrated this approach into the process-based land surface ET/heat fluxes algorithm (P-LSH, or P-LSHv1), developing an improved version, P-LSHv2. Using observations from 106 global flux towers, we calibrated biome- and climate-specific parameters and quantified moisture constraints across diverse climates and land cover types. P-LSHv2 achieves notable improvements in ET estimation, with a reduced Root Mean Square Error (RMSE) of 0.67 mm d⁻¹ and an increased Pearson correlation coefficient (R) of 0.81, indicating strong agreement with flux tower observations. As a result of these improvements, P-LSHv2 outperforming its predecessor, P-LSHv1, particularly in arid regions. Comparative analyses show that P-LSHv2 surpasses the Penman-Monteith-Leuning model and the Global Land Evaporation Amsterdam Model in capturing soil moisture anomalies' effects on ET, enhancing global ET accuracy. Employing the P-LSHv2 algorithm, we have produced a long-term global daily ET dataset spanning 1982–2023, providing a valuable resource for research on terrestrial water and energy cycles and climate change. The dataset is freely available at https://doi.org/10.11888/Terre.tpdc.301969 (Feng and Zhang, 2025).

1 Introduction

Evapotranspiration (ET) is considered the second largest component of the terrestrial water cycle after precipitation and also plays a critical role in atmospheric-terrestrial carbon and energy exchanges (Jung et al., 2010; Wang and Dickinson, 2012; Zhang et al., 2016). Therefore, accurate estimation of global ET is essential for improving our understanding of the surface energy budget, global water cycle, and climate change (Fisher et al., 2017; Oki and Kanae, 2006).

Among various ET estimation methods, remote sensing retrieval provides an effective means to estimate large-scale ET for its spatially continuous and temporally frequent measurements of surface biophysical variables (Fisher et al., 2008; Longo et al., 2019). Several remote sensing-based global ET datasets have been generated across different temporal spans, spatial scales and resolutions (Jung et al., 2011; Martens et al., 2017; Mu et al., 2011; Zhang et al., 2010, 2019). Microwave-based measurements of hydrology states inside the canopy or soil generally contribute to high temporal resolution (i.e., daily) of ET estimation but typically have coarse spatial resolution (Martens et al., 2017; Wang et al., 2019), which is unsuitable for evaluating crop water requirements over complex agricultural landscapes. Conversely, the Moderate Resolution Imaging Spectroradiometer (MODIS) represents a comprehensive (500 m, 8 d) multi-decadal global operational satellite ET record (Mu et al., 2011; Zhang et al., 2019), which has been used to evaluate more recent interannual variability and trends (Hall et al., 2023; Román et al., 2024). However, MODIS data do not cover the pre-2000 period and are of insufficient length to represent long-term interannual variability and trends, and attribution analysis of ET. Therefore, before MODIS-era, how to extend high-quality ET series in the case of Advanced Very High-Resolution Radiometer (AVHRR), and how to balance fine spatial and high temporal resolutions, are essential for extending the coverage of ET datasets.

An additional challenge lies in the reliability and robustness of ET datasets, particularly in complex global ecosystems and diverse climates. Although some global ET datasets have been developed, stark discrepancies are identified among various datasets in terms of magnitude and changing trends (Badgley et al., 2015; Li et al., 2022; Pan et al., 2020; Zhu et al., 2022). For instance, an investigation from the Global Soil Moisture Project-2 (GSWP-2) found that the global average annual ET from 15 models varied widely, ranging from 272 to 441 mm yr⁻¹, with the maximum value 1.5 times the minimum (Dirmeyer et al., 2006). Similarly, Jiménez et al. (2011) showed considerable disagreement among ET estimates from 12 global products, particularly in tropical rainforests. McCabe et al. (2016) highlighted that four widely-used remote sensing ET models exhibit significant variability across specific biomes and climate zones, with no model consistently outperforming the others. The Water Cycle Multi-Mission Observation Strategy-Evapotranspiration Project (WACMOS-ET) demonstrated low correlation and disagreement between various products in arid regions (Michel et al., 2016), and the MODIS product tends to underestimate tropical and subtropical ET fluxes, with obvious disagreement among various products in water-limited situations (Miralles et al., 2016). Chen et al. (2014) further found that the annual average ET across eight products in China ranged from 535 to 852 mm yr⁻¹, with higher uncertainties in tropical and humid regions. In summary, despite the numerous remote sensing ET models and datasets currently available, each has its own merits and disadvantages, and no consensus on which method is the best has been reached.

In ET estimation, soil moisture plays an essential role in water transport within vegetation. Due to the water potential gradient between leaf and air, water is transported from soil to vegetation roots, and leaves, and then dissipated into the atmosphere through stomata. In addition to this physical gradient, plant transpiration is fundamentally driven by biological needs such as nutrient transport, turgor maintenance, and leaf cooling. These physiological processes are tightly regulated and exert feedback control on stomatal conductance, thereby influencing transpiration dynamics. Importantly, root tips are capable of sensing subtle gradients in water potential (as low as 0.5 MPa mm⁻¹), exhibiting behaviors such as hydrotropism to actively seek water in heterogeneous soil profiles (Dexter, 1986; Gregory, 2006). Such sensory responses provide a physiological basis for root foraging behavior, which is particularly important for sustaining transpiration under drought conditions. Therefore, soil water content not only serves as the direct water pool for vegetation but also provides the essential hydraulic foundation that supports plant physiological functions. It ultimately regulates both the capacity and the demand side of transpiration through its dual role in water supply and physiological control. In previous remote sensing ET algorithms, the vapor pressure deficit (VPD) was generally adopted to reflect the associated moisture constraint on canopy conductance (Mu et al., 2011; Zhang et al., 2009, 2010, 2019). However, several studies have shown that water availability is the primary limiting factor for ET in arid areas (Yang et al., 2019; Zhang et al., 2015), namely the so-called water-limited regions, where soil moisture constraints on canopy conductance increase with climate dryness and drought severity (Han et al., 2020; Novick et al., 2016). In arid and semi-arid regions, even under high atmospheric demand (i.e., high VPD), the actual ET is often constrained by soil water availability and root uptake capacity. As the root–soil interface becomes hydraulically disconnected under drought, transpiration may decline despite strong evaporative demand. Although vapor pressure deficit and soil moisture are generally connected through land-atmosphere interactions (Liu et al., 2020a; Zhou et al., 2019), their anomalies may be temporally lagged (Zhang et al., 2022b), indicating a partial decoupling. Therefore, relying solely on VPD can misrepresent transpiration dynamics, especially under conditions of prolonged soil dryness or climate extremes. Incorporating explicit soil moisture constraints is essential for improving ET estimation at finer temporal scales (Fu et al., 2022a), and this need is expected to intensify under global warming scenarios (Liu et al., 2020b), where soil–atmosphere coupling may become even more unstable.

One of the main challenges of incorporating soil moisture constraints into remote sensing ET algorithms is the development of a robust and globally applicable moisture constraint scheme for quantifying water availability. Although soil moisture constraints have been applied to estimate ET at local (Xu et al., 2018; Zhu et al., 2014), regional, and global scales (Brust et al., 2021; Purdy et al., 2018), long-term available datasets are limited by the lack of universal constraint schemes and soil moisture datasets. Furthermore, existing soil moisture constraint schemes generally rely on soil hydraulic properties as thresholds (Martens et al., 2017; Zheng et al., 2022). However, global soil data present significant uncertainties, which substantially affect ET estimation (Dennis and Berbery, 2021). On the one hand, gridded datasets derived from field investigation are prone to great uncertainties, particularly in high-altitude and other harsh regions (Dai et al., 2019). On the other hand, the pedotransfer functions for generating gridded soil properties still lack proper extrapolation methods from local to global scales (Ma et al., 2021; Van Looy et al., 2017).

In response to these challenges, we integrated a new soil moisture constraint scheme into the previously developed Process-based Land Surface evapotranspiration/Heat fluxes algorithm (P-LSH, or P-LSHv1) (Zhang et al., 2010, 2015), creating the second version (P-LSHv2). The objectives of this study are to (1) develop a globally universal soil moisture constraint scheme, calibrate key parameters, and quantify moisture levels across diverse vegetation types and climate zones, (2) upgrade ET algorithm by incorporating soil moisture constraints and evaluate its robustness from sites to basins, and (3) generate a long-term, reliable global daily ET dataset spanning from 1982 to 2023.

2 Methodology

2.1 Baseline description of the P-LSH algorithm

The P-LSH algorithm separately estimates three components of the ground surface evapotranspiration: vegetation transpiration (λE_Canopy: W m⁻²), soil evaporation (λE_Soil: W m⁻²), and open water evaporation (λE_Water: W m⁻²). These components are calculated based on the pixel's land cover type. In vegetation pixels, the ET is partitioned into vegetation transpiration and soil evaporation by partitioning available energy according to the fractional vegetation cover. In barren and open water pixels, ET consists solely of soil evaporation and open water evaporation, respectively.

The Penman-Monteith equation is used to calculate vegetation transpiration:

$$\begin{array}{}\text{(1)}& \mathit{\lambda }{E}_{\mathrm{Canopy}}={\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{c}}+\mathit{\rho }{C}_p\mathrm{VPD}{g}_{\mathrm{ac}}}{\mathrm{\Delta }+\mathit{\gamma }(\mathrm{1}+g_{\mathrm{ac}}/g_{\mathrm{c}})}},\end{array}$$

where λ (J kg⁻¹) is the latent heat of vaporization, Δ (Pa K⁻¹) is the slope of the saturated water vapor pressure curve as a function of temperature; A_c (W m⁻²) is the available energy allocated to the canopy; ρ (kg m⁻³) is air density; C_p (J kg⁻¹ K⁻¹) is the specific heat capacity of air; VPD (Pa) is vapor pressure deficit; g_ac (m s⁻¹) is the aerodynamic conductance of the canopy; γ (–) is the psychrometric constant; g_c (m s⁻¹) is the canopy conductance calculated using the Jarvis-Stewart-type model:

$$\begin{array}{}\text{(2)}& {\displaystyle }{g}_{\mathrm{c}}=g_{\mathrm{0}}\left(\mathrm{NDVI}\right)\times m\left(T_{\mathrm{day}}\right)\times m\left(\mathrm{VPD}\right)\times m\left(C_{{\mathrm{CO}}_{\mathrm{2}}}\right),\text{(3)}& {\displaystyle }\begin{array}{rl}{g}_{\mathrm{0}}\left(\mathrm{NDVI}\right)& =\mathrm{1}/\left[b_{\mathrm{1}}+b_{\mathrm{2}}\times \exp \left(-b_{\mathrm{3}}\times \mathrm{NDVI}\right)\right]\\ & -\mathrm{1}/\left(b_{\mathrm{1}}+b_{\mathrm{2}}\right),\end{array}\end{array}$$

where g₀ (m s⁻¹) is the maximum value of g_c based on NDVI; m(T_day), m(VPD), and $m\left(C_{{\mathrm{CO}}_{\mathrm{2}}}\right)$ are stress factors associated with daylight temperature T_day (°C), VPD (Pa), and CO₂ concentration (ppm), respectively; b₁, b₂, and b₃ are biome-and-climate-specific parameters derived from flux towers measurements. In addition, T_opt (°) and β (°) in m(T_day) calculation represent optimal air temperature for photosynthesis and empirical parameters, respectively, which are taken as the biome-and-climate-specific parameters in this study. More details of stress factors are available at Zhang et al. (2010) and Feng et al. (2022).

The soil evaporation is reduced from its potential value using the moisture constraint:

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\lambda }{E}_{\mathrm{Soil}}=f\mathit{\lambda }{E}_{\mathrm{POT}},\text{(5)}& {\displaystyle }f={\mathrm{RH}}^{\left(\frac{\mathrm{VPD}}{k}\right)},\text{(6)}& {\displaystyle }\mathit{\lambda }{E}_{\mathrm{POT}}={\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{s}}+\mathit{\rho }{C}_p\mathrm{VPD}{g}_{\mathrm{as}}}{\mathrm{\Delta }+\mathit{\gamma }\times g_{\mathrm{as}}/g_{\mathrm{totc}}}},\end{array}$$

where f (–) is the moisture constraint estimated by the complementary relationship hypothesis here; RH (–) is the relative humidity; k (Pa) is a parameter to fit the complementary relationship and reflect the sensitivity of VPD, we regard it as a parameter for different land cover types and climates. The λE_POT (W m⁻²) is the potential soil evaporation; A_s (W m⁻²) is the available energy component allocated to the soil; g_as (m s⁻¹) is the aerodynamic conductance of the soil surface and is defined as the sum of the conductance associated with convective heat transfer (g_ch: m s⁻¹) and the conductance to radiative heat transfer (g_rh: m s⁻¹). The g_ch term is expressed in its resistance form r_c, a biome-and-climate-specific parameter, while the g_rh term is calculated by daylight temperature. The g_totc (m s⁻¹) is the corrected value of total aerodynamic conductance g_tot (m s⁻¹) and the correction is based on actual air temperature and pressure. In this study, the g_tot term is also expressed in the form of total aerodynamic resistance r_tot, which is regarded as a parameter determined by land cover types and climates.

The open water evaporation is calculated using a modified Penman equation that considers the effects of surface wind speed on aerodynamic conductance:

$$\begin{array}{}\text{(7)}& \mathit{\lambda }{E}_{\mathrm{Water}}={\displaystyle \frac{\mathrm{\Delta }A+\mathit{\rho }{C}_p\mathrm{VPD}{g}_{\mathrm{aw}}}{\mathrm{\Delta }+\mathit{\gamma }}},\end{array}$$

where A (W m⁻²) is the available energy for open water; g_aw (m s⁻¹) is the aerodynamic conductance of the open water estimated by wind speed:

$$\begin{array}{}\text{(8)}& g_{\mathrm{aw}}={\displaystyle \frac{\mathrm{1}+\mathrm{0.536}{U}_{\mathrm{2}}}{\mathrm{4.72}{\left[\ln (z_{\mathrm{m}}/z_{\mathrm{0}})\right]}^{\mathrm{2}}}},\end{array}$$

where U₂ (m s⁻¹) is the wind speed at 2 m height converted from measurement height z_m (m) using vertical wind speed function; z₀ (m) is the aerodynamic roughness of the water surface and set to 0.00137. Further details of the P-LSH algorithm are available in Zhang et al. (2010) and Feng et al. (2022).

2.2 Improvements on the P-LSH algorithm

In this study, a series of improvements are implemented to evolve P-LSHv1 into the P-LSHv2 algorithm, incorporating soil moisture constraints on both vegetation transpiration and soil evaporation, as well as a Bayesian and Sobol' uncertainty analysis framework. These improvements allow for a more explicit quantification of water supply limitations on ET and offer a structured approach to parameter uncertainty across diverse land covers and climates.

2.2.1 Soil moisture constraint scheme for vegetation transpiration

Previous studies (Liu et al., 2020b) have revealed that atmospheric moisture and soil moisture are generally decoupled in arid areas, which introduces significant uncertainty in quantifying soil water supply on ET. In this study, vegetation transpiration response to soil moisture is represented by the typical Jarvis-Stewart model (Jarvis, 1976; Stewart, 1988), in which a revised soil moisture constraint scheme is applied to constrain canopy conductance:

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{g}_{\mathrm{c}}& =g_{\mathrm{0}}\left(\mathrm{NDVI}\right)\times m\left(T_{\mathrm{day}}\right)\times m\left(\mathrm{VPD}\right)\\ & \times m\left(C_{{\mathrm{CO}}_{\mathrm{2}}}\right)\times m\left(\mathrm{SM}\right),\end{array}\text{(10)}& {\displaystyle }m\left(\mathrm{SM}\right)=\left\{\begin{array}{ll}\frac{\mathrm{SM}-{\mathrm{SM}}_{\min }}{{\mathrm{SM}}_{\mathrm{c}}-{\mathrm{SM}}_{\min }}& \mathrm{SM}\le {\mathrm{SM}}_{\mathrm{c}}\\ \mathrm{1}& \mathrm{SM}>{\mathrm{SM}}_{\mathrm{c}}\end{array}\right.,\text{(11)}& {\displaystyle }{\mathrm{SM}}_{\mathrm{c}}=F(\mathrm{SM},n),\end{array}$$

where m(SM) is the soil moisture stress factor, calculated by actual surface soil moisture (SM: m³ m⁻³), its critical value SM_c (m³ m⁻³), and its minimum value SM_min (m³ m⁻³) over a given period. Such a linear equation ensures the record of soil moisture dynamic characteristics as well as the simple normalization structure.

Different from other studies, the SM_c is calculated using the SM sequence, the quantile (n), and the quantile function (F) in this study, shown in Eq. (11). The SM_min and SM_c derived from quantiles and actual sequence are used as the threshold, which circumvents the use of soil hydraulic parameters that exhibit significant uncertainty. The parameter n represents the proportion of the period during which ET is subject to moisture constraints. Specifically, 100 indicates that ET is constantly constrained by soil moisture, while 0 indicates that ET is never constrained by soil moisture. Considering the potential impacts of different vegetation types and climates on constraints, we regard n as a biome-and-climate-specific parameter, which is first calibrated by the tower measurements, and subsequently applied to global ET estimation. That means, the quantile n takes the same value for specific land cover types and climates, but when mapped to various SM sequences in global grid estimation, the corresponding SM_c may vary. In addition to quantiles, this scheme also relies on cumulative distribution function of the SM sequence. The actual SM sequences used in this scheme are sufficient enough so that several heavy rainfall will not lead to significant changes in distribution function unless the local climate changes significantly. Even so, the variations in SM_c resulting from changes in the distribution function under climate change also reflect how our scheme responds to moisture constraints in a changing climate. This scheme effectively quantifies the moisture constraint level of the land cover types and climates, while also mitigating the uncertainty associated with soil hydraulic parameterization. The soil moisture stress function is illustrated in Fig. 1.

Figure 1 Example illustration of (a) SM_c estimation and (b) soil moisture stress function estimation.

2.2.2 Soil moisture constraint scheme for soil evaporation

For soil evaporation, the original algorithm uses atmospheric humidity as a proxy for soil moisture constraint, which introduces enormous uncertainty in estimating soil evaporation. In this study, a direct soil moisture constraint scheme following Feng et al. (2023) is implemented to replace the moisture constraint in Eq. (4) for soil evaporation:

$$\begin{array}{}\text{(12)}& f={\displaystyle \frac{\mathrm{SM}-{\mathrm{SM}}_{\min }}{{\mathrm{SM}}_{\max }-{\mathrm{SM}}_{\min }}},\end{array}$$

where SM_max (m³ m⁻³) and SM_min (m³ m⁻³) terms are the respective maximum and minimum SM sequences during a given period. Equation (12) is a specific version of Eq. (10) where the quantile n reaches its maximum of 100 and SM_c reaches SM_max. Here f is a normalized soil moisture index (or saturation fraction), assuming that the recording period of soil moisture observations can entirely cover the full range of soil moisture change at a given tower or pixel, which is suitably beneficial for long-term decades of ET estimation.

2.3 Parameters uncertainty analysis framework for the P-LSHv2 algorithm

With the structural upgrades to the ET algorithm, it is necessary to perform a parameter uncertainty analysis and optimization in the P-LSHv2 algorithm. The Bayesian and Sobol' uncertainty analysis framework is employed to quantify parameter characteristics.

Sobol' sensitivity analysis is based on variance decomposition, which can quantify both the single impacts of parameters on the model outputs and the interaction impacts among multiple parameters. The variance of model output is based on variance decomposition:

$$\begin{array}{}\text{(13)}& D\left(Y\right)={\sum }_{i=\mathrm{1}}^d{D}_i+{\sum }_{i<j}^d{D}_{ij}+\mathrm{\dots }+D_{\mathrm{12}\mathrm{\dots }d},\end{array}$$

where D(Y) is the variance of model output Y; d is the total number of parameters; D_i represents the partial variance for the first-order sensitivity of the ith parameter to the model output; D_ij denotes the partial variance for the second-order sensitivity corresponding to the interaction between the ith and jth parameters; D_12…d indicates the partial variance for the dth-order sensitivity involving all parameter interactions.

The first-order index and the total-order index can reflect the impacts of single parameters and multiple parameters' interaction on the model output:

$$\begin{array}{}\text{(14)}& {\displaystyle }{S}_i={\displaystyle \frac{D_i}{D}},\text{(15)}& {\displaystyle }{S}_{\mathrm{T}i}=\mathrm{1}-{\displaystyle \frac{D_{\sim i}}{D}},\end{array}$$

where S_i is the first-order index and represents the sensitivity from the effect of parameter X_i; S_Ti is the total-order index and represents the sensitivity from the effect of parameter X_i and its interactions with the remaining parameters; D_∼i is the amount of variance attributed to all other parameters after removing X_i.

After identifying sensitive parameters of the P-LSHv2 algorithm, it is necessary to quantify the uncertainty of the sensitive parameters. In recent years, the Bayesian theory has been widely used in parameter uncertainty analysis. Based on subjective prior knowledge and samples, the prior distribution of parameters can be updated to the posterior distribution, which significantly reduces the uncertainty of prior ranges:

$$\begin{array}{}\text{(16)}& \mathit{\pi }\left(X\mathrm{|}x\right)\propto f\left(\left.x\right|X\right)\mathit{\pi }\left(X\right),\end{array}$$

where x is the measurement sample; X is a group of parameter set; π(X) denotes the prior distribution, which is assumed to follow a uniform distribution within specified bounds; f(x|X) is the likelihood function.

Among various Bayesian-based methods, the Differential Evolution Markov Chain (DE-MC) method stands out due to its effective combination of the Differential Evolution algorithm and the Markov Chain Monte Carlo method. The DE-MC algorithm performs well in terms of convergence speed and computational efficiency. In the DE-MC algorithm, candidate proposals are initially generated using two random chains, and their difference is scaled and added to the current chain:

$$\begin{array}{}\text{(17)}& X_{\mathrm{p}}=X_i+\mathit{\gamma }(X_{\mathrm{R}\mathrm{1}}-X_{\mathrm{R}\mathrm{2}})+\mathit{\epsilon },\end{array}$$

where X_p is the proposed parameter set; X_i is the current parameter set; X_R1 and X_R2 are two random chains excluding X_i; γ is the scaling factor; ε is a term that reflects probability acceptance rules.

2.4 Parameters setting and multi-scale evaluation of the P-LSHv2 algorithm

Compared to P-LSHv1, the P-LSHv2 algorithm has been upgraded both in ET algorithm structure and biome-specific parameters. Considering the various responses of vegetation physiological processes, particularly with respect to soil moisture constraints under different climates, the P-LSHv2 algorithm takes distinct parameters for dry and wet zones within each land cover type. While the algorithm structure remains consistent, specific parameters are adjusted based on diverse zones.

To evaluate the P-LSHv2 algorithm and demonstrate its advances relative to the P-LSHv1 algorithm, a multi-scale evaluation procedure is developed. First, we quantify parameter uncertainty and calibrate optimal values based on global 106 flux towers containing daily ET measurements (Sect. 4.1). The performance of P-LSHv2 algorithm is then evaluated through both direct validation and leave-one-out cross-validation methods at 106 flux towers (Sect. 4.2). Second, we evaluate the P-LSHv2 algorithm at basins across the conterminous United States (CONUS) and other global regions, using the reconstructed values from the water balance method as the benchmark. The evaluation of basin scales is conducted on monthly, annual, and multiyear levels (Sect. 4.3). Finally, we generate global ET datasets using the P-LSHv2 algorithm and then compare them with those from P-LSHv1 algorithm, analyzing their diversity in spatial patterns and temporal trends (Sect. 4.4).

3 Data and materials

3.1 Eddy Covariance Flux Towers

In this study, climatic and flux measurements from eddy covariance towers that represent nine different land cover types, as well as their associated dry or wet climate zones, are used to calibrate and verify algorithm performance. These land cover types include evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), mixed forest (MF), open shrublands (OSH), woody savannas (WSA), savannas (SAV), grasslands (GRA), and croplands (CRO). The majority of towers come from the FLUXNET2015 datasets (Pastorello et al., 2020), with additional towers from the AmeriFlux network (https://ameriflux.lbl.gov/, last access: 15 September 2025) and the European Eddy Fluxes Database Cluster (European Flux) (http://www.europe-fluxdata.eu/, last access: 15 September 2025). To ensure robustness, we applied the following selection criteria: (a) each land cover type and climate zone should contain at least two towers, and two-year or longer measurements are available for each tower; (b) the land cover type from tower footprints should be consistent with the dominant type from coarse resolution remote sensing data; (c) the towers should have both global integrity of biome and climate zone and regional representation. As a result, 106 flux towers are selected for the ET algorithm development and evaluation (Fig. 2). More details involving flux towers are available in Table A1 of Appendix A.

Figure 2 Locations of 106 flux towers used for ET algorithm development and validation. The symbols represent land cover types and colored backgrounds represent climate zones.

The tower measurements including air temperature, VPD, wind speed, net radiation, air pressure, surface soil moisture, and latent flux are used to drive and evaluate ET algorithm. The 30 min measurements from AmeriFlux and European Flux are aggregated into daily averages when at least 80 % of the measurements are available.

3.2 Climatic and remote sensing forcing data

The global datasets used for ET generation include the satellite remote sensing inputs, daily surface meteorology, and climate zone database. The satellite inputs consist of NDVI, net radiation, and land cover, while daily surface meteorology and soil moisture data come from reanalysis datasets.

The arid index (AI) data come from the Global Aridity Index and Potential Evapotranspiration Database – Version 3 (Zomer et al., 2022). Based on the AI and the distribution of flux towers, the globe is divided into dry (AI < 0.65) and wet zones (AI ≥ 0.65). This separation acknowledges that soil moisture constraints and parameters differ substantially between dry and wet zones, even within the same biome.

The NDVI sequence is merged from three independent products of the AVHRR GIMMS NDVI (Tucker et al., 2005), University of Arizona Vegetation Index and Phenology Lab (VIP) NDVI (Didan, 2010), and MODIS NDVI (Didan, 2015) using the equidistant cumulative distribution function (EDCDF) matching method (Zhang et al., 2015). The adjusted GIMMS and VIP sequences are fused with the continuously updated MODIS sequence to produce a consistent long-term $\mathrm{1}/\mathrm{12}$° NDVI series. The temporal linear interpolation is applied to generate daily NDVI sequence from the original adjacent semi-monthly values, providing a simple but effective method to produce daily series of vegetation observations (Zhang et al., 2010).

The daily satellite-based net radiation term is derived from the NASA World Climate Research Programme/Global Energy and Water Cycle Experiment (WCRP/GEWEX) Surface Radiation Budget (SRB) Release-3.0 datasets and the Clouds and the Earth's Radiant Energy System (CERES) SYN1deg radiative fluxes. Since they are two independent datasets and neither dataset spans the entire period since 1982, the same merging method used for NDVI is applied to fuse SRB with the latest CERES. The 1° coarse resolution radiation inputs are further resampled to $\mathrm{1}/\mathrm{12}$° resolution using the widely-used bilinear interpolation method.

For land cover types, the MODIS Land Cover Dynamics product at 0.05° resolution (Sulla-Menashe et al., 2019) is selected. To match the spatial resolution of other inputs, the land cover type is resampled to $\mathrm{1}/\mathrm{12}$° resolution based on the dominant type within each pixel.

The daily meteorological data, including air temperature, vapor pressure, and wind speed, are obtained from the NCEP/DOE AMIP-II Reanalysis (NCEP2) (Kanamitsu et al., 2002). Despite its coarse spatial resolution of 1.9° × 1.875°, NCEP2 is chosen for its long-term consistency and generally high quality against ground observations (Zhang et al., 2015). The coarse-resolution meteorology is further resampled to $\mathrm{1}/\mathrm{12}$° resolution using bilinear interpolation. In addition, monthly CO₂ concentration records come from The Global Monitoring Laboratory of the National Oceanic and Atmospheric Administration (ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_gl.txt, last access: 15 September 2025).

The surface soil moisture data are derived from Global Land Data Assimilation System (GLDAS) Noah Land Surface Model (Rodell et al., 2004) for its long-term temporal consistency and global spatial coverage. To cover the entire period from 1982 to 2023, we use two versions: version 2.0 for data prior to 2000 and version 2.1 for data from 2000 onwards. For version 2.1, 3-hourly soil moisture data are averaged to produce daily values.

3.3 Other remote sensing evapotranspiration and associated datasets

To evaluate and compare our ET estimates with other products, two mainstream remote sensing ET datasets with high reputations are selected. The first is the Global Land Evaporation Amsterdam Model (GLEAM) Version 3.3a (Martens et al., 2017), which is a microwave-based dataset with a 0.25° daily resolution, covering a long time span since 1980. The second is the Penman-Monteith-Leuning (PML) Version 2 (Zhang et al., 2019), a diagnostic biophysical model that couples transpiration with gross primary productivity and offers a 500 m resolution and an 8 d interval since 2002.

Our new ET estimates are also validated against the reconstructed ET estimates using the terrestrial water balance method at the basin level. One dataset generated by Wan et al. (2015) provides monthly ET benchmarks of 592 subbasins across the conterminous United States (CONUS). The ET item is regarded as the residual of precipitation, runoff, and downscaled water storage changes. These subbasins, which range from 292 to 303 700 km², cover 73 % of the CONUS, making the dataset well-suited for evaluating remote sensing ET products. For our evaluation, we use the monthly sequence from 2003 to 2008 to evaluate ET estimations from P-LSHv2, P-LSHv1, GLEAM, and PML, as well as other remote sensing ET datasets presented in our previous study (Chao et al., 2021). Another dataset (Pan et al., 2012) provides monthly water budget estimates for 32 major global basins by merging numerous global datasets. Here we use its ET estimates, in which the imbalance error has been resolved by the constrained Kalman filter technique. The 32 basins cover all continents except Antarctica and represent a wide range of climates and land covers, spanning from 1984 to 2006.

4 Results

4.1 Uncertainty analysis and optimization of parameters

We used ET measurements from 106 global flux towers to conduct an uncertainty analysis and optimize the parameters of the P-LSHv2 algorithm. Details about these towers are provided in Appendix A. The 106 flux towers cover 9 main land cover types, which represent 97.1 % of the global land area, excluding permanent ice, snow, barren lands, and open water bodies. Each land cover type is further divided into dry and wet zones, where key parameters are identified and calibrated separately. Since the dry zone constitutes only 1.7 % of the EBF type, no further classification was made for EBF and there are 17 scenarios in total.

A Sobol' sensitivity analysis revealed significant variations in parameter sensitivity across multiple scenarios, as shown in Fig. 3. In general, parameters b₁ and b₃, which govern canopy conductance, demonstrate the highest sensitivity in the P-LSHv2 algorithm. The first-order indices for b₁ and b₃ are 26 % and 18 %, respectively, while their total-order indices are 55 % and 40 %. Next in sensitivity are β and T_opt, both of which are associated with temperature stress factors in canopy conductance estimation, with total-order indices of 18 % and 15 %, respectively. Other parameters such as r_tot, k, n, b₂, and r_c, which are involved in soil evaporation and canopy conductance, exhibit relatively low sensitivity, with total-order indices of 8 %, 5 %, 5 %, 2 %, and 1 % respectively. Parameters related to temperature and VPD thresholds (T_{close_min}, T_{close_max}, VPD_open, and VPD_close) are insensitive, with first-order and total-order indices close to 0. To summarize, among 13 parameters in the P-LSHv2 algorithm, b₁ and b₃ are the most sensitive, while β and T_opt have moderate sensitivity. The r_tot, k, n, b₂, and r_c exhibit weak sensitivity, and T_{close_min}, T_{close_max}, VPD_open, and VPD_close are found to be insensitive. Consequently, in subsequent parameter uncertainty analysis, these four insensitive parameters are removed, and the remaining nine sensitive parameters are retained for optimization.

Figure 3 The (a) first-order sensitivity index and (b) total-order sensitivity index of the parameters across diverse land cover types and climate zones. The average values of diverse zones are shown in (c).

The Differential Evolution Markov Chain (DE-MC) method is used for parameter calibration, with the root mean square error (RMSE) as the objective function, as well as the Nash-Sutcliffe efficiency coefficient (NSE) and coefficient of determination (R²) to quantify mismatches. 10 Markov chains were set to perform 20 000 iterations for parallel crossover operations, with a burn-in period of 5000 iterations. Each parameter can converge stably across all 17 scenarios. The posterior distributions of parameters, with 95 % high probability intervals, are shown in Fig. 4. Compared with the wider prior intervals, most parameters show significant updating and reduction in uncertainty.

Figure 4 Boxplots showing the 95 % highest posterior density intervals of parameters across diverse land cover types and climate zones. The red lines indicate the medians, gray dots represent outliers, and black dashed lines mark the bounds of the uniform prior distributions.

In general, all parameters exhibited considerable variability across scenarios, reflecting the diverse responses of different vegetation and soil to the ET process in diverse climates. Among all biome-and-climate-specific parameters, the parameter n, which quantifies the soil moisture constraint on vegetation transpiration, warrants further discussion. In all land cover types except EBF and CRO, n is noticeably higher in dry zones compared to wet zones, indicating stronger soil moisture constraints in dry regions. In almost all wet zones, except for OSH and SAV, n is close to 0, suggesting negligible soil moisture constraints. Even wet OSH and wet SAV represent weak constraint levels, with n of 23.1 and 19.8, respectively. The n values vary widely in dry zones, with the highest 92.9 of OSH, suggesting that only 7.1 % of the samples are not stressed by soil moisture. The lowest n in dry zones is 0 of CRO, suggesting that the cropland ecosystem is still not affected by soil moisture even in dry zones. This is because cropland ecosystem is generally subject to regular irrigation schedules, which ensure sufficient water supply for crop growth and ET dissipation, while the identification of climate zones only considers precipitation and potential ET without accounting for irrigation. Further discussion on uncertainties related to cropland moisture constraint is available in Sect. 4.5. Compared to shrublands and grasslands, the soil moisture constraint of forest towers is generally low. This may partially be attributed to the deep root system of forests, which extends beyond the monitoring coverage of the surface layer. Even so, the surface soil moisture still captures a certain level of moisture constraints on forest ET, with n values of 24.9, 7.8, and 16.0 for dry ENF, dry DBF, and dry MF, respectively.

In each specific scenario, the parametern is treated as a constant. Therefore, by combining global patterns of land cover types and climate zones, the spatial patterns of n are plotted in Fig. 5, revealing distinct gradients dominated by surface land cover and climate. The regions where n value exceeds 80, 50, and 20 account for 10.2 %, 36.0 %, and 44.8 % of the global vegetation area, respectively. The regions with strong constraints are primarily located in the Patagonian Plateau of South America, the western United States, central and southern Africa, central Asia, Oceania, and parts of the pan-Arctic where dry shrublands and grasslands are typically dominant. In contrast, regions with weak or no moisture constraints are concentrated in forests, croplands, and most other wet areas.

Figure 5 Global map of the parameter n on vegetation areas covering 9 land cover types.

The g₀ ∼ NDVI curves, determined by parameters b₁, b₂, and b₃, are plotted in Fig. 6. Despite sharing a similar sigmoid-type structure, these curves exhibit distinct biological and climatic differences, reflecting variations in leaf traits and physiology. In ENF, DBF, SAV, and GRA, the curves for dry and wet zones are largely consistent, indicating that vegetation stomata respond similarly to NDVI, regardless of climate. However, in some specific NDVI intervals, the corresponding g₀ in dry zones is notably higher than that in wet zones for MF, OSH, WSA, and CRO types. One explanation is that this is caused by the actual response mechanism of vegetation stomata under different climates, indicating that even in the same ecosystem, different climates have an essential impact on the maximum stomatal conductance. More importantly, some vegetation in dry zones is subject to strong moisture constraints from VPD and soil moisture simultaneously. In these cases, the linear structure of constraints in our scheme may reinforce the moisture constraint, which is then compensated for in the potential stomatal conductance response curve.

Figure 6 The sigmoid-type response curve between g₀ and NDVI. The colors represent climate zones and n is the parameter to control soil moisture constraints.

4.2 Performance related to flux towers

We estimated daily ET at 106 global flux towers using the optimized P-LSHv2 algorithm, driven by tower-based measurements of radiation, meteorology, soil moisture, and remote sensing-based NDVI. The estimated ET was then compared against flux tower measurements for evaluation. Due to the limited samples of flux towers in each scenario, the samples were not divided into calibration and validation groups. Instead, all samples were used for parameter calibration (i.e. calibration mode) to evaluate the algorithm performance. The ET estimates of P-LSHv2 and P-LSHv1 were evaluated at each flux tower, as illustrated in Fig. 7. Overall, the P-LSHv2 algorithm is significantly upgraded. The average RMSE across 106 flux towers drops from 0.79 to 0.67 mm d⁻¹, a decrease of 15.2 %. Similarly, the average NSE increases from 0.43 to 0.58, while the average R² sees an increase from 0.62 to 0.67. Among various land cover types, DBF, MF, OSH, WSA, and SAV experienced the most pronounced improvements, with RMSE reduction exceeding 0.15 mm d⁻¹, NSE rise level exceeding 0.15, and R² rise level ranging from 0.03 to 0.17.

Figure 7 Comparison of ET estimates from the P-LSHv2 algorithm (black circles) and the P-LSHv1 algorithm (gray circles) at 106 flux towers. Each circle represents a flux tower. The right panels summarize statistical metrics across all sites, as well as separately for dry and wet conditions. RMSE: root mean square error; NSE: Nash-Sutcliffe efficiency; R²: coefficient of determination.

The differences between the two algorithms are more distinct when comparing ET estimations in wet and dry zones. On one hand, the ET estimation in dry zones is generally worse than that in wet zones, the average RMSE of the P-LSHv1 algorithm is 0.89 and 0.70 mm d⁻¹ (dry zone versus wet zone), and the average RMSE of the P-LSHv2 algorithm is 0.75 and 0.60 mm d⁻¹ (dry zone versus wet zone). This demonstrates that the physical ET process in dry zones is generally more challenging to simulate than in wet zones. On the other hand, the rising level of the P-LSHv2 estimation in dry zones is higher than that in wet zones. Among all 106 towers, the dry towers show an RMSE decrease of 0.14 mm d⁻¹, an NSE increase of 0.19, and an R² increase of 0.07, compared to RMSE reductions of 0.10 mm d⁻¹, NSE improvements of 0.12, and R² gains of 0.04 in wet zones. We selected several typical dry towers to compare daily ET estimations between two algorithms, as shown in Fig. 8. At these towers, particularly US-Me6, US-Whs, and US-SRM, the ET estimations significantly improved during water-limited periods, with low values effectively constrained. The three towers represent ENF, OSH, and WSA, respectively, with high n values (24.9, 92.9, and 58.8) and heavy moisture constraints, which further demonstrates the reliability of moisture constraint scheme in P-LSHv2 algorithm. These findings suggest that although the new algorithm has been upgraded in both zones, the soil moisture constraint plays a more essential role in ET simulation in dry zones than in wet zones.

Figure 8 Time series of daily measured and modelled latent heat flux (LE: W m⁻²) using the P-LSHv2 and P-LSHv1 algorithms for typical dry towers.

Due to limited samples, all tower measurements for each scenario, covering the entire period of available records, were used for parameter calibration in the preceding parameter characteristic analysis and algorithm evaluation. To further illustrate the robustness of the P-LSHv2 algorithm and its parameter optimization scheme, we designed another controlled trial using the leave-one-out cross-validation method. In this trial, daily ET measurements from all flux towers within each scenario were randomly split into two groups – one for model calibration and the other serving as the reference truth value for validation using the optimized parameters. This process was repeated twice to ensure that both groups were used both for calibration and validation. The pink and dark symbols in Fig. 9 almost overlap, representing calibration and cross-validation modes respectively, indicating that the performance of the P-LSHv2 algorithm in cross-validation mode is nearly identical to that in the calibration mode. When rounded to two decimal places, no differences are observed compared to the calibration mode. The largest difference is observed in wet WSA types, where the RMSE difference is only 0.014 mm d⁻¹. These results confirm that the P-LSHv2 algorithm outperforms the P-LSHv1 algorithm in both the calibration and cross-validation modes, demonstrating its stability and robustness in the calibration mode across all climatic and land cover types.

Figure 9 Taylor diagram evaluating the P-LSHv2 and P-LSHv1 algorithms across all scenarios. Different symbols represent different land cover types. The blue, pink, and dark brown represent the ET estimation statistics from the P-LSHv1 algorithm, the P-LSHv2 algorithm in calibration mode, and the P-LSHv2 algorithm in cross-validation mode, respectively.

4.3 Multiscale comparison of the P-LSHv2 results with the other remote sensing products at the basin scale

The P-LSHv2 algorithm is further evaluated with other global remote sensing ET datasets. Since these datasets are raster-based and derived from remote sensing and reanalysis, their evaluation against flux towers typically encounters issues related to spatial scale inconsistency and upscaling accuracy. Therefore, the reconstructed values from the water balance method (ET_recon) at the basin scale are selected as the benchmark to assess the performance of various datasets at monthly, annual, and multiyear levels.

4.3.1 Comparison of ET remote sensing products across the CONUS

We evaluated gridded ET estimation of 592 basins across CONUS from the P-LSHv2 algorithm, using the reconstructed values from the water balance method as the benchmark. The results were also compared with the P-LSHv1 algorithm and two mainstream products: GLEAM and PML.

At the monthly level, the ET estimation from P-LSHv2 agrees well with ET_recon, accounting for 71 % of the variation, with an RMSE of 21.4 mm month⁻¹ and an NSE of 0.71 (Fig. 10). Compared to P-LSHv1, the RMSE value drops by 7.3 %. Moreover, the performance of P-LSHv2 is comparable to that of GLEAM and PML, with acceptable agreement considering the large number of basins and the relatively fine temporal resolution.

Furthermore, we aggregated various monthly ET estimations on an annual basis. On the annual level, the ET estimation from P-LSHv2 upgrades significantly from its predecessor, with an RMSE dropping from 114.4 to 104.4 mm yr⁻¹, a decrease of 8.7 % (Fig. 10). Particularly in arid areas (i.e., the basin average AI < 0.65), the P-LSHv2 results effectively correct the overestimation seen in P-LSHv1, reducing the deviation from −44.1 to 27.2 mm yr⁻¹. In the driest basins (AI < 0.3), the deviation is reduced even further, from −64.2 to 13.8 mm yr⁻¹. These improvements highlight the efficacy of the soil moisture constraints in P-LSHv2, especially in modelling water availability in arid regions. From the perspective of RMSE, NSE, and R², the ET estimation from P-LSHv2 is also better than that from GLEAM and PML (Fig. 10), with an RMSE of 104.4 mm yr⁻¹, an NSE of 0.77 and an R² of 0.80.

As for the multiyear average scale, the P-LSHv2 results have an RMSE value of 80.4 mm yr⁻¹, a decrease of 13.2 % compared to its predecessor. In regions with high ET density (around 400 mm yr⁻¹), both P-LSHv2 and GLEAM, which account for soil moisture constraints, deliver accurate ET estimates, while PML and P-LSHv1 tend to overestimate, which is an understandable result of their failure to account for soil moisture constraints. In mid-to-high ET intervals (500 to 1000 mm yr⁻¹), GLEAM shows higher uncertainty than P-LSHv2, as reflected in its lower NSE and R² values, indicating the significant impact of different soil moisture constraint schemes on ET estimation.

Figure 10 Comparison of monthly, yearly, and multiyear average ET derived from four remote sensing datasets (ET_RS) and ET_recon on 592 basins of CONUS. Colors are used to represent density, with red representing high density and blue representing low density.

The P-LSHv2 ET estimation also captures spatial gradients across CONUS, with a general relative deviation of −5.7 %, significantly lower than the 8.8 % deviation of P-LSHv1. About 95.5 %, 78.0 %, and 35.7 % of the area (weighted by size) have relative deviations within ±30 %, ±15 %, and ±5 %, respectively. Only about 5 % of the area experiences underestimation or overestimation beyond 30 %. Considerable underestimations are scattered across the Colorado basin, California basin, and parts of the northeast and northwest, while overestimations are concentrated in the West Gulf and North Central basins (Fig. 11). These deviations are likely influenced by reservoir control and interbasin water transfers, which are not accounted for in the model (Wan et al., 2015). Additionally, uncertainties in input data such as precipitation (Daly et al., 2008) and water storage (Landerer and Swenson, 2012) contribute to discrepancies in ET reconstruction. Despite these localized deviations, the general close spatial patterns and minimal differences indicate the high reliability and robustness of P-LSHv2 across CONUS.

Figure 11 Spatial patterns of multiyear average ET from ET_recon and P-LSHv2 and their difference (RE) on 592 basins of CONUS.

4.3.2 Comparison of ET remote sensing products across global large basins

Apart from CONUS, we also investigated the differences between various remote sensing ET datasets across 32 large basins around the world. Given that the ET_recon of these basins spans 1984 to 2006, the shorter PML dataset, which starts from 2002, is excluded.

On the monthly level, GLEAM, P-LSHv1, and P-LSHv2 all accurately capture seasonal ET variations across these basins (Fig. 12), with RMSE values below 15 mm month⁻¹ and R² values exceeding 0.89. Among them, P-LSHv2 achieves the lowest RMSE (9.8 mm month⁻¹), improving accuracy by 22.8 % over P-LSHv1 and aligning more closely with ET_recon than GLEAM.

At the annual scale, P-LSHv2 demonstrates a 34.6 % improvement in accuracy over P-LSHv1 and shows better agreement with ET_recon than GLEAM. Both GLEAM and P-LSHv1 tend to overestimate ET in these 32 basins (Fig. 12), with respective relative deviations of −23.4 % and −20.0 %, whereas P-LSHv2 reduces the deviation to −10.9 %, indicating that soil moisture constraint incorporated in our algorithm plays an essential role in correcting the overestimation. Compared to the other two datasets, P-LSHv2 particularly agrees well in the low-value intervals below 300 mm yr⁻¹ and the high-value intervals above 800 mm yr⁻¹.

The patterns of multiyear (from 1984 to 2006) average ET estimations of each dataset are similar to those of the annual basis, and P-LSHv2 has lower RMSE, higher NSE, and higher R² than the other datasets (Fig. 12). These comparisons confirm that P-LSHv2 performs well in capturing ET patterns across various climate zones and ecosystems around the world.

Figure 12 Comparison of monthly, yearly, and multiyear average ET derived from three remote sensing datasets (ET_RS) and ET_recon on 32 global large basins.

The spatial patterns of deviations from different datasets vary significantly across 32 basins (Fig. 13). GLEAM and P-LSHv1 show relative deviations exceeding 20 % in 19 and 20 basins, respectively, while P-LSHv2 has deviations above 20 % in only 10 basins. In basins such as the Amazon, Parana, Mississippi, Niger, Danube, Indus, and Murray-Darling, the deviations of the three datasets all fall within ±20 %, demonstrating their strong capability in representing ET budgets at large basin scales. However, some overestimation persists across datasets, notably in the Yukon, Nile, Ural, Indigirka, Amur, and Yangtze basins. This overestimation may be attributed to uncertainties in the water balance model. Since the benchmark data extend to the pre-GRACE era, Pan et al. (2012) drove the Variable Infiltration Capacity hydrological model to simulate water storage changes, introducing potential errors despite employing Kalman filter techniques to correct for balance errors. Despite some uncertainties, P-LSHv2 has better agreement with ET_recon in most basins and shows better accuracy than other remote sensing datasets.

Figure 13 Spatial patterns of relative deviations of multiyear average ET from three remote sensing datasets on 32 global large basins.

4.4 Global ET patterns

We calculated the global daily terrestrial ET at a spatial resolution of $\mathrm{1}/\mathrm{12}$° from 1982 to 2023 using the improved P-LSHv2 algorithm, excluding permanent snow and ice. The global multiyear average ET reveals distinct regional variations and latitudinal gradients (Fig. 14). The area-weighted global terrestrial multiyear average ET of P-LSHv2 is 524.8 mm yr⁻¹, lower than the 562.2 mm yr⁻¹ estimated by P-LSHv1. In dry and wet zones, the area-weighted multiyear averages of P-LSHv2 are 358.8 and 761.4 mm yr⁻¹, respectively, both lower than the corresponding estimates of 383.8 and 816.4 mm yr⁻¹ from P-LSHv1. The spatial patterns between the two datasets are generally similar, with deviations between P-LSHv2 and P-LSHv1 falling within ±100 mm yr⁻¹ for 56.6 % of area-weighted regions. In contrast, 29.5 % of the area shows a negative deviation below −100 mm yr⁻¹ and 14.0 % of the area shows a positive deviation above 100 mm yr⁻¹. The overestimations of P-LSHv1 are particularly concentrated in regions such as the Andes Mountains and Patagonia Plateau of South America, southern North America, and central and southern Africa. These areas are typically composed of OSH, WSA, and SAV biomes, where P-LSHv2 accounts for substantial soil moisture constraints. The higher estimates of P-LSHv2 are located in open water bodies (Great Lakes and the Caspian Sea), bare soils (Sahara Desert, Arabian Peninsula, and central Asia), and some grasslands in wet-dry boundary regions of the Southern hemisphere.

Figure 14 Global map of multiyear average ET from P-LSHv2 (1982–2023) and P-LSHv1 (1982–2015), as well as their differences and latitudinal profiles.

In terms of temporal trends, the P-LSHv2 dataset follows a similar trajectory to its predecessor, both exhibiting significant upward trends (Fig. 15). However, the linear upward trend of P-LSHv2 from 1982 to 2015 is 0.69 mm yr⁻² (p<0.001), slightly lower than that of P-LSHv1, with a trend of 0.87 mm yr⁻² (p<0.001). The trend of P-LSHv2 is comparable to PML (0.68 mm yr⁻²), and higher than GLEAM (0.38 mm yr⁻²). P-LSHv2 ET increased by 0.46 mm yr⁻² (p<0.001) from 1982 to 2023, although the rate of increase appears to have slowed in recent years.

Figure 15 Annual global ET estimates derived from P-LSHv1 and P-LSHv2.

5 Discussion

Although the P-LSHv2 dataset has good accuracy at multiple temporal and spatial levels around the world, there are still some limitations and uncertainties of the algorithm that are worth discussing.

Soil structural heterogeneity plays a crucial role in regulating root development and water availability, which further constrain ET. Looser soils facilitate root penetration, while denser soils can retain water at higher matric potentials, thus extending water availability during dry periods. These physical characteristics influence not only the spatial distribution of roots but also the efficiency of root water uptake and overall transpiration dynamics. Incorporating soil hydraulic properties (e.g., soil water retention curves and hydraulic conductivity) would ideally allow more accurate constraint of soil moisture availability in evapotranspiration estimation. However, global soil hydraulic datasets remain highly uncertain due to limited in situ observations, spatial variability, and differences in pedotransfer function assumptions. To address this challenge, we adopted a simplified quantile-based method to parameterize the soil moisture constraint function in our P-LSHv2 algorithm. Specifically, the parameter n was determined based on land cover and climate types, which serve as proxies for large-scale heterogeneity in both soil properties and vegetation root strategies.

In addition, the soil moisture data used in this study are derived from the GLDAS Noah land surface model rather than satellite remote sensing products. This choice is due to the challenges in obtaining over 40 years of continuous global remote sensing data. The GLDAS Noah model effectively captures seasonal variations and anomalies in soil moisture (Chen and Yuan, 2020; Spennemann et al., 2020) and performs well in ET moisture constraints (Feng et al., 2023). We used data from versions 2.0 and 2.1, covering the period from 1982 to 2023, without applying data fusion, as they originate from the same model and data source. In the P-LSHv2 algorithm, we used surface soil moisture as a moisture constraint surrogate for the entire soil layer. On the one hand, the current satellite remote sensing retrievals struggle to capture soil moisture characteristics of deep layers, and uncertainties associated with deep soil moisture from reanalysis or land surface models are generally greater than those from surface layers. On the other hand, the root zone depth of most tree species tends to be relatively shallow (Flo et al., 2022; Tumber-Dávila et al., 2022), and correlations between soil moisture thresholds at varying depths have been reported (Fu et al., 2022b). Regardless, further model refinement and testing of the P-LSHv2 algorithm could be conducted when high-quality root zone soil moisture datasets become available.

Although GLDAS data were used for soil moisture inputs, the other meteorological forcing variables – including air temperature, vapor pressure, and wind speed – were obtained from the NCEP2 reanalysis dataset. This choice ensures consistency with the original version of our algorithm (Zhang et al., 2015), enabling a fair comparison of model improvements independent of changes in input data. However, we acknowledge that reanalysis datasets have substantially improved in both spatial resolution and methodological sophistication in recent years, and NCEP2 may not capture fine-scale meteorological variability as accurately as newer products. On one hand, the relatively coarse resolution of NCEP2 may lead to an underestimation of spatial heterogeneity after interpolation to finer grids. Although the spatial variability of air temperature, vapor pressure, and wind speed is generally lower than that of precipitation or radiation, the coarse resolution could still smooth out meaningful local variations, introducing moderate uncertainty. On the other hand, biases inherent in NCEP2 – such as systematic overestimation or underestimation of specific variables – may propagate through the model and affect the magnitude and spatial pattern of ET simulations. Given these limitations, adopting higher-resolution and more recent reanalysis products represents a valuable direction for future studies, as it could improve the representation of local-scale meteorological variability and reduce uncertainties associated with coarse input data.Apart from the forcing data used for ET estimation, there are still certain uncertainties in other ET observations or reconstructions used for algorithm calibration and verification. For instance, the heat flux measurements from eddy covariance towers may suffer from energy imbalance due to the complex circumstances of wind, footprint, and sampling variability (Wilson et al., 2002). However, it has been shown that there is little difference between daily latent heat flux without the energy closure correction and after the correction (He et al., 2022), and the uncorrected latent heat flux is also widely used in the evaluation of ET models (Mu et al., 2011; Zhang et al., 2010). As for the benchmark data from 592 basins across CONUS and 32 global large basins, both introduced water storage change data from land surface models to either downscale or extend temporal coverage, which leads to some uncertainty in small basins or extended years. In the future, further enhancing the reliability of flux tower observations and water balance reconstructions will promote the development of ET algorithms and models.

Besides, the interaction between soil moisture and ET has not been explicitly considered in this study. We have used existing soil moisture data as an input to estimate ET, without simulating the feedback loop between them. This is primarily due to the focus of this study on assessing the impact of soil moisture on ET estimations. However, it is recognized that SM and ET interact dynamically, and this is a limitation of the current approach. Future work could incorporate the interaction between soil moisture and ET to improve the model's representation of soil water dynamics.

Interestingly, our analysis revealed that dry zones generally experience higher soil moisture constraints than wet zones across most land cover types, with cropland being a notable exception. Contrary to other scenarios where transpiration in dry zones is severely limited by soil moisture, flux tower data indicated that the parameter n of cropland ecosystem was 0 for both dry and wet zones, suggesting that neither was water-limited. Further investigation into five dry cropland towers showed that three of them (US-Ne2, US-Twt, and IT-CA2) were subject to irrigation exceeding 50 mm yr⁻¹, a factor not accounted for in our climate zone classification. According to the global irrigation dataset from Zhang et al. (2022a), climate zones, and land cover type data used in this study, approximately 56.5 % of global dry cropland is sufficiently irrigated, primarily located in eastern China, India, and central North America. Understanding how to accurately quantify the unique characteristics of irrigated areas within the soil moisture constraint framework presents a compelling avenue for future ET simulation research.

Moreover, the current climate zone classification based on the aridity index (AI) may misrepresent energy-limited ecosystems, such as partitioning energy-limited cold land areas, including northern taiga and tundra, into the dry climate category. This is because AI is calculated as the ratio of precipitation to potential evapotranspiration (P/PET), which may lead to low AI values even when surface soils remain moist during summer growing season. As a result, these high-latitude regions may exhibit apparent soil moisture constraints and model underestimation of ET, despite being primarily limited by energy rather than water. While AI offers a globally consistent framework, future work could benefit from incorporating growing season length, energy inputs, and phenology to better represent cold-region ecohydrology in large-scale ET modeling.

6 Data availability

All input data used in this study are freely available (see Sect. 3). The P-LSHv2 ET dataset is freely available at the National Tibetan Plateau Data Center (https://doi.org/10.11888/Terre.tpdc.301969, Feng and Zhang, 2025). Daily, monthly, and yearly datasets are provided. The dataset is published under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.

7 Conclusion

In this study, we incorporated a new global soil moisture constraint scheme into ET estimation and developed a new version of remote sensing ET algorithm (P-LSHv2). The biome- and climate-specific parameters and soil moisture constraint levels were calibrated using 106 flux towers around the world. The calibration and cross-validation modes at the flux towers showed that the new algorithm is robust, and the improvement in the towers' ET simulation for dry zones is higher than that for wet zones. The investigation of flux towers also quantified the moisture constraint levels across diverse land cover types and climate zones, revealing that constraints in dry zones generally surpass those in wet zones, particularly for dry OSH, GRA, and SAV biomes. To evaluate the P-LSHv2 algorithm, we employed two ET benchmarks reconstructed using water balance methods, assessing performance across multiple temporal scales at the basin level. The P-LSHv2 estimation effectively captures the impact of soil moisture anomalies on ET, outperforming P-LSHv1, GLEAM, and PML datasets. Consequently, the P-LSHv2 algorithm provides a more realistic representation of ET processes, and the overall quality of the dataset has been improved, especially in areas where available water is limited. Based on the P-LSHv2 algorithm, a reliable and continuous long-term global ET dataset spanning from 1982 to 2023 has been generated, which contributes to the global assessment of ET climatology and a better understanding of terrestrial water and energy cycle dynamics under climate change conditions.

Appendix A

Table A1 Details of 106 flux towers used for calibration and validation of ET algorithm.

^∗ The EBF is no longer categorized into various climate zones, as less than 2 % of the global EBF falls within dry zones.