Validation and Comparison of Seven Land Surface Evapotranspiration Products in the Haihe River Basin, China

Abstract

Evapotranspiration (ET) is an important part of the surface energy balance and water balance. Due to imperfect model parameterizations and forcing data, there are still great uncertainties concerning ET products. The validation of land surface ET products has a certain research significance. In this study, two direct validation methods, including the latent heat flux (LE) from the flux towers validation method and the water balance validation method, and one indirect validation method, the three-corned hat (TCH) uncertainty analysis, were used to validate and compare seven types of ET products in the Haihe River Basin in China. The products evaluated included six ET products based on remotely-sensed observations (surface energy balance based global land evapotranspiration [EB-ET], Moderate Resolution Imaging Spectroradiometer [MODIS] global terrestrial evapotranspiration product [MOD16], Penman–Monteith–Leuning Evapotranspiration version 2 [PML_V2], Global Land Surface Satellite [GLASS], global land evaporation Amsterdam model [GLEAM], and Zhangke evapotranspiration [ZK-ET]) and one ET product from atmospheric re-analysis data (Japanese 55-year re-analysis, JRA-55). The goals of this study were to provide a reference for research on ET in the Haihe River Basin. The results indicate the following: (1) The results of the six ET products have a higher accuracy when the flux towers validation method is used. Except for MOD16_ET and EB_ET, the Pearson correlation coefficients (R) were all greater than 0.6. The root mean square deviation (RMSD) values were all less than 40 W/m². The GLASS_ET data have the smallest average deviation (BIAS) value. Overall, the GLEAM_ET data have a higher accuracy. (2) When the validation of the water balance approach was used, the low values of the MOD16_ET were overestimated and the high values were underestimated. The values of the EB_ET, GLEAM_ET, JRA_ET, PML_ET, and ZK_ET were overestimated. According to the seasonal variations statistics, most of the ET products have higher R values in spring and lower R values in summer, and the RMSD values of most of the products were the highest in summer. (3) According to the results of the uncertainty quantification based on the TCH method, the average value of the relative uncertainties of the GLEAM_ET data were the lowest. The relative uncertainties of the JRA_ET and ZK_ET were higher in mountainous areas than in non-mountainous area, and the relative uncertainties of the PML_ET were lower in mountainous areas. The performances of the EB_ET, GLEAM_ET, and MOD16_ET in mountainous and non-mountainous areas were relatively equal. The relative uncertainties of the ET products were significantly higher in summer than in other periods, and they also varied in the different sub-basins.

1. Introduction

Terrestrial evapotranspiration (ET) consists of evaporation from soil and canopy interception, as well as transpiration from vegetation [1]. It has importance for surface energy and water balance, and has profound implications for climate change [2,3]. The types of ET products include physical process models [4,5,6], energy balance models [7,8], atmospheric re-analysis [9], and machine learning or integration models [10]. The rapid development of remote sensing satellite-based models in the past few decades has made it possible to capture the spatial and temporal variability of ET from regional to continental and even global scales. However, these remote sensing satellite-based models usually present different degrees of uncertainties depending on their theories, structural assumptions, and input parameterization, and those limitations are mainly governed by changes in landscape, climatological, and hydrological conditions [11]. RS-ET products are susceptible to spatial heterogeneity and meteorological complexity, so they have different applicability in different regions. Due to the quality and quantity of atmospheric data, atmospheric re-analysis products also have uncertainties. The uncertainties of the models and input data of the ET products affect the accuracy and application scope of the products. Therefore, the accuracies of ET products need to be evaluated.

Reliable information on ET is very important for various bio-geophysical applications, including those of the forestry, hydrological, agricultural, and meteorological disciplines [12]. The accurate quantification of ET is challenging because its variations are controlled by complex interactions among soil moisture availability, atmospheric feedback, and heterogeneous vegetation conditions [13]. The current land surface ET product accurate evaluation methods usually consist of direct validation and indirect validation [14]. Direct validation is based on the use of in situ measurements to obtain the ground truth values. This method can be used as the primary and reliable method of validating ET data and is usually employed at the pixel-scale and regional-scales [15,16]. In situ measurements, including the Bowen ratio, lysimeters, laser isotopes, and eddy co-variance, can be used as the relative true values in pixel-scale evaluations to evaluate ET products. Regional evaluations are based on the water balance method [17,18,19,20,21] and the multi-site scale expansion method [22]. In the absence of ground truth ET data, indirect validation becomes feasible, which includes the cross-checking method, multi-scale evaluation method based on high-resolution remote sensing data, and spatiotemporal change trend analysis method. Many ET products have good test results at specific observation stations, but since observation stations cannot cover all surface types, it is not clear whether they can be applied to other surface types [23]. The cross-validation method focuses on evaluating the relative accuracy of different surface evapotranspiration products and the consistency of temporal and spatial change trends. In particular, the cross-validation of multiple ET products has been further improved after the introduction of quantified methods, such as triple collocation (TC) [24] and three-cornered hat (TCH) [25,26,27,28,29,30]. TC and TCH are innovative and reliable approaches for estimating the error variance of various time-series products. The generalized TCH method allows for a relative comparison of at least three datasets based on their respective uncertainties without the need of a priori knowledge of their uncertainties [31].

