Spatio-temporal variations in terrestrial water storage and its controlling factors in the Eastern Qinghai-Tibet Plateau

The Qinghai-Tibet Plateau, the Earth's third pole, possesses water resources that supply billions of people downstream in China, India, Pakistan and other countries in the region (Shrestha et al. 2015; Pritchard 2017; Zhao et al. 2019). Temporal and spatial changes in its water storage capacity have a major impact on the social economic development and livelihood of the region. In the context of global climate change, the ablation of ice and snow in this region is accelerating (Brun et al. 2017), and the land water cycle has undergone significant changes, particularly in the Eastern Qinghai-Tibet Plateau (EQTP). In recent years, the EQTP has attracted increasing attention from various sectors due to the increased frequency and intensity of natural disasters, ecological degradation, and the decreasing social and economic development. Frequent natural disasters are mainly represented by large-scale and long-term droughts and regional floods, which are generally linked to abnormal terrestrial water dynamics (Thomas et al. 2014; Xu et al. 2019). However, the lack of in-situ observations limits our understanding of the water cycle in this region.

Fortunately, studies have reported that the terrestrial water storage (TWS) derived from Gravity Recovery and Climate Experiment (GRACE) satellite data shows both the potential for revealing changes in surface water storage and the ability to support large-scale hydrologic modeling (Klees et al. 2007; Luo et al. 2016; Yang et al. 2017; Zhao & Li 2017; Frappart & Ramillien 2018). As mentioned before, GRACE has been used in many river basins, and has displayed satisfactory results for TWS changes (Li et al. 2013, 2018b; Ni et al. 2014; Feng et al. 2017; Khaki et al. 2018; Wang et al. 2019). A promising application can also be found for demonstration of TWS variability in the Three Gorges Reservoir area, the Guanzhong Plain and some medium-scale basins.

There are few studies on characterizing TWS changes on the EQTP based on GRACE observations. The previous studies have mostly addressed the entire plateau and partly discussed the EQTP, for example, Deng et al. (2018) and Meng et al. (2019) inferred a declining trend from 2002 to 2016 in the EQTP, but neglected the spatial heterogeneity of TWS change. Moreover, existing studies have been mostly concentrated on the annual scales of TWS change; few studies have pointed out the comprehensive mechanisms of TWS change, especially in the region with large spatial relief and limited in-situ measurements like the EQTP.

Generally, the spatio-temporal characteristics of TWS and its change are the result of the combined effect of climate change and human activities (Hu et al. 2018). Precipitation, runoff, and evapotranspiration are the dominating processes of changes in TWS either in space or time (Tangdamrongsub et al. 2011; Frappart et al. 2013; Yang et al. 2015). As an active process in rivers, the depletion of glaciers and snow cover is an influential factor in TWS dynamics (Syed et al. 2008). Anthropogenic processes like irrigation, the abstraction of groundwater, land use, etc. also impact TWS in highly populated areas (de Beurs et al. 2015; Huang et al. 2015; Khandu et al. 2016).

Several attempts have been carried out to understand the possible causes of TWS changes on the EQTP and its surrounding region. According to a recent study, the dominant factor contributing to an increase in TWS in eastern India and southeastern Tibet was increased precipitation, while the decrease was mainly attributed to a decrease in precipitation and an increase in irrigation water (Rodell et al. 2018). Another study, however, pointed out that evapotranspiration dominated the TWS depletion in the Brahmaputra river and precipitation played an important role on the TWS accumulation in the upper Yangtze. Additionally, glacier mass loss was the most likely cause of TWS decrease in the Brahmaputra river (Meng et al. 2019). Two findings showed somewhat different dominating factors in TWS change in almost the same area (most of the EQTP), which leads to uncertainty of our understanding on the cause of TWS changes. Each region has its own particular local climates, relief conditions, and human conditions, and some of these factors affect the evaluation of TWS in large-scale regions (Deng et al. 2018); for instance, in the EQTP, human activities are less intense, and their impact on TWS can be ignored, but the widely distributed glaciers and snow in this area require us to consider the substantial impacts of snow and ice meltwater on TWS. Consequently, there is still a lack of consistent understanding on the causes of TWS change. Therefore, the controlling factors for variations in TWS should be further discussed with additional validated datasets and detailed geographic division in the EQTP.

