This study was conducted at the Haibei National Field Research Station, Qinghai, China (37°37′ N, 101°19′E), which is located on the Northeastern QTP at an elevation of 3,200 m MASL (Figure 1). This area is characterized by a plateau continental monsoon climate, with well-developed seasonally frozen ground. The average annual air temperature is −1.7°C, with the maximum monthly temperature in July (10.1°C) and the minimum monthly temperature in January (−15°C). The annual precipitation is ~580 mm, of which 80% falls in the growing season (i.e., from May to September), leading to high water content (close to field capacity) in the soil during the growing season, thus, the evapotranspiration during the growing season was not limited by soil water content (Dai et al., 2021). The average annual pan evaporation is ~1,191.4 mm (Zhang et al., 2018). The soil type around the lysimeter system is classified as Mat-Gryic Cambisol, which belongs to clay loam, and has a thickness of ~60–80 cm (Dai et al., 2019a), and the basic soil property is shown in Table 1. The grass crop is dominated by perennial sedge and graminoid species, including Kobresia humilis, Stipa aliena, and Elymus nutans, with ~8–15 cm in height, which together constitutes 60–80% of plant cover around the lysimeter system (Dai et al., 2021).

The actual evapotranspiration (ET_a) was measured using large-scale weighing lysimeters (height 2 m, diameter 1 m, and resolution 10 g) installed in the central alpine Kobresia meadow of the Haibei National Field Research Station. Lysimeters are vessels containing both monolites and repacked soils and are embedded completely in the soil at their top level with the soil surface. The topsoil, including the root, was composed of monolites, while deep soil, not including roots, comprised repacked soils. The soil was cut off from the parent soil at the base of the lysimeters, and the lower lysimeter boundary was exposed to atmospheric pressure. Measurements of ET_a were based on changes in weight, with measurements at 30-min intervals recorded by a data logger (CR1000, Campbell, USA) and converted to daily values. To ensure data quality, all negative or abnormal ET_a values caused by falling soil particles were discarded; outliers that differed more than three times from average ET_a, which yielded 393 days of data spanning June 2017–December 2018. According to in situ phenological observations of the dominant plant foliage, we defined the growing season as the period from May 1 to September 30, while the period from October 1 to April 30 of the following year was defined as the nongrowing season (Zhang et al., 2018).

Given that the soil moisture of the root region was more than the field capacity (Table 1), thus, the ET_a in this study was mainly not limited by surface moisture. Furthermore, the plant height was 8–12 cm, actively growing, and completely shading the ground, which satisfies the definition of reference evapotranspiration (defined as 8–15 cm high, uniform, actively growing, green grass that completely shades the soil and not limited by soil water availability), and we, thus, could assume that the ET_a in this study could be considered as reference crop evapotranspiration or potential evaporation. To verify the assumption that ET_a in this study was not limited by soil water availability, our previous studies showed that the annual Priestley–Taylor coefficient values (ratio of ET_a to ET_eq) ranged from 1.08 to 1.14 (ET_eq represents the minimum possible evapotranspiration from a moist surface and depends only on air temperature and available energy), and the annual decoupling coefficient Ω ranged from 0.43 to 0.48 (Zhang et al., 2018), which provided direct evidence that ET_a in this study occurred under strongly energy-limited conditions rather than soil water availability constraints.

