Extended global terrestrial evapotranspiration and gross primary production dataset from 1982 to near present

The Penman–Monteith–Leuning (PML) model is a widely recognized diagnostic framework for estimating coupled terrestrial evapotranspiration (ET) and gross primary production (GPP). To address the critical need for high-fidelity, long-term, and near-present eco-hydrological records, we developed the PML-V2.2 dataset, spanning from 1982 to 2025. Driven by observation-constrained Multi-Source Weighted-Ensemble Precipitation (MSWEP) and Multi-Source Weather (MSWX) meteorological variables, the dataset comprises three complementary products: (1) PML-V2.2a, an 8 d 500 m MODIS/VIIRS satellite-based product (2000–2024 and 2012–2025) optimized for near-present monitoring (updated annually); (2) PML-V2.2b, a half-month 0.1° AVHRR-based product (1982–2020) anchoring long-term climate attribution; and (3) PML-V2.2c, a consolidated half-month 0.1° record integrating the above products for seamless 44-year continuity (1982–2025). Our methodological framework features an expanded bottom-up calibration using 208 flux sites (∼ 1400 site-years) across various plant functional types (PFTs) and a refined parameterization that explicitly distinguishes between irrigated and rainfed croplands. This distinction effectively mitigated systematic biases in agricultural regions, reducing ET and GPP estimation errors by 8.7 % and 16.2 %, respectively. Performance evaluation reveals high accuracy across PFTs (cross-validation Nash-Sutcliffe Efficiency, NSE > 0.60, absolute bias < 5 %), while top-down water-balance validation across 56 large river basins during 1982–2016 and 152 basins during 2003–2020 confirms high reliability (NSE: 0.89–0.91) as compared with other products. The MODIS- and VIIRS-based PML-V2.2a datasets are internally consistent, and exhibit high agreement with PML-V2.2b during their overlapping period (NSE = 0.90 and 0.79 for annual ET and GPP anomalies), ensuring a seamless transition across satellite epochs. Based on the consolidated PML-V2.2c dataset, global terrestrial ET and GPP during 1982–2025 are estimated at 65.8 × 10³ km³ yr⁻¹ (with 58.2 % from transpiration) and 143.4 PgC yr⁻¹, respectively. Long-term analysis reveals significant (p < 0.05) increasing trends in GPP (0.343 PgC yr⁻²) and ET (0.019 × 10³ km³ yr⁻²) during 1982–2025, where vegetation greening impact on ET is partially offset by physiological water saving under rising atmospheric CO₂, consequently enhancing water use efficiency. By bridging the gap between satellite epochs, PML-V2.2 provides an internally consistent long-term global dataset for hydrology, ecology, and other Earth science studies. The dataset is freely accessible, with the 500 m resolution PML-V2.2a product hosted on Google Earth Engine, and all 0.1° PML-V2.2a/b/c versions archived at the National Tibetan Plateau Data Center under https://doi.org/10.11888/Terre.tpdc.303314 (Xu et al., 2026).

1 Introduction

Terrestrial evapotranspiration (ET), the primary pathway for water transfer from the land surface to the atmosphere, plays a pivotal role in regulating Earth's hydrological cycle and energy balance (Oki, 2006; Zhang et al., 2026). As a fundamental driver of surface-atmosphere interactions, ET dictates the partitioning of net radiation into latent and sensible heat fluxes, thereby modulating surface temperatures, boundary layer dynamics, and global atmospheric circulation. Beyond its physical impact, ET serves as a critical bridge between the water and carbon cycles via transpiration, which is coupled with gross primary productivity (GPP) through stomatal regulation. This mechanism inherently governs ecosystem water use efficiency (WUE) (Zhang et al., 2022; Li et al., 2023b; Yuan et al., 2025b). Furthermore, ET acts as a major moisture source for terrestrial precipitation (P) through moisture recycling processes – a feedback increasingly critical in a greening Earth (Hoek van Dijke et al., 2022; Cui et al., 2022; Yang et al., 2023; Yao et al., 2026). Given these multifaceted roles, the development of accurate, long-term, and spatially consistent ET and GPP datasets is imperative for addressing the escalating challenges in climate, hydrological, and ecological sciences amidst global warming, precipitation shifts, and vegetation changes (Zhang et al., 2016a; Zhao et al., 2022b; Xu et al., 2022; Zhang et al., 2023).

The increasing availability of satellite observations and the maturation of land surface modeling have fostered a diverse ecosystem of remote-sensing-based diagnostic ET and GPP products. Table 1 synthesizes current models, focusing on datasets with a spatial resolution of 0.25° or finer to ensure the spatial fidelity required for detailed regional and global assessments. Unlike prognostic land surface models (LSMs), these diagnostic frameworks directly assimilate satellite-derived vegetation dynamics (e.g., leaf area index, LAI) to circumvent the systematic biases often inherent in complex plant functional type (PFT) simulations. Broadly, these datasets can be categorized into three primary clusters based on their target variables (see Table 1 for references): (1) stand-alone ET algorithms rooted in Penman-Monteith (e.g., MOD16) or energy/water-balance principles (e.g., 3T, SiTH); (2) GPP-specific products driven largely by Light Use Efficiency (LUE) logic (e.g., MOD17, GLASS EC-LUE) or emerging proxies like SIF and NIRv (e.g., GOSIF); and (3) integrated frameworks that simultaneously estimate both fluxes to ensure physiological consistency (e.g., BEPS, BESS). Notably, within these categories, Machine Learning (ML) has emerged as a transformative force. This includes pure data-driven approaches that upscale flux observations directly (e.g., FLUXCOM-X-BASE, GloFlux) and hybrid strategies that embed ML into process-based schemes to refine specific parameters like stress factors (e.g., GLEAM4, CoSEB).

Table 1 Summary of major remote sensing-based diagnostic ET and GPP models and datasets.

Note: Datasets with a spatial resolution coarser than 0.25° are excluded from this summary. The literature search and data availability check were concluded on 15 July 2026. Access links for all products are provided in Table A1 in the Appendix. Time series of majority of datasets are shown in Fig. 13. In this study, “present” refers to data with a production delay of several days to one month (e.g., MODIS). Due to the latency in meteorological and remote sensing forcing data, “near present” refers to datasets updated annually through 2025.

Despite this methodological diversity, the current data landscape is characterized by a persistent trade-off between temporal span and biophysical consistency. This dichotomy is largely dictated by the evolution of satellite sensors: while higher-resolution diagnostic records are predominantly tethered to the post-2000 era of the Moderate Resolution Imaging Spectroradiometer (MODIS), and its continuity mission, the Visible Infrared Imaging Radiometer Suite (VIIRS), establishing pre-2000 baselines necessitates reliance on the Advanced Very High Resolution Radiometer (AVHRR). However, bridging these distinct epochs is difficult, and many long-term datasets covering earlier decades estimate ET and GPP in isolation. This challenge is further compounded by instabilities in meteorological driving data (e.g., solar radiation, vapor pressure deficit) derived from disparate reanalysis products, which collectively increase the uncertainty in detecting authentic climate-driven shifts (Kim et al., 2021; Xie et al., 2025). Beyond these consistency issues, a critical gap remains in data latency. As shown in Table 1, most biophysically consistent records suffer from significant production lags, precluding the inclusion of the near-present period (e.g., update annually to 2025). Consequently, few datasets are simultaneously long-term, consistent, and up-to-date, which hinders timely assessments of long-term changes and recent anomalies.

