A daily and 500 m coupled evapotranspiration and gross primary production product across China during 2000–2020

Accurate high-resolution actual evapotranspiration (ET) and gross primary production (GPP) information is essential for understanding the large-scale water and carbon dynamics. However, substantial uncertainties exist in the current ET and GPP datasets in China because of insufficient local ground measurements used for model constraint. This study utilizes a water–carbon coupled model, Penman–Monteith–Leuning Version 2 (PML-V2), to estimate 500 m ET and GPP at a daily scale. The parameters of PML-V2(China) were well calibrated against observations of 26 eddy covariance flux towers across nine plant functional types in China, indicated by a Nash–Sutcliffe efficiency (NSE) of 0.75 and a root mean square error (RMSE) of 0.69 mm d⁻¹ for daily ET, respectively, and a NSE of 0.82 and a RMSE of 1.71 g C m⁻² d⁻¹ for daily GPP. The model estimates get a small Bias of 6.28 % and a high NSE of 0.82 against water-balance annual ET estimates across 10 major river basins in China. Further evaluations suggest that the newly developed product is better than other typical products (MOD16A2, SEBAL, GLEAM, MOD17A2H, VPM, and EC-LUE) in estimating both ET and GPP. Moreover, PML-V2(China) accurately monitors the intra-annual variations in ET and GPP in the croplands with a dual-cropping system. The new data showed that, during 2001–2018, the annual GPP and water use efficiency experienced a significant (p<0.001) increase (8.99 g C m⁻² yr⁻² and 0.02 g C mm⁻¹ H₂O yr⁻¹, respectively), but annual ET showed a non-significant (p>0.05) increase (0.43 mm yr⁻²). This indicates that vegetation in China exhibits a huge potential for carbon sequestration with little cost in water resources. The PML-V2(China) product provides a great opportunity for academic communities and various agencies for scientific studies and applications, freely available at https://doi.org/10.11888/Terre.tpdc.272389 (Zhang and He, 2022).

1 Introduction

Terrestrial evapotranspiration (ET) and photosynthesis (or gross primary productivity, GPP) are indispensable processes in hydrological and carbon cycles at global and regional scales. Forming the second-largest water flux after precipitation in the terrestrial hydrological cycle, ET is the sum of plant transpiration (E_c), evaporation from the soil (E_s), and canopy evaporation from precipitation interception (E_i). Photosynthesis rate, tightly coupled to E_c by leaf stomata, is a key indicator of plant growth and provides food and fiber for human society.

In recent decades, numerous studies have been carried out to map ET at regional, continental, and global scales. Particular credit goes to remote sensing (RS)-based models that provide diagnostic ET estimates with a relatively high spatiotemporal continuity and reasonable biophysical significance. For example, Miralles et al. (2011b) and Martens et al. (2017) developed a global and daily ET product employing the Global Land Evaporation Amsterdam Model 3.0a (hereafter GLEAM) based on the Priestley–Taylor (P–T) equation. Although GLEAM does a good job in temporal resolution, it has a coarse spatial resolution of 0.25^∘. Another widely used semi-empirical formula is the Penman–Monteith (PM) equation. Mu et al. (2011) generated the MOD16A2 ET with 500 m and 8 d resolutions based on the PM equation as one global product of the Moderate-Resolution Imaging Spectroradiometer (MODIS). Leuning et al. (2008) developed the Penman–Monteith–Leuning (PML) model that described the physical characteristics of canopy–soil water loss by improving the surface conductance (G_s) formulations, and Zhang et al. (2010) estimated PML-based ET at 0.05^∘ and 8 d resolution across the Australian continent. Cheng et al. (2021) produced a 1 km and daily ET dataset across China by the Surface Energy Balance Algorithm for Land (SEBAL) but only evaluated the SEBAL product at eight eddy covariance (EC) sites of three land cover types.

Apart from the methods and products for the ET estimated above, numerous approaches have been used to estimate GPP, such as the enzyme kinetic (process-based) models (Houborg et al., 2013; Arain et al., 2006; Grant et al., 2005; Hanson et al., 2004; Medvigy et al., 2009), the light use efficiency (LUE) principle (Liu et al., 2003; Turner et al., 2003; Yuan et al., 2007; Running et al., 2015; Zhang et al., 2017a; Zheng et al., 2020), statistical methods (Potter et al., 1993; Hilker et al., 2008; Z. Zhang et al., 2020), and machine learning methods (Wolanin et al., 2019; Joiner and Yoshida, 2020; Huang et al., 2021). Among them, the LUE principle is well known because of its simple structure, strong portability, and relatively high temporal cover of inputs. Running et al. (2015) recently updated the global GPP product of MODIS (hereafter MOD17A2H) at 500 m and 8 d resolutions using the LUE principle. Zhang et al. (2017b, 2021) also mapped global GPP for 2000–2019, dubbed the Vegetation Photosynthesis Model GPP V20 (hereafter VPM), with the same spatiotemporal resolution as MOD17A2H using an improved LUE model adding the energy absorbed by chlorophyll. Zheng et al. (2020) generated a global GPP dataset at 0.05^∘ and 8 d intervals by a revised LUE model (hereafter EC-LUE) integrating the atmospheric CO₂ concentration.

Although significant efforts have been put into estimating ET and GPP, there are barely any coupled products available in China which meet the requirement of high temporal (≤ 1 d) and high spatial (≤ 500 m) resolutions simultaneously that is necessary to detect variations of the eco-hydrological cycle in diverse and large areas for the long term precisely (Table 1). For instance, products with low temporal resolutions are erratic in detecting subtle seasonal changes in areas seriously affected by human activities and in arid regions, such as irrigated farmland with a dry climate (Bodner et al., 2015) and an evergreen broadleaf Mediterranean forest during severe summer drought (Liu et al., 2015). On the other hand, products with high temporal resolutions like GLEAM can monitor the diurnal variability of ET, but their low spatial resolutions limit their effectiveness in fine-scale environment applications (Gevaert and García-Haro, 2015).

Table 1 Summary of typical ET and GPP products with high temporal and/or high spatial resolutions.

Secondly, the phenomenon of ignoring the water–carbon coupling process frequently appearing in the existing products has brought systematic errors. The photosynthesis and transpiration are coupled by the plant stomatal control on both water and carbon exchange between the land ecosystem and the atmosphere (Xiao et al., 2013; Y. Zhang et al., 2019). As indicated in Table 1, MODIS ET and MODIS GPP products are independent of each other and cannot ensure similar biophysical characteristics of vegetation in the same place. Furthermore, using ET and GPP from different products can lead to large uncertainty in analyzing the interaction between ET and GPP, such as the water use efficiency of an ecosystem (WUE is the ratio of GPP to ET). In that case, it is necessary to build a coupled ET and GPP model considering the water–carbon coupling process. Y. Zhang et al. (2019) developed the second generation of PML (i.e., the PML-V2 model) that estimates G_s using a water–carbon coupled model and mapped global ET and GPP at 500 m and 8 d resolutions in 2002–2017.

