Spatiotemporal Variation and Driving Analysis of Groundwater in the Tibetan Plateau Based on GRACE Downscaling Data

: The special geographical environment of the Tibetan Plateau makes ground observation of Ground Water Storage (GWS) changes difﬁcult, and the data obtained from the GRACE gravity satellites can effectively solve this problem. However, it is difﬁcult to investigate the detailed GWS changes because of the coarser spatial resolution of GRACE data. In this paper, we constructed a 0.1 ◦ resolution groundwater storage anomalies (GWSA) dataset on the Tibetan Plateau from 2002 to 2020 based on a phased statistical downscaling model and analyzed the spatiotemporal variation and driving factors of the GWSA in order to better study the changes of GWS on the Qinghai Tibet Plateau. The results show that: (1) In the Tibetan Plateau and 12 sub-basins, the GWSA before and after downscaling show a very high correlation in time series and relatively good performance in spatial consistency, and the downscaled GWSA indicate a consistent trend with the measured groundwater level. (2) The GWSA on the Tibetan Plateau shows a downward trend ( − 0.45 mm/yr) from 2002 to 2020, and the variation trend of the GWSA in the Tibetan Plateau shows signiﬁcant spatial heterogeneity. (3) The GWSA changes in the Tibetan Plateau are mainly dominated by natural factors, but the inﬂuence of human activities in individual sub-basins can not be ignored. Among the teleconnection factors, El Nino-Southern Oscillation Index (ENSO) has the greatest inﬂuence on the GWSA on the Tibetan Plateau.


Introduction
In recent years, the rapid growth of population, economy and technology has made human activities a non-negligible driver of climate change. The impact of climate change will be reflected in the water cycle, and the resulting water security issue is becoming a top policy issue, internationally [1]. Groundwater is an important source of water for human life, agriculture, and industry [2,3]. Groundwater is replenished in the wet season and released slowly in the dry season, which can alleviate the pressure on surface water resources and maintain the basic functions of rivers and wetlands. The traditional monitoring of Ground Water Storage (GWS) depends on the water level of wells, which is accurate but costly, and is difficult to implement over large-scale areas [4,5]. Therefore, effectively obtaining a high-precision, large-scale GWS dataset has become one of the urgent hydrology problems to be solved in recent years.
With the development of remote sensing technology, remote sensing hydrological monitoring can effectively fill the gap of ground well monitoring and obtain key information related to the changes of the GWS on a larger spatiotemporal scale [6]. The Gravity Recovery and Climate Experiment (GRACE) satellites and its follow-up mission, GRACE-FO, can obtain information on Terrestrial Water Storage Anomalies (TWSA) by monitoring the change of Earth's gravity field. GRACE data has the potential to provide changes in global

GRACE Data
GRACE data provide the basis for understanding the global terrestrial water cycle, ice cap and glacier mass balance, sea-level changes, and their link to global climate change [30]. The inversion of TWSA using GRACE data can be obtained in two forms: the Spherical Harmonic (SH) approach and the mass concentration computation (mascon) method. This study uses the more commonly used GRACE-FO RL06 mascon solution data provided by the University of Texas at Austin, Center for Space Research (UTCSR), http://www2.csr.utexas.edu/grace/RL06_mascons.html (accessed on 1 March 2021) [31]. The outstanding advantage of this data is that there is no need for additional post-processing, such as smoothing or striping, and the boundary between regions is clear [32]. The study period is from April 2002 to October 2020 and the missing data is filled with linear interpolation. GRACE data are usually provided in the form of TWSA anomalies relative to an average baseline from January 2004 to December 2009. The research period was from 2002 to 2020, and it may be appropriate to select the middle period of the study period (2007)(2008)(2009)(2010)(2011)(2012) as the base average baseline. In view of the relatively long research period of 2002-2020, the mean baseline is from January 2007 to December 2012 to reflect the groundwater change characteristics.

GLDAS Data
The Global Land Data Assimilation System (GLDAS) hydrological model is based on the principle of the water balance cycle and energy cycle to drive data from ground