The Penman–Monteith–Leuning (PML) model was initially conceived as a high-precision ET diagnostic framework (PML-V1) based on the Penman-Monteith equation coupled with Leuning's canopy conductance theory (Leuning et al., 2008; Zhang et al., 2008), which successfully established a long-term global ET record covering 1982–2012 (Zhang et al., 2016c). Another theoretical advance occurred with the development of the PML-V2 model, a coupled water and carbon flux model that estimates canopy conductance through the GPP process, thereby providing a more physically consistent and high-accuracy estimation of both ET and GPP (Zhang et al., 2019; Gan et al., 2018). Implemented on the Google Earth Engine (GEE) platform (Gorelick et al., 2017), this second generation leveraged 95 FLUXNET2015 sites of flux observations and MODIS data for parameterization to achieve a spatial resolution of 500 m at an 8 d interval, initially providing global coverage for 2003–2017 (Zhang et al., 2019), while an optimized regional daily version for China was subsequently released covering 2000–2020, integrating 26 local flux sites and local climate forcing (He et al., 2022). Consequently, despite its high accuracy, the original PML-V2 framework lacks a consistent methodology to bridge the gap between AVHRR and MODIS observations, leaving a critical void in our understanding of long-term ecosystem responses and recent climate-driven anomalies.

To provide an integrated, high-fidelity record that resolves the trade-offs between spatial resolution and temporal span, we introduce the PML-V2.2 dataset. Built upon a robust diagnostic framework, PML-V2.2 extends the coupled record to nearly 44 years (1982–2025, near present and updated annually). The dataset is strategically released as three complementary products: the “a” version provides 500 m near-present capability based on MODIS and VIIRS; the “b” version utilizes AVHRR to anchor historical estimates essential for climate attribution; and the “c” version systematically integrates their strengths into a continuous, globally consistent record. This design ensures that PML-V2.2 serves as a timely resource, supporting both high-resolution operational monitoring and robust long-term trend analysis in terrestrial hydro-climatic studies. Our new PML-V2.2 datasets have the following features:

1. Consistent forcing data: All simulations (1982–2025) share the same forcing data from bias-corrected MSWEP V2.8 and MSWX-Past datasets (Beck et al., 2019, 2022) to ensure internally consistent modelling results.

2. Expanded calibration: A significantly larger number of flux-tower sites for model calibration (208 versus 95 used in Zhang et al., 2019).

3. Refined cropland estimates: More accurate ET estimates for croplands that are separated into irrigated and non-irrigated.

4. Optimized parameterization: The new two-parameter function for vapor pressure deficit is employed, reducing the number of parameters by one compared to the previous version (Zhang et al., 2019; He et al., 2022).

A detailed description of these new features is provided in Sects. 2 and 3.

2 Methods and materials

2.1 Penman-Monteith-Leuning model

The PML-V2 model partitions terrestrial ET into three constitutive components: vegetation transpiration (E_c), soil evaporation (E_s), and canopy interception (E_i) (Zhang et al., 2019). Within a Penman-Monteith (PM) framework, the model utilizes LAI for canopy-soil energy partitioning. Crucially, E_c is estimated through a GPP-coupled canopy conductance (G_c) scheme to ensure biophysical consistency (Gan et al., 2018). E_s is estimated by determining the soil moisture stress by accumulated precipitation and soil equilibrium evaporation (Zhang et al., 2010), while E_i is quantified via the Gash rainfall interception model (van Dijk and Bruijnzeel, 2001).