The above methods have been applied in various basins. Previous studies have usually focused on the accuracy validation of only one ET model or product. For example, Zhou et al. [32] used the evapotranspiration data calculated using the Penman–Monteith (PM) formula to verify the SEBAL model in the Xilin River Basin, proving that the estimation results have a good accuracy. Liu et al. [33] used eddy correlation observations to verify the SEBAL model in the growing season in the Liaohe Delta region of China, and their validation results indicate that the SEBAL model is suitable for wetland evapotranspiration research. Studies have also been conducted on the applicability and comparison of several ET products. For example, Velpuri et al. [15] evaluated the accuracy of the MOD16 and operational simplified surface energy balance (SSEBop) products in the continental United States, and their results indicate that both products are suitable for most land cover types and that both ET products are reliable at the basin scale. Chao et al. [34] comprehensively evaluated the accuracy of five ET products (GLEAM, Penman–Monteith–Leuning (PML), MODIS, process-based land surface evapotranspiration/heat fluxes algorithm (P-LSH), and multi-tree ensemble (MTE)) based on ground-based and Gravity Recovery and Climate Experiment (GRACE) satellite observations over the continental United States. Their results showed that the P-LSH and GLEAM products were consistent with the recon data in the middle-value range. The MODIS and MTE had larger average deviation (BIAS) and root mean square deviation (RMSD) values on the yearly scale, and the MODIS and MTE datasets tended to underestimate and overestimate the ET values over the entire value range, respectively. Wartenburger et al. [35] evaluated the abilities of a range of monthly-scale ET products at the global scale and applied cluster analysis to identify the spatiotemporal differences between all of the datasets using factors such as the model choice, the meteorological forcing used to drive the assessed models, the data category, the ET scheme, and the number of soil layers in the model. Their results indicate that the model selection is mostly the dominant factor. Xie et al. [36] evaluated the applicability of three ET products using the water balance method in the Northwest Inland River of China, and their results revealed that the GLEAM data have the best fit in most of the alpine basins, and overestimation generally occurred in the ERA data in summer and autumn and in the MODIS data from October to January.

ET is affected by land use, soil moisture availability, and climate conditions [37], which makes it highly variable across heterogeneous landscapes [38]. All ET products have uncertainties [39,40], and it was suggested that combining artificial intelligence algorithms or data-driven algorithms and physical process algorithms will further improve the accuracy of ET estimation algorithms and the quality of ET datasets, as well as enhancing their capacity to be applied in different climate regions [34].

Although there is research that evaluated the ET products across China, the conclusion was that all the global ET products were unable to reasonably reproduce the ET time-series in most basins at the basin scale [41]. The Haihe River Basin is one of China’s main river basins, however, comprehensive comparative studies of multiple ET products are lacking. In this study, we mainly selected remotely-sensed ET (RS-ET) products for evaluation in this study. The principle of re-analysis data is different from that of remote sensing data, which combine ground station observations, satellite remote sensing data, and numerical simulation data to best reflect the atmospheric conditions. Therefore, a re-analysis data product was selected for comparison with the remote sensing products in this study. Direct and indirect evaluation methods were used to evaluate the applicability of ET products in the Haihe River Basin. The purposes of this study were: (i) to assess the performances of seven ET products in the Haihe River Basin using EC observation and the water balance method, and (ii) to systematically compare the relative uncertainties of the seven ET products.

2. Materials and Methods

2.1. Study Area

The Haihe River Basin is one of China’s main river basins. It is located in a semi-humid and semi-arid zone in the northern part of China (35–43°N, 112–120°E) (Figure 1). The Beijing–Tianjin–Hebei Economic Circle, which is one of China’s three major economic circles, is located in this basin. Geographically, the Haihe River Basin is bordered by Bohai Bay to the east, the Taihang Mountain to the west, the Yellow River to the south, and the Mongolian Plateau to the north. The total area of the basin is about 320,600 km². The Haihe River Basin has a temperate monsoon climate, with an annual average temperature of 8–12 °C and an annual average precipitation of 539 mm [42].

Using digital elevation model data (DEM) as the input data, the geospatial hydrologic modeling extension (HEC-GeoHMS) hydrological model was adopted to extract the sub-basins [43]. The Haihe River Basin was divided into 2605 small basins, and then, according to the data obtained from the Hydrological Yearbook of the People’s Republic of China, the small basins were merged into seven sub-basins: the Luan River Basin (R1), Chaobai River and Jiyun River Basin (R2), Inland River Basin (R3), Daqing River Basin (R4), Ziya River Basin (R5), South Canal Basin (R6), and Tuhai River and Majia River Basin (R7).

The land cover in the study area can be classified into eight classes, including Cultivated Land, Forest, Grassland, Artificial Surfaces, Water Bodies, Wetland, Shrubland, and Bare Land (Figure 2), and the area percentage of the above land-cover types accounted for 53.6%, 18.67%, 18.16%, 7.68%, 1.37%, 0.25%, 0.21%, and 0.06%, respectively.

2.2. Datasets

The data used in this study included ET product data, eddy correlation (EC) observation data from flux towers, DEM data, terrestrial water storage data, runoff data, and precipitation data.

2.2.1. The Evaluated ET Products

