2.2.1. GRACE/FO

The CSR GRACE/FO RL06 mascon dataset has been frequently used in the study of terrestrial water storage changes due to its highest spatial resolution at this stage compared with other products [40–42]. The spatial resolution was 0.25 degrees. We used the GRACE data for 237 consecutive months (from April 2002 to December 2021). The cubic spline interpolation method was used to interpolate the short-term missing data of GRACE in the time series. For the long-term data interruption caused by the replacement of the GRACE and GRACE-FO satellites, the reconstruction method based on precipitation, which was invented in previous research, was used, and more detailed information regarding this method can be found in previous studies [43,44]. Additionally, all the grids used in this study were relative to the 2004–2009 mean baseline. The TWS was recorded in terms of equivalent water height (EWH) in cm.

### 2.2.2. GLDAS

The Global Land Data Assimilation System (GLDAS) uses data assimilation technology to fuse satellite-based data and in situ ground observation data to generate surface state quantities and fluxes that are closest to the observation data [45,46]. GLDAS includes four land surface models, NOAH, VIC, CLM, and MOSAIC. Compared with the other models, the NOAH model has the advantages of being an advanced model, having a high spatial resolution, and having a stable driving field. Thus, the model selected in this study was the GLDAS Noah Land Surface Model L4 monthly 0.25 × 0.25-degree V2.1 [47,48]. It has been widely used worldwide and is a highly recognized data source [24,45,46,49]. The spatial resolution is 0.25 degrees and is consistent with GRACE.

### 2.2.3. Temperature and Precipitation

Temperature (Tmp) and precipitation (Pre) are the most widely studied and important meteorological factors. Previous studies also showed that Tmp and Pre have significant effects on the variation in GWS [23,50].

The Climatic Research Unit gridded Time Series (CRU TS) can provide 0.5-degreeresolution monthly data covering the global land surface [51]. The dataset of Tmp and Pre from the CRU TS from April 2002 to December 2021 was used in this study. It is worth noting that the spatial resolution was interpolated to 0.25 degrees using the bilinear interpolation method to match the resolution of GRACE.

### *2.3. Vegetation Response*

The vegetation response to the GWS changes was represented by the Normalized Difference Vegetation Index (NDVI) and Land Use and Cover Change (LUCC). The satellitebased NDVI, which has become the most popular index to reflect the state of regional climates and environments, is obtained by monitoring the vegetation growth status and estimating the vegetation coverage without damaging or altering the vegetation [52,53]. The satellite-based LUCC is one of the essential driving factors for regional climate variability, and its impact on hydrological cycle has become an important aspect of water resources [54,55]. The monthly data of the NDVI from April 2002 to December 2019 and the LUCC data of five periods in 2000, 2005, 2010, 2015, and 2020, were used in this study. The NDVI dataset is generated by calculating the maximum values of the first, middle, and last three ten days of each month. The LUCC types are divided into six first-level groups, which include cropland, forest, grassland, water, urban land, and unused land.

### *2.4. Method*

### 2.4.1. Groundwater Storage Anomalies (GWSA)

The Terrestrial Water Storage Anomaly (TWSA) consists of four parts, namely, groundwater storage anomalies (GWSA), soil storage anomalies (SMSA), snow water equivalent anomalies (SWEA), and canopy water storage anomalies (CWSA) [26]. The GWSA could be isolated using the following formula [26]:

$$\text{GWSA} = \text{TWSA} - \text{(SMSA} + \text{SWEA} + \text{CWSA)} \tag{1}$$

where TWSAs were derived from gravity anomalies observed by GRACE/FO and the other three components were provided by GLDAS.

It should be emphasized that all data are also relative to the 2004–2009 mean baseline, meaning the value of the GWSA can be positive, negative, or zero. If the GWSA calculated for a certain pixel in a certain month are positive, this indicates that the GWSA in this area increase compared with the 2004–2009 mean baseline; a negative value indicates a decrease, and zero indicates no change.

### 2.4.2. Trend Test and Significance Analysis

Temporally, unitary linear regression analysis was used to construct a linear regression equation, and the *k* was used to quantify the trends as follows [52,56]:

$$k = \frac{n \times \sum\_{i=1}^{n} \left(i \times x\_i\right) - \sum\_{i=1}^{n} i \times \sum\_{i=1}^{n} x\_i}{n \times \sum\_{i=1}^{n} i^2 - \left(\sum\_{i=1}^{n} i\right)^2} \tag{2}$$

where *n* is the number of months and *xi* is variable in month *i*. If *k* > 0, the variables show a positive trend; otherwise, they would show a negative trend. *k* = 0 means no change.

Spatially, the change trend for various factors in each pixel was analyzed using the modified Mann–Kendall (MMK) test and represented by *Slope* [52,56–59]. Additionally, whether this trend has statistical significance was determined by *t*-tests. The corresponding formulas were as follows:

$$Slope = \text{median}\left(\frac{\mathbf{x}\_j - \mathbf{x}\_k}{j - k}\right) \quad (1 \le k < j \le n) \tag{3}$$

$$\text{sgn}(\mathbf{x}\_{j} - \mathbf{x}\_{k}) = \begin{cases} 1 & \mathbf{x}\_{j} > \mathbf{x}\_{k} \\ 0 & \mathbf{x}\_{j} = \mathbf{x}\_{k} \\ -1 & \mathbf{x}\_{j} < \mathbf{x}\_{k} \end{cases} \tag{4}$$

$$S = \sum\_{k=1}^{n-1} \sum\_{j=k+1}^{n} \text{sgn}(\mathbf{x}\_j - \mathbf{x}\_k) \tag{5}$$

$$Var(S) = \frac{n(n-1)(2n+5)}{18} \tag{6}$$

$$\eta = 1 + \frac{2}{n(n-1)(n-2)} \times \sum\_{k=1}^{n-1} (n-k)(n-k-1)(n-k)\rho\_{\rm ACF}(k) \tag{7}$$

$$Z = \begin{cases} \begin{array}{c} \frac{S-1}{\sqrt{\eta \times Var(S)}} \quad \text{for } S > 0\\ 0 \qquad \text{for } S = 0\\ \frac{S+1}{\sqrt{\eta \times Var(S)}} \quad \text{for } S < 0 \end{array} \tag{8}$$

Here, *xj* is the parameter value at time *j*, and *xk* at time *k*. *Slope* > 0 represents an upward trend; *Slope* < 0 represents a downward trend; and *Slope* = 0 means no change. *ρACF*(*k*) is the autocorrelation function (ACF) of the ranks of the observations [58]. The *Z* value indicates the trend of factors. In this study, the given significance levels were 99%, 95%, and 90%, when *p* = 0.01, 0.05, and 0.1, and *|Z|* − *p/2* was equal to 2.58, 1.96, and 1.64, respectively [52,56]. Based on the factors trend, the nine levels were further divided and are presented in Table 1.

**Table 1.** Index classification.


### 2.4.3. Correlation Analysis

The correlation analysis method was used to analyze the relationship between two variables in each pixel [52]. The corresponding formula was as follows:

$$\text{CC}\_{xy} = \frac{\sum\_{i=1}^{n} \left[ (\mathbf{x}\_i - \overline{\mathbf{x}})(y\_i - \overline{\mathbf{y}}) \right]}{\sqrt{\sum\_{i=1}^{n} \left( \mathbf{x}\_i - \overline{\mathbf{x}} \right)^2 \sum\_{i=1}^{n} \left( y\_i - \overline{\mathbf{y}} \right)^2}} \tag{9}$$

where *CCxy* represents the correlation coefficient (*CC*) between variables *x* and *y*; *xi* and *yi* represent the value of the variables *x* and *y* in month *I*, respectively; and *x* and *y* are the mean value of the variables *x* and *y*, respectively. The *CC* ranged from −1 to 1. There was a positive correlation when the *CC* was greater than 0 and a negative correlation when it was less than 0. Additionally, the significance of the results was evaluated using a *t*-test, and it could be divided into four groups, extremely significant, significant, weakly significant, and insignificant. Thus, the relationships between the dependent variable and independent variable were further divided into nine levels, as detailed in Table 1.

### **3. Results and Discussion**