$$\begin{array}{}\text{(1)}& {\displaystyle }{E}_{\mathrm{c}}={\displaystyle \frac{\mathit{ϵ}{A}_{\mathrm{c}}+\left(\mathit{\rho }{c}_p/\mathit{\gamma }\right)D_{\mathrm{a}}{G}_{\mathrm{a}}}{\mathit{ϵ}+\mathrm{1}+G_{\mathrm{a}}/G_{\mathrm{c}}}}\text{(2)}& {\displaystyle }{E}_{\mathrm{s}}={\displaystyle \frac{f\mathit{ϵ}{A}_{\mathrm{s}}}{\mathit{ϵ}+\mathrm{1}}}\text{(3)}& {\displaystyle }{E}_{\mathrm{i}}=\left\{\begin{array}{ll}{f}_{\mathrm{v}}P,& P<P_{\mathrm{wet}}\\ f_{\mathrm{v}}{P}_{\mathrm{wet}}+f_{\mathrm{ER}}\left(P-P_{\mathrm{wet}}\right),& P\ge P_{\mathrm{wet}}\end{array}\right.\end{array}$$

where $\mathit{ϵ}=s/\mathit{\gamma }$, in which γ is the psychrometric constant (kPa °C⁻¹) and $s=\mathrm{d}{e}^{*}/\mathrm{d}T$ is the slope of the curve relating saturation water vapor pressure to temperature (kPa °C⁻¹). A_c and A_s are the available energy (net absorbed radiation minus soil heat flux, MJ m⁻² d⁻¹) separated by LAI for E_c and E_s. ρ is the density of air (g m⁻³). c_p is the specific heat of air at constant pressure (J g⁻¹ °C⁻¹). D_a is the water vapor pressure deficit of the air (kPa). G_a is the aerodynamic conductance (m s⁻¹). G_c is the canopy conductance (m s⁻¹). f is a dimensionless variable that determines the water availability for soil evaporation. P is the daily precipitation (mm d⁻¹). P_wet is the reference threshold rainfall amount if the canopy is wet (mm d⁻¹). f_ER is the ratio of average evaporation rate over average precipitation intensity storms (unitless), f_v is the fractional area covered by intercepting leaves (unitless).

The defining feature of PML-V2 is its coupled water-carbon mechanism. Unlike traditional PM-based models, it derives canopy conductance (G_c) directly from the photosynthesis process, constrained by vapor pressure deficit (D_a). Furthermore, it incorporates the impact of atmospheric CO₂ concentration (C_a) on carbon assimilation and G_c, providing a more physically robust simulation of ET and GPP under changing environmental conditions.

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{GPP}=A_{\mathrm{c},\mathrm{g}}\cdot f\left(D_{\mathrm{a}}\right)\text{(5)}& {\displaystyle }f\left(D_{\mathrm{a}}\right)=\left\{\begin{array}{ll}\mathrm{1},& D_{\mathrm{a}}\le D_{\min }\\ \frac{D_{\max }-D_{\mathrm{a}}}{D_{\max }-D_{\min }},& D_{\min }<D_{\mathrm{a}}<D_{\max }\\ \mathrm{0},& D_{\mathrm{a}}\ge D_{\max }\end{array}\right.\end{array}$$

where A_c,g is the gross assimilation rate without D_a constraint. f(D_a) is the D_a constraint function adopted from Wang et al. (2014) and Running et al. (2015).

In this model, the carbon assimilation and water exchange are subject to a unified constraint from D_a. Therefore, the same function f(D_a) can be used for deriving G_c:

$$\begin{array}{}\text{(6)}& G_{\mathrm{c}}=\mathrm{1.6}\cdot m\cdot A_{\mathrm{c},\mathrm{g}}\cdot f\left(D_{\mathrm{a}}\right)/C_{\mathrm{a}}\end{array}$$

where m is a dimensionless stomatal conductance coefficient. Note that 1.6 is the ratio of stomatal conductance to water vapor relative to that to carbon flux (Medlyn et al., 2011; Yebra et al., 2013; Gan et al., 2018).

More details on the PML modelling can be found on the supplement of Zhang et al. (2019). Compared with that, the D₀ parameter and the stress function of the original Ball-Berry-Leuning equation are replaced with the two-parameter function. One parameter can be reduced after this update.

2.2 Bottom-up parameterization and validation

The ten parameters of the PML-V2.2 model were optimized and validated using a comprehensive bottom-up framework based on the latest global eddy covariance (EC) flux observations, yielding the optimal parameter sets presented in Table 2. Crucially, versions a and b share this exact parameterization scheme, and any product discrepancies arise solely from their respective forcing datasets. The parameter optimization was performed exclusively at an 8 d temporal scale, which was chosen because the MODIS vegetation and albedo data employed in our framework provide an optimal balance between data quality and spatiotemporal resolution. The optimization process utilized these MODIS products alongside the highest-quality site-level meteorological observations, predominantly covering the period from 2000 to 2022.

Table 2 Optimized PML-V2.2 model parameters for the harmonized IGBP and CCI PFTs.

^a IGBP PFT are reported by data suppliers of the flux network. ^b Due to limited number of shrub sites for robust model calibration, the open shrublands and closed shrublands are combined. ^c ESA CCI PFT has four types of SHRUBS (NE, BE, ND, BD), which are also combined. ^d Croplands are split as rainfed (RFD) and irrigated (IRR) types based on the IGRIP30 dataset. The corresponding human-managed grass (GRASS-MAN) class is also split. Note that β is initial slope of the light response curve to assimilation rate (i.e. quantum efficiency), η is initial slope of the CO₂ response curve to assimilation rate (i.e. carboxylation efficiency), V_m is maximum catalytic capacity of Rubisco per unit leaf area at 25 °C. k_Q is extinction coefficient of PAR. k_A is extinction coefficient of available energy. D_minand D_max are thresholds below which there is no vapor pressure constraint and highest constraint, respectively. m is stomatal conductance coefficient. S_l is specific canopy rainfall storage capacity per unit leaf area. F₀ is specific ratio of average evaporation rate over average rainfall intensity during storms per unit canopy cover.

Flux measurements were compiled from six major international monitoring networks including FLUXNET2015 (Pastorello et al., 2020), the Australian and New Zealand flux tower network (OzFlux) (Beringer et al., 2016; Isaac et al., 2017; Beringer et al., 2022), AmeriFlux (Novick et al., 2018; Chu et al., 2023), the Integrated Carbon Observation System (ICOS) (Heiskanen et al., 2022; Warm Winter 2020 Team et al., 2022; ICOS RI et al., 2025), JapanFlux (Ueyama et al., 2025), ChinaFlux (Yu et al., 2006, 2024), and other independent data providers (Liu et al., 2018, 2013b; Ma et al., 2020; Zhou et al., 2023). For sites with missing data, gap-filling and flux partitioning were performed using the REddyProc R package (Wutzler et al., 2018). This processing pipeline aligns with the FLUXNET2015 data protocol (Pastorello et al., 2020), incorporating friction velocity (u^*) threshold filtering, Marginal Distribution Sampling (MDS) for gap-filling, and standardized carbon flux partitioning. GPP was calculated as the average of estimates derived from both daytime and nighttime partitioning methods (Reichstein et al., 2005; Lasslop et al., 2010). To ensure temporal continuity, missing meteorological forcing data were supplemented using Multi-Source Weighted-Ensemble Precipitation (MSWEP) V2.8 and Multi-Source Weather (MSWX) Past version (Beck et al., 2019, 2022), which integrate multi-source in-situ observations, satellite-based estimates and rescaled European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis V5 (ERA5) data. Following a protocol consistent with FLUXNET's use of ERA5-Interim, these seamless daily fields were subsequently used to drive the PML model. Finally, all observations were aggregated into an 8 d timescale, with a Quality Control (QC) flag assigned to represent the proportion of measured and high-quality gap-filled data. To maintain high data standards, only records with QC > 0.9 were retained for further calibration. Notably, rainy days were not filtered out, as precipitation noise is naturally smoothed at the 8 d scale. Moreover, removing them would severely disrupt data continuity, particularly in humid regions.

To ensure data reliability, a stringent quality control procedure was also applied to the selected sites, following Zhang et al. (2019). Only sites meeting the following criteria were retained: (i) an energy balance closure ratio exceeding 0.75, (ii) a spatially homogeneous flux footprint within a 500 m radius, and (iii) a continuous observational record longer than one year. After filtering, the final dataset comprised 208 sites (∼ 1400 site-years for ET), representing a substantial expansion compared to the 95 sites used in the earlier version 2.0 (Zhang et al., 2019).

Plant functional type (PFT) classifications were harmonized between the International Geosphere–Biosphere Programme (IGBP) and European Space Agency Climate Change Initiative (ESA CCI) schemes, which are intended to be used for the versions “a” and “b”, respectively. While IGBP classifications are directly provided by the flux networks as high-accuracy references, the ESA CCI classes were mapped from IGBP categories through a one-to-one correspondence at the selected sites. For the IGBP scheme, model parameters were calibrated for 12 vegetated PFTs using site-level observations (Fig. 1a). During this harmonization process, the four forest classes (ENF, EBF, DNF, and DBF) and other key types, including grasslands (GRA) and croplands (CRO), were mapped to their corresponding ESA CCI classes (e.g., GRASS-NAT and GRASS-MAN) to ensure a perfect match (Fig. 1b). To ensure robust cross-validation, the nine open and closed shrublands (CSH/OSH) were merged into a single shrubland (SH/SHRUBS) class. Furthermore, the 30 m 2015 Landsat-derived Global Rainfed and Irrigated-cropland Product (LGRIP30) (Teluguntla et al., 2023) was used to further split CRO and GRASS-MAN into rainfed (RFD) and irrigated (IRR) types. Finally, all vegetated PFTs were parameterized with EC data, while the four non-vegetated classes (water, ice, barren land, and urban areas) were assigned static values.

Figure 1 Global distribution of eddy covariance flux towers for model calibration and validation.

Consequently, we retained a subset of 172 sites (82.7 % of the 208 IGBP sites) that exhibited perfect one-to-one alignment between IGBP and CCI classification, enabling both schemes to share identical parameter sets (Table 2). This rigorous selection was necessary because initial investigations revealed that some ESA CCI labels did not align with the strict physical definitions of the classes. Specifically, the 36 sites where mixed forests (MF) were labeled as specific forest types, or where woody savannas (WSA), savannas (SAV), and wetlands (WET) were labeled as grasslands, were excluded to ensure classification integrity. During modelling, the CCI classification uses fractions of other PFTs to represent these IGBP PFTs.

Model calibration was conducted independently for each PFT, and cross-validation was performed using a leave-one-out strategy within each PFT (Zhang et al., 2019; He et al., 2022). This means that for each PFT, the parameters were optimized using all sites except one, and the calibrated model was then used to predict the omitted site. This process was iterated so that each site, in turn, was left out for prediction, ensuring that all sites within a PFT were ultimately predicted. Compared to version 2.0, PML-V2.2 employs an enhanced multi-objective error metric (F) (Viney et al., 2009; Zhang et al., 2016b), which places a heightened emphasis on bias control to improve simulation accuracy:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}F& =\mathrm{2}-({\mathrm{NSE}}_{\mathrm{ET}}+{\mathrm{NSE}}_{\mathrm{GPP}})+\mathrm{5}\left|\ln \right(\mathrm{1}+{\mathrm{Bias}}_{\mathrm{ET}}){|}^{\mathrm{2.5}}\\ & +\mathrm{5}\left|\ln \right(\mathrm{1}+{\mathrm{Bias}}_{\mathrm{GPP}}){|}^{\mathrm{2.5}}\end{array}\end{array}$$

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left|y_{\mathrm{sim},i}-y_{\mathrm{obs},i}\right|}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left|y_{\mathrm{obs},i}-\stackrel{\mathrm{‾}}{y_{\mathrm{obs}}}\right|}^{\mathrm{2}}}}\text{(9)}& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\sum _{i=\mathrm{1}}^n{y}_{\mathrm{sim},i}-\sum _{i=\mathrm{1}}^n{y}_{\mathrm{obs},i}}{\sum _{i=\mathrm{1}}^n{y}_{\mathrm{obs},i}}}\end{array}$$

where NSE is the Nash–Sutcliffe efficiency of the 8 d ET or GPP. Bias is the model bias. y_obs,i and y_sim,i are the observed and simulated values at time step i. $\stackrel{\mathrm{‾}}{y_{\mathrm{obs}}}$ is the mean of all observations. This objective function jointly optimizes the NSE and Bias for both ET and GPP. By simultaneously constraining temporal dynamics (via NSE) and long-term water and carbon balance (via Bias), this formulation improves the robustness and transferability of model parameters for global-scale applications.