Seven ET product datasets were comprehensively evaluated in this study: the surface energy balance based global land evapotranspiration (EB-ET in this paper) [2,44], GLEAM (GLEAM_ET in this paper) [4,45], MOD16 (MOD16_ET in this paper) [5], Penman–Monteith–Leuning Evapotranspiration version 2 (PML_V2) global evapotranspiration and gross primary production (PML_ET in this paper) [6], land evapotranspiration of Global Land Surface Satellite products (GLASS_ET in this paper) [46], a continuous satellite-derived global record of land surface evapotranspiration product created by Zhang Ke (ZK_ET in the paper) [47], and land evapotranspiration of the JRA_55 (JRA_ET in this paper) [9]. The ET products included a single-source model ET product (EB_ET), ET products based on the Penman–Monteith formula (MOD16_ET, PML_ET, and ZK_ET), an ET product based on the Priestley–Taylor (PT) formula (GLEAM_ET, ZK_ET), an integrated model product (GLASS_ET), and a re-analysis product (JRA_ET). A summary of these ET products used in the study are presented in

Table 1. Information about ET product datasets evaluated in this study.

Product Name	Approach	Temporal Resolution	Spatial Resolution	Units
EB_ET	SEBS, Re-analysis Data	Daily, monthly	0.05° × 0.05°	mm
GLEAM_ET	Modified Priestley–Taylor, soil stress factor	Daily, monthly	0.25° × 0.25°	mm
JRA_ET	Re-analysis	6 h, Monthly	0.5612° × 0.56162°	W/m2
MOD16_ET	Modified Penman–Monteith-Leuning, surface conductance model	8 day, monthly	1 km × 1 km	mm
PML_ET	Modified Penman–Monteith-Leuning	8 day	500 m × 500 m	mm
ZK_ET	Modified Penman–Monteith for canopy and soil, Priestley-Taylor for water	Monthly	0.1° × 0.1°	mm
GLASS_ET	Integrated	Daily every 8 days	1 km × 1 km	W/m2

.

2.2.2. Other Datasets