GRACE Data
GRACE data provide the basis for understanding the global terrestrial water cycle, ice cap and glacier mass balance, sea-level changes, and their link to global climate change [30]. The inversion of TWSA using GRACE data can be obtained in two forms: the Spherical Harmonic (SH) approach and the mass concentration computation (mascon) method. This study uses the more commonly used GRACE-FO RL06 mascon solution data provided by the University of Texas at Austin, Center for Space Research (UTCSR), http://www2.csr. utexas.edu/grace/RL06_mascons.html (accessed on 1 March 2021) [31]. The outstanding advantage of this data is that there is no need for additional post-processing, such as smoothing or striping, and the boundary between regions is clear [32]. The study period is from April 2002 to October 2020 and the missing data is filled with linear interpolation. GRACE data are usually provided in the form of TWSA anomalies relative to an average baseline from January 2004 to December 2009. The research period was from 2002 to 2020, and it may be appropriate to select the middle period of the study period (2007)(2008)(2009)(2010)(2011)(2012) as the base average baseline. In view of the relatively long research period of 2002-2020, the mean baseline is from January 2007 to December 2012 to reflect the groundwater change characteristics.

GLDAS Data
The Global Land Data Assimilation System (GLDAS) hydrological model is based on the principle of the water balance cycle and energy cycle to drive data from ground and satellite observations to generate land surface fluxes. The GLDAS model has been well validated in the field of hydrometeorology and can effectively fill the data gaps in areas lacking ground-based observations [33,34]. The GLDAS NOAH model is used to extract monthly runoff and evapotranspiration data for the Tibetan Plateau from April 2002 to October 2020. The GLDAS NOAH2.1 model monthly dataset with a 0.25 • spatial Water 2022, 14, 3302 4 of 23 resolution was released by the National Aeronautics and Space Administration Land Data Assimilation System (http://ldas.gsfc.nasa.gov/gldas/, accessed on 1 March 2021) [35]. As with the GRACE data, we calculated the mean values of GLDAS data from January 2007 to December 2012 to obtain the anomalies.

Precipitation Data
The Integrated Multi-satellite Retrievals for GPM (IMERG) combines infrared, microwave, and gauge observations from multiple satellites to provide relatively high spatial and temporal resolution precipitation estimates [36]. The IMERG precipitation products, which have been widely used, cover a wide range and can accurately detect trace precipitation [37]. We selected monthly GPM_IMERG V6 products from April 2002 to October 2020 with a spatial resolution of 0.1 • [38]. We also calculated the mean values of precipitation data from January 2007 to December 2012 to obtain the anomalies.

Groundwater Observation Data
Groundwater level observations from the Lhasa and Maoxian stations of the Chinese Ecosystem Research Network (CERN) are collected from 2005 to 2014 to verify the accuracy of the downscaled GWSA data, and the locations of the observation wells are shown in Figure 1. Because the data of these two observation wells are relatively complete, these data can verify the accuracy of the downscaled GWSA data. In view of the data limitation, we selected the data of these two locations. The anomalies of groundwater level observations are removed by the Normalized Median Absolute Deviation (NMAD) method [39]. After the abnormal values are removed, the time series of groundwater observation data can be obtained.

Teleconnection Data
Four large-scale climate factors from 2002 to 2020 are selected to study the relationship between the GWSA and teleconnection factors in the Tibetan Plateau, as shown in Table 1. The teleconnection data include AOI, ENSO, NAOI, and PDOI.

Method
Firstly, the rainfall, runoff, and evapotranspiration data with a resolution of 1 • are used to calculate the TWSA (water balance method). Secondly, the regression function is constructed between the TWSA (water balance method) and the GRACE-based TWSA (GTWSA). Thirdly, according to the regression function, the high-precision downscaled GRACE-based TWSA data is obtained. Finally, the downscaled ground water storage anomalies (GWSA) can be calculated. The flow chart representing the methodology is depicted in Figure 2.

Calculation of Groundwater Storage
The water balance method is the basis for the study of hydrological processes and the calculation of water storage [40]. Before the gravity satellite was launched and used, the water balance method was usually selected. The equation of water balance in the watershed system is expanded as follows: where TWSA is calculated based on the water balance method; PA is precipitation anomalies; QsA is runoff anomalies; EvapA is evapotranspiration anomalies. The GLDAS model provides monthly evapotranspiration and runoff data, and IMERG provides monthly precipitation data. In order to facilitate the calculation, the dimension is converted into a unified unit (mm). Because canopy water storage is small in the research area, the surface water in the study area is composed of snow depth, runoff, and soil moisture [41]. Therefore, the GWSA is calculated by the Formula (2): where SWEA is the snow depth anomalies, QsA is the surface runoff anomalies, and SMA is the soil moisture anomalies of 0-200 cm. The three data are all obtained from the GLDAS model, and GTWSA is obtained from the GRACE data.

Phased Statistical Downscaling Model
GTWSA is the terrestrial water storage anomalies obtained from the GRACE data, which is calculated according to GRACE-FO RL06 mascon solutions. Although the CSR GRACE mascon data are represented by a grid of 0.25°, they represent an equal area geodesic grid with a size of 1 × 1° at the equator, so the CSR-GRACE mascon data are

Calculation of Groundwater Storage
The water balance method is the basis for the study of hydrological processes and the calculation of water storage [40]. Before the gravity satellite was launched and used, the water balance method was usually selected. The equation of water balance in the watershed system is expanded as follows: where TWSA is calculated based on the water balance method; PA is precipitation anomalies; QsA is runoff anomalies; EvapA is evapotranspiration anomalies. The GLDAS model provides monthly evapotranspiration and runoff data, and IMERG provides monthly precipitation data. In order to facilitate the calculation, the dimension is converted into a unified unit (mm). Because canopy water storage is small in the research area, the surface water in the study area is composed of snow depth, runoff, and soil moisture [41]. Therefore, the GWSA is calculated by the Formula (2): where SWEA is the snow depth anomalies, QsA is the surface runoff anomalies, and SMA is the soil moisture anomalies of 0-200 cm. The three data are all obtained from the GLDAS model, and GTWSA is obtained from the GRACE data.

Phased Statistical Downscaling Model
GTWSA is the terrestrial water storage anomalies obtained from the GRACE data, which is calculated according to GRACE-FO RL06 mascon solutions. Although the CSR GRACE mascon data are represented by a grid of 0.25 • , they represent an equal area geodesic grid with a size of 1 × 1 • at the equator, so the CSR-GRACE mascon data are resampled to data with a resolution of 1 • . When the TWSA 1 • of the Tibetan Plateau is Water 2022, 14, 3302 6 of 23 moved backward one month relative to the GTWSA 1 • , the correlation coefficient increases from 0.33 to 0.46 using cross-correlation, so there is a one-month lag relationship between the TWSA 1 • and the GTWSA 1 • . Additionally, the Root Mean Squared Error (RMSE) is calculated between the GRACE-based TWSA and downscaled TWSA, in order to better understand the model's ability. The statistical downscaling steps are as follows: Step 1: according to Formula (1), using rainfall, runoff, and evapotranspiration data with a resolution of 1 • , TWSA 1 • is calculated; Step 2: the resolution of the CSR GRACE mascon data is resampled to 1 • to obtain the GTWSA 1 • ; Step 3: the regression function relationship between the TWSA 1 • and the GTWSA 1 • is constructed by Formula (3): Where G TWSA is the TWSA fitted by the linear regression model. The difference between the G'TWSA 1 • and the GTWSA 1 • is calculated as the residual of the downscaling system; Step 4: according to the Formula (1), using the monthly Evap, P and Qs data with a resolution of 0.1 • , the TWSA 0.1 • is calculated, and the G'TWSA 0.1 • is obtained according to the Formula (3); Step 5: deducting the system residual from step 3 in the G'TWSA 0 . 1 • , then the downscaled TWSA with a resolution of 0.1 • is obtained; Step 6: calculate the downscaled GWSA according to the Formula (2). Since GTWSA and TWSA time series in different periods are fluctuating, it is difficult to get a good regression function relationship in step 3. We divided the whole period (2002 to 2020) into several stages based on the Bernaola Galvan segmentation algorithm. The Bernaola Galvan segmentation algorithm is a suitable method for phase division of nonsmooth, nonlinear time series. Compared with traditional methods, the Bernaola Galvan segmentation algorithm can segment non-stationary sequences into multiple stationary subsequences with different means [42,43] 2020.10), as shown in Figure 3. After the phase division of the GRACE TWSA time series, the statistical regression relationship between the GTWSA (GRACE) and the TWSA (water balance method) is established by using a phased linear regression model. As shown in Figure 4, the original relationship between GTWSA and TWSA is low (R 2 = 0.212) for the whole period (2002-2020). After the segmentation, the relationship was enhanced obviously in three stages, A (2002.4-2007.4 with R 2 = 0.581), B (2007.5-2015.12 with R 2 = 0.68) and C (2016.1-2020.10 with R 2 = 0.516). A phased statistical downscale model is proposed to improve the effectiveness for obtaining downscaled GWSA.