### *3.1. Spatial–Temporal Patterns of GWSA*

The GWSA across the TP had obvious spatial heterogeneity, with a mean of 0.17 cm (Figure 2a) and a decreasing rate of −0.25 mm/a (Figure 2b). The higher values were mainly distributed in the QB, the IB, the source of the Yangtze River, and the southeastern edge of the TP. The difference in the GWSA reached 9.26 mm. The value of the GWSA in the QB was the highest, 3.42 cm, ranking first in the sub-regions. However, it had been decreasing with the highest decreasing rate of −1.24 mm/a (Figure 2b) and in a large area of 79.3% of the region (Figure 2c). On the contrary, the YZ-BRB had the lowest GWSA of −5.84 cm, with the highest increasing rate at 2.12 mm/a and a large increasing area of 80.2% of the basin. Furthermore, it revealed an increasing trend in 42.0% of the plateau area, while a decreasing trend was present in 58.0% of the plateau area (Figure 2c,d). The areas where the slope of the GWSA was positive were mainly HC, S-IRB, YZ-BRB, and south and southwest of the IB. In the HC and S-IRB, 53.9% and 33.3% of the regions showed a significant increasing trend, and a significant decrease was observed in <2% of these two basins, which indicated that the GWS levels in these sub-regions were increasing. In the YB and the YRB, 83.0% and 63.7% showed a decreasing trend, while 31.2% and 17.0% showed a significant decreasing trend, respectively.

The GWSA varied seasonally, with an increase in rainy seasons (May–October) and a decrease in dry seasons (November to April in the next year) (Figure 2e). The difference in the GWSA in the YB reached 1.06 cm, where the GWSA in rainy seasons and dry seasons were 1.15 cm and 0.09 cm, respectively. The variation in the GWSA in the YB was consistent with the variation in rainfall, indicating that the annual variation in groundwater may be affected by the recharge of rainfall. The GWSA of the IB in rainy seasons and dry seasons were 3.67 cm and 3.09 cm, respectively. As the amount of precipitation in this basin is still less in the rainy season, it is possible that permafrost degradation and glacier melting caused by increasing temperature have a stronger impact on the GWS than precipitation. The minimum negative anomaly of the interannual variation in the GWSA was −2.16 cm in 2009, while the maximum positive anomaly of 2.11 cm occurred in 2012 (Figure 2f). Although the GWS increased rapidly from 2009 to 2012, it was still declining overall from 2003 to 2021, which was similar to the result inferred from the monthly GWSA from April 2002 to December 2021 (Figure 2g). The latter was divided into five periods, April 2002 to June 2005, July 2005 to June 2010, June 2010 to May 2013, June 2013 to December 2015, and January 2016 to December 2021. Overall, the GWS in the TP was decreasing at an average rate of −0.89 mm/a from January 2003 to December 2021. Despite the lowest GWSA value of −2.66 cm occurring in October 2004, the GWS was still rising at a speed of 15.23 mm/a from April 2002 to June 2005. Then, there was a 60-month decrease at −7.31 mm/a from July 2005 to June 2010. The third period, from June 2010 to May 2013, was an increasing period. The increasing speed was 11.29 mm/a, and the GWSA reached the maximum value of 5.82 cm in May 2013. The largest decrease rate of −21.29 mm/a occurred from June 2013 to December 2015. It had been slowly rising since January 2016, with a speed of 1.47 mm/a. The results of the study by Liu et al. [22] were the same. Liu et al. [22] also came to a similar conclusion that GWS in the Qinghai province and the Tibet Autonomous Region reached a split point in 2016, and the groundwater storage changed from decreasing with a mean rate of −0.2 mm/a from 2003 to 2015 to increasing with a mean rate of 3.28 mm/a from 2016 to 2019. The GWSA in the TP showed an overall downward trend, but increases after 2016. This shows that the amount of groundwater storage is gradually recovering.