All meteorological variables needed to calculate ET₀ using various models were obtained or estimated from the weather station at Haibei Station, including precipitation (PPT), relative humidity (RH), wind speed (WS), net radiation (R_n), total radiation (R_s), soil heat flux (G), maximum air temperature (T_max), vapor pressure deficit (VPD), minimum air temperature (T_min), and mean air temperature (T_a). PPT was collected using a PPT gauge (52203, RM Young, USA) at 0.5 m height. Radiation was measured by four radiometers (CNR4, Kipp & Zonen, Netherlands) at 1.5 m height; RH, WS, and T were measured at 1.5 m height (HMP45C, Vaisala, Finland), and WS was converted to 2 m height (u2=uz4.87ln (67.82-5.42)) for calculating ET₀. The G was measured using three heat flux plates (HFT-3, Campbell, USA), which were buried 5 cm beneath the surface. Half-hourly means of meteorological data were stored using a data logger (9210 XLITE, Sutron, USA). The BREB method parameter was determined using an eddy covariance system installed near a lysimeter system, which included a three-dimensional ultrasonic anemometer (CSAT3, Campbell, USA) and an open-path infrared CO₂/H₂O gas analyzer (LI-7500A, LI-Cor, USA), and 30-min fluxes were calculated and logged with a SMARTFLUX system (LI-COR, USA) (Zhang et al., 2018). The soil water content was measured using the oven-drying method at depths of 0–10, 10–20, 20–30, 30–40, and 40–50 cm through a soil auger in July because the evapotranspiration reaches its maximum in July with its maximum water demand.

The reference evapotranspiration was defined as 8–15 cm high, uniform, actively growing, green grass that completely shades the soil and not limited by soil water availability (Doorenbos and Pruitt, 1977). A total of fourteen often-used ET₀ models were selected for comparison, including BREB, three combination models (1963 Penman, FAO-24 Penman, and FAO-56 PM), seven radiation models [Priestley–Taylor, DeBruin–Keijman, Makkink, Makkink (1957), Makkink (1967), IRMAK1, and IRMAK2], and three temperature-based models (Hargreaves, Hargreaves1, and Hargreaves2), to compare their performance using the lysimeter measurement. The specific equations and parameters of the models are listed in Table 2.

In this study, the ET_a measured by the large-scale weighing lysimeters and the performances of the ET₀ models were compared with these lysimeter system estimates on a daily basis. Pairwise comparisons were conducted using a general linear regression. For further comparison, the root means squared error (RMSE), percentage error of estimate (PE), mean absolute error (MAE), and coefficient of determination (R²) were used to evaluate the ET₀ models. The RMSE, PE, MAE, and R² are defined as follows:

where P_i is the predicted value; O_i is the observed value; P¯ and O¯ are the averages of P_i and O_i, respectively; and n is the total number of data points.

To achieve the best comparison between the models and measurements, we need to select the dominant meteorological factors affecting the measured ET_a. Given that it may not be appropriate to explore results based solely on the coefficient of independent variables in multiple regression analysis, owing to the strong collinearity and nonlinearities among meteorological factors, we adopted a boosted regression tree (BRT) model to quantitatively evaluate the relative influences of meteorological variables on measured ET_a. BRT is a machine-learning method based on the classification regression tree algorithm (CART). This method generates multiple regression trees through random selection and self-learning methods, which can improve model stability and prediction accuracy. During operation, some data are randomly selected many times to analyze the influence of independent variables on dependent variables and to quantitatively evaluate the relative effect of independent variables on dependent variables, while the remaining data are used to test the fitting results (Elith et al., 2008). Most importantly, the BRT can evaluate the relative influence of an independent variable on a dependent variable, without transformations, and can cope well with nonlinear relationships. Furthermore, the BRT exhibits good performance in dealing with stronger collinearity and nonlinearities. Thus, the BRT was adopted to evaluate the individual influences of the controlling factors on the measured evapotranspiration. All statistical analyses were performed using R software version 3.03 (R Development Core Team, 2006), and all figures were plotted using Origin 9.0.

The measured ET_a showed a clear seasonal pattern, and the growing season ET_a was significantly higher than that in the nongrowing season (P < 0.05) (Figure 2a). The average measured daily ET_a during the study period was 2.33 mm day⁻¹, with an average daily measured ET_a of 4.14 and 0.65 mm day⁻¹ during the growing season and nongrowing season, respectively (Figure 2a). Environmental variables showed a similar seasonal pattern, with maximum and minimum values in the growing and nongrowing seasons (Figure 2). The average daily R_n, R_s, T_a, VPD, and RH during the growing season were 9.53 MJ m⁻², 18.50 MJ m⁻², 10.15°C, 0.33 kPa, and 74.88%, respectively (Figure 2). The average daily R_n, R_s, T_a, VPD, and RH during the nongrowing season were 3.28 MJ m⁻², 12.65 MJ m⁻², −3.67°C, 0.23 kPa, and 57.16%, respectively.