More importantly, previous studies utilized sparse ground observations in China (as shown in Table 1) that covered few terrestrial ecosystems for model calibration and validation, resulting in improper input parameters and making it difficult to obtain more reliable estimates of ET and GPP in diverse land cover types (Heinsch et al., 2006; Anon, 2013). Although the EC flux sites have provided consecutive measurements of water and carbon fluxes since the early 1990s (Xiao et al., 2013; Wofsy et al., 1993; Baldocchi et al., 2001), the EC-observed data in China remain much sparser than those in North America and Europe, and most of them are not publicly available, impeding a national-scale constraint of RS-based models for improving ET and GPP estimates (Villarreal and Vargas, 2021; Chu et al., 2017). For instance, GLEAM, which only employed eight EC sites over China, overestimates ET on a large scale, especially for evergreen needleleaf forest, evergreen broadleaf forest, and mixed forest (Li et al., 2018). Therefore, the uncertainty in estimating ET and GPP is large, and it requires sufficient EC flux sites to calibrate and validate ET and GPP models for better local and regional applications.

In addition, understanding the spatial and temporal patterns of ET and GPP is particularly important for China, the largest contributor to the absolute growth of greenhouse gases that has directly induced global warming over the past decade (Minx et al., 2021). On the other hand, China has huge carbon sequestration potential of terrestrial ecosystems to slow accumulation of atmospheric carbon dioxide and mitigate climate change. Additionally, given China's water shortage, it is crucial to clearly understand water budgets and transportation (Ma et al., 2020). Therefore, it becomes vital to estimate GPP and ET accurately across China under changes in climate and land cover types.

Facing the above challenges, this study utilizes a water–carbon coupled and remote-sensing-based model, PML-V2, that is constrained against the most comprehensive observations in China (i.e., 26 EC observations across nine plant functional types – PFTs) to generate daily and 500 m ET and GPP gridded products from 2000 to 2020. We then test whether the newly developed product for China outperforms PML-V2(global) and other mainstream products (i.e., MOD16A2, SEBAL, GLEAM, MOD17A2H, VPM, and EC-LUE) and investigate annual change and spatial pattern in ET and its components: plant transpiration (E_c), evaporation from the soil (E_s), canopy evaporation from precipitation interception (E_i), GPP, and WUE in 2001–2018 across China. It is noted that the years of 2019–2020 have not been selected for trend analysis. This is because a different forcing dataset was used to drive PML-V2 (details are provided in Sect. 2.2.2).

The novelties of this study mainly include the following.

• i.

  Observation data from 26 EC flux stations in nine PFTs across China are employed to constrain the PML-V2 calibration for estimating ET and GPP.

• ii.

  The country-specific meteorological forcing, i.e., the China Meteorological Forcing Dataset (CMFD), is used to drive the PML-V2 in China, which is more accurate than those forcings extracted from global forcing products.

• iii.

  The PML-V2(China) product is generated with a daily resolution compared with the previous global product with a temporal resolution of 8 d.

• iv.

  It improves intra-annual ET and GPP dynamics for various ecosystems, particularly for the cropland ecosystem, which provides more accurate estimates and monitoring of agricultural water consumption compared with other mainstream products.

2 Materials and methods

2.1 Description of the PML-V2 model

The PML-V2 model is a water–carbon coupled diagnostic biophysical model. Compared to the old version that does not calculate GPP and the effect of CO₂ on evapotranspiration, PML-V2 couples a photosynthesis model (Thornley, 1976) and an improved canopy stomatal conductance model (Yu et al., 2004) with the PM equation to estimate GPP and E_c collectively (Gan et al., 2018). Y. Zhang et al. (2019) further improved PML-V2 by incorporating the vapor pressure deficit constraint into GPP that is then used to constrain canopy conductance and E_c. The detailed descriptions of PML-V2 are provided in the Supplement.

2.2 Model input data

2.2.1 Remote sensing data

MODIS collections with 500 m spatial resolutions from 26 February 2000 to 31 December 2020 (hereafter 2000–2020) are used, which includes the leaf area index (LAI) from MOD15A2H.006 (Myneni et al., 2015), the albedo from the MCD43A3.006 product (Schaaf and Wang, 2015), and the surface emissivity from MOD11A2.006 (Wan et al., 2015). The smoothed LAI inputs for PML-V2(China) utilized the weighted Whittaker smoother with dynamic lambdas instead of a constant lambda to eliminate missing or unreliable pixels due to noise contaminations by snow, shadow, cloud, etc., compared with PML-V2(Global) (Kong et al., 2019; Y. Zhang et al., 2019). The improved LAI can better express peaks and seasonal changes. The albedo and surface emissivity inputs were gap-filled by the linear interpolation of the nearest good-quality points. If there were not enough good-quality points close to the point with a missing value, it was filled by the historically averaged values for the same grid. In addition, the PML-V2 model needs land cover types to accurately estimate ET and GPP in different terrestrial ecosystems. Here we used the International Geosphere-Biosphere Program (IGBP) layer of the MCD12Q1.006 land cover product (Sulla-Menashe et al., 2019) during 2000–2020 since IGBP classification is annually continuous and has acceptable accuracy in China when compared with other land cover products (Feng and Bai, 2019).

2.2.2 Meteorological data

The meteorological inputs of PML-V2 include specific humidity, air pressure, air temperature, wind speed, precipitation, downward longwave radiation, downward shortwave radiation, and land surface temperature. In this study, the main meteorological data came from the CMFD from 2000 to 2018 (He et al., 2020), which has spatial and temporal resolutions of 0.1^∘ and 3 h, respectively. Generated through the fusion of five RS or reanalysis datasets and 753 China Meteorological Administration stations, the CMFD dataset shows the best accuracy among most available meteorological datasets and is widely employed in hydrological and land surface modeling of China (Z. Zhang et al., 2019; Wang et al., 2020). The Global Land Data Assimilation System Version 2.1 (GLDAS-2.1) forcing data with 0.25^∘ and 3 h resolutions (Beaudoing and Rodell, 2016) were used for 2019–2020, as the CMFD was not available during such a period. Since the CMFD and the GLDAS-2.1 are different meteorological forcing datasets for driving the PML-V2 model, it is necessary to correct the Bias of 2019–2020. Here, a widely used methodology, delta change (i.e., DC, also called the change factor), was selected for bias correction (Anandhi et al., 2011; Teutschbein and Seibert, 2012; Rasmussen et al., 2012; Hempel et al., 2013; Beck et al., 2018; Haro-Monteagudo et al., 2020). The underlying idea of the DC method is to use simulated future anomalies (i.e., GLDAS-2.1 in this study) for a perturbation of observed data (i.e., CMFD) rather than to use the simulations of future conditions directly. For each grid cell, we bias-corrected the daily meteorological data during 2019–2020 by monthly scaling factors. Specifically, a multiplicative approach was applied for precipitation (Eqs. 1 to 2).

$$\begin{array}{}\text{(1)}& P_{(i,j)}=\mathrm{\Delta }{P}_{\left(j\right)}\cdot P_{\left(i,j\right)}^{\ast };i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{31};j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{12},\end{array}$$

where P_(i,j) is the precipitation corrected by the relative change factor for day i and month j, and $P_{\left(i,j\right)}^{\ast }$ is the multiyear daily mean observed precipitation at the ith day, jth month in the historical or reference period (2000–2018 in this study). ΔP_(j) is the change factor that can be formulated as

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{P}_{\left(j\right)}=\stackrel{\mathrm{‾}}{P_{\mathrm{fut}}}/\stackrel{\mathrm{‾}}{P_{\mathrm{his}}};j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{12},\end{array}$$

