Complementary relationship estimation of actual evapotranspiration in extreme cold and arid areas: a case study of the Hotan River Basin, northwest China

Actual evapotranspiration (ET_a) is a key variable in the energy, water, and carbon cycles of terrestrial ecosystems (Vinukollu et al. 2011; Ma et al. 2015a; Crago et al. 2016; Aminzadeh & Or 2017; Szilagyi et al. 2017). Accurate estimates of ET_a are needed to address a wide range of problems such as global climate change, water balance computation, agricultural management, and hydrological modeling (Xu & Singh 2005; Crago et al. 2016). Evapotranspiration is the major water consumption pathway in the water budget, accounting for more than 60% of annual precipitation on a global scale (Oki & Kanae 2006) and more than 90% in extreme arid areas (Glenn et al. 2007). ET_a exhibits strong heterogeneity because of the wide spatial variability of precipitation, ground conditions, and vegetation types (Mallick et al. 2015). The harsh natural conditions in extreme cold and arid areas make it difficult to obtain sufficient hydrological and meteorological data in situ. Despite the urgent need for it, estimating ET_a rates and analyzing the temporal and spatial distribution features of ET_a with limited observation data in extremely cold and arid basins is a tremendous challenge.

Although ET_a can be directly observed by the eddy covariance method, weighing lysimeter method, energy balance Bowen ratio method, or the scintillometer method (Wang & Dickinson 2012; Ma et al. 2015b), these all require sufficient manpower, materials, and financial resources. Several approaches to the estimation of ET_a have been developed to date, some of which were derived from the land surface energy balance equation using remotely sensed surface temperatures (e.g., Bastiaanssen et al. 1998; Su 2002; Mariotto et al. 2011), while others, such as soil-vegetation-atmosphere transfer models, are physically based methods (e.g., Haverd et al. 2011) and some are derived from the Penman (1948) or Priestley & Taylor (1972) equations (e.g., Monteith 1965; Allen et al. 1998; Monteith & Unsworth 2013). The complementary relationship (CR) between ET_a and potential evaporation (ET_p) was first proposed by Bouchet (1963). Several models based on the CR have been suggested during the past several decades, such as the advection-aridity (AA) model developed by Brutsaert & Stricker (1979), the complementary relationship areal evapotranspiration (CRAE) model of Morton (1983), the Granger model described by Granger & Gray (1989), as well as some nonlinear normalized CR models (Han et al. 2012, 2014; Brutsaert 2015; Crago et al. 2016; Szilagyi et al. 2016, 2017). The CR models are intellectually appealing since they only use routine meteorological data and bypass detailed soil moisture, land cover type, and terrain and vegetation information (Xu & Singh 2005; Ma et al. 2015b; Szilagyi et al. 2016). The CR approaches have been widely applied over different underlying conditions, spatial scales, and temporal scales (e.g., Xu & Singh 2005; Huntington et al. 2011; Wang et al. 2011; Han et al. 2012, 2014; Crago & Qualls 2013; Mallick et al. 2015; Hobbins et al. 2016; McEvoy et al. 2016). However, one criticism of the CR is the lack of definitive derivation and physical constraints (Crago et al. 2016; Szilagyi et al. 2016), although some studies have attempted to rectify this (Szilagyi 2007; Han et al. 2012, 2014; Brutsaert 2015; Aminzadeh et al. 2016; Crago et al. 2016; Szilagyi et al. 2016, 2017). In addition, there is a strong argument between symmetric and asymmetric CR under different settings (Szilagyi 2007; Jaksa et al. 2013; Ma et al. 2015b; Crago et al. 2016). Han et al. (2012) proposed boundary conditions (BCs) in the CR formulation and a new nonlinear function satisfying these BCs. Through comparison of the nonlinear CR developed by Han et al. (2012), the CRAE, and the Penman–Monteith approaches, Ma et al. (2015b) suggested that the nonlinear CR performed well with limited observation data for local calibration of the model parameter values. Moreover, they found that this CR model had an obvious advantage in both theory and practice in that it did not require users to select a symmetric or asymmetric CR relationship. Brutsaert (2015) reformulated the CR as a general polynomial based on strictly physical considerations. Szilagyi et al. (2016) compared the original version of the generalized CR of Brutsaert (2015) with a revised version in which the BCs-(ii) of Brutsaert (2015) was replaced. Furthermore, Crago et al. (2016) elaborated why these BCs-(ii) needed to be replaced, and introduced a new concept, the maximum ET_p (denoted by E_pads in the literature), to normalize ET_p as x_min. They also proposed a new version of the CR based on rescaling, which not only reduced scatter of ET_a estimates, but led to formation of a proposed self-adjusting CR model. Szilagyi et al. (2017) redefined the maximum ET_p as an invariable result of adiabatic drying under near-neutral atmospheric conditions (denoted by E_pmax in the literature). E_pads was calculated by the mass-transfer-based method, while E_pmax was calculated by the Penman-based equation (Ma & Zhang 2017). How to define the maximum ET_p could be the key to the application of the modified Brutsaert (2015) generalized CR model by Crago et al. (2016) and Szilagyi et al. (2017).