A comparison of fourteen evapotranspiration equations against the lysimeter measurements is presented in Figure 3, Table 3, showing that the relationship between daily ET₀ calculated by the ET₀ models and lysimeter measurements ET_a was all significant (P < 0.01), with high coefficients of determination (R²) ranging from 0.59 to 0.86. For the combination models, FAO-24 Penman (FAO-24 Pen) yielded the highest correlation (0.77), followed by PM-56 (0.76) and Penman 63 (Pen-63) (0.75). For radiation-based models, PT obtained the highest correlation (0.80), followed by DK (0.79), Makkink (1967) (0.69), Makkink (1957) (0.69), IRMAK1 (0.66), Makkink (0.65), and IRMAK2 (0.62). For temperature-based models, HAR obtained the highest correlation (0.62), and HAR1(0.59) and HAR2 (0.59) obtained the same correlation. The daily estimates of combination models generally underestimated the ET_a values measured by lysimeter, with MAEs of −0.26 to −0.01 mm day⁻¹. However, the radiation-based models (except PT and IRMAK1) and temperature-based models generally overestimated the ET_a values, with MAEs ranging from −0.14 to 0.50 mm day⁻¹ for radiation-based models and 0.40–0.59 mm day⁻¹ for temperature-based models.

Throughout the study period, the RMSE of the BREB method was 0.98, and the RMSE of combination models ranged from 1.19 to 1.27 mm day⁻¹ and averaged 1.22 mm day⁻¹. Furthermore, the RMSE of FAO-56 PM increased from 1.22 to 1.29 mm day⁻¹ as aerodynamic resistance (r_s) changed from 20 to 60 s m⁻¹ (Figure 4). The RMSE for radiation-based models ranged from 1.02 to 1.42 mm day⁻¹ and averaged 1.27 mm day⁻¹, and the RMSE for temperature-based models ranged from 1.47 to 1.48 mm day⁻¹ and averaged 1.47 mm day⁻¹. Based on the RMSE, the performance of the fourteen evapotranspiration models decreased in the following order: BREB (0.98) > PT (1.02) > DK (1.03) > Pen-63 (1.19) > FAO-24 Pen (1.19) > PM-56 (1.27) > IRMAK1 (1.32) > Makkink (1957) (1.34) > Makkink (1967) (1.36) > Makkink (1.38) > IRMAK2 (1.42) > HAR (1.47) > HAR1 (1.47) > HAR2 (1.48). The best model was 34% and 30% more accurate than the poorest (HAR2) and the commonly used FAO-56 PM equation, respectively. Furthermore, Pen-63 and FAO-24 Pen demonstrated better performance than the commonly used PM-56 equation. Overall, for the entire study period, the BREB yielded the best performance, followed by the combination, radiation-based, and temperature-based models.

During the growing season, the daily ET₀ calculated by fourteen evapotranspiration models was significantly correlated with the lysimeter measurements ET_a (P < 0.01), with R² ranging from 0.32 to 0.64 (Figure 5, Table 4). Of the combination models, PM-56 had the highest R² (0.60), followed by FAO-24 Pen (0.59) and Pen-63 (0.58). Of the radiation-based models, PT and DK obtained the highest R² (0.59), followed by IRMAK1 (0.50), Makkink (1967) (0.47), Makkink (1957) (0.47), IRMAK2 (0.43), and Makkink (0.40). Notably, Makkink (1967) and Makkink (1957) had the same R² (0.47). Of the temperature-based models, HAR, HAR1, and HAR2 obtained the same R² values (0.32). Interestingly, all models (except HAR1 and HAR2) generally underestimated ET_a during the growing season, values of MAE ranged from −1.10 to −0.15 mm day⁻¹, with PM-56 having the largest underestimation (26.59%) and HAR the minimum underestimation (3.56%) (Table 4).

