Figure 2A shows the spatial pattern of the mean annual value of ET during 1981–2010. The high values (>1,000 mm) mainly existed in the Brazilian coast, the GulfofMexico and Atlantic coasts in America, the African coast, and the Oceania East coast. Specifically, the ET decreased from east (west) to west (east) across the North (South) America, and from southeast to northwest across China. By contrast, the spatial variability of ET CV was not line with the ET value: the high ET occurred in the Amazonian Plain and Brazilian plateau, whereas high ET CV occurred in South China (Figure 2B). Additionally, the ET tended to increase in the Eurasia and Brazilian plateau, while decreasing in the Amazonian Plain (Figure 2C). Overall, about 20% of basins showed the significant trends, and the significant increases were mainly in the Northwest China, Europe, and the midwest U.S, while the significant decreases were mainly in the Congo Basin and Amazonian Plain) (Figure 2D). In conclusion, ET regarding magnitude, temporal variability and trend showed the high spatio-temporal heterogeneity.

All ET products could reproduce the spatial distribution for climatological values of ET with high spatial R values ≥ 0.90 (Table 3). Among these products, the PML performed slightly better than the other products with the highest and R value of 0.96, though not with optimal β and γ . The β values for most ET products were consistently around the optimal value of 1.0, except for GLASS (1.27 for β ) and MERRA-2 (1.22 for β ), suggesting that the magnitudes of climatological values were well captured by most ET products. However, the spatial variabilities of ET tended to be underestimated by most ET products with 0.7< γ <1.0, especially for ERA5-Land with a γ of 0.74. The values of R, β and γ for capturing the ET by the products converged to the optimum value of 1.0, resulting in KGE values of integrative performances that outweighed 0.71. Notably, the model-based products had higher overall performance (KGE ≥ 0.81) in reproducing the climatological ET, compared to the reanalysis products (0.71 ≤ KGE < 0.81). Regarding the temporal variability, all ET products generally underestimated the CV (0.34≤ β ≤0.89), but evidently overestimated its spatial variability (1.21≤ γ ≤2.39). Moreover, the spatial distribution of ET CV were poorly captured be most ET products, with GLEAM-3.5a having the maximum R value of 0.25 among the eight ET products. Overall, the KGE values were mostly negative, ranging from −0.80 (MTE) to 0.04 (MERRA-2), indicating that most ET products had limited KGE-based ability to simulate ET temporal variability. In the view of the ET trend, the directions (i.e., upward and downward) could be hit by most products, with 59.29% for PML ≤ TPR + TNR≤65.66% for MERRA-2. However, the FPR, near to and even larger than the TNR, suggested that the negative trends would be misidentified as positive trends, especially for GLASS (39.68% versus 3.04%). The R values ranged from 0.09 (MTE) to 0.36 (MERRA-2) indicating that GLASS and MERRA-2 with values above 0.30 could capture the ET trends in space. Except for reanalysis products underestimating the ET trend, all others overestimated the ET trend, with β larger than 1.0. By contrast, all products underestimated the spatial variability of the ET trend, with −0.19 (MERRA-2)≤ γ ≤0.18 (GLEAM-3.5a). All KGE values were negative, indicating that these poor overall performance of these ET products in capturing the ET trend.