Additionally, the EQTP contains the source regions of the Yangtze, Lancang/Mekong, and Nujiang/Salween rivers. Their hydrologic dynamics are key to water resources in the downstream area. Huge fluctuations in downstream discharge have been focused on in recent years, especially in international rivers. Eliminating spatio-temporal patterns and the effect of climate features on the change of TWS may largely explain the water cycle process in source areas and the non-human impact on flow changes. Therefore, the present study is of great significance for local agriculture and animal husbandry, and ecological water demand predictions and management.

The overarching goal of this study is to characterize the spatial heterogeneity of the TWS change on the EQTP and explore a comprehensive understanding of possible causes of TWS change. We first obtain the GRACE-derived TWS in the EQTP and describe its variations by decomposing the total signal into inter-annual and seasonal variability. Second, we divide the EQTP into several parts according to its geographical conditions and present the dominant factors associated with the increase and decrease in TWS throughout the region and in its subbasins by using statistical methods and multi-source data. The main findings are presented in the Results and discussion section below.

We defined the EQTP as the region that covers the entirety of the three parallel river basins (TPRB) (90°E–101°E, 27°N–36°N) with an approximate area of 1,207,400 km² (TPRB 399,000 km²), including the sources of the Nujiang, Lancangjiang and Jinshajiang rivers, as well as their basins in the north of Yunnan province, and part of the Big Bend Region of the Brahmaputra (Figure 1). Part of the Brahmaputra is included because there is a glacierized area and the glacier/snow runoff of the area may affect notably on the TWS change in the context of glacier retreat in high mountain Asia. The topography of the region is complicated, with an average altitude exceeding 3000 m, and the terrain fluctuates substantially as it is generally high in the north and low in the south, high in the west and low in the east, and inclines downward from the northwest to the southeast. The study area is within the influence of the southwest and the southeast monsoon. It traverses the subtropical, plateau temperate and plateau subfrozen zones from the south to the north, with substantial differences in climate. Its ecological environment and meteorological conditions are affected not only by the uplift of the Qinghai-Tibet Plateau but also by the formation and evolution of the longitudinal range of the gorge. Dry-hot valleys surrounded by a relatively humid environment with high temperature are very common in this region, such as the TPRB region (Guo et al. 2006; Luo 2009). Complex terrain and special landform, coupled with varied climatic conditions, makes the spatial and temporal variability in water storage more prominent.

The region has been suffering from warming and drying in recent decades (Li et al. 2018a). A rise of 0.7 °C in temperature, and a reduction of 5.5 mm in area-averaged precipitation has been noted in the last 15 years. The arid climate has caused the region to suffer from frequent environmental problems, such as water shortage, grassland degradation, land desertification and soil erosion (Li et al. 2018a). Therefore, the temporal and spatial variations in TWS must be identified as an indicator for regional water resources and their responses to climatic events.

In this study, GRACE Level-2 (RL06) GSM (spherical harmonic coefficients of the geoid model) data produced by the Center for Space Research, University of Texas at Austin (UTCSR), were employed to calculate the TWS. The data were collected from April 2012 to August 2016, with a time resolution of one month, and the spherical harmonic coefficient of gravity was chosen to degree 60 and order 60. The GRACE monthly gravity field model had to be refined because of the error. Satellite laser ranging (SLR) C20 has higher precision and is more sensitive in reflecting seasonal characteristics than GRACE. As a result, the C20 of SLR was used to replace C20 of GRACE (Cheng & Tapley 2003; Chen & Wilson 2008). In addition, improved p3m9 decorrelation filtering was used to eliminate the north–south strip error, and 300 km fan filtering was used to reduce the high-order spherical harmonic coefficient error in the present study (Han et al. 2005; Zhen et al. 2018).