The RMSE of BREB was 1.31, the RMSE for combination models ranged from 1.38 to 1.58 mm day⁻¹ and averaged 1.47 mm day⁻¹, the RMSE for radiation-based models ranged from 1.19 to 1.55 mm day⁻¹ and averaged 1.40 mm day⁻¹, and the RMSE for temperature-based models ranged from 1.42 to 1.43 mm day⁻¹ and averaged 1.42 mm day⁻¹ (Table 4). Based on the RMSE values, the performance of the fourteen evapotranspiration models followed the order: DK (1.19) > PT (1.22)> Makkink (1967) (1.28) > BREB (1.31) > Pen-63 (1.38) > HAR1 (1.42) > HAR2 (1.42)> HAR (1.43) > FAO-24 Pen (1.46) > Makkink (1.47) > IRMAK1 (1.54) > Makkink (1957) (1.55) > IRMAK2 (1.55) > PM-56 (1.58). Evidently, the best DK was 25% more accurate than the poorest (FAO-56). Overall, for the growing season period, the BREB yielded the best performance, followed by the radiation-based, temperature-based, and combination models.

During the nongrowing season, the daily ET₀ calculated by the fourteen evapotranspiration equations was also significantly correlated with the lysimeter measurements ET_a (P < 0.01), but with lower coefficients of determination (R²) ranging from 0.17 to 0.64 (Figure 6, Table 5). Of the combination models, FAO-24 Pen obtained the highest R² (0.28), followed by Pen-63 (0.25) and PM-56 (0.21). Of the radiation-based models, PT and DK obtained the highest R² (0.34), followed by Makkink (1967) (0.27), Makkink (1957) (0.27), Makkink (0.26), IRMAK1 (0.26), and IRMAK2 (0.23). Among the temperature-based models, HAR had the highest R²-value (0.19). Interestingly, all models (except BREB) generally overestimated the ET_a values measured by lysimeter during the nongrowing season, with MBEs ranging from 0.40 to 1.20 mm day⁻¹ and averaging 0.78 mm day⁻¹; Makkink (1967) yielded the largest underestimate (by 185.83%) and PT the minimum underestimate (by 61.06%) (Table 5).

The RMSE of BREB was 0.52, the RMSE for combination models ranged from 0.86 to 1.00 mm day⁻¹ and averaged 0.92 mm day⁻¹, the RMSE for radiation-based models ranged from 0.80 to 1.44 mm day⁻¹ and averaged 1.12 mm day⁻¹, and the RMSE for temperature-based models ranged from 1.47 to 1.54 mm day⁻¹ and averaged 1.50 mm day⁻¹. Based on the RMSE, the ET₀ model performance decreased in the following order: BREB (0.52) > PT (0.80) > DK (0.84) > FAO-24 Pen (0.86) > PM-56 (0.90) > Pen-63 (1.00) > IRMAK1 (1.07) > Makkink (1957) (1.11) > Makkink (1.29) > IRMAK2 (1.30)> Makkink (1967) (1.44) > HAR (1.47)>HAR1 (1.51) > HAR2 (1.54). Evidently, the best model was 66.24% more accurate than the poorest model (HAR2). Overall, for the nongrowing season period, the BREB yielded the best performance, followed by the combination, radiation-based, and temperature-based models.

The fourteen evapotranspiration model estimations were consistent with the pattern observed in the lysimeter measurements (Figure 7), with a peak in July. As already noted, the combination models and BREB underestimated the measurements from May to September and overestimated those in the other months (Figures 7A,B). The radiation-based models underestimated the measurements from June to September and overestimated those in the other months (Figure 7C). However, the temperature-based models generally overestimated the measured evapotranspiration during most months (except for July and August) (Figure 7D). Overall, all models tended to underestimate the measured ET_a during the growing season (with larger evaporative demand) and overestimated ET_a during the nongrowing season (with reduced evaporative demand).