Spatiotemporal Analysis Method
In this study, the ESMD method is used to decompose the time series of the GWSA to study the variation period and trend of the GWSA in the Tibetan Plateau. The ex-

Extreme-Point Symmetric Mode Decomposition
In this study, the ESMD method is used to decompose the time series of the GWSA to study the variation period and trend of the GWSA in the Tibetan Plateau. The extremepoint symmetric mode decomposition (ESMD) method is an improvement of the traditional Hilbert-Huang transform based on the Fourier transform. ESMD has strong adaptability and local characteristics and is suitable for nonlinear and non-stationary time series analysis [44]. The ESMD method decomposes the original time series into a series of intrinsic mode function (IMF) components and a trend (R). The period of each IMF component can be obtained by a fast Fourier transform, and the trend term R represents the overall trend of the whole time series [45].

The Modified Mann-Kendall Trend Test
The MMK method is applied to reveal the change trend characteristics of the GWSA based on grid scale. The modified Mann-Kendall (MMK) trend test method can eliminate the autocorrelation component in the series and significantly improve the testability of the Mann-Kendall method, which is a nonparametric test method to effectively detect the changing trend of time series and is widely used in the field of hydro-meteorology [46,47].

Groundwater Storage Changes Attribution Analysis
The multiple linear regression model is usually used to study the relationship between changes in a dependent variable dependent on multiple independent variables and is one of the most widely used regression models, as shown in Equation (4) [48]. After standardizing all variables, the quantitative interpretation of the contribution of each variable can be given based on the standardized coefficients, as shown in Equation (5) [49].
where Y is the dependent variable, X i is the independent variable, a 0 is the intercept, b i is the standardized coefficient, and η i is the contribution rate of the independent variable. The evapotranspiration and runoff data of the GLDAS NOAH2.1 model and the precipitation data of the GPM_IMERG V6 precipitation product are selected as natural factor data, and the prolonged artificial nighttime-light dataset of China is selected as the human activity data [50]. Since the natural factor of glacier melting is not considered, the contribution of human activities may be overestimated. The contribution rate of natural factors and human activities to climate change is calculated by the standardized coefficients of the multiple linear regression model, and then the attribution of climate change is quantitatively analyzed.