The EC observation data from the flux towers (green dots in Figure 1) used for the ET product validation at the pixel-scale were obtained from the Observations Network of the Qinghai–Tibet Plateau Science Data Center (https://data.tpdc.ac.cn/zh-hans/, accessed on 17 January 2020) and the China Flux Observations Data Alliance Network (http://www.chinaflux.org/, accessed on 17 January 2020). The flux data included 30-min scale and daily-scale data.

Table 2. Information about the flux towers used in this study.

EC Site Name	Land Cover	Region	Period	Longitude (°E)	Latitude (°N)
Yucheng	Wheat/corn	Shandong	2003–2010	116.57	36.82
Huailai	Corn	Hebei	2013–2017	115.79	40.35
Guantao	Wheat/corn, cotton	Hebei	2008–2010	115.13	36.52
Miyun	Orchard	Beijing	2008–2010	117.32	40.63
Daxing	Wheat/corn, orchard	Beijing	2008–2010	116.43	39.62

presents the information about the five flux towers. The main land-cover type at the sites of the five flux towers was agricultural land. The latent heat flux of the EC data from the five flux towers were used for the direct validation of the ET products.

The SRTM DEM data from the United States Geological Survey (https://lpdaac.usgs.gov/, accessed on 22 February 2018) was used to divide the study area into mountainous and non-mountainous areas based on an elevation threshold of 200 m. The spatial resolution of the DEM data was 90 m × 90 m. The terrestrial water storage data, runoff data, and precipitation data were used in the ET product validation using the water balance method. The monthly-scale data from 2003 to 2012 were used in this study.

The Gravity Recovery and Climate Experiment (GRACE) was adopted to represent the terrestrial water storage change (TWSC) [48]. The values in each pixel represent the equivalent water height, which is an anomaly based on the mean value from January 2004 to December 2009. The RL06M version of the monthly-scale data was used (http://www2.csr.utexas.edu/grace/, accessed on 9 May 2021). The spatial resolution of the data was 0.5° × 0.5°.

The runoff data (https://geonetwork.nci.org.au, accessed on 27 May 2021) used in this study were a global gridded runoff product (Linear Optimal Runoff Aggregate, LORA) [49], with a spatial resolution of 0.5° × 0.5° at monthly timescales for the period of 1980–2012, and includes time-variant uncertainty. This product merges runoff estimates from hydrological models constrained with observational streamflow records. The product broadly agrees with published runoff estimates at many river basins and represents the seasonal runoff cycle for most of the globe. The LORA data were validated using the monthly observed discharge of five hydrological stations in the Haihe River basin in 2003 [50]. We needed to convert the observed discharge data into runoff depth data, which were consistent with the LORA data, according to the upstream drainage area of the hydrological station. The locations of the five hydrological stations were showed in Figure 1, including Shixiali station (SXL, with drainage area of 2.39 × 10⁴ km²) in Sanggan River, Xiangshuibao station (XSB, with drainage area of 1.45 × 10⁴ km²) in Yang River, Yanchi station (YC, with drainage area of 4.37 × 10⁴ km²) in Yongding River, Zhangjiafen station (ZJF, with drainage area of 0.85 × 10⁴ km²) in Bai River, and Xiahui station (XH, with drainage area of 0.53 × 10⁴ km²) in Chao River. The precipitation and temperature data used in this study were obtained from the Monthly Data of Surface Precipitation and Temperature in China dataset (http://data.cma.cn/, accessed on 16 January 2018), which was developed by the basic ground meteorological data construction project. The data have been updated since 1961, and this meteorological dataset is commonly used in many studies. Spatial kriging interpolation was performed on the data with missing pixels to ensure the integrity of the data. The spatial resolution was 1 km × 1 km.

The surface net radiation data used for the analysis of ET influencing factors is the GLASS-NR product (http://www.glass.umd.edu/NR/, accessed on 27 March 2021), with a spatial resolution of 0.05° × 0.05° and a temporal resolution of 1 day.

The land-cover type data used in this study was GlobeLand30 (http://www.globallandcover.com/, accessed on 8 August 2021) from 2010, which is the first global geo-information public product provided by China to the United Nations, with the spatial resolution of 0.5° × 0.5°.

2.3. Methods

Direct evaluation methods, including EC observation and the water balance method, were used to evaluate the applicability of the ET products in the Haihe River Basin. The applicability evaluation index adopts the BIAS, RMSD, Pearson correlation coefficient (R), and Taylor plots. Then, the TCH uncertainty quantification method, as an indirect evaluation method, was used to assess the performances of the ET products. In order to explore the influencing factor of ET, the spatial distribution of ET under different land-cover types and the cosine similarity between the temporal variation of ET products and meteorological data were adopted.

2.3.1. ET Product Validation Using Flux Tower EC Observations

Taking the multi-year observation data from five flux stations in the Haihe River Basin from 2003 to 2010 as the true values, the 30-min scale flux tower EC observation data with no more than 30% of the data missing were selected for the screening and daily-scale synthesis [51]. The ZK_ET product does not include daily-scale data, so it was not evaluated in this section. The MOD16_ET product provides the total evapotranspiration over a period of 8 days (5 or 6 days at the end of the year) in units of 0.1 mm/8 days, 0.1 mm/6 days, or 0.1 mm/5 days. In order to ensure the comparability of the data, the units were uniformly converted into the standard latent heat flux units, W/m². Through comparison of the pixel values of the EC products and the values of the EC observations, the six types of ET products (EB-ET, GLEAM_ET, JRA_ET, MOD16_ET, PML_ET, and GLASS-ET,) were verified. The applicability evaluation index adopts the BIAS, RMSD, Pearson correlation coefficient (R), and normalized central root mean square deviation (RMSD’) for the Taylor plots. Taylor plots are used to present an accuracy assessment, using metrics such as the Pearson correlation coefficient (R), standard deviation (STD), and RMSD’. The above evaluation indexes are calculated as follows:

$$BIAS=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N\left(X_i-Y_i\right),$$

(1)

$$RMSD={\left(\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(X_i-Y_i\right)}^2\right)}^{\frac{1}{2}},$$

(2)

$$R=\frac{{{\displaystyle \sum }}_{i=1}^N\left(X_i-X_{ave}\right)\left(Y_i-Y_{ave}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left(X_i-X_{ave}\right)}^2 {{\displaystyle \sum }}_{i=1}^N{\left(Y_i-Y_{ave}\right)}^2}},$$

(3)

$$RMSD’=\sqrt{RMS{D}^2-BIA{S}^2},$$

(4)

$$STD={\left(\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(X_i-X_{ave}\right)}^2\right)}^{\frac{1}{2}},$$

(5)

where N is the number of sample groups, $X_i$ and $Y_i$ are the evapotranspiration value to be verified and the true value of the relative evapotranspiration in group $i$, respectively, and $X_{ave}$ and $Y_{ave}$ are the mean values of $X_i$ and $Y_i$ in the N groups of data, respectively.

2.3.2. ET Product Validation Using Water Balance Method

The water balance principle is a method of calculating the actual evapotranspiration in a closed watershed [52]. The calculation formula is:

$$WB\_ET\left(i\right)=P\left(i\right)- R\left(i\right) - TWSC\left(i\right),$$

(6)

where $P\left(i\right)$ is the average precipitation in the basin (mm/month), $R\left(i\right)$ is the average runoff depth in the basin (mm/month), and $TWSC\left(i\right)$ is the monthly water storage change of the underlying surface of the basin (mm/month). The $TWSC\left(i\right)$ data were calculated from the GRACE data using the following equation:

$$TWSC\left(i\right)\approx \frac{TWSA\left(i+1\right)-TWSA\left(i-1\right)}{2∆i},$$

(7)

where $TWSA\left(i\right)$ is the anomalous value of the water storage in a certain month, $i$ denotes the month, and the $∆i$ is the monthly unit. This formula is helpful in reducing the noise.

Since the GLASS_ET is composed of single-day evapotranspiration data on an 8-day scale, there are large errors in the monthly-scale synthesis. In this study, the other six ET products were compared and analyzed on the monthly scale using the water balance method.

The period of the data used was 2003–2012. Among them, the monthly scale standard products of the GLEAM_ET, JRA_ET, MOD16_ET, and ZK_ET were directly used. The daily-scale data of the EB_ET and PML_ET were synthesized into monthly-scale data. In this study, the water balance method was used to analyze the product accuracies of these six ET products on the annual scale and in different seasons.

2.3.3. ET Product Comparison Using TCH Uncertainty Quantification Method

The three-cornered hat method (TCH) is a method for evaluating multiple sets of data sequences, which was developed in the last century. Its methodology was proposed and further developed by Premoli and Tavella [26,27]. Assuming that there are N groups of different observation sequences {$x_i$} (i = 1,…, $N$), each group of $x_i$ is the result of adding the true value and the measurement error, that is, $x_i$= $x_{true}$+ $e_i$. The traditional three-cornered hat method is limited to the case of $N$ = 3 [25], while the generalized three-cornered hat method can be used for the case of N > 3.

First, any observation sequence is selected as the reference values, and then the difference sequence $y_i$ between the remaining observation sequences and the reference values is obtained:

$$y_i=x_i-x_r=e_i-e_r \left(i=1,\dots ,N-1\right).$$

(8)

An interpolation sequence matrix $Y$ and difference sequence co-variance matrix S are constructed:

$$Y=\left[\begin{array}{cccc}{y}_1& y_2& \cdots & y_{\left(N-1\right)}\end{array}\right].$$

(9)

$$S=cov \left(Y\right).$$

(10)

An N × N noise co-variance symmetric matrix R and auxiliary matrix J are introduced:

$$J_{N-1,N}=\left[\begin{array}{cccc}1& 0& ⋮& -1\\ 0& 1& \cdots & -1\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & -1\end{array}\right].$$

(11)

$$R=\left[\begin{array}{cccc}{r}_{11}& r_{12}& ⋮& r_{1N}\\ r_{12}& r_{22}& \cdots & r_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ r_{1N}& r_{2N}& \cdots & r_{NN}\end{array}\right].$$

(12)

$$S=J·R·J^T.$$

(13)

By combining these matrices, a set of equations $r_{ij}=s_{ij}-r_{NN}+r_{iN}+r_{jN} (i,j<N)$ can be obtained. At this time, the number of equations is less than the number of unknowns, so the equations cannot be solved [26]. In order to ensure the positive definiteness of the co-variance matrix, constraint condition H is used as the constraint condition [26]:

$$H_2(r_{1N},\dots ,r_{NN})=-\frac{H_1(r_{1N},\dots ,r_{NN})}{K}<0,$$

(14)

$$K=\sqrt[N-1]{\left|S\right|},$$

(15)

$$H_1(r_{1N},\dots ,r_{NN})=\frac{\left|R\right|}{\left|S\right|}=r_{nn}-\left[r_{1n}-r_{nn},\dots r_{\left(n-1\right)n}-r_{nn}\right]S^{-1}{\left[r_{1n}-r_{nn},\dots r_{\left(n-1\right)n}-r_{nn}\right]}^T.$$

(16)

Introducing the squared and mean squared values of the off-diagonal elements of the noise matrix $\sqrt{\frac{1}{n}{\displaystyle \sum }_{i<j}^{n-1}{r}_{ij}{}^2}$ ensures that the global correlation of all of the datasets satisfies $\left|R\right|>0$ under the smallest premise. Thus, the following equation can be obtained:

$$F(r_{1N},\dots ,r_{NN})=\frac{1}{K^2}·{\displaystyle \sum }_{i<j}^N{r}_{ij}{}^2.$$

(17)

In order to ensure that the initial value is within the constraint conditions, the initial value of the above equation iteration is as follows:

$$r_n^0=0 for i<N,$$

(18)

$$r_{nn}^0=\frac{1}{2s^{*}},s^{ *}=\left[1,\dots ,1\right] S^{-1}{\left[1,\dots ,1\right]}^T,$$

(19)

Equation (17) is minimized under the constraints in order to obtain the optimal solution. The mean squared value of the diagonal elements corresponding to the symmetric matrix R is the uncertainty of the corresponding input dataset.

In this paper, the TCH method is used to calculate the relative uncertainties of the monthly-scale ET products, and the period is 2003–2012. The relative uncertainties of six products, including EB_ET, GLEAM_ET, JRA_ET, MOD16_ET, PML_ET, and ZK_ET, are calculated.

2.3.4. Cosine Similarity

Cosine similarity [53] is a measure of similarity between two datasets. The cosine of two sets can be derived by the Euclidean dot product formula:

$$cos\left(A,B\right)=\frac{{{\displaystyle \sum }}_{i=1}^n{A}_i\times B_i}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(A_i\right)}^2}\times \sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(B_i\right)}^2}}$$