where $\stackrel{\mathrm{‾}}{P_{\mathrm{his}}}$ is the daily mean precipitation of the jth month in the historic simulation (i.e., 2000–2018) and $\stackrel{\mathrm{‾}}{P_{\mathrm{fut}}}$ is similar but is the average of the future simulation (i.e., 2019 or 2020). The multiplicative approach is also used for the variables radiation, air pressure, and wind speed, whereas an additive approach was used to adjust temperature and specific humidity (Eqs. 3 to 4).

$$\begin{array}{}\text{(3)}& X_{(i,j)}=\mathrm{\Delta }{X}_{\left(j\right)}+X_{\left(i,j\right)}^{\ast };i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{31};j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{12},\end{array}$$

in which ΔX_(j) is calculated as

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{X}_{\left(j\right)}=\stackrel{\mathrm{‾}}{X_{\mathrm{fut}}}-\stackrel{\mathrm{‾}}{X_{\mathrm{his}}};j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{12}.\end{array}$$

The land surface temperature during 2000–2020 was from the Land component of the fifth generation of European ReAnalysis (i.e., ERA5-Land; Muñoz-Sabater et al., 2021) with spatial and temporal resolutions of 0.1^∘ and 1 h, respectively. Note that the above meteorological data were first aggregated into the daily scale, followed by resampling into 500 m by the bilinear interpolation method (Y. Zhang et al., 2019). The atmospheric CO₂ concentration data came from the National Oceanic and Atmospheric Administration (NOAA, 2021). The preprocessing of remote sensing data and meteorological data for model inputs is summarized in Fig. 1.

Figure 1 Flowchart of EC flux and AWS data pre-processing and PML-V2 model processing which is used to convert RS images and meteorological forcing images into GPP, E_c, E_i, E_s, E_w, and ET. For the pre-processing part: NEE (net ecosystem exchange, µmol m⁻² s⁻¹), LE (latent heat, W m⁻²), Rg (incoming radiation, W m⁻²), rH (relative humidity, %), T_a (air temperature, ^∘C), and QC (quality control). For the PML-V2 model part: T_max (daily maximum temperature, ^∘C), T_min (daily minimum temperature, ^∘C), T_avg (daily mean temperature, ^∘C), Pa (atmosphere pressure, kPa), U (wind speed at 10 m height, m s⁻¹), q (specific humidity, kg kg⁻¹), Prcp (precipitation, mm d⁻¹), Rl (inward longwave solar radiation, W m⁻²), Rs (inward shortwave solar radiation, W m⁻²), Pi (the difference of Prcp and Ei, mm d⁻¹), Es_eq (equilibrium evaporation, mm d⁻¹), ET_w (evaporation from water body, snow, and ice, mm d⁻¹) and GEE (© Google Earth Engine).

2.3 ET and GPP from eddy covariance observations

We collated EC flux towers and automatic weather station (AWS) data from 26 sites across China (Fig. 2 and Table 2) and generated the high-quality ET and GPP observed for calibration and validation of PML-V2. These data came from the following sources: FLUXNET2015 (Pastorello et al., 2020), the National Tibetan Plateau Data Center (Ma et al., 2020), the Heihe Integrated Observatory Network (Liu et al., 2011, 2018), and the Chinese Terrestrial Ecosystem Flux Research Network (ChinaFLUX) (Yu et al., 2006). These 26 sites encompass nine major PFTs in China, including two in evergreen needleleaf forests, one in evergreen broadleaf forests, one in mixed forests, one in open shrublands, one in savannas, eight in grasslands, three in wetlands, seven in croplands, and two in barren sparse vegetation. The observation variables, including air temperature, relative humidity, incoming shortwave radiation, latent heat flux, sensible heat flux, and net ecosystem exchange, were collected from the interval of 0.5 h or 1 h to prepare for gap filling and flux partitioning (Reichstein et al., 2005; Wutzler et al., 2018). Considering that certain gaps exist in the original half-hourly or hourly latent heat (LE) and net ecosystem exchange (NEE) flux data, we employed the marginal distribution sampling method (Reichstein et al., 2005) to fill these gaps using the station-observed air temperature, relative humidity, and solar radiation data. Subsequently, we partitioned NEE into gross carbon uptake (GPP) and respiration of the ecosystem according to the nighttime-based method of Reichstein et al. (2005). Because any gap filling of EC data may introduce extra uncertainties, we only used the days during which the percentage of the original observed and good-quality gap-filled data was no less than 60 % in the present study. Note that the data of the sites from the ChinaFLUX (i.e., CN-CBF, CF-HBG_S01, CF-HBG_W01, CF-NMG, CF-QYF, and CF-YCA) and the FLUXNET2015 (i.e., CN-Cng, CN-Du2, and CN-HaM) have already been gap-filled by the original data providers. Therefore, they were used directly in this study. Note that while energy imbalance does exist at many EC sites, correcting such a problem may also introduce more uncertainties (Foken, 2008). Therefore, we used the observed LE directly in the present study.

Figure 2 Geographical locations of 26 EC flux towers for nine major IGBP PFTs, the main rivers, and the 10 major river basins in China. Overlain are 20-year mean annual aridity index (AI) values during 2001–2020 using GLDAS-2.1, that is, the ratio of annual precipitation to Penman potential evapotranspiration. PFTs shown in the legend are ENF (evergreen needleleaf forests), EBF (evergreen broadleaf forests), MF (mixed forests), OSH (open shrublands), SAV (savannas), GRA (grasslands), WET (permanent wetlands), CRO (croplands), and BSV (barren sparse vegetation).

Table 2 Details of 26 EC flux towers employed in this study. Note that AP indicates mean annual precipitation and AT refers to mean annual temperature in its observed period.

2.4 Basin-scale water-balance-based evapotranspiration data

The water-balance method is generally regarded as a simple and accurate approach for calculating land evapotranspiration at the basin scale (Liu et al., 2016). Here we used the water-balance-based evapotranspiration (ET_wb, mm) of 10 major basins across China to evaluate PML-V2 ET estimates at the basin scale, that is,

$$\begin{array}{}\text{(5)}& {\mathrm{ET}}_{\mathrm{wb}}=\mathrm{Prcp}-Q-\mathrm{TWSC},\end{array}$$

where Prcp, Q, and TWSC (all with a unit of millimeters) are basin-wide precipitation, runoff, and change in terrestrial water storage at an annual scale, respectively. Among them, Prcp and Q are the annual values of 10 major river basins in China from 2003 to 2013, including the Hai, Huai, Liao, Northwest, Pearl, Songhua, Southeast, Southwest, Yangtze, and Yellow (Fig. 2), from the National Water Resources Bulletin (2021), which is extensively used in water resource calculation (Miao et al., 2022) and assessment (Yang et al., 2004; Xie et al., 2018). TWSC was quantified using three Gravity Recovery and Climate Experiment (GRACE) products (Landerer and Swenson, 2012; Landerer, 2021), including the NASA Jet Propulsion Laboratory, the GeoForschungsZentrum Potsdam, and the Center for Space Research, which have been available since April 2002 at a monthly scale. To reduce the uncertainties, this study used the mean values of these three products. We regarded the differences in the terrestrial water storage anomaly between two consecutive Decembers as the annual TWSC. Note that the ET_wb was only calculated from 2003 to 2013 in the present study since the December values from GRACE were not available in 2014.

2.5 Model calibration and model validation