In recent years, the rapid development of satellite remote sensing has led to remote sensing data providing great convenience for spatio-temporal samplings (Singh et al. 2008; Liu et al. 2009). However, reanalysis-based methods (Mesinger et al. 2006; Yang et al. 2010; Lavers et al. 2012) have unique advantages, such as less sensitivity to cloud cover, longer temporal coverage, and better representation of actual conditions for spatially distributed data (Szilagyi et al. 2017). The China Meteorological Forcing Dataset is a reanalysis dataset of surface meteorological and environmental factors developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences. Using Princeton reanalysis forcing data, Global Land Data Assimilation System data, Global Energy and Water Exchanges-Surface Radiation Budget data, and Tropical Rainfall Measuring Mission data (3B42) as the background, this dataset also combines China Meteorological Administration station data from 1979 to 2015. The spatial resolution of this dataset is 0.1° and its temporal resolution is 3 h (He & Yang 2011).

The Hotan River Basin (HRB) is located in the hinterland of the Taklimakan Desert of northwest China (the largest desert in China) and at the north side of the Kunlun Mountains (part of the Tibetan Plateau). There are limited observation data available for extreme cold and arid areas of the river basin. In this study, the modified generalized CR framework of Crago et al. (2016), with the maximum ET_p defined by Crago et al. (2016) and Szilagyi et al. (2017), was applied to estimate the daily ET_a rates in the HRB from 2006 to 2013. The China Meteorological Forcing Dataset was used to calculate the CR components and analyze regional water balance, while Moderate Resolution Imaging Spectroradiometer (MODIS) data were used to calculate the net solar radiation (R_n). The daily precipitation data and monthly runoff data were collected to calculate a regional water budget. The objectives of this study were: (i) to analyze the improvement effects of different modified methods and select an appropriate definition of the maximum ET_p according to conditions in the study basin; (ii) to calibrate the local parameter values of the Priestley–Taylor coefficient and the adjustable parameter to estimate daily ET_a rates and validate the annual ET_a by using a regional water balance; and (iii) to evaluate the temporal and spatial distribution variations of ET_a and the water budget in the HRB during 2006–2014.

Utilizing the Surface Energy Balance System method of Su (2002), R_n was calculated using remote sensing data as an instantaneous value. The daily R_n values were estimated by implementing the sinusoidal model given by Bisht et al. (2005). The other climate-related parameters, G,⁠, ⁠, ⁠, ⁠, ⁠, and ⁠, were calculated using the method recommended by the Food and Agriculture Organization (Allen et al. 1998).

The HRB is located in northwest China (77.40°E–81.59°E, 34.84°N–40.44°N), covering an area of 88,753 km² (draining to the Xiaota hydrological station) (Figure 1). The HRB is one of the main headwaters of the Tarim River, which is the largest inland river in China. The altitude of the HRB varies between 1,014 m and 6,858 m above sea level (asl) (Figure 1). The source of the Hotan River is located at the north side of the Kunlun Mountains, and the lower reaches pass through the hinterland of the Taklimakan Desert. The mountainous regions (above 2,000 m asl) are covered with ice and snow, alpine grasslands, and bare rock (Figure 2), and are characterized by a frigid arid climate (Zheng et al. 2013). The plains regions (below 2,000 m asl) are mainly covered with croplands, desert steppes, and desert (Figure 2). The oasis-desert systems are characterized by a warm temperate arid climate (Zheng et al. 2013). Because of the colossal irrigation system in the HRB, drainage collects downstream of the irrigated area, where an oasis-desert transition area is formed. The Kalakashi River and the Yulongkashi River are two branches of the Hotan River that cross the Hotan oasis. There are only two meteorological stations in the Hotan oases and five hydrological stations on the Hotan River (Figure 1). All of these stations are at less than 2,000 m asl, and there are no stations in the mountainous regions or desert regions of the HRB. Two mountain-pass hydrological stations, Wuluwati and Tonguziluoke, are at around 2,000 m asl. The Aigeliya and Tuzhiluke hydrological stations are at the convergence of the two branches. The Xiaota hydrological station is at 319 km downstream of the convergence, which is the single basin export to the Tarim River (Zhao et al. 2009; Lyu et al. 2015).