The BRT model indicated that R_n was the dominant factor controlling the seasonal variation in measured ET_a throughout the study period, accounting for 69.02% of the total variability, followed by VPD (7.13%), T (6.75%), RH (5.73%), R_a (5.15%), R_s (4.33%), and WS (1.85%) (Figure 8a). During the growing season, R_n remained the main control of seasonal variation in measured ET_a (Figure 8b), accounting for 44.30% of the total variability, followed by T (14.12%), R_s (12.56%), R_a (8.34%), RH (7.80%), VPD (6.87%), and WS (6.02%). However, the seasonal variation in the measured ET_a in the nongrowing season was dominated by RH (Figure 8c), accounting for 27.99% of the total variability, followed by R_n (20.99%), R_a (12.96%), VPD (12.56%), T (10.18%), R_s (9.39), and WS (5.94%) (Figure 8).

Given that the evapotranspiration in our study was not limited by soil water conditions and the plant height was 8–12 cm, actively growing, and completely shading the ground, which is close to the definition of ET₀, we thus assumed that ET_a was comparable to ET₀ during the growing season, with the aim of selecting the best fit model over a humid alpine meadow on the northeastern QTP to estimate ET₀. Previous studies have shown that the Penman family models are generally the most accurate when evaluating ET₀ across various climate scenarios and regions (Liu et al., 2017). Of the Penman models for ET₀, the PM-56 has been considered the standard equation for estimating evapotranspiration (Allen et al., 1998). For instance, Yoder et al. (2005) found that PM-56 displayed the best performance in the humid southeast United States. López-Urrea et al. (2006) tested seven evapotranspiration equations using lysimeter observations in a semiarid climate and found that the PM-56 equation was the most precise method compared with other evapotranspiration equations. Unlike previous studies, the PM-56 model was not the best in our study; we found that Pen-63 and FAO-24 Pen were more accurate (Pen-63 and FAO-24 Pen had smaller RMSE than PM-56 throughout the study period). Similar results have been reported in other studies (Berengena and Gavilán, 2005; Martel et al., 2018). A more recent study also reported the poor performance of PM-56 when compared with data from 20 FLUXNET towers (Ershadi et al., 2014). These results suggest that PM-56 might not be the only standard model for evaluating ET₀ because it did not yield better accuracy than the other Penman models. Given the better performance of Pen-63 and FAO-24 than PM-56 in this study, we may apply old Penman family models to our study region, especially considering that the PM-56 requires many meteorological inputs, which limits its use in areas with sparse data, especially in harsh environments (Tabari et al., 2012). Overall, the poor performance of PM-56 may be attributed to the higher aerodynamic resistance (r_s), and there is increasing evidence indicating that PM-56 underestimation is related to the fixed r_s = 70 s m⁻¹ in the equation, which may be too large (Liu et al., 2017). This result was also confirmed by our results that the values of RMSE increased from 1.22 to 1.29 mm day⁻¹ as r_s changed from 20 to 60 s m⁻¹, and the RMSE was nearly unchanged (from 1.12 to 1.13 mm day⁻¹) when r_s varied between 0 and 20 s m⁻¹ (Figure 4). Therefore, reducing the value of r_s from 70 to 0–20 s m⁻¹ can improve daily PM-56 estimates. Other studies also found that r_s should be a variable rather than a fixed value (Allen et al., 2006; Liu et al., 2017). For instance, r_s should be smaller when the result is underestimated and should be larger when it is overestimated (Ventura et al., 1999).