To provide greater transparency and evaluate parameter behavior across different PFTs, a variance-based global sensitivity analysis was conducted using the Sobol method. We calculated the Sobol sensitivity indices for all calibrated parameters to quantify their relative contributions to the variance of simulated ET, GPP, and their underlying components. As illustrated in Fig. A1 in the Appendix, the parameters β, η, and V_m strongly influence photosynthesis and GPP. Parameter m, which regulates the stomatal conductance process, dominates E_c and WUE, making it the most influential parameter for total ET. Additionally, soil evaporation (E_s) is primarily driven by k_A, which controls canopy-soil energy partitioning, while interception evaporation (E_i) depends on the canopy interception capacity (S_l) and evaporation rate (F₀). This sensitivity analysis enhances the physical interpretability of the PML model and further justifies our PFT-specific parameterization strategy.

2.3 Model versions and consolidation

To meet diverse research needs ranging from high-resolution monitoring to long-term climate attribution, this study develops three primary versions of the PML-V2 model. PML-V2.2a and PML-V2.2b serve as the foundational products for the MODIS/VIIRS and AVHRR satellite eras, respectively. These versions provide the basis for the final consolidated dataset, PML-V2.2c, which bridges the AVHRR, MODIS, and VIIRS records to offer a seamless global dataset from 1982 to the near-present. While the core physics of the PML-V2.2 model remain consistent across all versions, these datasets are distinguished by their specific satellite sensors, spatial resolutions, and Plant Functional Type (PFT) inputs to ensure maximum utility across different spatiotemporal scales (Fig. 2).

Figure 2 Methodological framework and temporal coverage of the PML-V2.2 dataset.

A critical enhancement in this framework is the adoption of MSWEP V2.8 and MSWX-Past as meteorological forcing (Beck et al., 2019, 2022), replacing the Global Land Data Assimilation System (GLDAS) V2.1 used in previous versions (Beaudoing and Rodell, 2020). The primary driver for this transition was the temporal limitation of GLDAS V2.1, which is restricted to the post-2000 era. In contrast, MSWEP and MSWX offer a seamless, high-accuracy record from 1982 to 2025, enabling the consistent detection of decadal eco-hydrological trends. Beyond temporal coverage, these datasets utilize a unique multi-source fusion strategy that optimally combines gauge, satellite, and reanalysis data. Unlike purely model-based products, this integration significantly reduces biases in precipitation and meteorological variables, providing a climatic baseline with superior accuracy. Furthermore, the native 0.1° resolution of MSWEP/MSWX minimizes the footprint mismatch caused by the coarser 0.25° GLDAS, ensuring better spatial alignment with satellite-derived inputs. For implementation, all forcing variables (e.g., precipitation, radiation, temperature, vapor pressure, wind speed) were spatially aligned using bilinear interpolation (Fig. 2). Additionally, the global monthly atmospheric CO₂ concentration (C_a) data were sourced from Mauna Loa Observatory administered by the National Oceanic and Atmospheric Administration (NOAA). These monthly concentrations were applied as static values for the 8 d and half-monthly simulation steps within each corresponding month. With these consistent climate forcing and parameter sets, the three datasets are produced as follows.

PML-V2.2a serves as the primary moderate-resolution (500 m, 8 d) version of the framework from 2000 to the present. While it is predominantly driven by MODIS Collection 6.1 inputs (including MOD15A2H LAI, MOD43A3 albedo, MOD11A2 emissivity, and MOD12Q1 IGBP PFT data) (Myneni et al., 2021; Schaaf and Wang, 2021; Wan et al., 2021; Friedl and Sulla-Menashe, 2022), we also incorporated subsidiary inputs from the corresponding VIIRS V002 products (including VNP15A2H LAI, VNP43MA3 albedo, VNP21A1D emissivity) (Myneni, 2023; Schaaf et al., 2025; Hulley and Hook, 2025). Integrating VIIRS ensures the future continuity of the dataset as the aging MODIS sensors are gradually decommissioned. While it maintains consistency with the physical structure of earlier versions, PML-V2.2a incorporates the latest model parameters calibrated against an expanded dataset of 208 flux tower sites. A key technical advancement is the integration of the weighted Whittaker LAI smoothing algorithm V0.1.6 (Kong et al., 2019). This upgraded version enhances outlier detection and specifically avoids the common pitfall of over-smoothing. Consequently, it accurately preserves the signals of bimodal seasonality in double-cropping croplands and mitigates the chronic underestimation of LAI in tropical regions caused by frequent cloud cover and rainfall (He et al., 2022). During the stable overlapping period of 2013–2022, the MODIS and VIIRS LAI and albedo exhibited highly consistent relative changes (NSE = ∼ 0.80, Figs. A2 and A3), which provides robust support for the consolidation of these two datasets. However, it is important to note that severe sensor degradation and orbital drift have been observed in the post-2022 MODIS records (Twedt et al., 2023; Wu et al., 2025). These uncompensated radiometric anomalies are known to propagate systematic, non-physical biases into downstream products such as cloud droplet number concentration (Liu et al., 2026). Specifically, these uncompensated shifts caused artificial declines in LAI retrievals particularly over tropical regions (Fig. A2d), even though their current impact on albedo retrievals seems marginal. Consequently, the MODIS-based version a ceases updates after the year 2024. Future extensions will rely solely on VIIRS until the MODIS Bidirectional Reflectance Distribution Function (BRDF) calibration issues are officially resolved. These sensor transitions and algorithmic refinements contribute to a theoretically higher accuracy and better representation of vegetation heterogeneity.

PML-V2.2b is the AVHRR-based version designed for historical eco-hydrological analysis, providing half-month 0.1° outputs from 1982 to 2020. This version utilizes the GIMMS LAI4g dataset (Cao et al., 2023; Zhu et al., 2013) as its primary vegetation input. Unlike its predecessors, GIMMS LAI4g employs a deep-learning-based consolidation strategy that effectively eliminates satellite orbital drift and sensor degradation, ensuring consistency across the pre- and post-2000 eras. For model parameterization, these vegetation dynamics are paired with GLASS albedo and emissivity products (Liu et al., 2013a; Cheng et al., 2016; Liang et al., 2021) and ESA CCI PFT data (Harper et al., 2023). To match the model grid, the high-resolution land cover data were aggregated to 0.1° resolution by calculating the fractional coverage of five major PFTs within each grid. Although restricted by the coarser spatial resolution of AVHRR compared to MODIS, PML-V2.2b provides the necessary decadal continuity for studying terrestrial water and carbon cycles over the past 39 years. By employing the same climate forcing, model and parameter sets as version a, it ensures that the foundational diagnostic mechanisms remain consistent across different satellite epochs.

PML-V2.2c is the final consolidated dataset (half-month 0.1°) generated through a rigorous, pixel-scale bi-directional consolidation that bridges three satellite eras: AVHRR, MODIS, and VIIRS, providing a continuous record from 1982 to the present (Fig. 2). This consolidation was performed in two successive stages. First, to bridge the pre- and post-2000 eras, we treated the MODIS-based estimates (version a) as the primary benchmark and harmonized the historical AVHRR-based estimates (version b) using the 2001–2003 overlapping period. Second, to extend the dataset reliably into the post-MODIS era, a similar bidirectional consolidation was applied to link the MODIS and VIIRS records, utilizing their highly stable 2020–2022 overlapping window. As detailed in the description of version a, this specific period was carefully selected to avoid the post-2022 MODIS sensor drift, which renders subsequent years unsuitable for benchmark calibration. For most pixels, the primary variables from the subsidiary sensors (AVHRR or VIIRS) were scaled to match the mean values of the benchmark sensor to ensure statistical alignment. In extreme scenarios, such as desert or ice-covered regions or instances where the benchmark detected vegetation (> 0) but the subsidiary sensor remained at zero (precluding standard scaling), a reverse scaling or masking to zero was applied to maintain absolute consistency. This harmonization effectively eliminates systematic jumps caused by sensor transitions and PFT classification differences, ensuring that the 44-year record preserves both accurate mean states and genuine long-term trends.