The 11 parameters of the PML-V2 model for each PFT were calibrated and cross-validated against 26 EC sites by a global optimization method – the genetic algorithm (GA). The GA generates a randomly initialized population and then evaluates the fitness of solutions according to its objective function. As generations iterate, the population includes more appropriate solutions, and eventually it will converge (Holland, 1992; Konak et al., 2006). Specifically, we applied the GA algorithm with a population size of 1000 and number generations of 50. All EC-observed ET and GPP data within a PFT are used to minimize the following objective function (F_opt):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{F}_{\mathrm{opt}}& =\mathrm{2}-{\mathrm{NSE}}_{\mathrm{ET}}-{\mathrm{NSE}}_{\mathrm{GPP}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left({\mathrm{ET}}_{\mathrm{est}}-{\mathrm{ET}}_{\mathrm{obs}}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left({\mathrm{ET}}_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{{\mathrm{ET}}_{\mathrm{obs}}}\right)}^{\mathrm{2}}}}\\ & +{\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left({\mathrm{GPP}}_{\mathrm{est}}-{\mathrm{GPP}}_{\mathrm{obs}}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left({\mathrm{GPP}}_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{{\mathrm{GPP}}_{\mathrm{obs}}}\right)}^{\mathrm{2}}}},\end{array}\end{array}$$

where NSE_ET and NSE_GPP are the Nash–Sutcliffe efficiency of the daily ET and the daily GPP, respectively. The subscripts est and obs stand for estimated and observed, respectively. In this way, each of the nine PFTs gained a unique set with 11 calibrated parameter values, illustrated in Table S1 in the Supplement.

The “leave-one-out” cross-validation method was utilized to evaluate the robustness of the PML-V2 model. For each PFT, the data from one “ungauged” observation were excluded from the optimization, while the data from all other observations at the same PFT were used for model calibration to obtain the simulated data at the “ungauged” position. All nine PFTs were actualized in this way. Note that the PFT including EBF, MF, OSH, and SAV only has one ground site (Table 2). Therefore, it is appropriate to divide the data at each of the four sites into two sub-groups for cross-validation. The CF-CBF and the CF-HBG_S01 covering 2003 to 2010 were divided into two sub-groups, each of which had 4 years: 2003–2006 and 2007–2010, while both the BNXJL and YJGRHG only covered 1 year and were divided into two sub-groups by a 2 d time step separately. After that, the daily estimates in the cross-validation mode were against the daily observation from the 26 stations to explore the model transferability from known observations to any location.

2.6 Model performance metrics

We assessed the performance of calibration and cross-validation of PML-V2 (and another seven mainstream ET and GPP products) against the observed sites or water-balance basins utilizing the following four metrics:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathrm{NSE}}_X=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left(X_{\mathrm{est}}-X_{\mathrm{obs}}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left(X_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{X_{\mathrm{obs}}}\right)}^{\mathrm{2}}}},\text{(8)}& {\displaystyle }{R}_X={\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(X_{\mathrm{est}}-\stackrel{\mathrm{‾}}{X_{\mathrm{est}}}\right)\left(X_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{X_{\mathrm{obs}}}\right)}{\sqrt{\sum _{i=\mathrm{1}}^N{\left(X_{\mathrm{est}}-\stackrel{\mathrm{‾}}{X_{\mathrm{est}}}\right)}^{\mathrm{2}}\times \sum _{i=\mathrm{1}}^N{\left(X_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{X_{\mathrm{obs}}}\right)}^{\mathrm{2}}}}},\text{(9)}& {\displaystyle }{\mathrm{RMSE}}_X=\sqrt{{\displaystyle \frac{\sum _{i=\mathrm{1}}^N(X_{\mathrm{est}}-X_{\mathrm{obs}}{)}^{\mathrm{2}}}{N}}},\text{(10)}& {\displaystyle }{\mathrm{Bias}}_X={\displaystyle \frac{\sum _{i=\mathrm{1}}^N(X_{\mathrm{est}}-X_{\mathrm{obs}})}{N\times \stackrel{\mathrm{‾}}{X_{\mathrm{obs}}}}},\end{array}$$

where NSE, R, RMSE, and Bias are the Nash–Sutcliffe efficiency, the correlation coefficient, the root mean square error, and the ratio of the difference between the estimated and observed values to the observed average. The subscript X represents ET or GPP; the subscripts est and obs stand for estimated and observed, respectively.

3 Results

3.1 Model calibration and model validation

The simulated ET and GPP from the calibrated PML-V2(China) were first evaluated against EC measurements of 26 flux sites at a daily scale (Fig. 3). Overall, PML-V2(China) shows an excellent performance in estimating daily ET and daily GPP, as evidenced by the NSE (0.75 and 0.82, respectively), R (0.88 and 0.90, respectively), RMSE (0.69 mm d⁻¹ and 1.71 g C m⁻² d⁻¹, respectively), and Bias (−5.81 % and −2.3 %, respectively). For the mean values of each site, the simulated daily ET and daily GPP show higher NSE (≥ 0.87) and R (≥ 0.93) values (Fig. 3).

Figure 3 Scatter plots between the observed ET and GPP against PML-V2(China) simulations in the calibration and cross-validation modes: daily comparisons in the left panels and site mean comparisons in the right panels.

PML-V2(China) is only slightly degraded from calibration to cross-validation, indicated by a slightly declined performance in ET and GPP (Fig. 3). For daily ET, the NSE and R values decreased by 0.06 and 0.04, respectively, from the calibration mode to the cross-validation mode. Correspondingly, the RMSE and Bias of ET in the cross-validation mode increase by 0.08 mm d⁻¹ and 3.5 %, respectively. For daily GPP, the NSE and R values in the cross-validation mode reduce by 0.11 and 0.06, respectively; the RMSE and Bias increase by 0.45 g C m⁻² d⁻¹ and 1.79 %, respectively. A similarly slight degradation is applied to their site means. These results demonstrate that PML-V2(China) is robust for estimating daily ET and daily GPP across large regions and suitable for generating good-quality daily ET and daily GPP data for China.

Figure 4 further summarizes PML-V2(China) performance at 26 flux sites across nine PFTs. The estimates of ET and GPP from the model calibration show high consistency with the EC-observed values in all terrestrial biomes. For daily ET (Fig. 4a), the NSE values vary in the range 0.36∼ 0.82, the RMSE 0.39 ∼ 0.88 mm d⁻¹, and Bias −10.09 % ∼ −0.21 %. For daily GPP (Fig. 4b), the ranges of statistical metrics become 0.41–0.91 for NSE, 0.3–3.19 g C m⁻² d⁻¹ for RMSE, and −10.52 %–3.26 % for Bias. In terms of cross-validation, nine PFTs all showed slight declines in the statistical metrics when compared to those in the calibration mode. For daily ET, the declines in NSE values are less than 0.14 in most PFTs except barren sparse vegetation (BSV) and ENF, whose NSE values decreased by 0.36 and 0.33, respectively. As expected, RMSE values all increased to some extent in all PFTs (ranging from 0.002 to 0.305 mm d⁻¹) when compared with those in calibration mode. The Bias values in the cross-validation mode were almost identical to those in the calibration mode for most PFTs except WET and ENF, of which the absolute value of Bias increased by 10.59 % and 17.42 %, respectively (Fig. 4a). From calibration to cross-validation, the degradation of BSV, ENF, and WET is more serious than that for the remaining PFTs, which is mainly caused by the small samples (2, 2, and 3, respectively) for ET estimates. Regarding daily GPP, the NSE values all degraded by less than 0.04 for most PFTs except BSV, GRA, and WET, where there exists 0.21–0.32 NSE degradation. In the meantime, the declines in R values are all within 0.19. Regarding RMSE, the increases are particularly marginal for most PFTs except WET, with an increase of 1.58 g C m⁻² d⁻¹ (Fig. 4b). The above PFT tests suggest that the present PML-V2(China), with parameter values being calibrated against 26 EC flux stations, does perform satisfactorily in estimating both ET and GPP across different PFTs in China.