Figure 3 shows the metrics of β , γ , R, and KGE for 1,381 basins. The majority of ET products overestimated ET at more than 50% of basins, especially GLASS and MERRA-2 which overestimated ET at above 92% of basins (Figure 3A). However, PML, GLDAS2.0-Noah and MTE underestimated the ET at more than 50% of basins. Spatially, the β values displayed evident spatial differences, with most of ET products greatly overestimated the ET in China, Europe, and North America. Considering the metric of γ (Figure 3B), all ET products tended to underestimate the ET temporal variabilities at over 70% of basins. When γ <0.2, ET products, especially MTE and GLASS, underestimated the ET temporal variabilities at around 30% of basins worldwide. Additionally, the overestimates of ET temporal variabilities tended to be at American Midwest. About the spatial patterns of β and γ , it is worth noting that the higher ET magnitude estimates were accompanied by lower ET variability estimates, since the ratio of basins having β >1.0 outweighed the ratio of basins having γ <1.0 for most ET products except PML, GLDAS2.0-Noah, and MTE (Figures 3A,B). Regarding temporal fluctuation (Figure 3C), positive R values were observed for 59.30% (MERRA-2) to 84.50% (ERA5-Land) of basins, especially MERRA-2 with R > 0.6 at nearly 20% of basins, indicating that ET products had a broad R-based ability to simulate ET temporal fluctuation. High R values (around 0.8) mainly appeared in the Midwest United States, South Africa, Western Australia. However, the average R values for all ET products were slightly low, ranging from 0.06 for GLDAS2.0-Noah–0.24 for ERA5-Land. Based on KGE (Figure 3D), negative KGE values were found in 47.65% (E2O-En) to 71.76% (MTE) of basins, with general negative basin-averaged KGE values (−0.14 (GLASS) to 0.03 (E2O-En)), indicating that all ET products had the limited overall performance for temporal scale. Relatively, the KGE > 0.2 mainly existed at Australia and Midwest America.

In terms of the climatological values of ET under dry and wet conditions (Figure 4), except GLASS under all conditions and MERRA-2 under wet condition, the ET products could reproduce the magnitudes of ET with 0.86 for PML≤ β ≤1.17 for ERA5-Land under dry condition and with 0.96 for MTE≤ β ≤1.07 for ERA5-Land under wet condition, which was consistent with the results presented in Section 3.1 (Table 3). In particular, most of the ET products underestimated the water balance ET above 1,200 mm (Figure 4), which mainly occurred in Amazonian Plain and Brazilian Plateau (Figure 2). As for γ , most of the ET products could generally detect the spatial variability for the climatological values of ET under dry and wet conditions, corresponding to a range of 0.69 (ERA5-Land)≤ γ ≤1.32 (GLASS) and 0.72 (ERA5-Land)≤ γ ≤1.07 (GLASS), respectively. Broadly, the spatial variability estimates of ET under dry condition tended to be higher than those under wet condition (represented as γ _ d r y > γ _ w e t ) except ERA5-Land and MTE. Regarding R, the ET products had a high R-based ability to simulate spatial distribution of ET with R > 0.8 under dry and wet conditions. Meanwhile, the ET products could better represent the spatial distribution of climatological ET (except for GLASS, ERA5-Land, and MTE) under dry basins than wet basins (represented as R_dry > R_wet). As for KGE, ET products exhibited the high overall performance on climatological ET conditioned by aridity, especially generating the highest KGE values for GLDAS2.0-Noah (0.89) under dry condition and PML (0.94) under wet condition.

With β at ∼ 1.0, the magnitude of temporal variability of ET tended to be more easily simulated under dry condition, compared with wet condition (Figure 5). As for γ , the spatial variability of ET temporal variability was generally overestimated by ET products under all aridity conditions, with 0.87 for GLDAS2.0-Noah≤ γ ≤1.47 for ERA5-Land under dry condition, and 0.79 for GLASS≤ γ ≤1.95 for PML under wet condition. As for R, the ET products could detect the spatial distribution of ET CV under dry basins, of which the highest R value was 0.76 for MERRA-2, followed by 0.69 for E2O-En. However, under wet condition, the ET products presented a contrasting performance, compared with dry condition, with R values ranging from −0.07 to 0.12. Considering KGE, similar to R, the ET products could not simulate the ET CV under wet condition, whereas, under dry condition, ERA5-land, MERRA-2, GLEAM3.5a, E2O-En and GLDAS2.0-Noah showed better overall performances, generating a KGE above 0.40.