In the present study, 0.1° grid-based meteorological data of the China Meteorological Forcing Dataset in the HRB during 2006–2014 were used, including mean daily air temperature, mean specific humidity, mean wind speed, precipitation, and 3-hourly mean air temperature. As mentioned in the Introduction, the China Meteorological Forcing Dataset was produced by merging a variety of data sources (He & Yang 2011). This study used 3-hourly mean temperature from this dataset to facilitate the preparation of daily maximum and minimum temperatures. Note that the wind speed in this dataset was at a height of 10 m (U₁₀); therefore, U₁₀ values were transformed to U₂ values via U₂=0.748 U₁₀ (Allen et al. 1998). We also used MODIS data at a spatial resolution of 1 km, including MOD11A2 (emissivity, land surface temperature), MCD43B3 (albedo), MCD15A2 (leaf area index), and MCD12Q1 (land cover type). Additionally, digital elevation model (DEM) data of approximately 90 m resolution were collected. The forcing datasets and DEM data were used with the cubic convolution resampling technique to attain the 1 km spatial resolution matching the remotely sensed data. MODIS data collected at 8-day time steps was linearly interpolated to obtain daily values. Spatio-temporal gap filling of missing data during cloudy periods was accomplished through compositing time series data (Cleugh et al. 2007). Moreover, the monthly runoff data from the five hydrological stations were collected to analyze regional water balance.

In the above description, was replaced by ⁠, and in Equation (3) was evaluated at T_w rather than at T_a. In this section, the improvement effects of different modified methods were analyzed. In addition, the two definitions of ET_pmax of Crago et al. (2016) and Szilagyi et al. (2017) described above, were compared and selected according to conditions in the study basin. Analysis of annual precipitation data from 2006 to 2014 revealed that the precipitation in 2011 was in the median value. Therefore, 2011 was chosen as an example to represent the climate of the study region.

As shown in Figure 3(a) and 3(b), the spatial distribution of the daily mean wind function derived from Equation (4) was compared with that derived from Equation (5) ⁠. The mean value of (3.11 mm day⁻¹ kPa⁻¹) was smaller than the mean value of (4.54 mm day⁻¹ kPa⁻¹). The spatial distribution of the absolute percentage difference (APD) between and (denoted by APD_f (U)), as well as the spatial distribution of the APD between ET_p derived from and that derived from (denoted by APD_ET_p) can be seen in Figure 3(c) and 3(d). The mean values of APD_f (U) and APD_ET_p were 0.32 and 0.23, respectively, with larger values in the plains regions than in the mountainous regions. Moreover, the values of APD_ET_p in desert regions were larger than those in oasis regions. These results indicate that the correction of the wind function had a more significant reduced-ET_p effect in the plains regions than in the mountainous regions, while the correction of the wind function to the estimate of ET_p had the most significant influence in desert regions of the plains regions, followed by influence in oasis regions of the plains regions; meanwhile, it had the least influence in the mountainous regions.

Figure 4(a) and 4(b) display the spatial distribution of the daily mean T_a and the daily mean T_w, respectively. The spatial distribution of the mean absolute error (MAE) between T_a and T_w (denoted by MAE (T_a−T_w)), and the MAE between ET_rad derived at T_a and that derived at T_w (denoted by MAE (ET_rad)) can be seen in Figure 4(c) and 4(d). The values of MAE (T_a−T_w) in the plains regions were larger than those in the mountainous regions, while those in oasis-desert transition regions were larger than those values in desert and/or oasis regions. The spatial distribution of MAE (ET_rad) was similar to the distribution of MAE (T_a−T_w). The correction of T_w to the estimate of ET_rad had a relatively small influence, with the maximum MAE (ET_rad) value of 32.87 mm (about 3.5% of the ET_rad derived at T_a).