2.4 Top-down water balance validation

We use top-down water balance ET (ET_WB) estimates to independently validate our new PML ET estimates. At the annual scale, water balance ET is estimated by solving the basin-scale water balance equation using precipitation (P), streamflow (Q), and changes in terrestrial water storage ($\mathrm{d}S/\mathrm{d}t$):

$$\begin{array}{}\text{(10)}& {\mathrm{ET}}_{\mathrm{WB}}=P-Q-\mathrm{d}S/\mathrm{d}t\end{array}$$

Figure 3 Multi-source water balance ET and its associated uncertainty.

To account for inherent uncertainties in individual components of the water balance, the validation integrates multiple observation-based data sources. Similar to Ma et al. (2024), precipitation is obtained from several widely used products, including the MSWEP V2.8 (Beck et al., 2019), the Global Precipitation Climatology Project (GPCP) V3.3 (Huffman, 2024), the Global Precipitation Climatology Centre dataset (GPCC) version 2022 (Schneider et al., 2022), and the Climatic Research Unit Time Series (CRU-TS) V4.09 (Harris et al., 2020). Streamflow data are compiled from in situ observations provided by the Global Runoff Data Centre (GRDC), the United States Geological Survey (USGS), the Australian Bureau of Meteorology (BOM), and the Ministry of Water Resources (MWR) of China. Terrestrial water storage anomalies (TWSA) are obtained from the Gravity Recovery and Climate Experiment and its Follow-On mission (GRACE/FO), including solutions from the Jet Propulsion Laboratory (JPL) (Save et al., 2016; Save, 2022), the Center for Space Research (CSR) (Watkins et al., 2015; Wiese et al., 2016, 2018), and the Goddard Space Flight Center (GSFC) (Loomis et al., 2019) during 2003–2020. The 11 gap months during 2017–2018 and other isolated gaps are filled by two independent methods, including a climate adjustment scheme using singular spectrum analysis (climSSA) (Zhang et al., 2025) and a transformer-based deep learning approach (Wang and Zhang, 2024). Following Ma et al. (2024), we use the GRACE reconstruction dataset (GRACE-REC) to extend TWSA records during 1982–2002 (Humphrey and Gudmundsson, 2019). $\mathrm{d}S/\mathrm{d}t$ is computed as a mean of the differences between consecutive Januarys and Decembers at an annual scale.

By combining these independent data sources, the analysis explicitly accounts for observational uncertainties, which can vary substantially across regions. Basin-scale water-balance ET estimates (ET_WB) are obtained using an arithmetic mean method, which allows for the estimation of a consensus ET signal while minimizing the influence of individual data-source errors. The associated uncertainty (σ) is quantified as the standard deviation across all feasible combinations of P, Q, and $\mathrm{d}S/\mathrm{d}t$ datasets.

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_{\mathrm{ET}}=\sqrt{{{\mathit{\sigma }}_P}^{\mathrm{2}}+{{\mathit{\sigma }}_Q}^{\mathrm{2}}+{{\mathit{\sigma }}_{\mathrm{d}S/\mathrm{d}t}}^{\mathrm{2}}}\end{array}$$

where σ are computed as standard deviation among different data sources and processing methods. Following Ma et al. (2024), a relative uncertainty of 5 % is assumed for Q.

As shown in Fig. 3, this validation framework is applied across multiple spatial and temporal scales. Specifically, it includes an analysis of 56 large river basins (Ma et al., 2024), covering approximately 29 % of the global land surface for the period 1982–2016, and an expanded assessment of 152 large basins, representing about 42 % of global land area, for the period 2003–2020 (streamflow data extended from Zhang et al., 2023).

2.5 Statistical analysis

The predictive accuracy of the model is evaluated using a suite of performance metrics, including NSE, Bias, the Pearson correlation coefficient (R), and Root Mean Square Error (RMSE). These metrics are selected to provide a balanced assessment, as they compensate for each other's specific limitations. While R measures the strength of the linear relationship and NSE evaluates the overall fit relative to the observed mean, both can be sensitive to extreme outliers or peak values. To address this, RMSE is used to quantify the absolute error in the same units as the data, and Bias is employed to detect any systematic over- or under-estimation (volume error). Collectively, this multi-metric approach ensures that the model is validated for its ability to capture both the timing and the magnitude of the observed data.

To interpret long-term changes, the non-parametric Mann-Kendall (MK) test and Sen's Slope estimator are applied to quantify the significance and magnitude of trends. These methods are robust against outliers and non-normal distributions, ensuring that the results represent long-term trends rather than impact of anomalies.

Figure 4 Model calibration and cross-validation result under IGBP classification of the total 208 sites.

Figure 5 Calibration and cross-validation result over different PFTs under IGBP classification. Note that 8 of the 12 IGBP PFTs have the same parameters but with different names under CCI classification (Table 2).

3 Results

3.1 Model calibration and validation

The PML model demonstrates high proficiency in estimating 8 d ET and GPP across the 208 eddy covariance sites (Fig. 4). Overall, the model performs robustly in capturing flux variability, indicated by ET statistical metrics: NSE = 0.72, R = 0.86, RMSE = 0.64 mm d⁻¹, Bias = −1.7 %, and GPP statistical metrics: NSE = 0.76, R = 0.87, RMSE = 1.97 g C m⁻² d⁻¹, Bias = 0.2 %. The corresponding Kling-Gupta Efficiency (KGE) values for ET and GPP are 0.86 and 0.82. Similarly, it has robust model performance in cross-validation for both ET and GPP, with only slight degradation in overall NSE and other metrics. The validation of WUE (GPP per unit of ET consumption) at the site scale demonstrates good model performance across both multi-year and annual scales, with cross-validation NSE values of 0.52 and 0.35, respectively, and absolute biases consistently maintained within 5 % (Fig. A4). The clustering of points along the 1 : 1 line across both calibration and validation phases demonstrates the model's robustness in capturing flux variability across diverse site conditions.

Building upon these global site-level assessments, the model's performance remains highly stable when analyzed across different PFTs (Fig. 5). For ET performance, accuracy is consistent across most biomes, with PFT-specific NSE values typically ranging from 0.45 to 0.78 and R values between 0.74 and 0.92. PFTs such as grasslands (GRA), deciduous broadleaf forest (DBF), and wetlands (WET) exhibit the highest precision in water flux estimation, with RMSE generally staying between 0.5 and 0.95 mm d⁻¹. Bias remains within ± 5 % for the majority of categories.

Regarding GPP performance, the model achieves even higher metrics across the PFT spectrum, with NSE reaching up to 0.88 and R values exceeding 0.90 for categories such as WET and deciduous needleleaf forest (DNF). While agricultural types – specifically irrigated (CRO-IRR) and rainfed (CRO-RFD) croplands – display higher absolute errors (RMSE) compared to natural vegetation, their biases remain tightly clustered around zero. This indicates that the model effectively captures the magnitude of carbon productivity with minimal systematic error across different vegetation structures.

Overall, the model demonstrates high robustness and regionalization capability, with high consistency between calibration and cross-validation metrics across biomes. This robustness is particularly evident in GPP estimates, which show greater stability than ET during validation. Despite this reliability, some PFT-specific sensitivities persist: EBF remains challenging to simulate due to low seasonal variability limiting the signal-to-noise ratio, while DNF shows a sharper cross-validation drop in ET NSE likely due to the limited representation (only four sites). Nonetheless, the high consistency across most PFTs confirms the model's reliability for capturing complex eco-hydrological processes under varying conditions.

