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

The nature of the heterogeneity of terrestrial water storage (TWS) in the Eastern Qinghai-Tibet Plateau (EQTP) is poorly understood because of the lack of validated datasets and the complex topographical conditions. In this study, monthly GRACE Level 2 Release 6 (RL06) products were employed to characterize TWS changes between April 2002 and August 2016 in the EQTP. Based on the observations and hydrological model output, the dominant factors contributing to the changes in TWS in sub-basins, and areas of TWS decrease and increase were analyzed systematically. We concluded that the TWS in the EQTP showed a slight decreasing trend from 2002 to 2016 with obvious spatial heterogeneity. The decrease in TWS may be attributed to the increase in evapotranspiration, which explains approximately 59% of the variations. In the region where a substantial decrease in TWS was observed, the trend primarily depended on evapotranspiration, and was certainly affected by glacial ablation. Moreover, the expansion of lakes supplemented by glaciers was the main cause of TWS change in the areas where TWS increased. A decrease in TWS mainly occurred in summer and was mainly due to the increase in evapotranspiration because of warming, an increase in wind speed, and a decrease in relative humidity.

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. ), 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. ; Xu et al. ). However, the lack of in-situ observations limits our understanding of the water cycle in this region. 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. ). 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. ).
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. ); 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.

STUDY AREA
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 2 (TPRB 399,000 km 2 ), including the sources of the Nujiang, Luo ). 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. a). 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. a). Therefore, the temporal and spatial variations in TWS must be identified as an indicator for regional water resources and their responses to climatic events.

GRACE RL06
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 ; Chen & Wilson ). 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. ; Zhen et al. ). 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.

Hydrometeorological data
Precipitation, evapotranspiration and runoff data. 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. ). As shown in Figure 2, GLDAS precipitation displayed good correlation with the station data (R 2 ¼ 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 ), TRMM.3b43v7

Measured data
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,

SEASONAL-TREND DECOMPOSITION PROCEDURE BASED ON LOESS (STL)
Since the data used in this study are time series data, they can be decomposed to obtain details and to clarify the 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. 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.

WATER BALANCE
In the vertical direction, regional TWS includes soil moist- changes is expressed as T variable =T ΔTWS :

RESULTS AND DISCUSSION
The overall trend of TWS in the EQTP 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 2 (hereafter R 2 ; here R 2 ¼ 0.84 before removing the seasonal cycle, but R 2 ¼ 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  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.

Spatial differences in TWS in the EQTP
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 According to Zhang et al. (a), 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.

Attributes of spatial heterogeneity in TWS
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 2 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 2 . 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. (), 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 2 , 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. ). 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. b; Peng ).
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

Factors affecting the decrease in TWS in summer
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 ), 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. ), and thus we used the Indian Monsoon Index (IMI) (Wang & Fan ) 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   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. 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.

CONCLUSIONS
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. T ΔSW þ T ΔG 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.