(20)

where n is the number of observations, Σ is the summation symbol, A_i is the A value for observation i, and B_i is the B value for observation i. In this paper, cosine similarity is used to study the similarity between the temporal variation of evapotranspiration products and meteorological data.

3. Results and Analysis

3.1. Results of ET Product Validation Using Flux Tower EC Observations

The pixel scale ET product validation depends on the Latent Heat Flux (LE) data observed by EC of the flux station. Taking Yucheng flux station as an example, Figure 3 shows the scatter diagram of the daily scale LE data from 2003 to 2010. The data shows a relatively regular temporal change trend. In a year, the LE gradually increased from January to June, then decreased since June and began to rise around the beginning of July, and reached a peak in August, then decreased gradually from August to December. The farmland where Yucheng station was located adopts the Winter Wheat–Summer Maize Rotation System. The data information of “key points for the integrated management of winter wheat and summer corn rotation” issued by the Ministry of agriculture and rural affairs of the People’s Republic of China shows that the winter wheat sown in the previous year will be harvested in mid-June, and the summer maize will be sown in late June. Therefore, the emergence of the sudden drop of LE was caused by crop harvesting.

The ET product pixel value at the location of the flux tower was extracted and compared with the latent heat flux observations. According to the results of the flux tower validation (Figure 4), the ET products generally overestimated the low value and underestimated the high values. There is a variation of 0.48–0.77 in R. Except for MOD16_ET and EB_ET, the R values were all greater than 0.6. According to the calculated BIAS values, with a variation 2.93–14.27mm, the GLASS data were less biased. The RMSD ranged from 29.8mm to 39.3mm, with the lowest value coming from PML_ET.