Taking the ET trend into consideration (Figure 6), more than 50% of the total number of basins were located in the first and third quadrants, with 46.89% for ERA5-Land ≤ TPR + TNR≤78.00% for MERRA-2 under dry condition and 54.20% for GLASS ≤ TPR + TNR≤63.26% for ERA5-Land under wet condition. This indicates that most of the ET products can capture the ET trend directions. Despite that, it is worth noting that FPRs outweighed the TNRs under wet condition. This suggested that under the wet condition, these products tended to change the negative ET trends to the positive ET trends. Based on β , under wet condition, most of the ET products (except ERA5-Land and MERRA-2) tended to underestimate the magnitude of ET trend, with −19.63 for GLASS≤ β ≤-4.93 for GLEAM-3.5a. By contrast, under dry condition, the underestimations of the ET trend got relieved, with −2.06 for ERA5-Land≤ β ≤4.76 for GLASS, except that general underestimations still existed in dry condition. In addition, ET products underestimated the extreme ET trends over the wet basins (<−5 and >5 mm yr⁻¹), which mainly occurred in the Amazonian Plain and Brazilian Plateau (Figure 2). As for the spatial variability of ET trend, the γ values were around zero for all ET products under wet condition, ranging from −0.09 (GLEAM-3.5a) to 0.18 (ERA5-Land), while the γ values were more deviated from optimal value (1.0) for most ET products under dry condition, ranging from −23.21 for GLEAM-3.5a to 23.53 for E2O-En. Overall, all ET products exhibited limited R-based ability to simulate spatial distribution of the ET trend, with 0.09 for GLASS ≤ R ≤ 0.65 for MERRA-2 under dry condition and 0.05 for MTE ≤ R ≤ 0.33 for GLASS under wet condition. Furthermore, the overall performance for each ET product under all aridity conditions was poor with −23.25 (GLEAM-3.5a)≤KGE ≤ 0.44 (MERRA-2) under dry condition and −19.66 (GLASS)≤KGE ≤ −1.05 (ERA5-Land) under wet condition. Notably, the overall performances were generally worse for the latter.

Temporally, as shown in Figure 7A, regarding β , except PML, GLDAS2.0-Noah, and MTE, the magnitudes of ET were overestimated at 51.33% for GLEAM-3.5a to 98.89% for GLASS of dry basins, and at 51.66% for GLEAM-3.5a to 96.46% for MERRA-2 of wet basins. The basin-averaged β values for the ET products (except GLASS and MERRA-2) were both near to 1.0 under dry and wet conditions (Figure 7B). The basin-averaged γ values under dry condition were also close to 1.0 for most products, while under wet condition the values were overwhelmingly low, ranging from 0.17 for MTE to 0.72 for MERRA-2. As for R (Figure 7C), more than 50% of basins exhibited a value over zero for most products under all conditions. Despite that, each ET product showed a higher R-based ability to simulate ET temporal fluctuation under dry condition than wet condition, with average R values ranging from 0.15 for GLDAS2.0-Noah to 0.46 for MERRA-2 under dry condition and −0.03 for MERRA-2 to 0.17 for ERA5-Land under wet condition. As for KGE (Figure 7D), compared with the R-based ability, the overall performance of ET products under wet conditions worsened. For example, 64.12% (ERA5-Land) to 93.88% (MTE) of basins showed negative KGE values under wet conditions, whereas 63.11% (GLDAS2.0-Noah) to 87.33% (E2O-En) of basins exhibited positive KGE values under dry conditions. Furthermore, all ET products showed a negative average KGE value under wet conditions, while most products, except for GLASS, showed a positive average KGE value under dry conditions.

From perspective of climatological ET, the magnitude and spatial variability of ET could be represented by most of the ET products across all vegetation conditions (Figure 8), with both β and γ around 1.0. However, most of the ET products (excluding GLASS) also underestimated the ET values above 1,200 mm under LAI-III condition, which mainly exist in Amazonian Plain and Brazilian Plateau (Figure 2). Concerning R, the capacity to simulate the spatial distribution of climatological ET increased first, and then decreased as vegetation became greener for most ET products except GLDAS2.0-Noah, In terms of KGE, most ET products show good KGE-based performance. In addition, GLASS, ERA5-Land, MERRA-2, E2O-En, PML, and MTE showed that the KGE-based performance was the best under LAI-II condition.