**Figure 2.** Characteristics of GWS in the TP and ten sub-regions. (**a**) Spatial distribution of monthly mean GWSA from January 2003 to December 2021; (**b**) spatial distribution of slope of annual mean GWSA from 2003 to 2021; (**c**) spatial distribution of significant changes in annual mean GWSA from 2003 to 2021; (**d**) significance analysis of changes in annual mean GWSA from 2003 to 2021; (**e**) variation in monthly mean GWSA; (**f**) variation in annual mean GWSA; (**g**) temporal variations in monthly GWSA from April 2002 to December 2021.

### *3.2. The Influence of Climate Change on GWS in the TP*

### 3.2.1. Fluctuation Characteristics of Regional Climate

The changes in temperature and precipitation were calculated, and the spatial characteristics of climate change are shown in Figure 3.

**Figure 3.** Spatial characteristics of climate change from 2003 to 2021. (**a**) Spatial distribution of slope of annual mean temperature; (**b**) spatial distribution of significant changes in annual mean temperature; (**c**) spatial distribution of slope of annual precipitation; (**d**) spatial distribution of significant changes in annual precipitation.

As shown in Figure 3a,b, the changes indicate that the temperature increased in 64.1% of the TP and decreased in 35.1% of the TP. The areas where the temperature increased were mainly distributed in the YRB, the L-MRB, the N-SRB, the S-IRB, the east of the YB, the east of the YZ-BRB, and the south of the IB. The increase rate of the temperature in the YRB was 0.016 ◦C/a. Additionally, 89.0% of the YRB showed an increasing trend. The decrease rate of the temperature in the QB, in which 88.3% showed a decreasing trend, was −0.015 ◦C/a.

The precipitation showed a decreasing trend at −0.04 mm/a and in 49.6% of the TP from 2003 to 2021, mainly over the YZ-BRB, the N-SRB, the TB, the south of the YB, the south of the L-MRB, the center of the IB, and the northwest of the QB (Figure 3c,d). The Brahmaputra River Basin had the highest decreasing rate of −3.55 mm/a and a decreasing area of 90.5%. Meanwhile, the largest increase in precipitation was shown to be in the Indus River Basin, with an increasing rate of 4.57 mm/a and an increasing area of 92.4%. These results indicate that there have been different or even opposite trends in the rate of climate change in different sub-regions of the TP in recent years.

### 3.2.2. Relationship between Climate Change and GWS

The correlation coefficients between the GWSA and Tmp and Pre were calculated to analyze the relationships between the GWS and climate change. Then, significance tests were carried out. The results are shown in Figure 4.

Overall, the GWSA were positively correlated with Tmp, and the mean correlation coefficient was 0.0492 (Figure 4a). The correlation coefficient in the YB was 0.143, which was the highest among the ten sub-regions. The HC, the QB, the YRB, the YB, the IB, and the S-IRB were the main regions that showed positive correlations. This shows that the areas that had positive correlations accounted for about 74.94% of the study area. In these areas, 15.23% showed extremely significant positive correlations. There were mostly insignificant

positive correlations in 47.51% of the study area, which were mainly the HC, the QB, and the YRB (Figure 4b). The negatively correlated area covered 25.06% of the study area, with insignificant negative correlation occupying 13.66% of the area (Figure 4a,b). In the L-MRB, the N-SRB, the YZ-BRB, and the TB, the GWSA were insignificantly negatively correlated with the temperature in 25.87% of the area.

**Figure 4.** Relationship between climate change and GWS. (**a**) The correlation coefficients between GWSA and temperature; (**b**) the significance test between GWSA and temperature; (**c**) the correlation coefficients between GWSA and precipitation; (**d**) the significance test between GWSA and precipitation.

As shown in Figure 4c, the mean correlation coefficient was 0.054, which was greater than 0, indicating a positive correlation between the GWSA and Pre in the TP. The variables were positively correlated in 70.81% of the study area and negatively correlated in 29.19% of the study area. The mean correlation coefficient was 0.148 in the YB, which was the highest among the ten sub-regions. The L-MRB, the N-SRB, the YZ-BRB, the IB, and the TB were the main regions showing negative correlations, accounting for 68.24% of the negatively correlated areas in total. The results of the significance tests in Figure 4d indicate that the GWSA were extremely significantly positively correlated with precipitation in 13.96% of the TP, mainly in the YB and south of the IB.