Cross Wavelet Transform
The cross wavelet transform can reveal the common power and relative phase in timefrequency space between two sequences, and the region with large transform coefficient represents that the two signals have strong correlation. In this study, the cross wavelet transform method is adopted to clarify the influence of teleconnection factors on the change of the GWSA on the Tibetan plateau.  Figure 5. We use monthly cell average value (2002.4-2020.10) representing the whole region or sub-basin. Overall, the correlations between the downscaled data and GRACE data in three time series are high (r > 0.997) and the root mean square error (RSME < 2.75 mm) is low. Among them, the downscaled data in section B has the lowest RMSE (1.77 mm), and the downscaled data in section C has the highest correlation coefficient (r = 0.9992). Therefore, the time series data before and after downscaling are highly consistent. Spatially, the Tibetan Plateau can be subdivided into 12 sub-basins, and the correlation coefficients of the data before and after the downscaling of each sub-basin fluctuate between 0.9878-0.9997, and the overall correlation coefficient can reach 0.9967, and all of them are significantly correlated. The RSME of the 12 sub-basins is relatively small, and the maximum RSME of the Ganges River Basin is 24.34 mm. The smallest RSME in the Qaidam Basin is 1.91 mm, and the RSME in the Tibetan Plateau as a whole is 2.34 mm. The downscaling data almost completely maintain the time series variation characteristics of GRACE data, and the conservation of mass in sub-basins also existed objectively.
Water 2022, 14, x FOR PEER REVIEW 9 of 23 reach 0.9967, and all of them are significantly correlated. The RSME of the 12 sub-basins is relatively small, and the maximum RSME of the Ganges River Basin is 24.34 mm. The smallest RSME in the Qaidam Basin is 1.91 mm, and the RSME in the Tibetan Plateau as a whole is 2.34 mm. The downscaling data almost completely maintain the time series variation characteristics of GRACE data, and the conservation of mass in sub-basins also existed objectively. The spatial comparison of the GTWSA and downscaled GTWSA data in January 2003 and September 2012 were tested in the 12 sub-basins and the Tibetan Plateau. As shown in Figure 6, the downscaled data had finer spatial variation characteristics than the original GRACE data, especially in regions with larger absolute values. In Figure 7, by using 200 random sampling points, corresponding cell value were used to verify the reliability of the downscaled results in the 12 sub-basins and the Tibetan Plateau. Overall, the correlations are high (R 2 > 0.77) between the downscaled GTWSA and GTWSA in these two periods. Spatially, the correlation coefficients in 12 sub-basins are various. Most of them are higher than 0.7, and the correlation coefficients in Qaidam are 0.866 and 0.929 in two periods. However, the correlation coefficients in both periods are very low in the Hexi Corridor. The spatial comparison of the GTWSA and downscaled GTWSA data in January 2003 and September 2012 were tested in the 12 sub-basins and the Tibetan Plateau. As shown in Figure 6, the downscaled data had finer spatial variation characteristics than the original GRACE data, especially in regions with larger absolute values. In Figure 7, by using 200 random sampling points, corresponding cell value were used to verify the reliability of the downscaled results in the 12 sub-basins and the Tibetan Plateau. Overall, the correlations are high (R 2 > 0.77) between the downscaled GTWSA and GTWSA in these two periods. Spatially, the correlation coefficients in 12 sub-basins are various. Most of them are higher than 0.7, and the correlation coefficients in Qaidam are 0.866 and 0.929 in two periods. However, the correlation coefficients in both periods are very low in the Hexi Corridor.         Figure 11a shows the monthly average GWSA on the Tibetan Plateau during the study period. The water balance on the Tibetan Plateau shows a large intra-annual oscillation, with the GWSA reaching its lowest value in July and August and its highest value in March and June. According to geological knowledge, the plateau climate is dominated by monsoon winds, and tropical tropospheric easterly winds, subtropical westerly winds and stratospheric easterly winds can bring a lot of precipitation to the plateau, but in summer, precipitation only accumulates on the surface of the plateau, and it takes some time for groundwater to produce flow. After October, the plateau climate enters a dry and cold period, when the rivers dry up and freeze, and groundwater replenishment is impeded, while in the spring, groundwater rises significantly when the temperature rises and large amounts of snow and ice melt. The frequency analysis of seasonal maxima is carried out for each sub-basin to better understand the spatial characteristics of GWSA in different seasons. Spring on the Qinghai-Tibet Plateau is March, April and May, and the specific months of summer, autumn and winter are obtained in turn. The number of months containing the maximum value in each season is counted, and the season with the highest number represents the seasonal color of the whole basin, as shown in Figure 11b Figure 11a shows the monthly average GWSA on the Tibetan Plateau during the study period. The water balance on the Tibetan Plateau shows a large intra-annual oscillation, with the GWSA reaching its lowest value in July and August and its highest