The evaluation results of the ET validation conducted using the EC observations are presented as a Taylor diagram (Figure 5). The red origin represents the actual observation data. The horizontal and vertical axes represent the STD. The blue radial line represents the R value. The green dashed line represents the RMSD’ values. The results indicate that GLEAM_ET had the highest R value and the lowest RMSD’.

Due to the uneven distribution of the observation sites in China, the data were very limited; and because the observation footprint of the flux site itself was affected by the meteorological environment, the observation range was irregular, and the uncertainty caused by it cannot be overcome. Therefore, there are still errors in the evaluation results, which need to be combined with other validation methods.

3.2. Results of ET Product Validation Using Water Balance Method

The LORA verification result was shown in Figure 6. The Nash–Sutcliffe Efficiency (NSE), BIAS, and R value are 0.35, 0.35, and 0.76 respectively, which indicate that the accuracy of the LORA runoff dataset in the Haihe River Basin is not very good, while the BIAS and R value were acceptable.

In order to show the changes in the water-balance-derived ET (WB_ET in this paper) and the six types of ET products in the different months, the multi-year average monthly evapotranspiration of the data was calculated. The statistical results are shown in Figure 7. The results indicate that the six ET products had an obvious seasonal pattern. The evapotranspiration values were slightly lower in winter and autumn and were relatively high in spring and summer. During most parts of the study period, all six ET products captured the seasonal ET estimations and corresponded well with the WB_ET. The evapotranspiration values of the JRA_ET, PML_ET, and ZK_ET were significantly higher than those of the WB_ET. The highest value of the WB_ET occurred in July, while the highest values of the ZK_ET and MOD16_ET occurred in August. The EB_ET and GLEAM_ET were in good agreement with the WB_ET.

The monthly runoff depth data were usually less than 6 mm, and this part of the runoff in LORA has been overestimated, which may lead to the low ET value that was estimated by the water balance method.

Figure 8 presents the validation scatter plots for the comparison of the six ET products with the WB_ET, and Figure 9 presents the results as a Taylor diagram. The results show that the six ET products were strongly correlated with the WB_ET. The low values of the EB_ET and MOD16_ET were overestimated, while the high values were underestimated in the Haihe River Basin. The validation points were divided by season. The results indicate that the six ET products exhibited obvious seasonal variation characteristics. As is shown in Figure 6, the JRA_ET and GLEAM_ET had higher R and lower RMSD’ values.

As is shown in

Table 3. Accuracy evaluation index values for the ET products.

Index	Season	EB_ET	GLEAM_ET	JRA_ET	MOD16_ET	PML_ET	ZK_ET
R	Spr.	0.77	0.8	0.84	0.82	0.81	0.79
	Sum.	0.47	0.46	0.6	0.31	0.23	0.3
	Fall	0.76	0.88	0.91	0.89	0.89	0.87
	Win.	0.48	0.52	0.46	0.81	0.56	0.33
	Yearly	0.85	0.91	0.92	0.82	0.90	0.88
RMSD (mm)	Spr.	19.83	16.36	19.68	9.82	27.3	13.11
	Sum.	14.53	22.72	27.73	25.8	53.28	30.22
	Fall	15.76	12.02	12.36	13.89	12.47	14.77
	Win.	3.85	7.74	12.81	16.43	7.19	5.68
	Yearly	16.88	14.3	17.56	16.01	28.21	17.08
BIAS (mm)	Spr.	17.32	13.29	17.15	−0.73	22.84	8.42
	Sum.	9.68	13.44	20.36	8.83	47.12	20.82
	Fall	−8.56	2.36	4.16	−3.49	5.07	8.87
	Win.	1.09	0.8	0.84	0.82	0.81	0.79
	Yearly	7.06	9.14	13.72	4.39	20.02	11.16
RMSD’ (mm)	Spr.	9.66	9.54	9.65	9.79	14.95	10.05
	Sum.	10.84	18.32	18.83	24.24	24.87	21.90
	Fall	13.23	11.79	11.64	13.44	11.39	11.81
	Win.	3.69	7.70	12.78	16.41	7.14	5.62
	Yearly	15.33	11.00	10.96	15.40	19.87	12.93