The Global Land Data Assimilation System (GLDAS) datasets (version 2.1), established by the Goddard Space Flight Center (NASA) and the National Centers of Environmental Prediction (NCEP) were used in the present study, which is proved to be an efficient tool for examining meteorology and hydrology (Xia et al. 2014; Yang et al. 2018; Wu et al. 2019). The GLDAS, which includes three land surface models (NOAH, CLM, and Mosaic) and one hydrological model (VIC), provides datasets (global precipitation, soil temperature, runoff, radiation flux and other hydrometeorological data) that are simulated based on multiple observation data, atmospheric assimilation products and reanalysis data for use by various researchers and relevant institutions (Rodell et al. 2004). The four-layer soil moisture data of NOAH have been widely used in the inversion of GRACE-derived TWS (Li et al. 2018b; Zhen et al. 2018; Zhong et al. 2018). The data derived from NOAH from 2002 to 2016 were employed in this study, with a spatial resolution of 1 × 1° and a temporal resolution of one month, and they were mainly used in this study to calculate the scale factor and characterize the soil moisture in the TWS.

We employed GLDAS precipitation, evapotranspiration, and runoff to depict the water balance to minimize the error caused by different data sources. As the largest input to the TWS, precipitation data are crucial for identifying the regional water balance process. Thus, we evaluated the accuracy of GLDAS precipitation here. The accuracy was validated by comparing it with the station data using R-square, bias, root mean square error (RMSE) and relative RMSE (rRMSE) (Junzhi et al. 2012). As shown in Figure 2, GLDAS precipitation displayed good correlation with the station data (R² = 0.45), which generally reflects the characteristics of regional precipitation. However, from the perspective of RMSE and rRMSE, GLDAS was not close to the actual observed value. In addition, GLDAS slightly underestimates the precipitation. Considering that we cannot adequately quantify the errors owing to different data sources, and the GLDAS-derived TWS using the water balance method which was quite consistent with GRACE derived TWS (Figure 3), we suggest that GLDAS precipitation basically satisfies the requirements for analysis, although we also tested precipitation from CN501 (Jia & Xuejie 2013), TRMM.3b43v7 (Huffman et al. 2007), and GPCP (Adler et al. 2018), which exhibited better results than GLDAS in the EQTP.

Observed wells data can accurately reflect groundwater reserves, but the number of wells is small and data acquisition is difficult. In recent years, data from the WaterGAP Global Hydrology Model (WGHM) simulated by IPG (University of Frankfurt) have better reflected changes in groundwater reserves. In addition, the model also reflects the changes in all forms of water, except glaciers. Currently, the results from this model have been widely used (Eicker et al. 2014; Feng et al. 2018). The groundwater data from 2002 to 2013 were used in the present study, with a spatial resolution of 0.5 × 0.5° and a temporal resolution of one month (Döll et al. 2015).