Periodic Characteristics of GWSA
As shown in Figure 12, the original GWSA time series is decomposed into five IMF components (IMF1-IMF5) and one trend term (R). The sum of these IMF components and trend term completely coincides with the original GWSA sequence, which indicates that the ESMD method is complete and the decomposition result is reliable. The GWSA on the Tibetan Plateau shows a trend of increasing, then decreasing, and finally stabilizing from  Table 2. The variance contribution rate of 4.65 years for IMF5 is the largest, with a value of 37.58% and a very obvious oscillation signal, followed by the variance

Periodic Characteristics of GWSA
As shown in Figure 12, the original GWSA time series is decomposed into five IMF components (IMF1-IMF5) and one trend term (R). The sum of these IMF components and trend term completely coincides with the original GWSA sequence, which indicates that the ESMD method is complete and the decomposition result is reliable. The GWSA on the Tibetan Plateau shows a trend of increasing, then decreasing, and finally stabilizing from  Table 2. The variance contribution rate of 4.65 years for IMF5 is the largest, with a value of 37.58% and a very obvious oscillation signal, followed by the variance contribution rate of 0.43 years period for IMF1, with a value of 23.14%. Therefore, the variability of the GWSA time series in the Tibetan Plateau is mainly determined by IMF5 and IMF1, with the first and second principal periods of 4.65 and 0.43 years, respectively. variability of the GWSA time series in the Tibetan Plateau is mainly determined by IMF5 and IMF1, with the first and second principal periods of 4.65 and 0.43 years, respectively.