Figure 4 Comparison of ET and GPP between PML-V2(China) model calibration and validation across 10 PFTs.

To investigate the model performance at each EC site, this study also compares the variations in daily ET and GPP between PML-V2(China) in calibration mode and the EC observations (Figs. 5 and 6 and Table S2). Overall, the estimates from PML-V2(China) show a similar amplitude and phase to the EC observations, indicating that it performs well in capturing the seasonal phase of ET and GPP at most flux sites. From Fig. 5, the PML-V2(China) ET reveals a single peak in each annual cycle for most flux sites, except those with double-cropping systems such as the CF-YCA, DXZ, and GTZ cropland sites. In particular, the ET values are consistent with observed data at the desert site (e.g., HZZHMZ and QZ-NAMORS), with NSE ranging from 0.41 to 0.48, indicating that the model performs well in the sparse vegetation area. The measurements fluctuate higher than the PML-V2(China) ET in peak areas for most sites. For the site mean ET, the difference values between PML-V2(China) and EC observations range from −0.46 to 0.21, with the minimum at QZ-NAMORS and the maximum at CN-HaM.

Figure 5 The daily ET simulated by PML-V2(China) in calibration mode and the observed daily ET variation in time series from 26 EC sites (see Fig. 2) across China. “ALL” represents the site mean value for each EC site.

Figure 6 The daily GPP simulated by PML-V2(China) in calibration mode and the observed daily GPP variation in time series from 26 EC (see Fig. 2) sites over China. “ALL” represents the site mean value for each EC site.

For daily GPP, the model also performs well in depicting the seasonal variation. However, at certain stations (e.g., DMCJZ) the peak GPP values within a year appear to be underestimated. In terms of the cropland flux sites, the GPP also shows double peaks within a year because of the double-cropping system (e.g., winter wheat and maize rotation), which is similar to ET. This is especially apparent for the GTZ, DXZ, and CF-YCA flux sites located in the North China Plain. For the site mean GPP, the discrepancies between the model and EC observations mainly range from −0.66 to 0.96 g C m⁻² d⁻¹ for most flux sites except CF-YCA, MYZ, and DXZ, where the differences exceed 1.3 g C m⁻² d⁻¹.

3.2 Comparing with other products

3.2.1 Comparisons at a plot scale using EC-observed data

To explore whether the ET and GPP of PML-V2(China) simulations are more accurate than the previous products, we also evaluated ET and GPP accuracy from its global version, MOD16A2, MOD17A2H, and another five widely available ET or GPP products against EC observations from the 26 sites at a daily or 8 d resolution. In this study, PML-V2(China) uses its cross-validated simulations to compare with other products instead of the calibration results to avoid introducing a priori knowledge. Additionally, to compare at a consistent time resolution, PML-V2(China) estimates in cross-validation mode need to be upscaled to an 8 d average or remain at a daily scale, depending on the temporal resolution of the comparison products. Specifically, PML-V2(China) is compared with GLEAM and SEBAL at the daily scale compared with PML-V2(Global), MOD16A2, MOD17A2H, EC-LUE, and VPM at the 8 d scale.

Table 3 provides a direct comparison of model performances among varieties of ET or GPP products against observations overall from 26 ground stations. It is evident that PML-V2(China) excels over other state-of-the-art ET or GPP products, presented by NSE being 0.12–7.76 higher for ET and 0.07–0.79 higher for GPP, R being 0.07–0.68 higher for ET and 0.0–0.51 higher for GPP, and RMSE being 0.15–3.62 mm d⁻¹ lower for ET and 0.24–1.98 g C m⁻² d⁻¹ lower for GPP. Specifically, at a daily scale, PML-V2(China) ET exhibits the highest NSE value of 0.66, followed by GLEAM (0.44) and SEBAL (−7.10). PML-V2(China) daily ET achieves the highest R (0.84), followed by GLEAM (0.69) and SEBAL (0.16); correspondingly it obtains the smallest RMSE (0.33 mm d⁻¹), followed by GLEAM (1.04 mm d⁻¹) and SEBAL (3.94 mm d⁻¹). SEBAL is the worst performer, although its Bias is closest to 0 (Table 3) because it is far away from the observed values, significantly yielding a Bias value of over 50 % or less than −50 % among the five PFTs (Fig. 7). At the 8 d scale, PML-V2(China) outperforms PML-V2(Global) and MOD16A2 for estimating ET, with the highest values of NSE (0.74) and R (0.87) and the lowest RMSE (0.66 mm d⁻¹). Moreover, PML-V2(China) also has the best performance in estimating 8-daily GPP, followed by PML-V2(Global), MOD17A2H, VPM, and EC-LUE, indicated by three statistics: NSE, R, and RMSE. In summary, PML-V2(China) performs well when compared with other mainstream ET or GPP products in China.

Figure 7 Statistical indicators of PML-V2(China) and other models for estimating ET at each PFT. Open and solid dots represent PML-V2(China)-estimated ET in cross-validation mode and other models. For the daily temporal resolution of GLEAM and SEBAL, PML-V2(China) is also daily scale, while for 8 d resolution of PML-V2(Global) and MOD16A2, the referred PML-V2(China) is upscaled to 8 d. Note that “(+)” indicates that the model's simulation statistics dot is more than the upper bound, while “(–)” indicates that the model's simulation statistics dot is less than the lower bound.

Table 3 Statistical indicators of PML-V2(China) and other models for simulating ET and GPP at 26 EC flux towers. NSE and R values are unitless. The unit of RMSE for ET is mm d⁻¹, while it is g C m⁻² d⁻¹ for GPP. The unit of Bias is %.

Figure 7 displays the performance comparison of four ET products with PML-V2(China) under nine PFTs. The simulated ET by PML-V2(China) has greater NSE and R and fewer RMSE values than the other four products in most PFTs, especially in EBF, SAV, and WET. All the models have poor performance, with NSE being lower than 0 except for PML-V2(China) in EBF and SAV. However, PML-V2(China) is not the best in both ENF and BSV. In BSV, most models perform poorly, rendered by NSE being lower than 0 except for GLEAM (0.50) and PML-V2(China) (0.07 for the daily scale, 0.11 for the 8 d scale). Only SEBAL achieves worse results than PML-V2(China) in ENF. As shown in Fig. 8, PML-V2(China) performs significantly better than other advanced products in simulating the GPP of CRO, MF, ENF, EBF, SAV, and BSV, producing higher NSE, R, and lower RMSE and Bias. While PML-V2(China) ranked second in GRA (following MOD17A2H), OSH (following PML-V2(Global)), and WET (following PML-V2(Global)). Synthetically, PML-V2(China) successfully captures the sites' seasonality in most PFTs compared to the high-resolution ET/GPP datasets currently available.

