Long-term variations in actual evapotranspiration over the Tibetan Plateau

Actual terrestrial evapotranspiration (ETa) is a key parameter controlling land–atmosphere interaction processes and water cycle. However, spatial distribution and temporal changes in ETa over the Tibetan Plateau (TP) remain very uncertain. Here we estimate the multiyear (2001–2018) monthly ETa and its spatial distribution on the TP by a combination of meteorological data and satellite products. Validation against data from six eddy-covariance monitoring sites yielded root-mean-square errors ranging from 9.3 to 14.5 mm per month and correlation coefficients exceeding 0.9. The domain mean of annual ETa on the TP decreased slightly (−1.45 mm yr−1, p < 0.05) from 2001 to 2018. The annual ETa increased significantly at a rate of 2.62 mm yr−1 (p < 0.05) in the eastern sector of the TP (long> 90 E) but decreased significantly at a rate of −5.52 mm yr−1 (p < 0.05) in the western sector of the TP (long< 90 E). In addition, the decreases in annual ETa were pronounced in the spring and summer seasons, while almost no trends were detected in the autumn and winter seasons. The mean annual ETa during 2001–2018 and over the whole TP was 496± 23 mm. Thus, the total evapotranspiration from the terrestrial surface of the TP was 1238.3±57.6 km3 yr−1. The estimated ETa product presented in this study is useful for an improved understanding of changes in energy and water cycle on the TP. The dataset is freely available at the Science Data Bank (https://doi.org/10.11922/sciencedb.t00000.00010; Han et al., 2020b) and at the National Tibetan Plateau Data Center (https://doi.org/10.11888/Hydro.tpdc.270995, Han et al., 2020a). Published by Copernicus Publications. 3514 C. Han et al.: Monthly evapotranspiration on the Tibetan Plateau from 2001 to 2018


Introduction
As the birthplace of Asia's major rivers, the Tibetan Plateau (TP), famous as the "Water Tower of Asia", is essential to the Asian energy and water cycles (Immerzeel et al., 2010;Yao et al., 2012). Along with increasing air temperature, evidence from the changes in precipitation, runoff, and soil moisture indicates that the hydrological cycle of the TP has been intensified during the past century . Consuming around two-thirds of global terrestrial precipitation, evapotranspiration (ET) is a crucial component that affects the exchange of water and energy between the land surface and the atmosphere (Oki and Kanae, 2006;Fisher et al., 2017). ET is also a key factor modulating regional and global weather and climate. As one essential connecting component between the energy budget and the water cycle in the terrestrial ecosystems (Xu and Singh, 2005), ET and its variations have been drawing more attention worldwide (Xu and Singh, 2005;Li et al., 2014;T. Zhang et al., 2018;Yao et al., 2019;G. Wang et al., 2020). Total evaporation from large lakes of the TP has been quantitatively estimated recently (B. ; however, the terrestrial ET on the TP and its spatial and temporal changes remain very uncertain. Many studies have tried to evaluate ET's temporal and spatial variability across the TP using various methods. The pan evaporation (E pan ), which represents the amount of water evaporated from an open circular pan, is the most popular observational data source of ET. Long time series of E pan are often available with good comparability among various regional measurements. Thus, it has been widely used in various disciplines, e.g., meteorology, hydrology, and ecology. Several studies have revealed the trend of E pan on the TP (Zhang et al., 2007;Liu et al., 2011;Shi et al., 2017;C. Zhang et al., 2018;Yao et al., 2019). Although E pan and potential ET suggest the long-term variability in ET according to the complementary relationship (CR) between E pan and actual ET (ET a ) (Zhang et al., 2007), these measures cannot precisely depict the spatial pattern of trends in ET a . Recently, several studies applied revised models based on the CR of ET to estimate ET a on the TP (T. Ma et al., 2019;G. Wang et al., 2020). Employing only routine meteorological observations without requiring any vegetation and soil information is the most significant advantage of CR models (Szilagyi et al., 2017). However, numerous assumptions and requirements of validations of key parameters limit the application and performance of CR models over different climate conditions. The application of eddycovariance (EC) technologies in the past decade has dramatically advanced our understanding of the terrestrial energy balance and ET a over various ecosystems across the TP. However, the fetch of the EC observation is on the order of hundreds of meters, thus impeding the ability to capture the plateau-scale variations in ET a . Therefore, finding an effective way to advance the estimation of ET a on the TP is of great importance.
Satellite remote sensing (RS) provides temporally frequent and spatially contiguous measurements of land surface characteristics that affect ET, for example, land surface temperature, albedo, and vegetation index. Satellite RS also offers the opportunity to retrieve ET over a heterogeneous surface (Zhang et al., 2010). Multiple RS-based algorithms have been proposed. Among these algorithms, the surface energy balance system (SEBS) proposed by Su (2002) has been widely applied to retrieve land surface turbulent fluxes on the TP (Chen et al., 2013b;Ma et al., 2014;Zou et al., 2018;Zhong et al., 2019;Han et al., , 2017. Chen et al. (2013b) improved the roughness length parameterization scheme for heat transfer in SEBS to expand its modeling applicability over bare ground, sparse canopy, dense canopy, and snow surfaces in the TP. An algorithm for effective aerodynamic roughness length had been introduced into the SEBS model to parameterize subgrid-scale topographical form drag (Han et al., 2015(Han et al., , 2017. This scheme improved the skill of the SEBS model in estimating the surface energy budget over mountainous regions of the TP. A recent advance by Chen et al. (2019) optimized five critical parameters in SEBS using observations collected from 27 sites globally, including 6 sites on the TP, and suggested that the overestimation of the global ET was substantially improved with the use of optimal parameters.
While the spatial and temporal pattern of the ET a in the TP had been investigated in many studies (T. G. Wang et al., 2020;Zhang et al., 2007), considerable inconsistencies for both trends and magnitudes of ET a exist due to uncertainties in forcing and parameters used by various models. Thus, in this study, with full consideration of the recent developments in the SEBS model over the TP, we aim to (1) develop an 18-year (2001-2018) ET a product of the TP, along with independent validations against EC observations; (2) quantify the spatiotemporal variability in the ET a in the TP; and (3) uncover the main factors dominating the changes in ET a , using the estimated product.

Model description
The SEBS model (Su, 2002) was used to derive land surface energy flux components in the present study. The remotely sensed land surface energy balance equation is given by (1) R n is net radiation flux (W m −2 ), H is sensible heat flux (W m −2 ), LE is latent heat flux (W m −2 ), and G 0 is ground heat flux (W m −2 ). Note that this equation neglected energy stored in the canopy, energy consumption related to freezethaw processes of permafrost and glaciers, etc. This equation is not applicable to any condition where a phase change in water occurs, except the liquid-to-vapor phase change.
The land surface net radiation flux was computed as where α is the land surface albedo derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) products. Downward shortwave (SWD) and longwave (LWD) radiation was obtained from the China Meteorological Forcing Dataset (CMFD). Land surface temperature (T s ) and emissivity (ε) values were also obtained from MODIS products. In vegetated areas the soil heat flux, G 0 , was calculated from the net radiation flux and vegetation cover: where r s and r c are ratios of ground heat flux and net radiation for surfaces with bare soil and full vegetation, respectively. Fractional vegetation cover (f c ) was derived from the normalized difference vegetation index (NDVI). Over water surfaces (NDVI < 0 and α < 0.47), G 0 = 0.5R n was used (Chen et al., 2013a;Gao et al., 2011). On glaciers, G 0 is negligible (Yang et al., 2011) and G 0 = 0.05R n .
In the atmospheric surface layer, sensible heat flux and friction velocity were calculated based on the Monin-Obukhov similarity (Stull, 1988), U is the horizontal wind velocity at a reference height z (m) above the ground surface; θ 0 is the potential temperature at the land surface (K); θ a is the potential temperature (K) at the reference height z; d 0 is the zero-plane displacement height (m); ρ is the air density (kg m −3 ); C p is the specific heat for moist air (J kg −1 • C −1 ); κ = 0.4 is the von Kármán's constant; u * is the friction velocity; L is the Monin-Obukhov length (m); θ v is the potential virtual temperature (K) at the reference height z; ψ m and ψ h are the stability correction functions for momentum and sensible heat transfer, respectively; and g is the gravity acceleration (m s −2 ). To account for the form drag caused by subgrid-scale topographical obstacles, effective roughness lengths for momentum (z eff 0 m , m) and sensible heat (z eff 0 h , m) transfer were introduced into the SEBS model by Han et al. (2017). These modifications are parameterized as follows (Grant and Mason, 1990;Han et al., 2015): ln h 2z eff where h is the average height of the subgrid-scale roughness obstacles; λ is the average density of the subgrid-scale roughness elements calculated from digital elevation models; D is the form drag coefficient, and D = 0.4 is used for the mountainous areas of the TP as suggested by Han et al. (2015); and z 0 m and z 0 h are the local-scale roughness lengths for momentum and heat transfer (m), respectively. Detailed calculations can be found in Su (2002). A revised algorithm for z 0 h developed by Chen et al. (2013b) was applied as this algorithm outperforms the original scheme of the SEBS model on the TP. To constraint the actual evapotranspiration, the evaporative fraction was applied in the SEBS model, which is determined by taking energy balance considerations at dry and wet limiting cases. Under the dry-limit condition, the evaporation becomes zero due to the limited supply of available soil moisture, while water vapor evaporates at the potential rate under the wet-limit condition (Su, 2002). The evaporative fraction ( ) is defined as After calculating evaporative fraction based on the assumption of dry and wet limits, latent heat flux was calculated by inverting Eq. (9). Finally, latent heat flux was converted to ET a . Details are available in Su (2002) and Chen et al. (2013a). Note that the dry-wet limit assumption did not apply to frozen soil, water, snow, and ice surfaces. The latent heat flux was obtained as the residual of the surface energy balance Eq. (1) after calculating net radiation, sensible heat flux, and ground heat flux when the dry-wet limit assumption is not applicable.

Data
In situ observations, satellite-based products, and meteorological forcing data were used in this study to estimate monthly ET a over the TP area. The CMFD, which was developed based on the released China Meteorological Administration (CMA) data (He et al., 2020), was used as model input. The CMFD covers the whole landmass of China at a spatial resolution of 0.1 • and a temporal resolution of 3 h. The CMFD dataset was established through the fusion of in situ observations, remote sensing products, and reanalysis datasets. In particular, the dataset benefits from the merging of the observations at about 700 CMA weather stations and by using the Global Energy and Water Cycle Experiment -Surface Radiation Budget (GEWEX-SRB) shortwave radiation dataset (Pinker and Laszlo, 1992). The GEWEX-SRB data have not been used in any other reanalysis dataset. In addition, independent datasets observed in western China, where weather stations are scarce, were used to evaluate the CMFD. This includes data collected through the Heihe Watershed Allied Telemetry Experimental Research (Hi-WATER) (Li et al., 2013) and the Coordinated Enhanced Observing Period (CEOP) Asia-Australia Monsoon Project (CAMP) (Ma et al., 2003). The CMFD dataset has been validated against in situ meteorological observations and compared with other reanalysis datasets on the TP, demonstrating that it is one of the best meteorological forcing datasets over the TP area (Zhou et al., 2016;Xie et al., 2017;. Therefore, it is suitable for this study to drive the SEBS model. Detailed information for the CMFD dataset is listed in Table 1. MODIS monthly land surface products, including land surface temperature and emissivity, land surface albedo, and vegetation index, provide land surface conditions for the SEBS model. Detailed information on MODIS land surface variables is listed in Table 1. The values of land surface variables in the MODIS monthly products are derived by compositing and averaging the values from the corresponding month of MODIS daily files. Validations of MODIS land surface temperature and albedo against in situ observations on the TP suggest a high quality of MODIS land surface products with low biases and small root-mean-square errors Wang et al., 2004;Ma et al., 2011).
In situ EC data observed at six flux stations on the TP were used to validate model results. Locations of the six observation sites are illustrated in Fig. 1, and detailed descriptions for these six sites are shown in Table 2. The instrumental setup at each site consists of an EC system comprising a sonic anemometer (CSAT3,Campbell Scientific Inc) and an open-path gas analyzer (LI-7500, LI-COR); a fourcomponent radiation flux system (CNR-1, Kipp & Zonen), installed at a height of 1.5 m; a soil heat flux plate (Hukseflux, HFP01), buried in the soil at a depth of 0.1 m; and soil moisture and temperature probes, buried at a depth of 0.05, 0.10, and 0.15 m, respectively (Han et al., 2017). The EC data were processed with the EC software package TK3 (Mauder and Foken, 2015). The main post-processing procedures of the EC raw data were as follows: spike detection, coordinate rotation, spectral loss correction, frequency response corrections (Moore, 1986), and corrections for density fluctuations (Webb et al., 1980). The ground heat flux was obtained by summing the flux value observed by the heat flux plate and the energy storage in the layer above the heat flux plate . A more comprehensive dataset including the EC data used in this work has been published and is freely available .
The 3-hourly CMFD data were averaged into daily and then into monthly data to be consistent with MODIS products in terms of temporal resolution. Daily land surface albedo has been averaged into a monthly variable. MODIS land surface products and canopy height data were remapped onto the CMFD's grid. Monthly EC data and in situ meteorological observations, which are used for model validation, were generated from half-hourly variables.

Model evaluation metrics and data analysis methods
The model performance was assessed using the Pearson correlation coefficient (R), the root-mean-square error (RMSE), and the mean bias (MB) between the estimated and observed monthly ET a at the six stations on the TP.
The least-square regression technique was used to detect the long-term linear annual trends in ET a values. The linear model to simulate ET a values (Y t ) against time (t) is defined as below, and the slope of the linear equation (b) is taken as the changing trend, The Student's t test, having an n − 2 degree of freedom (n is the number of samples), was used to evaluated the statistical significance of the linear trends, and only tests with a p value less than 0.05 were selected as having passed the significance test.

Validation against flux tower observations
The SEBS-estimated ET a was validated against EC observations at six flux stations on the TP at a monthly scale (Fig. 2). The SEBS model is capable of capturing both the magnitude and seasonal variation in the monthly ET a signal at all the six stations. The correlation coefficients are all larger than 0.9 and have passed the significance test at the p = 0.01 level. Specifically, the SEBS model performed particularly well at the short-grass sites (BJ and NAMORS), with correlation coefficients as high as 0.98 and MB values below 5.0 mm per month. At the high-grass site (SETORS) and the gravel site (QOMS), the SEBS model is capable of reproducing the EC-observed monthly ET a with RMSE values of 14.5 and 13.2 mm per month, respectively. In addition, the underestimates of ET a by SEBS are mostly in the dry season, when the canopy is withered. The validation at the site scale indicates that the SEBS model used in this work can be applied to a wide range of ecosystems over the TP.

Spatial distribution
There was a clear spatial pattern to the multiyear average of annual ET a between 2001 and 2018 (Fig. 3). In general, the SEBS-estimated ET a decreases from the southeast to the northwest of the TP, with the maximum value above 1200 mm in southeastern Tibet and the Hengduan Mountains and the minimum value less than 100 mm on the northwestern edge of the TP. In the central TP, where there are several lakes, ET a was typically from 500 to 1000 mm. ET a was lower than 200 mm over the high snow-and ice-bound mountainous areas, for example, over the northern slopes of the Himalaya, the Nyenchen Tanglha Mountains, and the eastern section of the Tanggula Mountains. The reason is that these snow-and ice-bound mountainous areas have a higher ability to reflect downward shortwave radiation and hence have less available energy to evaporate. On the whole, the domainaveraged multiyear mean annual ET a over the TP is 496 ± 23 mm. The total amount of water evapotranspired from the terrestrial surface of the TP is around 1238.3±57.6 km 3 yr −1 , considering the area of the TP to be 2.5 × 10 6 km 2 . Figure 4 shows the multiyear-average spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February) ET a on the TP. Generally, the distribution pattern of seasonal ET a was comparable with that of the annual ET a . Both seasonal and annual ET a show a decreasing trend from the southeastern TP to the northwestern TP. Note that the spatial contrast of ET a almost faded out in the winter season owing to a minimum in available energy and precipitation (Fig. 4d). The ET a in spring is higher than that in autumn, except for some high mountainous areas (e.g., mountain ranges of the Himalaya and Hengduan Mountains). The spring ET a ranges from 50 to 450 mm, while autumn ET a ranges from 50 to 250 mm. In summer, the ET a is larger than 250 mm in most of the TP, while the ET a is still below 100 mm in large areas of the northwestern TP. The multiyear seasonal ET a averaged over the whole TP is 140 ± 10, 256±12, 84±5, and 34±4 mm for spring, summer, autumn, and winter, respectively.

Trend analysis
The trend of annual ET a during 2001-2018 is shown in Fig. 5. Overall, an increasing trend of SEBS-simulated ET a is dominant in the eastern TP (long > 90 • E), while a decreasing trend is dominant in the western TP (long < 90 • E). The trends pass the t test (p < 0.05) in most parts of the areas. The decreasing trend in the western TP is pronounced and passes the t test (p < 0.05). This trend is larger than −7.5 mm yr −1 in most parts of the area and even larger than −10 mm yr −1 in a few parts. In the eastern TP, the increasing trend is mostly between 5 and 10 mm yr −1 and passes the t test (p < 0.05). The ET a trend tends to be greater along the marginal region of the northern, eastern, and southeastern TP. Along the marginal region of the southwestern TP and in the western section of the Himalaya Mountains this trend weakens.  The trends of seasonal ET a between 2001 and 2018 are spatially heterogeneous over the TP (Fig. 6). Decreasing trends in spring and summer are generally at a rate between −2.5 and −7.5 mm yr −1 , and increasing trends are generally at a rate below 5.0 and 7.5 mm yr −1 in spring and summer, respectively. Areas showing decreasing ET a tend to become larger in autumn and winter seasons. Both the decreasing and increasing trends are subdued in autumn and winter compared with that in the spring and summer seasons. Decreasing rates of ET a in autumn and winter are generally below −2.5 mm yr −1 , and only a few areas have a rate larger than −2.5 mm yr −1 .
Due to the contrast in the trends in the eastern and western halves of the TP, we divided the TP into two regions: the eastern TP (long > 90 • E) and the western TP (long < 90 • E). Trends of the ET a anomaly averaged over the entire TP, the western TP, and the eastern TP are shown in Fig. 7a. The domain means of ET a on the TP as a whole and in the western TP decreased at rates of −1.45 and −5.52 mm yr −1 , respectively. However, the ET a in the eastern TP increased at a rate of 2.62 mm yr −1 . The decreasing rate of ET a in the entire TP is influenced mainly by the significant decrease in ET a in the western TP. Seasonally, the rates of change in ET a over the whole TP are −0.82 mm yr −1 (p < 0.05) and −0.79 mm yr −1 (p < 0.05) in spring and summer, respectively (Fig. 7b). However, in autumn and winter the ET a changes at a rate of 0.10 and 0.06 mm yr −1 , respectively, and does not pass the t test (p < 0.05). ET a in the spring and summer seasons accounts for 75.7 % of the annual ET a . The variation in amplitude and changing rates in these two seasons are much larger than in the other seasons. Moreover, spatial distributions of spring and summer ET a trends are close to that of the annual ET a trend (Fig. 6). Thus, changes in ET a in the spring and summer dominate the variations in ET a in the whole year.
The decrease in ET a over the whole TP and in the western TP during 2001-2018 can be explained by the decrease in R n in the same time period (Fig. 8a). From 2001 to 2012,   ET a averaged over the entire TP increased slightly and then decreased dramatically from 2012, reaching a minimum in 2014. The significant decrease in ET a between 2012 and 2014 was due to the rapid decline in the R n (Fig. 8a). In the eastern TP, ET a increased during 2001-2018, while R n decreased in the same period. Thus, R n was not the dominant factor controlling the annual variations in ET a . However, the increasing trends of both precipitation and air temperature can explain the increase in ET a in the eastern TP during the period 2001-2018 ( Fig. 8b and c). The increasing precipitation increased the water resource available for ET a . Moreover, the increasing air temperature accelerated the melting of permafrost and glaciers on the TP. Hence, the melting water replenished the ecosystem and increased the ET a of the eastern TP.
Although the domain-averaged trend in ET a has been decreasing across the entire TP from 2001 to 2018, ET a values in some areas have increased. Moreover, the changing rates also depend on the time series of ET a . For example, the ET a increased slightly from 2001 to 2012, while it decreased from 2001 to 2018. This demonstrates the necessity to evaluate the spatial distribution of changing trends in ET a and utilize long time series to investigate the trends in ET a over the TP.

Data availability
The dataset presented and analyzed in this article has been released and is available for free download from the Science Data Bank (https://doi.org/10.11922/sciencedb.t00000.00010, Han et al., 2020b) and from the National Tibetan Plateau Data Center (https://doi.org/10.11888/Hydro.tpdc.270995, Han et al., 2020a). The dataset is published under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.

Summary and conclusions
The SEBS-estimated ET a is at a resolution of around 10 km, while the footprint of EC-observed ET a values ranges from a few dozen meters to a few hundred meters. SEBS-estimated ET a compares very well with observations at the six flux towers, showing low RMSE and MB values. These estimates were able to capture annual and seasonal variations in ET a , despite these two datasets being mismatched in their spatial representation. Note that the energy consumption related to freeze-thaw processes and sublimation is neglected. Thus, the dataset is likely to be less reliable over the glacier and permafrost and in the winter season.
Heterogeneous land surface characteristics and nonlinear changes in atmospheric conditions resulted in hetero-geneities in spatial distributions of ET a and changes in ET a . The SEBS-estimated multiyear (2001-2018) mean annual ET a on the TP was 515 ± 22 mm, resulting in approximately 1287.5 ± 55.0 km 3 yr −1 of total water evapotranspiration from the terrestrial surface. Annual ET a generally decreased from the southeast to the northwest of the TP. The maximum was over 1200 mm, in southeastern Tibet and the Hengduan Mountains, while the minimum was less than 100 mm in the northwestern marginal area of the TP. Moreover, ET a was typically lower than 200 mm over snow-and ice-bound mountainous areas as there was limited available energy to evaporate the water.
Averaged over the entire TP, annual ET a increased slightly from 2001 to 2012 but decreased significantly after 2012 and reached a minimum in 2014. Generally, there was a slight decreasing trend in the domain mean annual ET a on the TP at the rate of −1.45 mm yr −1 (p < 0.05) from 2001 to 2018. However, trends of annual ET a were opposite in the western and eastern TP. The annual ET a decreased significantly in the western TP at a rate of −5.52 mm yr −1 (p < 0.05) from 2001 to 2018, while annual ET a in the eastern TP increased at a rate of 2.62 mm yr −1 (p < 0.05) in the same period.
The spatial distributions of seasonal ET a trends were also noticeably heterogeneous during 2001-2018. The spatial patterns of ET a trend in spring and summer were similar to the annual changes in ET a . ET a decreased as well in the spring and summer seasons but at slower rates compared with the annual ET a ; however, only very weak trends were found in the autumn and winter seasons.
Author contributions. CH and YM conceived the study and analyzed simulation output. CH ran simulations, post-processed output, and generated the figures. CH and YM wrote the manuscript with assistance from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement.
This article is part of the special issue "Extreme environment datasets for the three poles". It is not associated with a conference.