Trend Characteristics of GWSA
The GWSA trend characteristics of the Tibetan Plateau from 2002 to 2020 based on the MMK gridded trend test method are shown in Figure 13. The GWSA trend of the Tibetan Plateau exhibits significant spatial heterogeneity at the 0.05 level, which also reflects the unique topographic and climatic conditions of different regions of the Tibetan

Trend Characteristics of GWSA
The GWSA trend characteristics of the Tibetan Plateau from 2002 to 2020 based on the MMK gridded trend test method are shown in Figure 13. The GWSA trend of the Tibetan Plateau exhibits significant spatial heterogeneity at the 0.05 level, which also reflects the unique topographic and climatic conditions of different regions of the Tibetan Plateau.
Specifically, the GWSA in the western part of the Tibetan Plateau and the eastern part except the Yangtze basin indicates a decreasing trend, with the Brahmaputra basin showing the most obvious decreasing trend, reaching −13.84 mm/yr. However, the GWSA in the Inner, Tarim, and Yangtze basins show an upward trend, in which the rising trend of Inner is the most obvious, reaching 4.40 mm/yr, which may be affected by glacier melting caused by climate change in recent years [51]. Specifically, the GWSA of the Lhasa station shows a downward trend, while the GWSA of the Maoxian station shows an upward trend, which is consistent with the results obtained in Section 4.1.2.
Water 2022, 14, x FOR PEER REVIEW 16 of 23 Plateau. Specifically, the GWSA in the western part of the Tibetan Plateau and the eastern part except the Yangtze basin indicates a decreasing trend, with the Brahmaputra basin showing the most obvious decreasing trend, reaching −13.84 mm/yr. However, the GWSA in the Inner, Tarim, and Yangtze basins show an upward trend, in which the rising trend of Inner is the most obvious, reaching 4.40 mm/yr, which may be affected by glacier melting caused by climate change in recent years [51]. Specifically, the GWSA of the Lhasa station shows a downward trend, while the GWSA of the Maoxian station shows an upward trend, which is consistent with the results obtained in Section 4.1.2.

GWSA Change Attribution Analysis
The contribution of natural factors and human activities to the GWSA changes on the Tibetan Plateau and its sub-basins are quantified based on the multiple linear regression, as shown in Table 3. The goodness-of-fit of the multiple linear regression models on the sub-basins of the Tibetan Plateau ranged from 0.21 to 0.89, with most of the basins having a goodness-of-fit of more than 0.6, indicating the reliability of the multiple linear regression model results. The GWSA changes of the Tibetan Plateau are dominated by natural factors, which is consistent with the sparsely populated nature of the Tibetan Plateau. However, there are some individual basins where human activities dominate. Due to the climate warming caused by human activities in recent years, there exists a non-negligible accelerated glacier melting on the Tibetan Plateau, which in turn has an impact on the spatial and temporal changes of groundwater. However, due to the difficulty in obtaining accurate glacier melt data, the natural factor of glacier melt was not considered in this study, which may have overestimated the impact of human activities.