In terms of the ET CV (Figure 9), most ET products (except GLASS and MTE) reasonably estimated ET magnitude under LAI-I condition, with 0.86 for PML≤ β ≤1.34 for ERA5-Land. However, the β values were limited for other vegetation conditions, with 0.22 for MTE≤ β ≤0.62 for PML under LAI-II condition and 0.14 for PML≤ β ≤0.67 for MERRA-2 under LAI-III condition. The β values for the ET temporal variability decreased as the vegetation turned green for each ET product. And the γ values for the spatial variability of ET temporal variability tended to be overestimated under all vegetation conditions. For R, all the ET products (except PML) had the limited R-based ability to simulate the spatial distribution of ET temporal variability, with vegetation greening. For example, the R values under LAI-I, LAI-II, LAI-III conditions ranged from 0.24 to 0.69, 0.24 to 0.38, and −0.09 to 0.13, respectively. Similar trends were occurred to KGE, except that the overall performance of KGE was even worse than that of R-capacity.

In the view of ET trend (Figure 10), its condition is similar to the aridity condition. The ET products could hit the ET trend directions, with 50.73% for PML ≤ TPR + TNR79.88% for MERRA-2 under LAI-I condition, 52.58% for GLASS ≤ TPR + TNR≤66.39% for PML under LAI-II condition, and 52.09% for PML ≤ TPR + TNR≤66.85% for ERA5-Land under LAI-III condition. Additionally, FPRs outweighed the TNRs for ET products (except ERA5-Land and MERRA-2) under LAI-II and LAI-III conditions, for example, for GLDAS2.0-Noah, FPR versus TNR was 39.27% versus 8.82% under LAI-II condition, and 40.96% versus 7.32% under LAI-III condition, indicating that the ET products tended to misidentify the negative ET trends as positive ET trends. Based on β , except ERA5-Land and MERRA-2, the ET products tended to seriously underestimate the magnitudes of ET under LAI-III condition, with −8.21 for GLASS≤ β ≤-3.15 for MTE. And the values of β were much larger than 1.0 under LAI-II condition (excluding ERA5-Land and MERRA-2), suggesting that the overestimation occurred in LAI-II condition. As for γ , all the ET products underestimated the spatial variability of the trends (excluding MERRA-2, GLEAM-3.5a, E2O-En and MTE for LAI-I condition). As for R values, the ET products (except GLASS, MTE and PML) showed lower correlations with the greening of vegetation. Interestingly, the ET trends were remarkably overestimated by most products in LAI-II condition, and slightly underestimated under LAI-I and LAI-III conditions. As for KGE, most of the ET products had bad performance with negative values under all conditions. Especially under LAI-II and LAI-III conditions, they had almost no simulability.

Temporally, the basin-averaged β values were around the 1.0 for all vegetation conditions (Figure 11A), though the temporal magnitudes of ET were either overestimated or underestimated by the ET products. Considering γ (Figure 11B), the basin-averaged values for all the ET products significantly decreased with the vegetation turning green, and were overestimated under LAI-I condition, but underestimated under the other vegetation conditions. It is worth noting that as vegetation was getting greener, the R-based ability for all the ET products was significantly constrained (Figure 11C). Specifically, all the ET products consistently performed, and the average R value and the basin percentages of the R values over zero decreased with vegetation greening. Figure 11D clearly shows that, like R-based ability, the basin-averaged overall performances of all the ET products decreased, as the vegetation was getting greener, except that the KGE values were lower than R values.