For the performance comparison of radiation models against lysimeter measurements, we found that the PT models yielded the best performance of the radiation-based models, which was consistent with a previous study conducted in humid areas where the PT method yielded a good accuracy estimate for ET₀ (Ershadi et al., 2014). There is increasing evidence indicating that the input parameters are the dominant factors affecting their performance (Lang et al., 2017), and we conclude that the better performance of the PT models might be associated with the most important meteorological factors affecting ET, such as R_n, which is supported by our results (Figure 8). Compared with the PT method, other radiation models that use R_s as the main driving variable may overestimate ET₀ due to the reduction in R_s through atmospheric reflection in this region (Zhang et al., 2018). Furthermore, each model was developed based on its specific underlying surface and climatic conditions. For instance, the PT model was established in a humid climate condition, which is suitable for a humid alpine meadow in the northeastern QTP. Most importantly, the PT models require fewer meteorological inputs than the combination models. However, given that the ET_a was mainly limited by available radiation (i.e., R_n-G), the structures of the PT and DK models both included available radiation items. Combining these factors, we recommend the PT and DK models for use in humid alpine meadows in the northeastern QTP, especially when considering the difficulty in obtaining evapotranspiration in this harsh climate.

For the performance comparison of temperature models against lysimeter measurements, a previous study reported that the Hargreaves equation is one of the simplest empirical methods used for ET₀ estimation because of its lower meteorological data input, including some meteorological data required in the standard PM-56 model (Jensen et al., 1990). To further select the best Hargreaves version equation, we compared the performance of the original (HAR) and two modified versions (HAR1 and HAR2) and found that the original HAR model had the lowest error (RMSE = 1.47 mm day⁻¹, MAE = 0.40 mm day⁻¹, and PE = 17.37%), which was consistent with previous studies conducted in humid regions (Tabari, 2010) but contrasted to studies conducted in arid regions in which the modified Hargreaves equation displayed a more accurate estimation of evapotranspiration (Ravazzani et al., 2012). Overall, the temperature models displayed poor performance compared with radiation models because the Hargreaves method was established in semiarid areas (Tabari, 2010). Furthermore, the structure of temperature models missing the most important parameter (i.e., R_n) results in poor performance compared with radiation models. Therefore, a local calibration is required to improve the Hargreaves method accuracy in nonarid regions.

By comparing the four model types, we found that the BREB yielded the best performance, followed by the combination, radiation-based, and temperature-based models (Table 3). Overall, most radiation-based models underestimated the measured ET_a throughout the study period, whereas the temperature-based models tended to overestimate ET_a. This was consistent with previous studies in which the Makkink and PT models generally underestimated ET_a (Priestley and Taylor, 1972; Xu and Singh, 2002), while the Hargreaves equations often overestimated ET_a in cold-humid conditions and required local calibration (Berti et al., 2014). Given that our study region was a humid alpine meadow, ET_a tended to be overestimated. An alternative explanation for the poor Hargreaves model performance in humid regions may be that the Hargreaves method was established in semiarid areas (Tabari, 2010), and the R_a parameter used in the Hargreaves model structure was based on the maximum possible radiation value and does not consider the atmospheric transmissivity. However, the atmospheric transmissivity in humid regions is affected by many factors, such as atmospheric moisture; thus, the solar radiation reaching the surface is significantly reduced because of the high atmospheric moisture content (Temesgen et al., 1999), resulting in the overestimation of solar radiation, ultimately leading to an overestimation of evapotranspiration using the Hargreaves method.

Furthermore, there were common features in all four groups of models. All the models tended to underestimate the measured ET_a during the growing season (with larger evaporative demand) and overestimated ET_a during the nongrowing season (with reduced evaporative demand), which was consistent with a previous study conducted in a semiarid region (Liu et al., 2017). The underestimated measured ET_a during the growing season may be related to the ET_a in alpine meadows under strongly energy-limited conditions rather than soil water content during the growing season; thus, higher solar radiation could lead to a higher ET_a during the growing season (Zhang et al., 2018). Therefore, both the Hargreaves equations and other models require further local or regional calibration before being applied to a given region (Xu and Singh, 2002). It should also be noted that the data used in this study were obtained from 1 year and a single weather station, which may be insufficient to represent the whole humid climate or the alpine ecosystem but represent only a humid alpine meadow in the northeastern QTP. Thus, a longer period and more lysimeter systems should be used in the alpine ecosystem in the future to obtain more accurate estimates of evapotranspiration over humid alpine meadows in the northeastern QTP.