GWSA Change Attribution Analysis
The contribution of natural factors and human activities to the GWSA changes on the Tibetan Plateau and its sub-basins are quantified based on the multiple linear regression, as shown in Table 3. The goodness-of-fit of the multiple linear regression models on the sub-basins of the Tibetan Plateau ranged from 0.21 to 0.89, with most of the basins having a goodness-of-fit of more than 0.6, indicating the reliability of the multiple linear regression model results. The GWSA changes of the Tibetan Plateau are dominated by natural factors, which is consistent with the sparsely populated nature of the Tibetan Plateau. However, there are some individual basins where human activities dominate. Due to the climate warming caused by human activities in recent years, there exists a non-negligible accelerated glacier melting on the Tibetan Plateau, which in turn has an impact on the spatial and temporal changes of groundwater. However, due to the difficulty in obtaining accurate glacier melt data, the natural factor of glacier melt was not considered in this study, which may have overestimated the impact of human activities.

Teleconnection Driving Forces on GWSA
Groundwater changes are not only affected by natural factors and human activities, but also highly influenced by the teleconnection factors. The influence of teleconnection factors (AOI, ENSO, NAOI, PDOI) on the change of the GWSA on the Tibetan plateau is revealed by cross wavelet transform, as shown in Figure 14. There are three significant resonance periods between the GWSA and AOI, which  Figure 14a. Figure 14b shows that there is a long-term significant resonance period between the GWSA and ENSO that almost covers the study period and is largely positively correlated. Figure 14c  connection factors on the Tibetan Plateau. The teleconnection factors have a significant statistical correlation with the GWSA on the Tibetan Plateau, indicating that they play an important role in the change of the GWSA on the Tibetan Plateau. It can be seen from the cross wavelet transform that the area of significant correlation between ENSO and GWSA is the largest among the teleconnection factors. In general, ENSO has the greatest influence on the GWSA on the Tibetan Plateau, while AOI, NAOI, and PDOI have similar and smaller effects on the GWSA on the Tibetan Plateau.