In this study, the simulations of ET derived from the eight methods were evaluated by the water balance ET of global 1,381 basins under various water, energy, and vegetation conditions. Since water, energy, and vegetation are crucial for accurately simulating ET, the lack of sufficient their information, caused by the lack of ET algorithm, forcing data and calibration methods, affects the performance of ET simulation (Xu et al., 2019; Elnashar et al., 2021; Li et al., 2022; Yu et al., 2022). As is shown, the comprehensive performance of ET products (Figures 7, 11) and the capture of ET variance (Figures 5, 9) regularly decrease, with the humidity and vegetation greenness increasing. These phenomena imply that the accuracy of the ET simulations may decrease, when the regional climate is wetting and the global vegetation is greening (Mankin et al., 2017; Lian et al., 2021; Zhang et al., 2022). Additionally, the ET products tend to misidentify the negative trends as the positive trends, especially under wet and LAI-III conditions, implying that the estimates of ET trends may be overestimated across the globe or in wet and LAI-III conditions (Figures 6, 10). These issues will be further discussed in the following.

In terms of the impact of water and energy denoted by AI, ET process in dry or wet regions can be conceptualized as a water- or energy-limited process, respectively: ET under dry conditions is water-limited, in that it is constrained by the soil moisture available for ET, while ET under wet conditions is energy limited, since there is sufficient soil moisture available for ET. Therefore, the maximum rate and temporal variations of ET proceeds are determined by atmospheric water demand (potential evapotranspiration) rather than soil moisture (Draper et al., 2018). All the ET products could better capture the mean annual value of all aridity conditions. However, the ET CV in wet basins tend to be more remarkably underestimated than in dry basins, by the ET products except GLASS and PML (Figures 5, 7). Indeed, wet zones have more active land-atmosphere coupling than dry zones, in that the inevitable ET algorithm errors or data forcing errors magnify the uncertainties under wet zones. For instance, Penman-Monteith method (GLDAS2.0-Noah, MERRA-2 and PML) is primarily driven by net radiation (R_n) under wet zones using a linearized approximate solution (Gao, 1988; Grignon, 1992; Leca et al., 2011), which is sensitive to low vapor pressure deficit (VPD) and may induce considerable problems in the extreme conditions (such as the water balance ET higher 1,200 mm (Figure 4) and extreme ET trends (Figure 6) and the soil evaporative term (Bai and Liu, 2018; Blatchford et al., 2020). More importantly, the presence or absence of ET products TWSC components in simulating ET under dry and wet areas cannot be ignored. However, most ET methods do not have an aquifer storage component, and LSMs lack a good representation of groundwater withdrawal for agricultural depletion, such as irrigation (Liu et al., 2016; Zeng and Cai, 2018). Additionally, the errors in the ET estimates and differences among the ET products are also mainly dependent on various inputs (Li et al., 2018).

The surface variables also control the ET process, especially vegetation (Wang et al., 2022; Zheng et al., 2022). Similarly, the response of ET products to vegetation was investigated. Regarding the mean annual value, we found that the simulability of datasets first increased and then decreased, with the increase of vegetation density (Figure 8), in line with the Lu et al. (2021). In addition, we also confirmed that the comprehensive performance (KGE) of ET products decreases as the vegetation is getting greener (Figures 8, 10). The first reason for this is that whether ET algorithms take the LAI or vegetation dynamics into consideration. For example, GLEAM-3.5a model lacks vegetation-related information, though it considers the vegetation optical depth, which may result in lower accuracy in high vegetation regions (Martens et al., 2017; Xu et al., 2019; Qiu et al., 2022). Another reason is that the ET algorithm do not comprehensively consider the vegetation process in hydrology or energy cycle. MERRA-2 overestimates the interception loss fraction defined as the fraction of rainfall, i.e., rainfall intercepted by the canopy and reevaporating back into the atmosphere without infiltrating into the soil or causing surface runoff (Reichle et al., 2011; Bosilovich et al., 2017; Gelaro et al., 2017; Reichle et al., 2017; Hinkelman, 2019), which could explain why the MERRA-2 generally has the highest β under various LAI conditions among the eight ET products (Figures 8, 11, and Lv et al. (2020)). The last easily neglected issue is related to the model forcing data. One aspect of the issue is the forcing data errors. The accuracy of LAI dataset is impacted by the leaf shadowing (Mehrez et al., 1992), especially tall and dense vegetations. Besides, shaded leaves are not light-saturated, leading to diffuse sunlight conditions and then having a higher fraction of FAPAR (the fraction of photosynthetically active radiation absorbed by the canopy) (Jimenez et al., 2011; He et al., 2013; Xu et al., 2019). Another aspect of the issue is the forcing data settings, for instance, MTE ET product was generated from machine learning method by compiling the 253 globally distributed flux towers data and remote sensing data, including vegetation information (FAPAR). We speculated that the varying performance of MTE product with various LAI conditions was probably driven by data settings. For example, the vegetation was used to do split not regression, which results in inadequate vegetation information (Jung et al., 2010). Or ERA5-Land was used to generate land elements data including ET, by using a static monthly climatology of a fixed land use and leaf area index (LAI) (Muñoz-Sabater et al., 2021). And GLDAS2.0-Noah also uses a static land use, though with high spatial resolution (Rodell et al., 2004). Therefore, they ignored the change of land cover and cities, and lost more frequent LAI anomalies during the reanalysis period (Muñoz-Sabater et al., 2021).

The model calibration methods also have a significant impact on the performance of ET simulation. One problem concerning the methods is that the ET simulations are often calibrated with the mean annual value not the variance and trend of actual ET, though considering multiple calibration metrics. Another problem is that the data used for calibration are often EC site data that are not representative of the regional scale (Bai and Liu, 2018; Xu et al., 2019). In addition, as far as we know, the ET products except GLASS and MTE are accompanied by component data such as soil evaporation, vegetation evapotranspiration and water surface evaporation, but these data are not calibrated with sufficient actual measurements (Swanson, 1994; Brunel et al., 1997; Chen et al., 2014).

The uncertainties in Pre and TWSA products are the largest source of uncertainties in assessing the ET products (Liu et al., 2016). According to the water balance budget, the assessment of global-scale ET products needs to rely on grid-scale Pre and TWSA products, although the global-scale observatory data is difficult to collect. As for the three Pre products selected in this study (GPCC, CPC-Unified, and CRU TS4.05), the uncertainties derive from the number of stations used, the time homogeneity and the quality control procedures (Trenberth et al., 2014; Sun et al., 2018). However, these products are interpolated from an unprecedented number of station data and are the most reliable precipitation products currently available (Sun et al., 2018). Regarding the TWSA data (GRACE-REC and GRID-CSR-GRACE-REC), the uncertainties arise mainly from the models used for the reconstruction (pre-2002) and the driving data (Gyawali et al., 2022). However, the correlation of GRACE-REC with yearly streamflow anomalies have median value of around 0.60 over 1981–2010 (Humphrey and Gudmundsson, 2019); the GRID-CSR-GRACE-REC has high correlation with Global Mean Sea level with R of 0.91 (Li et al., 2021). We further investigated the uncertainties in water balance evapotranspiration defined as the CV of the six Pre-TWSC-Q combinations, and found that most of the basins with uncertainties of <0.10 and uncertainties above 0.10 were located mainly in the Midwest USA and Southwest China (rainfall gauges are more sparsely distributed in high mountain areas) and the Arctic (Figure 12). In addition, Q data may be affected by human harvesting of deep groundwater and inter-basin water transfers (Liu et al., 2016). However, TWSC can reasonably take into account the impact of human activities on Q. Moreover, in validating the model, only the terrestrial water balance (not the atmospheric water balance) is considered, and the measured evapotranspiration values lack the cross-validation to further reduce the error with the true values (Li et al., 2019). The generalizability of our results to other regions of the world may be subject to additional uncertainty, as the basins included in this study do not cover the entire globe. However, it is important to note that the performance of evapotranspiration products varies with dryness and vegetation greenness, and it is necessary to ensure that all types of dryness and vegetation greenness are covered (Figure 1). To minimize errors caused by different spatial resolutions, all ET products were re-interpolated linearly to 0.5° before evaluation. Furthermore, our analysis is based on observed ET using the water balance method, which represents the average ET of watersheds controlled by hydrological stations, reducing the uncertainty caused by a single grid point to some extent. The scale effect on ET product performance related to aridity and vegetation greenness response needs further exploration in future research.