Increasing precipitation will lead to more groundwater recharge, and thus, an increase in the GWS. Precipitation in the TP has been slightly increasing since the 1960s and is projected to further increase in the future in most areas of the TP [60]. The direct manifestation of permafrost degradation caused by increased temperature is the increase in the number and area of thermal karst lakes, which leads to the reduction in GWS in the form of runoff or evaporation [61,62]. However, when a large number of glaciers are distributed in the basin, the increase in glacial meltwater caused by the increase in temperature will recharge the groundwater, weaken the decrease in GWS, and even cause the increase in GWS. This is the reason why the GWS in the north and south of the IB showed different trends. There are many glaciers in the southern part of the IB and almost no glaciers in the northern part of the IB. Similar results were obtained in previous studies in the central Qiangtang Nature Reserve and the Upper Indus Basin [16]. The GWS in the L-MRB, N-SRB, and YZ-BRB increased in the context of decreased precipitation and increased temperature. As decreased precipitation and increased glacial melt due to temperature increases may cause decreased GWS, these three regions may require more groundwater recharge. The

decreased GWS in the YRB and YB, where the temperature and precipitation increased, may indirectly support this hypothesis that the groundwater in the YRB and YB may seep to the L-MRB, N-SRB, and YZ-BRB [16].

### *3.3. Ecological Effect of Groundwater Storage*

### 3.3.1. Vegetation Change

As shown in Figure 5, grassland has the largest area of 84.285 × <sup>10</sup><sup>4</sup> km<sup>2</sup> of the TP (Table 2), mostly distributed in the IB, YB, YRB, and YZ-BRB. The area of unused land has the second largest area, which is 59.736 × <sup>10</sup><sup>4</sup> km2, mainly covering the IB and QB. The forest land that is mainly distributed in the YB and YZ-BRB covers 23.457 × <sup>10</sup><sup>4</sup> km2. The area of water, mainly lakes and glaciers, is 9.577 × 104 km2. The cropland and urban land types cover 1.994 × <sup>10</sup><sup>4</sup> km2 and 0.210 × 104 km2, respectively, which means there is a low level of human activity on the TP. However, these areas have increased by 0.19 × <sup>10</sup><sup>4</sup> km<sup>2</sup> and 0.10 × <sup>10</sup><sup>4</sup> km2 from 2000 to 2020, respectively, indicating that human activities cannot be ignored.

The highest decrease and increase in land cover changes from 2000 to 2020 were grassland and unused land, with changed areas of −17.19 × 104 km2 and 12.18 × <sup>10</sup><sup>4</sup> km2, respectively. This means that large areas of grassland have degraded to sand, desert, and other kinds of unused land, especially in the IB, where the area of grassland decreased to 14.926 × <sup>10</sup><sup>4</sup> km2 and the unused land increased to 12.284 × 104 km2. It should be noted that different trends occurred in the LUCC of the YRB and YB. The area of unused land decreased while the area of forest and grassland increased, which proved the effectiveness of ecological protection and construction projects, such as eco-migration, grazing bans, forest and wetland reservations, etc. [16,63].

As LUCC is discrete data, it is difficult to quantify the vegetation response to GWS. Therefore, the influences of groundwater change on vegetation in time and space were discussed using the NDVI, an indicator of vegetation change.

**Figure 5.** Land use map of the TP in 2020.


**Table 2.** Area of different land use types in 2020 and changed area from 2000 to 2020 (104 km2).

### 3.3.2. Vegetation Responses to GWS Changes

Overall, the NDVI increased with the mean rate of 0.0035 per ten years, but the increase was lower in the west and faster in the east from 2003 to 2019 (Figure 6a). The YB showed the fastest increasing rate of 0.0144 per ten years, while the lowest was in the QB, with a rate of −0.0071 of per ten years. The areas where the NDVI increased significantly are mainly distributed over the east edge of the TP (Figure 6b). In Figure 6c, it can be seen that the GWS and the NDVI were positively correlated, and the correlation coefficient was 0.068. Additionally, they had a positive correlation relationship in 73.17% of the study area. The highest correlation coefficient was 0.113 in the YB, while the lowest was −0.016 in the L-MRB. The weakly to extremely significant positive correlations were in 36.91% of the aforementioned areas, mostly in the QB, YB, and YZ-BRB (Figure 6d). The regions where there was a negative correlation between the GWS and the NDVI covered 26.83% of the TP, of which only 23.59% were significantly or extremely significantly negatively correlated (Figure 6d).