Figure 8 Statistical indicators of PML-V2(China) and other models for estimating GPP at each PFT. Open and solid dots represent PML-V2(China)-estimated GPP in cross-validation mode and other models. PML-V2(China) is upscaled to 8 d to compare with the 8 d resolution of PML-V2(Global), MOD17A2H, EC-LUE, and VPM. Note that “(+)” indicates that the model's simulation statistics dot is more than the upper bound, while “(–)” indicates that the model's simulation statistics dot is less than the lower bound.

3.2.2 Comparisons at the basin scale using ET_wb

In addition to testing the model at a plot scale, Fig. 9a–e presents the ET validations from PML-V2(China) and four ET products based on RS against the annual ET_wb in the 10 major river basins of China during 2003–2013. It illustrates that PML-V2(China) shows the best performance among them, as indicated by the highest NSE (0.82) and lowest RMSE (69.59 mm yr⁻¹) and Bias (6.28 %) values. This is closely followed by GLEAM and PML-V2(global) with NSE values of 0.36 and 0.26, respectively. However, MOD16A2 and SEBAL tend to overestimate ET in the majority of basins, with much smaller NSE values of −0.21 and 0.02, respectively, which are consistent with the performance evaluations using the EC observation shown in Fig. 7. The above basin-wide evaluations, together with the plot-scale validations against EC observation data, demonstrate PML-V2(China) overall performing best among the tested products selected in this study.

Figure 9f illustrates an inter-basin comparison of the 11-year mean ET of 2003–2013 within the five ET products. PML-V2(China) performs well in most basins (Biases within ±15 %) except for the Northwest and Southwest river basins, where ET is overestimated by 24.96 % and 61.57 %, respectively. Even so, PML-V2(China) performs best in the Southeast River Basin, with a Bias of −3.80 %, which is still better than the five ET products selected herewith, Bias ranging from 15.93 % to 62.42 %. Although PML-V2(China) overestimates ET in the Northwest basin to a large extent, it performs best relative to its global version (84.88 %) and another ET product, MOD16A2 (152.70 %).

3.3 Spatial patterns and annual variations of ET, E_c, E_i, E_s, GPP, and WUE

Figure 10 illustrates the spatial distribution of the multiyear (2001–2018) mean annual ET and three components (i.e., E_c, E_i, and E_s) from PML-V2(China) across China. In general, the ET shows an increasing gradient from the northwest to the southeast (Fig. 10a). High annual ET (> 900 mm yr⁻¹) is mainly located in the water bodies, Hainan, and western Taiwan, while most parts of the Northwest River Basin exist with low annual ET (< 100 mm yr⁻¹), especially in western Inner Mongolia and Gansu and southern Xinjiang. Annual ET experiences a statistically insignificant increasing trend during 2001–2018 with a tendency of 0.43 mm yr⁻² (p>0.05). On the whole, the mean annual ET over China is 392.12 ± 10.67 mm yr⁻¹ (mean ± standard deviation) over the last 18 years. For the three components, E_c and E_i products display a similar spatial distribution to annual ET, while high E_s values (> 400 mm yr⁻¹) are mainly scattered in higher soil moisture content areas, including the Tibet Plateau, Pearl River Delta, and Yangtze River Delta (Fig. 10b–d). High E_c (> 600 mm yr⁻¹) and E_i (> 80 mm yr⁻¹) values overall occur in the tropical and subtropical forests (e.g., Southwest River and Pearl River basins) but low E_c (< 50 mm yr⁻¹) and E_i (< 5 mm yr⁻¹) values in the Northwest River Basin except for the Tianshan, Altai, and Qilian Mountains. In particular, low E_i values also appear in the cropland areas, such as the Northeast Plain, North China Plain, and Sichuan Basin (Fig. 10c). For annual variation over China, E_c and E_i increase significantly during 2001–2018 with rates of 0.91 and 0.16 mm yr⁻² (p<0.001), respectively. However, annual E_s shows a declining trend with an insignificant rate of −0.69 mm yr⁻² (p<0.05).

Figure 9 Annual evapotranspiration (ET) of (a) PML-V2(China), (b) PML-V2(Global), (c) MOD16A2, (d) GLEAM, and (e) SEBAL plotted against the water-balanced derived ET (ET_wb) values for 10 major river basins over China during 2003–2013. The boxplot in panel (f) shows a multiyear mean of the five ET products above per river basin.

Figure 10 Spatial pattern of mean annual ET, E_c, E_i, E_s, and their annual variation during 2001–2018. In all the insets, the shaded areas represent the 95 % confidence interval based on the linear regression modeling. The number in the parentheses of each inset is mean ± standard deviation of the annual simulated variables during the 18 years.

The mean annual GPP shows similar spatial patterns compared to the mean annual ET, as indicated by Figs. 10a and 11a. A high annual GPP (> 2000 g C m⁻² yr⁻¹) mainly occurs in the tropical and subtropical forests and the North China Plain, where there exists a double-cropping system but low annual GPP (< 100 g C m⁻² yr⁻¹) in the arid zones such as the Northwest River Basin. On average, the multiyear GPP over China is 721.62 ± 51.83 g C m⁻² yr⁻¹, and interannual change displays a steady rising trend with a rate of 8.99 g C m⁻² yr⁻² (p<0.001) since 2001. Using the coupled estimation of the PML-V2(China) model, we study WUE (GPP divided by ET) during 2001–2018 across China (Fig. 11b). This result indicates the high annual WUE (> 3 g C mm⁻¹ H₂O) occurring in the forests and cropland, particularly in Northeast China and the North China Plain. The annual variation of WUE is similar to that of GPP, with a significant increasing trend (slope = 0.02 g C mm⁻¹ H₂O yr⁻¹, p<0.001).

Figure 11 Spatial pattern of mean annual GPP, WUE and their annual variation during 2001–2018. In all the insets, the shaded areas represent the 95 % confidence interval based on the linear regression modeling. The number in the parentheses of each inset is mean ± standard deviation of the annual simulated variables during the 18 years.

4 Discussion

4.1 Magnitude and trend in annual ET and GPP over China

For annual ET over China, the multiyear (2001–2018) mean annual ET from PML-V2(China) is 392.12 ± 10.67 mm yr⁻¹ (Fig. 10a). This result is overall consistent with the country-wide averaged annual ET estimated by the machine learning method (Yin et al., 2021: 397.65 mm yr⁻¹ for 2000–2018) and land surface models (N. Ma et al., 2019: 395.34 mm yr⁻¹ for 2001–2012) and slightly higher than MOD16A2 ET at about 359.61 ± 59.52 mm yr⁻¹ for 2001–2018 (Cheng et al., 2021). However, they are all less than the annual ET of about 482.27 ± 192.31 mm yr⁻¹ from SEBAL for 2001–2018 (Cheng et al., 2021). Furthermore, previous studies by Ren et al. (2015) and G. Q. Wang et al. (2012) show that the long-term mean precipitation and runoff in China are about 720 mm yr⁻¹ and 280 mm yr⁻¹, respectively. Hence, it is believed that an annual ET of less than 440 mm could be reasonable in China (N. Ma et al., 2019). The annual mean GPP over China from our results is 721.62 ± 51.83 g C m⁻² yr⁻¹ during 2001–2018 (Fig. 11a), which is lower than that of Jia et al. (2020) (771 g C m⁻² yr⁻¹) and higher than those of Yao et al. (2018) and J. Ma et al. (2019) (690 and 710 g C m⁻² yr⁻¹, respectively). These differences may be associated with the distinctions in the time window and data sources (Jia et al., 2020).