To improve the validity and accuracy of the results, we derived the precipitation, evapotranspiration, wind speed, wind direction, temperature, humidity and other observation data recorded by 36 stations in the study area, obtained from the national meteorological center (http://data.cma.cn/) (Figure 1). After preprocessing the observation data, the time series of the data recorded by the 36 stations from 2002 to 2016 were obtained. The measured data from each station were mainly used to verify the accuracy of the data obtained from the model and evaluate the effects of meteorological factors on TWS.

In the present study, STL is mainly used to decompose time series data, such as TWS, and to analyze the characteristics of each component and extract effective information.

One of the main explanations for the change in the spherical gravity field is the dynamic changes in the terrestrial water mass (Tapley et al. 2004), which is expressed in the form of equivalent water thickness (EWT) using the equation described by Wahr et al. (1998). More information about TWS inversion is described in the first section of the Supplementary material. Errors may occur in a series of postprocessing operations of the GRACE gravity field data (Klees et al. 2007), such as spherical harmonic expansion and spatial filtering (decorrelation filtering, fan filtering, etc.), which cause leakage of the signal (Feng et al. 2017). In this study, we used double scale factors to reduce the leakage error and acquired satisfied results. A detailed description of the processes can be found in the second section of the Supplementary material.

Using the TWS inversion methods described above under ‘Inversion of TWS’, the GRACE data were processed to obtain the water storage equivalent water thickness of the EQTP. Figure 3(a) shows the temporal variations in the regional averaged monthly TWS anomalies. Two TWS results are highly correlated, which indicates the strength of the linear relationship between the GRACE-derived TWS and GLDAS-derived TWS; the coefficient of determination R² (hereafter R²; here R² = 0.84 before removing the seasonal cycle, but R² = 0.52 after removing the seasonal cycle) shows we were able to calculate a non-seasonal variation in TWS of approximately 52% in the EQTP using GLDAS. Therefore, we could evaluate controlling factors of TWS change by analyzing the changes of TWS components derived from GLDAS.

As shown in Figure 3, the TWS showed a decreasing trend from 2002 to 2016, with an approximate rate of –0.15 mm/month, but it showed an increasing trend from 2002 to 2005, with an approximate rate of 0.89 mm/month. The time series was decomposed to the trend, seasonal and residual items with STL to refine the characteristics of TWS. The seasonal variation in TWS was distinguished from the decomposed signal, which showed substantial annual losses in winter and spring and the largest annual surpluses in summer and autumn. Its amplitude was maintained at ±50 mm, and no significant inter-annual variation in the seasonal characteristics was observed. The characteristics of the changes in the trend term were basically consistent with the overall change in TWS, indicating that the change in TWS mainly depended on the trend term. The residual term of TWS reflected the characteristics of TWS fluctuation and was the most powerful evidence for drought or flood events in the region. No obvious pattern of fluctuation in TWS was observed, but some obviously abnormal signals were detected, such as in the spring of 2003 and 2004, and in the winter of 2006 and 2015, etc. We counted the annual fluctuations in TWS to capture the surplus and deficits in TWS, and the results are shown in Table 1. In general, except for certain years, the loss of water storage capacity is greater than 35 mm, resulting in a drying condition; however, the surplus of water reserves is maintained at greater than 30 mm. These findings are consistent with the characteristics of a dry-hot valley. Regarding the balance of TWS, losses were observed for seven of the years between 2003 and 2015, particularly in 2004, 2006 and 2011. According to the statistics from the China Meteorological Disaster Yearbook issued by the China Meteorological Administration, the drought in northern Yunnan and southeastern Tibet was severe in these years, and agriculture and animal husbandry were substantially affected. Therefore, we confirmed that the losses of TWS in these years were related to droughts. From the peak of TWS surplus and deficit, except for 2003, 2004, 2008, 2010 and 2015, the largest losses in other years appeared in the summer and autumn, results that were completely opposite to the seasonal variation in TWS. We speculated that, in these years, the study area experienced a certain degree of drought. More interestingly, in 2012 and 2014, the maximum deficit and maximum surplus appeared only one month apart, which may have been caused by the lagged effect of heavy precipitation on regional soil moisture supplementation.

To understand the spatial differentiation of TWS in the EQTP, we next performed a grid-to-grid process to obtain the long-term trend, which is defined as the monthly average rate of TWS change. As shown in Figure 4, overall, the spatial variation rate of TWS ranged from –0.94 to 0.91 mm/month (the loss in local areas exceeded 1.4 mm/month). The spatial differences were significant. The southwest direction showed a large decreasing trend (defined as R1, the blue line frame in Figure 4), and this direction is the main flow of the Brahmaputra river, and the north direction showed an increasing trend (defined as R2, the yellow line frame in Figure 4), particularly in the area recognized as the source of the Jinshajiang Basin. TWS in the central and eastern parts of the study area (defined as R3, between the blue line and the yellow line in Figure 4) changed slightly and generally maintained a balance. In addition, TWS in the middle and lower reaches of the TPRB experienced a significant loss that was largely related to the unique topography of a typical dry-hot valley. This loss has led to frequent large-scale drought events in recent years. In the TPRB, the Jinshajiang Basin showed an increasing trend overall, with a rate of 0.19 mm/month, and a weak downward trend was observed downstream of the basin. The headstream of the Lancangjiang Basin showed a slight increase in TWS; nevertheless, the TWS of the whole basin generally decreased at a rate of 0.28 mm/month. The Nujiang Basin experienced a vast deficit in TWS and its negative change rate reached 0.49 mm/month.

We used soil moisture (SM), snow water equivalent (SWE) and total amount of water within canopy (CW) from the GLDAS-NOAH model and groundwater (GW) from WGHM model to further explore the characteristics of the variations in each component of TWS. As shown in Figure 5, the regional GW does not display an obvious change trend, consistent with the results reported by Feng et al. (2018), who showed that GW was not the main component contributing to the regional change in TWS. CW shows little contribution to the TWS change and tendency. The change in the SWE tends to be 0, which is mainly because the snow and ice is mainly concentrated in the upstream areas of the basins, and the area covered by ice and snow is much smaller than the total area of the EQTP. Regarding the whole area, the contribution of the SWE to changes in TWS is small, but the amount of change is generally greater than 0 (approximately 0–2 mm), particularly after 2011, indicating that the snow melt in basins is increasing. A continuous time series data representation for the change in glacier meltwater was not available. However, according to previous studies, the glacier in the EQTP has shown a negative mass balance in recent years, namely, the glacier melt water is increasing (Kääb et al. 2015; Berthier et al. 2016; Zhou et al. 2018). SM is the main component of TWS, which is a necessary condition for the survival and growth of regional plants. Drought occurs when SM is unable to meet the needs of plant growth. In 2009–2014, southwestern China, particularly the Yunnan-Guizhou Plateau, continued to experience a drought (Rong et al. 2018), and the drought was the most severe in 2010. As shown in Figure 5, the extent to which the SM decreased during this period was large. Consequently, we inferred that the drought during this period was mainly caused by insufficient soil moisture. Notably, in the southwestern drought in 2010, the reduction in TWS was large, while the SM did not fluctuate substantially, and a rapid decrease in SM occurred in early 2011. This difference may have been caused by the lag effects of P and ET on soil moisture. The decrease in the soil water content in 2011–2015 was greater than TWS. According to Zhang et al. (2015a), an increase in the surface temperature enhanced the regional ET, so the drought was strengthened during the period. At the same time, due to the increase in temperature, the melting of the glaciers was accelerated. Thus, we speculate that the acceleration of glacial ablation and the increase in ET partially contributed to the changes in TWS in the southwest of the region.

In this section, we analyze the causes of spatial differences in water reserves based on changes in precipitation (P), evapotranspiration (ET), and runoff (R). The TWS in R3 remained stable and will not be discussed here. We mainly analyzed the entire study area, zones R1 and R2.

Based on the results described above, in the study area, the variability in ET, P and R was calculated as 0.0180, −0.0107, and −0.0326 mm/month, respectively (Figure 6(a)). The general increase in ET and the decrease in P agreed with the negative TWS change with high correlation coefficient of 0.77 (p < 0.0001) and 0.66 (p < 0.0001), respectively. In terms of the variability rates, ET was nearly twice as high as P, which showed that ET more sensitive to TWS change. Additionally, the R² between ET and TWS was approximately 0.59, which indicated that ET could explain 59% of the variations of TWS. R integrated the effects of P and ET and it was mainly supplied by glacier melting water in this area, and thus it made a limited contribution to the changes in TWS, although the rate of R was even larger than that of ET. Thus, the increase in ET was a primary driver for water storage dynamics across the EQTP during 2002–2016.

R1 displayed a severe TWS deficit, with a loss rate of −0.71 mm/month. As shown in Figure 6(b), the P and ET exhibit higher variability in R1 than in the study region as a whole, and the P shows a positive tendency. This positive trend in P means that the loss of TWS was mostly caused by a decrease in ET. Similarly, we noted that 54% of the variations in TWS could be explained by ET by calculating the R². Figure 6(b) shows that the regional R displays a decreasing trend, which is inconsistent with the variation of P. Thus, we can infer that the water supply from glaciers and snow meltwater was decreasing. Based on the results reported by Brun et al. (2017), the mass balance of representative glaciers in the study area was −0.62 ± 0.23 m w.e. yr^–1, we concluded that the glacial melt variability in R1 was approximately −0.64 mm/month (the glacier area is 4,368.836 km², and the water equivalent of the glacial melt is 620 mm/year). The effect of glacier change on TWS change is indeed not negligible, although we cannot quantify the contribution of glacier meltwater to R in the present study. Therefore, our results suggest that glacier mass loss exerts a significant role in the variation of TWS apart from ET in the R1.

R2 showed an increase in TWS at a rate of 0.49 mm/month. As shown in Figure 6(c), the regional P, R, and ET show changes that are inconsistent with TWS, although the extent of these changes is small, which means that P and ET are not the key factors that affect change in TWS. Figure 1 shows the presence of many large lakes in R2. According to the relevant research, the water volume of lakes in the Qinghai-Tibet Plateau has increased significantly in recent years, and many new lakes have appeared (Li et al. 2019). Therefore, the contribution of the increase in lake water storage to the change in TWS must be considered. Currently, detailed data regarding the amount of change in the water storage capacity of these lakes are not available, thus we can only speculate that the increase in TWS observed in this area is closely related to lake expansion. The recharge source of these lakes is mainly the glaciers in R2 and those located outside the study area, which indirectly shows that glacial changes contribute to the increase in TWS. In addition, in recent years, the depth of the frozen soil active layer on the Tibetan Plateau has significantly increased, and the freezing depth has decreased significantly (Zhang et al. 2015b; Peng 2018). From this perspective, the transfer of mass from the ablation of frozen soil around the study area may also affect the changes in TWS.

To obtain a better understanding of the effect of lake expansion and glacial melt on TWS, we calculated the contribution of each variable. As shown in Figure 7(a), CW and SWE exhibit a slight influence on TWS change and SM dominates the TWS change before 2012. was estimated at 0.35 mm/month from 2002 to 2013 and showed a contribution of approximately 54% to TWS change, which verifies the speculation above. It is noteworthy that GLDAS-derived TWS was quite different from GRACE-derived TWS after 2012 (Figure 7(b)) and the SM showed a continuous downward trend after 2012. The decrease in SM was consistent with the increase in ET and the decrease in R. This suggests that the contribution of to the change in TWS was increasing.

Based on the analysis described above, the increase in ET was the main cause of the decrease in TWS in the EQTP during the study period. Additionally, the decrease in P also contributed to the decrease to some extent. Generally, of all seasons, the highest levels of ET were observed in the northern regions during the summer (defined as June, July and August, and represented by Q3) (Su & Feng 2015), when the study area was also significantly affected by the Indian summer monsoon. Therefore, in this section, we discuss the changes in TWS in the northern region during the summer and the factors that substantially influence ET. Namely, whether the variations in air temperature, air humidity and wind speed converge with the characteristics of change in ET. We will also analyze the variations in the Indian summer monsoon and its feedback effect on P to analyze the factors that affect the change in TWS in the northern region during the summer.

As shown in Figure 8, it was observed that the water deficit in summer dominates the TWS annual trend in the northern region. Meanwhile, P shows a decreasing trend that converges with the varying trends of TWS. Although the rate of decrease in P is small, as the main component of the regional water reserves this small decreasing rate slows the loss of TWS. P varies substantially in regions with complex topographic conditions, but the long-term tendencies of P are mainly affected by atmospheric conditions, and topographic and geomorphic features are amplifiers of atmospheric conditions. During summer, the study area is mainly affected by the Indian summer monsoon (Yao et al. 2017), and thus we used the Indian Monsoon Index (IMI) (Wang & Fan 1999) to measure the relationship between the Indian summer monsoon and P during the study period. In Figure 9, the IMI exhibits a decreasing trend, indicating that the IMI is weakening and the moisture brought by the Indian summer monsoon is decreasing during this period, which is one reason for the overall decrease in regional P. After the summer of 2005, the Indian summer monsoon continues to weaken, but the P only decreases slightly with a larger interannual fluctuation in general, which is particularly large from 2011 to 2015. Thus, the weakening of the Indian summer monsoon is not the main cause of the decrease in P in the EQTP. IMI is more consistent with the change in TWS than the change in P. The correlation coefficient between the TWS change and IMI is 0.53, which indicates that the TWS loss in summer is partly affected by the Indian summer monsoon. To quantify the relationship, we calculated the amplitude of effect of the Indian summer monsoon on TWS using a method that serves to characterize the association between ENSO and TWS (Phillips et al. 2012) (see ‘Quantification of the relationship between TWS change and Indian Monsoon Index’ in the Supplementary Material). The amplitude of 0.018 mm/(IMI*year) suggested that the Indian summer monsoon exerts a small influence on the yearly change in TWS. The nonconformity between the weakness of the Indian summer monsoon and the decrease in P requires further research and we will not delve into this topic in the present study. However, we speculate that the increases in summer temperature and ET (Figure 9), together with the changes in underlying surface in the region (the decrease in vegetation coverage caused by human activities and the decrease in glacier cover area caused by warming in the region, etc.), have resulted in a decrease in the atmospheric water vapor content, which may be another explanation for the regional decrease in P.

As shown in Figure 9, ET is the main factor contributing to the decrease in TWS, and it increases significantly in summer, showing a higher trend than is observed throughout the year. Based on these findings, the decrease in TWS that is mainly observed in summer is generally controlled by the increase in ET. To clarify the reasons for the increasing ET in summer, we selected the meteorological elements that had a greater impact on ET in the region and analyzed their changes. Generally, during the study period, the air temperature increase, the wind speed increase, and the variation in air humidity are consistent with the variation in ET in the region. In terms of the rate of change, the change in relative humidity has the greatest contribution to ET, followed by the wind speed. However, the relative humidity is also affected by the changes in temperature and air pressure. Therefore, the extent of the contributions of the three factors to ET is difficult to determine. Nevertheless, in terms of fluctuation, we determined that changes in temperature were the most consistent with changes in ET, indicating that changes in temperature are the factor with the greatest direct effect on ET. As shown in Figure 9, the change in regional ET is consistent with the change in TWS, but the temperature is consistent with TWS after a lag. Regarding the points at which ET and temperature changed in the time series data, asynchronous correlations were observed between the two variables as well, indicating that the temperature change exerted a certain lag effect on ET.

Overall, from 2002 to 2016, the TWS in the EQTP shows a slight decreasing trend, and the spatial differences were significant. In the years when droughts occurred, TWS fluctuated abnormally. Soil moisture (SM) was the main component of TWS. The drought in the EQTP was mainly due to the continuous decrease in SM. Groundwater (GW) fluctuated substantially during the dry season, particularly in the years of persistent droughts.

The increase in evapotranspiration (ET) was the dominant factor contributing to the decrease in TWS in the EQTP (approximately 59% of the variations in TWS can be explained), while the decrease in precipitation (P) had a certain effect on TWS change as well. In the area R1, 54% of the variation in regional TWS could be explained by the variation in evapotranspiration. In addition, glacier ablation played a significant role in the loss of TWS, although its contribution could not adequately be quantified. In the area where water reserves were increasing (R2), lake expansion and glacial melt played a leading role in changes in TWS. Furthermore, frozen soil around the study area may also have affected the variation in TWS. exhibited a contribution of approximately 54% to the increase in TWS from 2002 to 2013, and its effect was increasing.

The decrease in TWS mainly occurred in summer. The increase in ET in summer, which was the result of an increase in temperature, an increase in wind speed and a decrease in relative humidity, was the main cause of the decrease in TWS. The Indian monsoon exerted a certain effect on the decrease in P in summer, but understanding the dominant factor that contributed to the decrease requires further study.

The authors thank members of the Institute of International Rivers and Eco-Security, Yunnan University, plus Fuming Xie, Yongpeng Gao, Wenfei Miao, Xinxin Qiang and Shimei Duan for their input and stimulating discussions over past few months. We also thank Dr Wei Feng at Institute of Geodesy and Geophysics, Chinese Academy of Sciences, for helpful comments on this manuscript. This work is supported by the National Natural Science Foundation of China (Nos. 41761144075 and 41661144044), the Second Tibetan Plateau Scientific Expedition and Research Program (STEP, No. 2019QZKK0208), the Research Funds for Introducing Talents of Yunnan University (No. YJRC3201702) and Innovation Fund Designated for Graduate Students of Yunnan University (Grant No. 2018Z099).

Data cannot be made publicly available; readers should contact the corresponding author for details.