The GWS was positively correlated with the NDVI, which indicated that the vegetation was dependent on the groundwater condition. In the IB, large areas of grassland have degenerated into unused land because of the decline in the groundwater level in the context of climate change [64,65]. This is analogous to the earlier findings reported in the TP. Peng et al. [66] also found that the decline of the groundwater level in the source area of the Yellow River was closely related to the deterioration of the ecological environment in the permafrost-degraded area. This may have been due to the decrease in the groundwater level, which would have caused a great loss to occur in the shallow soil moisture, which

would have meant that the root system of some vegetation in cold areas could not effectively use soil water, eventually leading to the degradation or disappearance of vegetation [67–69].

**Figure 6.** Relationship between NDVI and GWS. (**a**) Spatial distribution of slope of annual mean NDVI; (**b**) spatial distribution of significant changes in NDVI; (**c**) the correlation coefficients between GWS and NDVI; (**d**) the significance test between GWS and NDVI.

### **4. Conclusions**

Based on the linear regression and modified Mann–Kendall test analysis methods, the temporal and spatial variations in GWS were analyzed using the GRACE and GLDAS data from April 2002 to December 2021. Then, the driving factors and ecological effects of the changes in GWS in the Tibetan Plateau were discussed, using the correlation analysis model and multi-source remote sensing data. This study of groundwater storage changes provides a useful reference for ecological protection and the high-quality development of the TP. The conclusions are as follows:


and eventually led to the degradation of 17.19 × <sup>10</sup><sup>4</sup> km2 of grassland to sand, desert, or other kinds of unused land on the TP.

The results of this study could offer a new opportunity to reveal the groundwater changes in a cryosphere region and to assess the impact of the changes in hydrological conditions on ecology. However, this study only provided a preliminary investigation of the changes to GWS, and more in-depth studies should focus on the following specific aspects: (1) in situ groundwater monitoring with high temporal and spatial resolution at the regional scale to gain insight into the accurate variation mechanism of groundwater storage on the TP; (2) the integration of multivariate satellite data and improving modeling capacity to deal with the quality change signals caused by rapid structural uplift; (3) the quantitative analysis of the response of groundwater to climate change and the evaluation of the effects of groundwater changes on vegetation evolution under climate change.

**Author Contributions:** W.R., conceptualization, methodology, formal analysis, and writing—original draft; Y.G., conceptualization and supervision; H.Q., conceptualization, supervision, resources, and writing—review and editing; Y.M., Z.S. and W.M., supervision and resources; Y.L. and P.X., methodology. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (Grant No. 2019QZKK0103), the National Natural Science Foundation of China (41931285, 42007184), the Fundamental Research Funds for the Central Universities, CHD (Grant No. 300102292901), and the Program of Introducing Talents of Discipline to Universities (Grant No. B08039).

**Data Availability Statement:** The CSR GRACE/FO RL06 mascon dataset is provided by the Center for Space Research at the University of Texas, Austin and downloaded at http://www2.csr.utexas. edu/grace (accessed on 27 March 2022). The GLDAS Noah model is developed by the National Aeronautics and Space Administration and the National Oceanic and Atmospheric Administration and downloaded at https://disc.gsfc.nasa.gov/datasets (accessed on 27 March 2022). The CRU TS data is produced by the UK's National Centre for Atmospheric Science (NCAS) and downloaded at https://crudata.uea.ac.uk/ (accessed on 21 May 2022). The NDVI and LUCC data are both obtained from the Resources and Environment Science and Data Center, Chinese Academy of Sciences at http://www.resdc.cn/ (accessed on 1 May 2022).

**Acknowledgments:** The authors gratefully acknowledge the supports of various foundations. The authors are grateful to the editor and anonymous reviewers whose very insightful comments have contributed to improving the quality of this paper.

**Conflicts of Interest:** The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