Figure 5(a) and 5(b) show the spatial distribution of the annual ET_pds estimated from Equation (11) and the annual ET_pdry estimated from Equation (12), respectively. Figure 5(c) and 5(d) display the spatial distribution of the difference between ET_pds and ET_p (denoted by ET_pds−ET_p) and the difference between ET_pdry and ET_p (denoted by ET_pdry−ET_p), respectively. The values of ET_pds−ET_p in the plains regions were greater than 0, while those in the mountainous regions were less than 0. These results indicated that ET_pds was less than ET_p in the mountainous regions. According to the basic assumptions of Crago et al. (2016), ⁠, and the defined by Crago et al. (2016), i.e., ET_pds, was not applicable in the study basin, especially in the mountainous regions. However, the values of ET_pdry−ET_p were greater than 0 in the whole basin, with the minimum value of 170.2 mm yr⁻¹ in the mountainous regions and the maximum value of 660.2 mm yr⁻¹ near Xiaota station in the plains regions. Therefore, defined by Szilagyi et al. (2017), i.e., ET_pdry, should be adopted in the present study.

Priestley & Taylor (1972) and subsequent studies have suggested a value of the Priestley–Taylor coefficient in Equation (3) near 1.26 with a narrow range (e.g., Eichinger et al. 1996; Crago et al. 2010). Hobbins et al. (2001) obtained a value of = 1.3177 using data from 92 basins of the USA. Brutsaert (2005) suggested that for saturated surfaces typically falls between 1.20 and 1.30. Xu & Singh (2005) determined values for three study regions at 1.18, 1.04, and 1.00. Yang et al. (2008) obtained an average = 1.17 with a range of 0.87–1.48 from 108 basins of China, whereas Gao et al. (2012) suggested falls between 1 and 1.23 for nine sub-basins of the Haihe River Basin. Similar seasonal variability for in the range of 1.1–1.4 was reported for the Asian monsoon (Yang et al. 2013). Brutsaert (2015) and Ma et al. (2015b) used near 1.13. The value was also reported to be a function of the water content of drying soil surfaces, reaching values as low as 0.25 for relatively dry soil surfaces (Aminzadeh & Or 2014). To accommodate the observed range of ⁠, we assume ⁠. Szilagyi et al. (2017) proposed a novel approach to obtain the value of that applied to a calibration-free formulation of the CR. However, their approach to obtain the value of is unsuitable in extreme cold and arid basins. Due to the lack of measured data, a value of cannot be accurately determined. However, the fact that ET_w ≤ ET_p can estimate the possible maximum value of ⁠, thus, for the HRB. We assume that a fixed value of is applicable in the whole basin, even though we agree with the view that changes over time and underlying surface. Hence, we assume that for the HRB in the present study. In addition, we assume that ET_a = 0 when T_a is less than the dew point temperature.

The adjustable parameter in Equation (16) was treated as the only calibration parameter. Note that one can deduce from Equation (16) that ET_a rates increase monotonously as the value increases. Because measured data for the HRB in extreme cold and arid areas were limited, the calculated results were not verified at daily and monthly time scales. Therefore, the calibration of and the validation of ET_a used the regional water balance method at annual time scales.

Calibration was conducted as follows. First, the coefficient in desert or urban regions of the plains regions was determined since precipitation does not produce runoff and is almost completely evaporated (Ma et al. 2014). Second, the coefficient in vegetated areas of the plains regions was determined according to the regional water balance method. Third, the coefficient in rock bare regions of the mountainous regions was determined. Considering that the daily maximum temperature was less than 0 °C above 5,000 m asl during most of the year, ET_a rates above 5,000 m asl were quite low. As mentioned above, ET_a rates increased monotonously with increasing ⁠; therefore, we assumed that the values of decreased monotonously with increasing elevation. For conveniently calibrating, elevation of the mountainous regions was divided into three parts, and we assumed that the relationship between and elevation was linear in each elevation interval. Note that the relationship between and elevation might be actually nonlinear. Finally, the coefficient in vegetated areas of the mountainous regions was determined using a method similar to that employed for the rock bare regions. The optimal values of used for Equation (16) were found (Table 1) by minimizing the sum of absolute errors (MSAE) between ET_wb and ET_a. The optimal values of were applied to estimate ET_a in the validation period, and the results are listed in Table 2. Note that Zhao et al. (2009) obtained mean annual evapotranspiration (198.78 mm) from year 1954 to 2003 in middle regions. In this study, the mean annual evapotranspiration in middle regions was 209.2 mm which is very close to the results of Zhao et al. (2009).

In general, the mean annual ET_a estimated using the optimal values of matched well with the mean annual ET_wb, which indicates that the MSAE in the validation period was −12.8 mm yr⁻¹. However, this does not necessarily indicate that the annual MSAE was also small. Indeed, there was a certain deviation between the annual ET_a and annual ET_wb with the maximum MSAE value of 143.65 mm in 2010.