, most of the ET products had higher R values in fall and lower R values in summer. As for the RMSD and BIAS values, the difference between the two was quite large for the different products. Except for the EB_ET, most of the products had the highest RMSD values in summer. According to the BIAS value results, most of the BIAS values were positive, and the negative BIAS values between the EB_ET and MOD16_ET mainly occurred in autumn, indicating that the product values were generally higher than the WB_ET values.

3.3. Results of ET Product Comparison Using TCH Uncertainty Quantification

The relative uncertainties of the six products were calculated using the TCH method. The relative uncertainty is the ratio of the uncertainty to the mean ET value, and it can describe the model skill by neglecting the ET magnitude. The average values of the relative uncertainties of the EB_ET, GLEAM_ET, JRA_ET, MOD16_ET, PML_ET, and ZK_ET were 12.50%, 6.04%, 10.40%, 17.20%, 16.30%, and 13.50%, respectively. It can be seen that the relative uncertainties of the GLEAM data were significantly lower than those of the other products. According to the results presented in Figure 10, the higher relative uncertainty values of the different products occurred at different locations. For example, the higher values of the PML_ET data were mainly located in the southeast. The JRA_ET had a low spatial resolution, and its higher values were concentrated in the northwest. According to the division index of the basic landform morphology [54], the Haihe River Basin was divided into mountainous and non-mountainous areas. The results are presented in

Table 4. Relative uncertainties (%) of the six products in mountainous and non-mountainous areas.

	EB_ET	GLEAM_ET	JRA55_ET	MOD16_ET	PML_ET	ZK_ET
Non-mountainous area	12.39	6.69	9.51	17.48	19.15	12.04
Mountainous area	12.81	5.58	11.18	17.12	14.03	14.89

. According to the distribution of the relative uncertainties, the uncertainties of the EB_ET, GLEAM_ET, and MOD16_ET were relatively consistent in the mountainous and non-mountainous areas. The JRA_ET and ZK_ET had higher relative uncertainties in the mountainous areas; while the PML_ET had higher relative uncertainties in the non-mountainous areas.

The relative uncertainties in each sub-basin in the different months are shown in Figure 11. The ET products had slightly higher relative uncertainties in summer. The climate influencing factors of evapotranspiration vary greatly in summer, as well as with high temperature and precipitation simultaneously, which lead to a large range of evapotranspiration in summer, thereby resulting in the increased uncertainty of ET.

In terms of the sub-basins, there was not much difference between the different sub-basins. Except for the JRA_ET data, the highest or second-highest values of the other products all occurred in the Tuhai River and Majia River Basin (R7).

4. Discussion

4.1. Spatial Distribution of ET under Different Land-Cover Types

The spatial distributions of the mean annual ET from 2003 to 2012 are shown in Figure 12. The ranges of different ET products vary greatly. Generally, the coarser the spatial resolution of the ET products, the smaller the range of the ET values. According to the spatial distribution of ET, the evapotranspiration in the northwest of the Haihe River Basin is lower than that in the southeast.

The areas of wetland, shrubland, and bare land were relatively less in the study area, so the mean annual ET of the other five land-cover types from 2003 to 2012 were calculated in Figure 13. Generally, the ET value is the highest in water bodies and the lowest in artificial surfaces. The ET values of forest and grassland are higher than those of bare land and artificial surfaces. Forestlands have greater abilities at retaining incoming precipitation due to a deeper root zone and thus higher soil moisture storage, resulting in greater occurrences of ET compared to other land uses [55]. The ET value in cultivated land is affected by the growth of crops, which also leads to fluctuations in evapotranspiration throughout the year. Most of these features were reflected by GLEAM_ET, EB_ET, and JRA_ET, while the ET of all types of land cover were relatively high in PML_ET and low in MOD16_ET.

4.2. Analysis of the Influence of Climatic Factors on ET

A time-series of evapotranspiration data, precipitation data (P), air temperature data (T), and surface net radiation data (RN) at a monthly scale from 2003 to 2012 are displayed in Figure 14. It can be seen that the fluctuation of ET is similar to those of P, T, and RN, and the cosine similarity between climatic factors and ET products shows a very high similarity (

Table 5. Cosine similarity between climatic factors and ET products.

	EB_ET	GLEAM_ET	JRA_ET	MOD16_ET	PML_ET	ZK_ET
T	0.946	0.940	0.928	0.852	0.952	0.942
RN	0.976	0.971	0.975	0.889	0.963	0.937
P	0.932	0.944	0.938	0.935	0.950	0.948

), which showed that all of the ET products were closely related to the above three climatic factors. All the cosine similarity values were above 0.930, except for those between MOD16_ET and temperature/net radiation. Under the temperate monsoon climate of the study area, the abundant precipitation in summer increases the soil water storage in the region and the higher temperatures increase the soil water evaporation capacity. Therefore, the phenomenon of high summer evapotranspiration is consistent with this theory.

4.3. Sampling Resolution and the Mechanism of the ET Products