3.2 Improved parameterization of croplands

The implementation of a split parameterization scheme for rainfed (CRO-RFD) and irrigated (CRO-IRR) croplands led to significant performance improvements (Fig. 6). For rainfed systems, the model showed a notable increase in the NSE by 0.07 and a reduction in RMSE by 0.07 mm d⁻¹, effectively mitigating systematic biases with a decrease in ET Bias of 10.0 %. While the irrigated ET precision metrics remained relatively stable, the Bias shifted by 7.4 % toward the zero-baseline, suggesting a more accurate representation of water availability and management in irrigated agricultural lands compared to the traditional mixed approach.

Figure 6 Improvement in cross-validation by split parameterization across rainfed and irrigated croplands compared with mixed parameterization. The evaluation was conducted across 29 flux sites, including 14 irrigated and 15 rainfed sites.

Substantial enhancements were also observed in GPP simulations, where the split scheme addressed systematic errors more effectively than the mixed model. GPP RMSE decreased by 0.17 and 0.24 g C m⁻² d⁻¹ for RFD and IRR, respectively, while Bias was substantially reduced by 18.9 % and 13.6 %. Overall, this parameterization scheme reduced the magnitude of systematic errors across all agricultural lands by averages of 8.7 % for ET and 16.2 % for GPP. This suggests that the conceptual simplification of treating irrigated and rainfed crops as distinct entities enables a more physically consistent representation of agricultural cycles. By separately capturing high-productivity irrigated regimes and moisture-limited rainfed growth, this approach effectively achieves a better closure of both water and carbon balances.

Figure 7 Comparison of simulated ET from three versions of PML-V2.2 datasets with water-balance-derived ET (ET_WB) averaged for two basin datasets over different historical periods (1982–2016 and 2003–2020). Note that ET_WBdata are shown in Fig. 3. The four evaluation metrics here are computed without area weights, while an additional figure considering area weights is shown in Fig. A5.

3.3 Basin-scale water balance validation

Basin-scale validation against water-balance-derived ET (ET_WB) confirms the robust performance of all three PML-V2.2 products across different historical periods (Fig. 7). Across both the 2003–2020 (152 basins) and 1982–2016 (56 basins) evaluation windows, the products consistently yield high accuracy with NSE values of 0.89–0.91 and R ≥ 0.95. Subtle discrepancies exist between versions a and b, primarily driven by differences in LAI, albedo, and PFT inputs. While the point-wise statistics indicate a slight underestimation for most products, the primary biases for all versions remain well within ± 6 %. Notably, supplementary analysis (Fig. A5) reveals a marginal bias when results are aggregated by basin area, highlighting the relative effects inherent in basin-scale water balance validation. While ET_WB estimates are not entirely bias-free (Fig. 3), the validation results remain within a reasonable and acceptable range, particularly given the robust NSE values across basins.

3.4 Consistency of consolidation

To evaluate the consistency and reliability of the data integration process, we first assessed the transition between the recent satellite epochs. Supported by highly consistent relative changes in LAI and Albedo between 2020–2022 and 2013–2015 (NSE = 0.80 and 0.79, Figs. A2 and A3), the corresponding shifts in ET and GPP based on MODIS and VIIRS exhibit high spatial agreement (NSE = 0.99 and 0.94, Fig. A6). Consequently, the primary challenge for long-term consistency lies in bridging the AVHRR-based product (PML-V2.2b) with the MODIS-based version (PML-V2.2a).

Figure 8 Consistency of interannual variability and decadal trends between AVHRR-based PML-V2.2b and MODIS-based PML-V2.2a during the overlapping period of 2001–2020. Panels (a) and (b) compare the annual anomalies for ET and GPP, respectively, indicating the agreement in interannual variability. Panels (c) and (d) compare the decadal trends (see Fig. A7 for the corresponding spatial trend maps). The metrics are calculated based on aggregated grids at a 0.5° resolution to demonstrate the robustness of the consolidation process.

To address this, their decadal trends and annual anomalies were compared for the overlapping period 2001–2020 (Fig. 8). The annual anomalies exhibit excellent agreement, with high consistency for ET (NSE = 0.90) and GPP (NSE = 0.79), indicating that the reconstructed dataset robustly captures interannual hydro-climatic fluctuations. The decadal trends of both ET and GPP also show a high degree of consistency across the global land surface. For ET, the two versions exhibit a strong correlation (R = 0.92) and a high NSE = 0.83 (slightly lower than NSE = 0.90 for ET anomalies due to inherent uncertainty in trend estimation). Similarly, GPP trends demonstrate significant agreement with an R value of 0.81 and an NSE of 0.64. These statistical results are further corroborated by the spatially coherent global patterns shown in Fig. A7. While slight regional variations exist due to remote-sensing forcing data inputs, the comparison of their mean climatology highlights that the PML-V2.2b product yields higher ET and GPP estimates across most non-tropical rainforest regions, primarily driven by the inherently higher LAI values in these areas (Fig. A8). Nevertheless, the overall alignment along the 1 : 1 line confirms that the consolidation process preserves the core long-term hydro-climatic signals of the original datasets.

Figure 9 Global patterns of terrestrial evapotranspiration (ET) and its components of the PML-V2.2c dataset during 1982–2025.

Figure 10 Global patterns of terrestrial gross primary production (GPP) and water use efficiency (WUE) of the PML-V2.2c dataset during 1982–2025. WUE is computed by GPP/ET.

3.5 Global patterns and trends during 1982–2025

The global spatial distribution of long-term mean ET during 1982–2025 exhibits strong latitudinal gradients, reflecting the fundamental role of energy and water availability (Fig. 9). The highest ET fluxes, exceeding 1000 mm yr⁻¹, are concentrated in tropical rainforests across the Amazon Basin, Congo Basin, and Southeast Asia. In these humid regions, vegetation transpiration and interception loss are the primary contributors to the total water flux. Conversely, ET values remain minimal in arid and semi-arid regions such as the Sahara and Central Australia, where water limitation severely constrains surface evaporation. The PML-V2.2 product effectively captures these hydro-climatic constraints, showing that transpiration dominates the global ET budget over most vegetated land surfaces. Overall, the PML-V2.2c product estimates global ET of 65.8 × 10³ km³ yr⁻¹ during 1982–2025. This corresponds to a mean flux of 509 mm yr⁻¹ averaged over a land area of 1.292 × 10⁸ km² (excluding inland water bodies, permanent snow and glaciers). The partitioning reveals that 58.2 %, 30.2 %, and 11.6 % of this flux are contributed by transpiration, soil evaporation, and interception evaporation, respectively.

Figure 11 Interannual variations and long-term trends of global water and carbon fluxes of the PML-V2.2c dataset during 1982–2025. All global calculations exclude water bodies to focus on terrestrial surfaces. The blue solid lines indicate the linear trends, where α denotes the Sen's Slope estimates (10³ km³ yr⁻², PgC yr⁻², or gC mm⁻¹ H₂O yr⁻¹ for different variables) and the p-value represents the significance level determined by the MK test.

Terrestrial productivity and its coupling with water consumption are characterized by the spatial distributions of GPP and WUE (Fig. 10). The global total GPP in the PML-V2.2c dataset is approximately 143.4 PgC yr⁻¹, peaking in tropical biomes (above 2500 gC m⁻² yr⁻¹), where photosynthesis is maximized by abundant sunlight and moisture. However, the WUE pattern reveals a different ecological strategy: the highest WUE values are found in temperate forests and certain semi-arid shrublands, indicating high carbon assimilation per unit of water lost. These spatial patterns highlight the regional differences in how ecosystems balance their carbon-water trade-offs, with the PML-V2.2 dataset providing a consistent representation of these global productivity gradients.