Figure 6 shows the spatial distribution of annual ET_a, and Figure 7 shows the spatial distribution of mean annual precipitation, ET_a, and P−ET_a difference during 2006–2014 in the HRB. The mean annual precipitation ranged from 46.0 mm to 342.0 mm, while the mean annual ET_a ranged from 2.3 mm to 800 mm. The higher precipitation values were in the mountainous regions, indicating that the mean annual precipitation was higher at high elevation. ET_a rates in the plains regions were higher than those in the mountainous regions, with the highest values being observed in croplands and shrublands along the riverside. Perhaps the irrigation or groundwater plays a key role in this issue. Relatively lower ET_a occurred in ice and snow regions, as well as in barren land of the mountainous regions. The difference (P−ET_a) reflects water surplus or deficits (Figure 7(c)). The mean annual P−ET_a difference in the mountainous regions was positive (except in alpine grasslands), while it was negative in the plains regions. These findings indicate that a great deal of precipitation and runoff dissipated in the form of ET_a in the plains regions.

Figure 8 shows the monthly mean P, ET_a, and P−ET_a difference from 2006 to 2014 in the upper, middle, and lower regions. As shown in Figure 8(a), most of the P was concentrated in the months of May to September, accounting for 65% (128 mm) of annual P in the upper regions, 61% (66 mm) in the middle regions, and 70% (61 mm) in the lower regions. As shown in Figure 8(b), the variations in ET_a in different regions were similar, following a single peak normal distribution. Most of ET_a was concentrated in the months of May to September, accounting for 89% of the mean annual ET_a (42 mm) in the upper regions, 76% (159 mm) in the middle regions, and 80% (105 mm) in the lower regions. The inter-annual changes in ET_a in winter months were lower than in summer months. As shown in Figure 8(c), a water deficit occurred in the middle and lower regions, and a water surplus occurred in the upper regions. Because the inter-annual changes in P were relatively large, the inter-annual changes in water budget were relatively large as well. Water deficits occurred during most months throughout the year in the middle regions (except January and February). In the lower regions water deficits occurred during May to September. From October to April, water budget was a basic equilibrium, because the river was drying up during these months. Additionally, the annual deficit was estimated to be 100 mm in the middle regions and 44 mm in the lower regions. Furthermore, the amount of available water was sufficient to satisfy ET_a in most months in the upper regions, where the mean annual water surplus was 150 mm.

There is great debate regarding whether proportionality parameter b is constant. Huntington et al. (2011) applied a modified AA model with b = 1 to estimate monthly and annual ET_a from arid shrublands of the southwestern USA and obtained quite satisfactory prediction accuracy. Ma et al. (2015b) applied a modified CR-based AA model to obtain daily ET_a in the alpine steppe of the Tibetan Plateau and believed that a symmetric CR (b = 1) contradicted previous research that used default parameter values to claim an asymmetric CR in arid and semiarid regions of the Tibetan Plateau. Sugita et al. (2001) demonstrated that b equals unity only when the underlying surface is smooth enough and soil moisture is sufficient. Furthermore, Ma et al. (2015b) suggested that a strictly symmetric CR is difficult to achieve because the actual wet surface is too small or too large, or there is the influence of some additional heat transfer though the bottom or side of pans. Kahler & Brutsaert (2006) also suggested that the asymmetry with ET_p given by pan evaporation was due to the effect of a significant bottom or sidewall effect of pans. However, Szilagyi (2007) and Brutsaert (2015) demonstrated that the asymmetry is the inherent nature rather than the cases. In the present study, Equation (14), proposed by Crago et al. (2016), was used to estimate the ET_a. The present CR model is a self-adjusting CR since b could be obtained for each daily and each grid value as b= (1 −x_min)/x_min (Crago et al. 2016). The CR is only linear (symmetric) when x_min= 0.5. The definition of variable X defined by Crago et al. (2016) is essentially the same as the definition of Szilagyi et al. (2017). However, the relationship between y and x or X is actually still unspecified. Brutsaert (2015) and Crago et al. (2016) recommended the simple polynomial to add some flexibility in CR formulations, so that its parameters can be calibrated easily. Hence, Equation (16) was adopted in this study.