Discussion
Currently, groundwater storage estimation based on GRACE data has been widely used in many areas [52][53][54]. Recognizing how GWS responds to climate change is essential for future water availability and water resources management in the Tibetan Plateau. Compared with other areas, glaciers contribute to the most GWS variations in the Tibetan Plateau. However, the dominant factor controlling long-term GWS variation and its responses to climate change are not well understood. Famous for the Asian Water Tower, the Tibetan Plateau plays a critical role in providing water resources for more than two billion people living downstream [55]. A better understanding of GWS changes is essential for water resources management to ensure water security and food production [56].
Furthermore, GRACE data can be used to assess regional groundwater storages in deserts and remote mountainous areas [57,58]. Related research shows that GRACE data can make up for the shortcomings of groundwater monitoring networks and show great potential in tracking changes in groundwater storages [59]. However, GRACE data have the problem of low spatial resolution. Normally the spatial resolution of GRACE data is about 350 km, which is higher than most groundwater research areas scale [60]. Chambers [61] performed first-order polynomial fitting of all even (odd) orders of the same order to improve the resolution of GRACE data when calculating changes in seawater mass. Gao et al. [62] combined the distribution of regional hydrogeological parameters to refine the GRACE data (downscaling). Our results showed that the trend of GWSA in the Tibetan Plateau exhibited significant spatial heterogeneity, which was consistent with the findings of Xiang et al. [63]. The GWSA in the Inner, Tarim, and Yangtze basins showed an increasing trend, which was comparable to Zhang et al. [64]. The western part of the Tibetan Plateau and the eastern region, except the Yangtze basin, showed a decreasing trend in GWSA, which seemed to be different from Zou et al. [65]. Zou et al. concluded that the GWS increased during 2003-2016 in all subdivisions of the Tibetan Plateau, except for the Inner Plateau. However, GRACE-based GWS change estimates still contain large uncertainties [66]. In contrast, this paper used the linear regression method to downscale the terrestrial water storage, constructing a GWSA dataset based on a phased statistical downscaling model. The downscaled GWSA is reliable in both time and space, and is consistent with the groundwater level trend of ground observation wells. Therefore, the phased statistical downscaling model proposed in this paper is reliable in the Tibetan Plateau. Furthermore, this study provides a comprehensive and detailed analysis of the spatial and temporal characteristics, periodic patterns, spatiotemporal trends, and driving factors of groundwater storage changes on the Tibetan Plateau, which provides a reference for groundwater resource management on the Tibetan Plateau. It is useful and convenient to use this method to study other areas by considering temporal variation over a long period.
GWSA change attribution analysis results indicate that the GWSA changes of the Tibetan Plateau are dominated by natural factors (Table 3). Noticeably, the average contributions of natural factors and human factors are 62% and 38%, respectively. Among the natural factors, compared with precipitation and runoff, evaporation is the main factor affecting the change of GWSA, with an average value of 27%. In total, the attributions of 10 basins are more than 20%. Furthermore, the human factors are mainly reflected by the contribution of night lights, which ranges from 20% (in the AmuDayra and Hexi basins) to 61% (in the Qaidam basin). The contribution of human activities to the change of groundwater storage in the Qaidam basin is 61%, which may be related to the construction of hydraulic facilities for agricultural development in the Qaidam basin in recent years and the extensive ecological restoration works in the area [67]. As for the AmuDayra and Hexi watersheds, where the contribution of natural factors is high, the former is an extensive uninhabited area and the latter is adjacent to the Qilian Mountains National Nature Reserve, where development has been restricted in recent years.
Among teleconnection factors, ENSO, a strong signal of global climate change, is a climatic phenomenon that occurs when the ocean and atmosphere interact and lose their balance, often arousing the atmospheric circulation change, leading to global climate anomalies [68]. Lu et al. [69] found that North China had the strongest response to ENSO by analyzing the correlation between the affected area ratio and the disaster area ratio of drought in various regions of China from 1949 to 1990. Zhang et al. [70] showed that ENSO events have a greater impact on the drought and flood changes in the Qinling Mountains. In general, ENSO has the greatest influence on the GWSA on the Tibetan Plateau. Zhu et al. [71] have also shown that ENSO has a certain impact on the change of groundwater storage in the southeast of the Tibetan Plateau, and the impact of ENSO is relatively strong in areas with significant loss of water storage.
However, there are still shortcomings in the current study. First, there are several institutions providing GRACE data products, and this paper only selects one of them for calculation. However, some studies have shown that the arithmetic averaging of GRACE data products from different institutions gives more accurate data on the TWSA [72,73]. Filling in the missing GRACE data with linear interpolation may also bring about uncertainties [74]. Secondly, the statistical downscaling method relies on the correlation between model data and GRACE data, and the method needs to be carried out in the edge-regularized area and then clipped according to the basin boundary, which may introduce more edge leakage errors. Thirdly, limited by the availability of data, the number of observation wells used for verifying the accuracy of the downscaled GWSA is a little small. Only the time series data of the two observation stations are relatively complete, which is one of the limitations of this paper. Finally, the GWSA is obtained by subtracting precipitation, evapotranspiration, and snow from the GRACE TWSA, but does not take into account the influence of lakes and reservoirs, so the accuracy of downscaled GWSA can be further improved. In the future, groundwater level fluctuations can also be simulated and predicted with a desirable accuracy using artificial neural network and intelligence models [75,76].

Conclusions
To improve the low spatial resolution of gravity satellite data products, this study successfully downscaled the GWSA on the Tibetan Plateau based on a phased statistical downscaling model and tested the reliability of the downscaling results at the spatial and temporal scales. Second, the spatiotemporal variation of groundwater storage from 2002 to 2020 is analyzed by using the downscaled GWSA data. Finally, a quantitative attribution analysis of the drivers of the GWSA changes on the Tibetan Plateau is conducted based on a multiple linear regression model. Then, the cross wavelet transform is used to analyze the relationship between groundwater and teleconnection factors. The main conclusions are as follows: (1) The correlation coefficients of the data before and after downscaling exceed 0.99 in the all-time series, and the RSME is low, while the data before and after downscaling in the 12 sub-basins also show high correlation. The groundwater levels measured at the Lhasa and Maoxian stations and the GWSA after downscaling of the corresponding grid show consistent trends, and are significantly correlated at the 0.01 level. The downscaling of GRACE and the study of the temporal and spatial variation of groundwater in the Tibetan Plateau still have certain challenges. Accurate GWSA data are helpful to land use, urban planning, and water resources management. This study provides high-resolution GWSA data products on the Tibetan Plateau, which can effectively provide regional groundwater storage change dynamics and provide data support for managers to formulate policies, such as balancing the relationship between regional irrigation water use and groundwater consumption, and can help evaluate regional groundwater consumption more effectively.