The annual ET displays a statistically insignificant increasing trend from 2001 to 2018, which is consistent with the calculated ET using the Budyko equation (Feng et al., 2018; Su et al., 2022). In terms of annual GPP, we found that there is a significant (p<0.001) increasing trend with a rate of 8.99 g C m⁻² yr⁻² during 2001–2018, in line with some other studies (J. Ma et al., 2019, 2018; Yao et al., 2018). The most likely reason for the remarkable rise in GPP is the effect of ecological restoration projects in China (Tong et al., 2018). In fact, a large number of ecological restoration projects have been conducted since the 1990s, such as the Grain for Green Project (Cao et al., 2009). These findings also confirmed that a significant increase in vegetation growth occurred in China over the past years, which agreed well with Ma et al. (2018).

4.2 Advantages of this new dataset

The multiscale testing using EC observations and water balance showed that the accuracy in ET and GPP by the present PML-V2(China) is better than the global product of PML-V2 and other mainstream ET or GPP models (Table 3). The reasons may be twofold. The first is that the water–carbon coupled process is particularly important for estimating ET and GPP since the water and carbon processes are highly coupled by the stomatal aperture at the leaf level. This result was also supported by Xiao et al. (2013) and Y. Zhang et al. (2019) in their recent studies. Second, this study employed 26 EC observations to calibrate the PML-V2 in China, which shows better accuracy than the previous global-scale ET and GPP estimates that were obtained using few EC observations to constrain the parameters. This indicates that more local observations will facilitate the improvement of ET and GPP estimates at regional and national scales. In fact, although the EC sites of MOD16A2 (72 EC sites), GLEAM (91 EC sites), and PML-V2(Global) (95 EC sites) are more than in this study, there are only 0, 8, and 8 sites in China, respectively (Table 1). In particular, the SEBAL model only used eight EC sites for three PFTs (i.e., forests, cropland, and grassland). In addition to the advantages of the overall accuracy in ET and GPP in the present study, PML-V2(China) showed its strong ability to reveal the characteristics of the water consumption from the croplands that have a double-cropping system. The GTZ in Hebei, DXZ in Beijing, and CF-YCA in Shandong are the only three observed sites with winter wheat–summer maize rotation cropping systems. We compared the intra-annual variation of the simulated ET and GPP between PML-V2(China) and other products against the EC-observed values at the three cropland sites (Fig. 12). In theory, when the winter wheat is harvested, the ET or GPP should decline to their valley values in June, which often occur between the two peaks (i.e., the reproductive growth stage of winter wheat and summer maize, respectively) within a given year. With this in mind, it can be seen that PML-V2(China) has improved its ability when compared to its global version, as has been indicated by its better performance in capturing the time when the lowest ET and GPP values emerged. This is mainly because an improved weighted Whittaker smoother was carried out to get better quality of LAI, as described in Sect. 2.2.1. While GLEAM is also able to detect the time when the valley values appeared, it underestimates ET evidently during the wheat growing season. In terms of SEBAL and MOD16A2, both have much poorer performances in detecting such intra-annual variations in ET. Regarding the other GPP models, only MOD17A2H can catch the time when the valley values appear. However, it substantially underestimates GPP in both the winter wheat and summer maize growing seasons. Moreover, this study also estimated the ability of the simulated ET to identify the crop phenology at the regional scale (Fig. S1). We extracted the cropland with peaks and identified the dates of peaks appearing within a year at each pixel by a faster peak detection algorithm (Liu et al., 2020). Taking the typical double-cropping system as an example, we quantified the cropping intensity in croplands (Fig. S1a1) and identified the dates of the first peak and the second peak appearing in 1 year (Fig. S1a2, a3). To verify the reliability of the results, we mapped the double-cropping cropland areas of winter wheat and summer maize rotations (Fig. S1b1) and the heading date distribution of winter wheat and summer maize (Fig. S1b2, b3) in 2015 based on the crop phenological dataset (Luo et al., 2020). The croplands with double cropping show similar spatial patterns, as indicated by Fig. S1a1 and b1. In particular, we also compared the first ET peak date (Fig. S1a2) with the heading date of winter wheat (Fig. S1b2) in 2015. The first ET peak date (i.e., day of year – DOY) is mainly between DOY 120 and 150, occurring after the heading date of winter wheat at about DOY 100 to 130. During the entire growth period of winter wheat, the ET intensity was highest from the heading date to the filling date, although the critical periods of water demand for winter wheat are the “jointing date-heading date-filling date” (He et al., 2022). This is consistent with our result that the ET peak appears slightly later than the crop heading date. Similarly, the heading date of summer maize also occurs earlier than the second ET peak date (Fig. S1a3, b3), which further proves the ability of the simulated ET in evaluating crop phenology.

Figure 12 The intra-annual variation of (a) ET at three crop-rotation stations between those observed and simulated by PML-V2(China) in validation mode, PML-V2(Global), SEBAL, GLEAM, and MOD16A2, respectively; (b) GPP at three crop-rotation stations between those observed and simulated by PML-V2(China) in validation mode, PML-V2(Global), EC-LUE, VPM, and MOD17A2H, respectively. The blue dotted lines above pass through the lowest values between the two peaks of the observed ET or GPP per year. Note that all variables are averages every 8 d, although units are per day.

4.3 Implications of PML-V2(China)

Based on the substantial advantages discussed above, PML-V2(China) has great implications and application prospects. For instance, daily outputs from PML-V2(China) can be better used by the agricultural and water sectors for operational applications. Timely access to daily data at the regional or national scale helps the Ministry of Agriculture and Water Resources to develop better policies. Indeed, there is a remarkable relationship between soil water content and ET (Graf et al., 2014; Brust et al., 2021), so getting daily ET information accurately is of great significance for soil water depletion assessment, irrigation system design, and water resource management in agricultural areas, such as in the North China Plain. On the other hand, this dataset has better simulations of carbon consequences and water use efficiency, which is important for carbon–neutron policy (Yang et al., 2022). Specifically, for 2001–2018, the annual GPP and water use efficiency experienced a significant increase (8.99 g C m⁻² yr⁻² and 0.02 g C mm⁻¹ H₂O yr⁻¹, respectively), but annual ET showed a non-significant increase (0.43 mm yr⁻²). This indicates that vegetation in China exhibits a huge potential for carbon sequestration with little cost in water resources, which plays an important role in the global carbon cycle.

4.4 Uncertainties

4.4.1 Eddy covariance method and water-balance formula