The long-term dynamics from 1982 to 2025 reveal an intensifying global water and carbon cycle, characterized by widespread increasing trends (Fig. 11 and Fig. A9). Global GPP has shown a robust and significant increase (α = 0.343 PgC yr⁻², p < 0.01), particularly in greening hotspots like China, India, and Europe. Parallel to carbon gains, global ET also exhibits a significant upward trend (α = 0.019 × 10³ km³ yr⁻², p < 0.05), primarily driven by the increasing transpiration and interception. While most regions show positive trajectories, notable regional heterogeneity exists, with declining ET trends in parts of the western U.S., South America and Southern Africa. Overall, the global-scale consistent rise in GPP, ET, and WUE (Fig. A9) underscores the significant impact of CO₂ fertilization and climate change on terrestrial ecosystem functions over the past 44 years.

4 Discussion

The transition to the PML-V2.2 series represents a significant advancement in model stability and accuracy. By implementing a refined optimization framework and updated PFT schemes, the absolute bias remains within 5 % across diverse biomes, marking a substantial improvement over previous versions where ET bias could reach 10 % and GPP bias could be exceeding 25 % during cross-validation (Fig. 5 in Zhang et al., 2019). A critical improvement is observed in cropland simulations. While remote sensing-based LAI partially captures the reduced water limitation in agricultural areas, our results suggest that differences in stomatal conductance are another factor in determining accuracy. By incorporating two basic irrigation types into the parameterization, the model better reflects the enhanced productivity of managed lands. Despite these gains, substantial potential remains for further refinement. Future iterations should focus on moving beyond static PFT-based parameters toward a scheme that integrates dynamic vegetation attributes, such as leaf traits and hydraulic properties, directly into the parameterization framework to better capture physiological responses to environmental change (Yan et al., 2025).

Figure 12 Consistency in interannual variations and long-term trends in global water and carbon fluxes across different PML-V2 versions. Note that the higher variability in the initial PML-V2.0 is likely attributable to unstable climate forcing from GLDAS V2.1. The 2025 MODIS-based PML-V2.2a estimates (dashed line) are shown to illustrate satellite-drift induced biases but will not be available in the public data repository.

Figure 13 Comparison of global mean and interannual variability in ET and GPP from mainstream diagnostic products listed in Table 1. The mean values are calculated for the period 2003–2016, during which data records are available for most products. White circles at the base of the bars denote products starting from 1982 or earlier, while white squares indicate datasets updated to the near-present of 2025.

The temporal consistency of the PML-V2.2 product from 1982 to 2025 (Fig. 12) demonstrates its robustness for long-term climate studies. This stability is primarily attributed to the use of unified meteorological forcing and parameter sets across different satellite eras. During the overlapping period, the historical AVHRR-based PML-V2.2b yields moderately higher global magnitudes (approximately +2.6 % for ET and +3.1 % for transpiration fraction) compared to the MODIS/VIIRS-based versions. This difference is driven by residual uncertainties in the consolidated historical inputs (specifically GIMMS LAI4g and GLASS albedo) (Cao et al., 2023), resulting in inherently higher LAI and GPP values across non-tropical regions (Fig. A8). Furthermore, the global and regional interannual trends show remarkable coherence (Figs. 12 and A7). To maintain this physical consistency for long-term trend analysis, we deliberately applied the unified parameter set across all epochs (Fig. 2), avoiding the scale-mismatch issues associated with separately calibrating the 0.1° AVHRR model against site-level flux towers with footprint ranging from several to hundreds of meters.

By fully integrating VIIRS to replace the aging MODIS, this framework successfully ensures high-resolution monitoring continuity into the near present and future eras (Yan et al., 2021; Endsley et al., 2025). This sensor transition is critical, as the severely degraded MODIS LAI in 2025 introduced an artificial, sharp decline in both ET and GPP estimates, as clearly captured by the v2.2a (MODIS) anomaly stream in 2025 (dashed blue line in Fig. 12), justifying our decision to discontinue the MODIS-based updates (Figs. A2 and A3). Under the current consolidation structure, PML-V2 can be continuously updated to provide a reliable, long-term record for critical assessments, such as the State of Global Water Resources Report compiled by World Meteorological Organization (2025), where high-fidelity, multi-decadal data are required to monitor global hydro-climatic shifts.

Differences in global means between the current and earlier PML versions are likely driven by the shift from GLDAS 2.1 to MSWEP/MSWX climate forcing (Kim et al., 2021; Xie et al., 2025), alongside updated parameterization, updated LAI data and smoothing algorithms. In version 2.0, the model predicted an aggressive increase in ET and GPP during 2003–2010 but followed by sharp declines. Our analysis reveals that these spurious trends may result from the high volatility in the solar radiation terms of the GLDAS forcing used in that version (Fig. 12). Comparatively, the new global mean GPP and WUE exhibit a continuing upward trend, a more physically plausible result in a greening world. Overall, while both versions have similar range of global patterns and total fluxes, the new version is recommended for users as it covers a longer period and offers higher quality as validated from eddy flux observations and water balance.

Our global ET estimate of 65.8 × 10³ km³ yr⁻¹ during 1982–2025 is slightly lower than recent estimates by Zhang et al. (2026). In their updated assessment of global water cycle variables, they report a mean annual land ET of 69.0 ± 6.2×10³ km³ yr⁻¹ for 1981–2014, excluding Greenland and Antarctica. Our estimate also excludes these regions, but additionally omits evaporation from surface inland water bodies, which is approximately 1.5 ± 0.15 × 10³ km³ yr⁻¹ (Zhao et al., 2022a). When this water body evaporation is included, our adjusted global ET becomes very similar to or 2.5 % less than the mean value reported by Zhang et al. (2026) (Fig. 13). Our updated parameter sets differ from the PML-V2.0 version (Zhang et al., 2019) due to three key refinements: (i) an expanded calibration using 208 global flux sites from the original 95 sites for improved robustness; (ii) a new distinction between irrigated and rainfed croplands that reduced their internal ET biases by 8.7 %; and (iii) the removal of a parameter to enhance model parsimony without losing accuracy. Compared to 17 mainstream diagnostic products during basin-scale water balance validation (Figs. A10 and A11), the PML-V2.2 products consistently demonstrate top-tier performance, ranking as the best or second-best among all products evaluated. Specifically, PML-V2.2 achieves the highest NSE (0.89–0.91) and the lowest RMSE (77.0 to 96.3 mm yr⁻¹).

Our global GPP estimate of 143.4 PgC yr⁻¹ during 1982–2025 is in closer alignment with the evolving consensus in the carbon cycle community. Traditional optical remote sensing-driven products, specifically LUE-based upscaling (∼ 115 PgC yr⁻¹), machine learning-based upscaling (∼ 124 PgC yr⁻¹), and process models (∼ 131 PgC yr⁻¹), tend to converge within a lower range. Our result is significantly higher but aligns more closely with contemporary estimates derived from novel proxies and ensemble assessments, such as the soil respiration constraints (∼ 149 ± 25 PgC yr⁻¹), the ¹⁸O signature of atmospheric CO₂ (150–175 PgC yr⁻¹), and the latest carbonyl sulfide (OCS)-based results by Lai et al. (2024) (157 ± 8.5 PgC yr⁻¹). This position is clearly illustrated in our multi-product intercomparison (Fig. 13c), where PML-V2.2c yields a higher magnitude than 11 of the 15 mainstream diagnostic GPP datasets. Although global GPP magnitude remains highly debated (Tian et al., 2026; Lai et al., 2026), this higher estimate reflects a physically constrained output rather than an unguided inflation. Notably, among the very few diagnostic products that simultaneously provide both water and carbon fluxes, PML-V2.2c is the only one that delivers internally consistent ET and GPP records continuously extending to the near-present of 2025 (Fig. 13).