Currently, various remote sensing ET products are constantly being developed, greatly facilitating academic research, and the evaluation of these products is helpful in ensuring the appropriate application of these products. An applicability evaluation is an important method used to evaluate the errors of ET products. In the applicability evaluation of ET products, the flux tower validation method and the water balance validation method are the most commonly used methods. However, in terms of applicability evaluation, the flux tower validation’s errors cannot be ignored. This is due to the fact that: (1) Since the spatial scale of the EC observation data is one hundred meters and the observation range is irregular due to factors such as the wind speed, a more accurate method of calculating the flux footprint needs to be applied in future studies. (2) The spatial resolutions of the ET products evaluated in this paper range from 500 m to 50 km approximately, this mismatch increases the error in the flux tower validation method and weakens the comparability of the evaluations. (3) There are not many flux site data for the study area. The representativeness of the flux site data has also become a major problem in the flux tower validation method. At present, five agricultural surface-type flux tower EC observation datasets are available in the Haihe River Basin, but it is difficult for these five datasets to represent the entire river basin. It is difficult to use observational data to characterize the overall accuracy of the ET products.

As for the water balance method, the applicability of auxiliary data, such as water storage, runoff, and precipitation data, will affect the accuracy of the entire water balance formula. The spatial resolution of each product is different, and it is inevitable that some information will be lost during resampling. In addition, the different evaluation methods also have their own limitations. The EB_ET adopts the improved SEBS algorithm [44,56], which was found to have the smallest deviation in the water balance method evaluation conducted in this study, but it differs from the measured EC data. The results of the overestimation and underestimation of MOD16_ET were similar to previous research results for the Hanjiang River Basin [57] and the Northwest Inland River in China [36]. The GLEAM_ET, JRA_ET, PML_ET, and ZK_ET produced overestimation. The overestimations of the PML_ET, ZK_ET and JRA55 have also been confirmed in the Hanjiang River Basin in China [58].

The MOD16_ET, PML_ET, and ZK_ET are calculated based on the modified PM formula. The PM formula has a solid theory and high accuracy. It is a commonly used ET estimation method. According to relevant research conclusions [59], the PM formula has a high correlation with climate data. Therefore, the accuracy of the climate data determines the accuracy of the PM formula.

The PT formula is a revised version of the PM formula, and it omits the aerodynamic processes. Compared with the PM formula, the accuracy of the PT formula is somewhat lower; however, the results of this study demonstrate that the applicability of the PT model product (GLEAM_ET) in the Haihe River Basin is much higher than that of the PM data products, and the good applicability of GLEAM_ET in the Haihe River Basin has also been confirmed in other studies [60].

The JRA-55 is re-analysis data, and it was created using a complete observation system from 1958 to the present. However, the data from the observation system are not uniformly distributed geographically, so some of the characteristics of the Haihe River Basin cannot be observed.

The GLASS_ET is a Bayesian method integrated product, which integrates five traditional latent heat flux algorithms: the MODIS algorithm, improved PM model and the Priestley–Taylor Jet Propulsion Laboratory (PT-JPL), modified satellite-based Priestley–Taylor (MS-PT), and semi-empirical Penman algorithms. In this study, only the flux tower validation of the GLASS_ET was carried out, and due to the lack of monthly-scale data, the water balance method and TCH uncertainty validations were not applied to the GLASS_ET. The results indicate that this product has large errors in the high ET value seasons, especially in summer.

5. Conclusions

Seven ET products, including the EB_ET, GLEAM_ET, JRA_ET, MOD16_ET, PML_ET, GLASS_ET, and ZK_ET, were evaluated in this study using three methods: pixel-scale flux tower validation, regional-scale water balance method validation, and TCH uncertainty analysis. The applicability of the above data products in the Haihe River Basin was evaluated. The results of this study provide data support for the promotion of ET products and water resource evaluation in the Haihe River Basin. The main conclusions of this study are as follows:

• For validation based on flux tower EC observations, the results indicate that, except for MOD16_ET and EB_ET, the R values were all greater than 0.6. The BIAS values of GLASS_ET were the lowest (2.93 w/m²). GLEAM_ET had the highest R (0.77) value and the lowest RMSD’ (27.3 w/m²). Overall, the GLEAM data fitted the EC measured data well, and the RMSD values were relatively low;
• Based on the validation using the water balance method, the EB_ET, GLEAM_ET, JRA_ET, PML_ET, and ZK_ET overestimated the values in the Haihe River Basin. The multi-year averaged monthly evapotranspiration of the WB_ET and the ET products showed that the EB_ET and GLEAM_ET were in good agreement with the WB_ET. Overall, the JRA_ET (R = 0.92, BIAS = 13.72 mm/month, and RMSD′ = 10.96 mm/month) and GLEAM_ET (R = 0.91, BIAS = 9.14 mm/month, and RMSD′ = 11.00 mm/month) data have higher accuracies;
• For the uncertainty analysis based on the TCH method, the average relative uncertainties of the GLEAM data were significantly lower than those of the other ET products. The relative uncertainties of the JRA_ET and ZK_ET were relatively high in the mountainous areas. The PML_ET had higher relative uncertainties in the non-mountainous areas. The performances of the EB_ET, GLEAM_ET, and MOD16_ET in the mountainous and non-mountainous areas were relatively equal.

Overall, among the seven products, GLEAM_ET shows the best consistency with the point EC observations, better consistency with the basin-scale benchmark data, and the lowest relative uncertainties. We thus recommend GLEAM_ET as the preferred choice for application among the seven products in the Haihe River Basin. This study offers useful information for product users to choose the appropriate ET product(s) to conduct their specific studies and also helps product developers to improve the accuracy of their ET products.