Although PML-V2(China) showed relatively good performance when compared to the EC sites and water-balance-based evapotranspiration (ET_wb), there exist several uncertainties related to the observed data (i.e., flux sites and ET_wb). First, the EC technique, considered a standard way to measure surface fluxes (Aubinet et al., 1999; Liu et al., 2011; Baldocchi, 2014), meets some issues and includes corrections when processing the turbulence data from the energy non-closure problem. The corrections, such as the spike detection, lag correction of H₂O and CO₂ based on the vertical wind, coordinating rotation, corrections for density fluctuation, and frequency response correction, had been pre-processed before the investigators shared data. Nevertheless, there is evidence that diverse data processing designs may lead to errors of 10 %–15 % (Mauder et al., 2007). In addition, systematic bias in the device, the loss from the contribution of low-frequency eddies to energy transmission, and the ability to capture larger eddies and the secondary circulations could cause the energy-imbalance problem (Liu et al., 2011). Hence, there are usually two schemes to deal with the energy non-closure issue for EC users, that is, to perform energy closure correction (Cheng et al., 2021) or to maintain the original four-component data including latent heat (LE), sensible heat (H), soil heat flux (G), and net radiation (R_n) (Y. Zhang et al., 2019; Ma and Zhang, 2022). We chose the second method in this study, considering that (i) forcing energy closure will introduce new errors artificially and (ii) most Chinese EC observation towers lack G and R_n. The observed ET calculated from the latent heat flux of the site without energy closure correction will be slightly less than the real value, resulting in a smaller ET simulated by the model determined by calibration using the sites. This phenomenon is not fully reflected in comparisons with the basin ET based on the water-balance calculation that only in the Pearl and Southeast river basins does PML-V2(China) underestimate multiyear mean ET (Fig. 9f), because the ET_wb used for assessing the simulation performance of models also needs to be explored in its accuracy.

Second, the inconsistency of the grid cell and EC footprint could also result in uncertainty when compared rudely to the measurements. Generally, the EC towers have a footprint of 100–1000 m², which is usually decided by the tower height and heterogeneity of the underlying surface (Liu et al., 2016; Xu et al., 2017). For example, the footprint of the forest sites is larger than the grassland and wetland sites (Chen et al., 2012). In this study, the PML-V2(China) model is first calibrated at a 500 m grid cell, avoiding the inconsistency issue to some extent. However, there still exists a mismatch between the grid cell center and EC sites. In this regard, higher spatial resolution products or additional ground observations at relevant scales would be beneficial for the cross-validation of the modeling grid cell (N. Ma et al., 2019). Note also that the limited flux sites for some PFTs may introduce extra uncertainties for model parameters since only one site was available for the OSH, SAV, EBF, and MF in this study. In fact, we artificially construct multiple site samples in the above PFT sites, utilizing the characteristics of the long time series of these sites. Nevertheless, the terrestrial biome of the sole flux site may not be typical in other climate zones for the same PFTs (Cheng et al., 2021). Therefore, more flux sites for the same PFTs are necessary to calibrate the model.

Third, ET derived from water balance over China was invested as a reference at the basin scale, although the results may be affected by some sources of uncertainty. For instance, the applicability of the water balance relies on its formula composition. In this paper, the water storage change (in Eq. 5) from GRACE was included for the purpose of reducing the uncertainties in estimating annual ET_wb. Those studies not using TWSC in ET_wb may not explain the decreasing groundwater such as in the Hai River Basin due to the human water extraction, which is not conducive to the credibility of the verification results (Cheng et al., 2021). In addition, the precision of its input data also affects the reliability of ET_wb (Mao and Wang, 2017). Precipitation, as the main source of ET, impacts ET_wb to varying extents. Nevertheless, precipitation data are derived from observations of the field rain gauge network, and their usability relies on the intensive and high-quality ground observations, which makes the Prcp estimates from statistics of stations worse in the less-populated remote regions or areas with highly various topography, particularly in western China (Immerzeel et al., 2015; Tang et al., 2016; Zhong et al., 2019). Zhong et al. (2019) evaluated three precipitation products in China and found that there was a slight overestimation in the west of China and an obvious underestimation in the western Tibet Plateau. Accordingly, it could be stated that ET_wb was overestimated in the Pearl River Basin and the Southeast River Basin and underestimated in the Southwest River Basin using Eq. (7).

4.4.2 Input data

While the daily ET and GPP of the PML-V2(China) product (in the calibration mode) simulated well against 26 flux sites overall (Fig. 3) and in most PFTs (Fig. 4), PML-V2(China) in calibration is degraded compared to its cross-validation (Fig. 4), such as in GRA. In this study, we got one parameter set for the GRA type by employing eight sites, including the sparse grassland (QZ-BJ) and the dense grassland (QZ-NAMORS). Although it may be appropriate to use diverse parameter values for estimating ET and GPP by further dividing the grassland type into finer land types, this comes at the expense of ignoring the possible interannual changes in land types because of few land cover maps with fine classifications and annual resolutions simultaneously (Ma and Zhang, 2022).

PML-V2(China) mainly used the remote sensing and meteorological data (e.g., MODIS, CMFD, GLDAS-2.1, and ERA5) as the inputs (see Sect. 2.2). However, there are still some uncertainties in these data. For example, we used the land cover datasets (MCD12Q1.006) as the PFT data across China. However, there exist misclassification issues for MCD12Q1 because of spectral confusion (e.g., savannas and grasslands) and coarse resolution (e.g., the mixed pixel of cropland and natural vegetation) (Y. Zhang et al., 2019; Liang et al., 2015; Adzhar et al., 2022). Moreover, LAI is a critical variable describing vegetation growth, and its temporal changes affect stomatal conductance and further affect the transpiration rate from the vegetation canopy. Using LAI, the PML-V2(China) model simulated E_c and E_i components at the canopy scale. To avoid noise issues caused by clouds, shadows, snow, and so on, MODIS LAI in this study has been smoothed by the weighted Whittaker smoother which can deal with underestimation and inefficiency issues (Kong et al., 2019). However, there are still underestimates in the sparse vegetation areas. This may explain why the ET and GPP estimates are poor in BSV (Figs. 7 and 8).

Additionally, downscaling uncertainties could also be introduced by the bilinear interpolation method which has been applied to minimize the footprint impact of coarse meteorological inputs, such as CMFD (Fig. 1). This approach depends only on nearby grid cells to downscale, which could neglect the other relative supports. For instance, precipitation is not only related to the surrounding precipitation, but also location and terrain (e.g., elevation and aspect) (Yue et al., 2020). Chao et al. (2018) found that gridded precipitation products in the high-altitude regions are far below what is inversely inferred by glacier mass balances. Consequently, the geographically weighted regression method coupled with a weighting function could work well to interpolate meteorological data (Chao et al., 2018).

5 Data availability

The product named PML-V2(China) with daily and 500 m resolutions from 26 February 2000 to 31 December 2020 is freely available at the National Tibetan Plateau Data Center (https://doi.org/10.11888/Terre.tpdc.272389, Zhang and He, 2022).

6 Code availability

The core source code of PML-V2(China) is available on request.

7 Conclusions

This study developed a daily, 500 m ET and GPP data product (PML-V2(China)) using the locally calibrated water–carbon coupled model, PML-V2. The model has been well-calibrated against observations at 26 flux sites across nine plant functional types, and it performs satisfactorily in the cross-validation mode. More importantly, the plot- and basin-scale evaluations suggest that the newly developed product outperforms not only the global version of PML-V2, but also other mainstream RS-based ET and/or GPP products. With such a new product, we investigated spatial patterns and trends in ET and its components (E_c, E_i, E_s), GPP, and WUE from 2001 to 2018 across China. In short, the present PML-V2(China) product has the following advantages. (i) The water output is tightly constrained by carbon flux. (ii) It has high spatial and temporal resolutions simultaneously. (iii) It obtains the improved accuracy in ET and GPP across different plant functional types because of the optimal parameter sets for China by exploiting 26 EC sites. (iv) It shows the better ability to reveal the ET and GPP for the croplands with the double-cropping system. In summary, we provide a novel daily and 500 m resolution ET and GPP product across China, which can be used by research communities and various water and ecological departments for operational applications.