The long-term trends of PML-V2.2c offer a nuanced view of the global water cycle, showing a more conservative ET growth of approximately 1.25 % during 1982–2025. In terms of linear trends, our estimate of 0.24 mm yr⁻² during 1982–2011 is significantly more moderate compared to the ensemble mean of six diagnostic products (0.66 ± 0.38 mm yr⁻²) (Yang et al., 2023). Similarly, when benchmarked against a broader ensemble of 11 ET products (1982–2016) with a median trend of 0.28 mm yr⁻² (Ma et al., 2021), our estimate of 0.18 mm yr⁻² for the same period remains on the lower end of the spectrum. This global trajectory results from the offsetting dynamics of internal components. This multi-product divergence is evident in Fig. 13. While global ET anomalies exhibit highly scattered trajectories due to their heavy reliance on varying meteorological or radiation forcings, global GPP anomalies display a remarkably coherent upward trend across nearly all datasets, especially during the MODIS era. This high GPP consistency is fundamentally driven by unified satellite-derived greenness inputs (e.g., LAI or NDVI) that carry a robust global greening signal. From PML-V2.2c estimates, the marginal increase in transpiration is likely driven by the interplay between vegetation greening and physiological water saving from rising atmospheric CO₂ (Zhang et al., 2021; Wei et al., 2024; Lesk et al., 2025; Wei et al., 2025). Simultaneously, interception evaporation shows a significant upward trend following LAI increases, while soil evaporation has slightly decreased due to the reduced available energy at the soil surface caused by canopy shading. The magnitude of global GPP remains consistent with other independent data sources (Yang et al., 2022; Lai et al., 2024; Wang et al., 2024), reflecting a strong increase in WUE due to vegetation greening (Fig. 12). While uncertainties remain, the internal consistency of the PML-V2 framework makes these results a valuable tool for multi-model assessments and climate attribution research, providing a stable benchmark for understanding how the terrestrial biosphere regulates water and carbon fluxes.

5 Data availability

The original PML-V2.2a dataset at 8 d 500 m resolution is freely accessible via Google Earth Engine at https://developers.google.com/earth-engine/datasets/publisher/pml_evapotranspiration (last access: 15 July 2026). The aggregated 0.1° versions of PML-V2.2a (MODIS, 2000–2024; VIIRS, 2012–2025), along with the PML-V2.2b (1982–2020) and PML-V2.2c (1982–2025) datasets, are available at the National Tibetan Plateau Data Center (TPDC) under https://doi.org/10.11888/Terre.tpdc.303314 (Xu et al., 2026).

The active PML-V2.2a and V2.2c datasets will be updated annually with a six-month processing latency. Users are encouraged to select the dataset version that best balances spatiotemporal resolution and temporal length based on their needs. To support open science, the PML-V2.2 model code (R and C$++$ language), calibration data for 208 flux sites, and figure data for global patterns and trends are available at Zenodo under https://doi.org/10.5281/zenodo.18385275 (Xu, 2026).

6 Conclusions

The development of PML-V2.2 establishes an internally consistent, 44-year record (from 1982 to the near present of 2025) of the global terrestrial water and carbon cycles. The high statistical and spatial agreement between the MODIS/VIIRS-based (V2.2a) and AVHRR-based (V2.2b) products validates the multi-tier architecture, ensuring that the consolidated PML-V2.2c dataset provides a reliable multi-decadal benchmark for climate attribution. By leveraging an expanded calibration of 208 flux sites and distinguishing between irrigated and rainfed systems, the dataset achieves high fidelity, yielding global mean ET and GPP estimates of 65.8 × 10³ km³ yr⁻¹ (with 58.2 % from transpiration) and 143.4 PgC yr⁻¹, respectively. The successful top-down water-balance validation (NSE = 0.89–0.91) across 152 large basins during 2003–2020 and 56 basins during 1982–2016 further confirms that this framework provides a physically plausible representation of global eco-hydrological states.

Beyond data production, this long-term record offers critical insights into the intensifying global water-carbon cycle. Specifically, global GPP has significantly increased at a rate of 0.343 PgC yr⁻² during 1982–2025 (p<0.01), while global ET exhibits a significant but more constrained rise of 0.019 × 10³ km³ yr⁻² (p<0.05). These diverging trends highlight a major eco-hydrological tradeoff, where the increased transpiration driven by vegetation greening is largely offset by the physiological water-saving effect from rising atmospheric CO₂.

With the successful integration of VIIRS securing data continuity beyond the aging MODIS, the consolidation strategies employed in PML-V2.2 provide a robust template for mitigating sensor degradation. Building on this stable multidecadal foundation, future refinements will focus on integrating dynamic vegetation traits and improved water stress mechanisms, while leveraging next-generation remote sensing data to deliver these physiological updates at higher spatial resolutions. As an evolving, near-present, and high-accuracy record, PML-V2.2 is positioned to support critical international assessments, such as the WMO State of Global Water Resources reports, providing the foundational data continuity required for climate adaptation, trend attribution, and global water management.

Appendix A

Figure A1 Global sensitivity analysis of PML-V2 model parameters across different Plant Functional Types (PFTs) and all sites combined (ALL) using the Sobol index. For each specific PFT and the ALL category, 5000 parameter sets were generated using Latin Hypercube Sampling (LHS). To maintain visual clarity and emphasize the dominant drivers, numerical ranks (where 1 represents the most sensitive parameter) are displayed only for parameters with a Sobol index greater than 0.05. Meaning of all parameters can be found in Table 2.

Figure A2 Divergence between MODIS and VIIRS LAI records highlighting post-2022 MODIS sensor degradation. Changes are computed between means of 2020–2022 and 2013–2015, as large uncertainties exist in trend estimates over short periods.

Figure A3 Consistency between MODIS and VIIRS albedo records even after post-2022 MODIS sensor degradation.

Figure A4 Calibration and cross-validation of WUE at annual scale. Annual WUE is derived from ET and GPP averages with >80 % 8 d availability, while multi-year WUE necessitates at least five years of valid data.

Figure A5 Comparison of simulated ET with water-balance-derived ET averaged for two basin datasets over different historical periods (1982–2016 and 2003–2020). Note that compared with Fig. 7, ET considers basin area coverage and is expressed as km³ yr⁻¹.

Figure A6 Spatial consistency of ET and GPP changes between MODIS and VIIRS-based PML-V2.2a datasets. Changes are computed between means of 2020–2022 and 2013–2015, as large uncertainties exist in trend estimates over short periods.

Figure A7 Spatial consistency of ET and GPP trends between the baseline MODIS-based PML-V2.2a and the consolidated AVHRR-based PML-V2.2b during the overlapping period (2001–2020). Panels (a) and (b) show the spatial distribution of ET trends from V2.2a and V2.2b, respectively. Panels (c) and (d) show the comparison for GPP trends. The highly similar spatial patterns verify that PML-V2.2b successfully captures the climate-driven trends observed in the MODIS-based PML-V2.2a. Scatter plots of trends at grid scale are provided in Fig. 8.

Figure A8 Differences between the PML-V2.2a and V2.2b products are mainly attributed to the discrepancies in LAI between MODIS and AVHRR (GIMMS LAI4g, though consolidated with MODIS). Across most non-tropical rainforest regions, the GIMMS LAI4g dataset exhibits higher LAI values, consequently leading to higher ET and GPP estimates. Despite these differences in mean magnitudes, their long-term trends remain highly consistent (see Fig. A7).

Figure A9 Global trends in ET, GPP, and WUE during 1982–2025 from the consolidated PML-V2.2c dataset.

Figure A10 Comparison of simulated ET with water-balance-derived ET averaged for 152 large basins over 2003–2020. To ensure a fair comparison, minor spatial gaps over barren land and water bodies in some products were gap-filled using PML data.

Figure A11 Comparison of simulated ET with water-balance-derived ET averaged for 56 large basins over 1982–2016. To ensure a fair comparison, minor spatial gaps over barren land and water bodies in some products were gap-filled using PML data.

Table A1 Data links of summarised diagnostic products in Table 1.