In the current study, the HRB was divided into two terrain conditions: the plains regions (below 2,000 m asl) and the mountainous regions (above 2,000 m asl), while the HRB was divided into two underlying surface conditions: vegetated regions and unvegetated regions. Hence, to precisely describe the spatial distribution of ET_a, the of each condition was recalculated according to the regional water balance and the basic evaporation laws, as mentioned above (Table 1). The underlying surface of the mountainous regions changed with elevation; therefore, a linear function of elevation was used to illustrate the change of in the mountainous regions. Note that calibrate varies significantly. This may show that this CR method might not be used in extreme cold and arid areas.

Obviously, estimated ET_a has some uncertainties. For example, extreme cold and arid areas often have precipitation saved in the form of snow and ice until they melt in subsequent years. This means that has a greater effect at smaller time scales, but this should not be assumed to be 0. However, in the present study, it had to be assumed to be 0, as explained above. Moreover, over the past few years, the data obtained from the Gravity Recovery and Climate Experiment (GRACE) have notably improved our understanding of water storage change at monthly time scales in large basins globally (Scanlon et al. 2015; Save et al. 2016). However, the area of the HRB makes it unsuitable for analysis because it uses GRACE data. Secondly, there are uncertainties regarding the function relations of and elevation in the mountainous regions. A linear function was selected since it was easy to understand and be calibrated; however, this does not mean that linear function is the best relationship between and elevation. In addition, the values of changed with elevation rather than time. The values of were also reported to be a function of surface meteorological and environmental factors (Aminzadeh & Or 2014). Obviously, the setting of values of can make a significant difference in the calibration parameter ⁠. As noted by Szilagyi et al. (2017), any CR-based method is quite sensitive to the parameter value of ⁠. Finally, the reanalysis dataset and runoff data still had significant uncertainties in extreme cold and arid areas because of the extremely sparse meteorological stations and harsh natural conditions. Although ET_a estimated in the mountainous regions remains uncertain and could not be verified at daily and monthly time scales, this study not only provides a new concept for calibration, but also a potential solution for different underlying surfaces and time scales.

In this study, we used the generalized CR to estimate ET_a rates during 2006–2014 at daily time scales in extreme cold and arid areas. The HRB served as a prime example. The China Meteorological Forcing Dataset and MODIS data were used to calculate CR components and to analyze regional water balance. Comparison and analysis revealed that correction of the wind function had more significant effects than correction of T_w. ⁠, defined by Crago et al. (2016), was not applicable in the study basin, especially in the mountainous regions. Hence, we adopted ⁠, defined by Szilagyi et al. (2017). Since available data in the HRB were limited, the validity and accuracy of the calculated ET_a were verified using a regional water balance method at annual time scales. The coefficient was calibrated based on the elevation and underlying surface types. Mean annual ET_a estimated using optimal values of matched the mean annual ET_wb well, while there was deviation between annual ET_a and annual ET_wb. Based on the model estimations, spatial and temporal distributions of ET_a and the water budget in the basin from 2006 to 2014 could offer a reference at monthly or even daily time scales. The mean annual ET_a ranged from 2.3 mm to 800 mm during 2006–2014. ET_a rates in the plains regions were higher than those in the mountainous regions, with the highest values being observed in croplands and shrublands along the riverside. Most of ET_a was concentrated in the months of May to September, accounting for 89% of the mean annual ET_a (42 mm) in the upper regions, 76% (159 mm) in the middle regions, and 80% (105 mm) in the lower regions. A water deficit occurred in the middle and lower regions, and a water surplus occurred in the upper regions. The annual deficit was estimated to be 100 mm in the middle regions and 44 mm in the lower regions. The amount of available water was sufficient to satisfy ET_a in most months in the upper regions, where the mean annual water surplus was 150 mm.

If meteorological stations and in situ flux stations could be established in the mountainous and desert regions in the future, the accuracy of estimated ET_a and water management efficiency could be improved. Moreover, a model-independent calculation of considering underlying surface conditions and environmental conditions has important implications for a calibration-free CR model. That CR model could be easily applied to all kinds of conditions, including extreme cold and arid areas. However, additional studies should be conducted to further improve the model.

The forcing datasets used in this study were developed by the Data Assimilation and Modeling Center for Tibetan Multi-spheres, Institute of Tibetan Plateau Research, and Chinese Academy of Sciences (http://westdc.westgis.ac.cn/data/). MODIS data and DEM data were provided by the Geospatial Data Cloud site, Computer Network Information Center, and Chinese Academy of Sciences (http://www.gscloud.cn). This work was supported by the National Natural Science Foundation of China (51509202 and 51679185). We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript. We are grateful to the four anonymous reviewers and the editor whose valuable comments and suggestions led to important revisions.