Analysis of Groundwater Storage Changes and Influencing Factors in China Based on GRACE Data

Abstract

Groundwater is a primary freshwater resource for human consumption and an essential source for industry and agriculture. Therefore, understanding its spatial and temporal trends and drivers is crucial for governments to take appropriate measures to manage water resources. This paper uses Gravity Recovery and Climate Experiment (GRACE) satellite data and the Global Land Data Assimilation System (GLDAS) to derive groundwater storage anomalies (GWSAs) and to analyze the spatial and temporal trends of GWSA in different regions of China (Xinjiang, Tibet, Inner Mongolia, North China Plain, South China, and Northeast China). It used groundwater-level observation data to verify the accuracy of GWSA estimates and analyzed the drivers of regional GWSA changes. The results showed that: (1) GWSA in South China increased at a rate of 4.79 mm/a from 2003 to 2016, and GWSA in other regions in China showed a decreasing trend. Among them, the decline rates of GWSA in Xinjiang, Tibet, Inner Mongolia, North China Plain, and Northeast China were −6.24 mm/a, −3.33 mm/a, −3.17 mm/a, −7.35 mm/a, and −0.75 mm/a, respectively. (2) The accuracy of the annual-scale GWSA estimates was improved after deducting gravity losses due to raw coal quality, and the correlation coefficient between GWSA and groundwater levels monitored by observation wells increased. (3) In Xinjiang, the annual water consumed by raw coal mining, industrial, and agricultural activities had a greater impact on GWSA than rainfall and temperature, so these human activities might be the main drivers of the continued GWSA decline in Xinjiang. Water consumption by raw coal mining and industry might be the main drivers of the continued decline in GWSA in Inner Mongolia and the North China Plain. The increase in groundwater storage in South China was mainly due to the recharge of rainfall.

1. Introduction

Groundwater is an essential part of the land water cycle and is one of the critical water sources for agricultural irrigation, industrial water consumption, and urban water consumption [1]. In recent years, many countries and regions have faced a severe lack of groundwater resources due to climate change, population growth, and over-exploitation [2,3]. As a country with a large population, China’s water resource demand has increased in the last several decades due to the rapid development of urbanization, agriculture, and industrialization. At the same time, as a country rich in coal resources, the impact of coal-mining activities on water resources cannot be ignored [4]. Coal mining can cause a loss of gravity by removing coal from Earth’s crust. Meanwhile, much groundwater must be artificially discharged from underground mines during the coal-mining process. Both of these aspects can cause a great gravity loss in a large area. Tang et al. [5] found that coal transport caused mass loss in the west of China at a rate of −3.13 mm yr⁻¹ from 2003 to 2011. The groundwater in many regions in China is extracted to meet people’s demand for water resources. Improper exploitation of water resources has caused many serious problems, such as land subsidence, seawater intrusion, and soil salinization, which have greatly affected the development process of society [6,7]. Therefore, it is vital to investigate the change in groundwater storage and the influencing factors of such changes.

Traditional groundwater monitoring methods have relied on groundwater-level observation wells. However, this method is costly and unevenly distributed. Gravity Recovery and Climate Experiment (GRACE) gravity satellite observations have the advantage of uniform scale and even spatial distribution, solving the problems of uneven spatial distribution and missing data of hydrological and meteorological stations. This technique allows one to measure Earth’s gravity field to obtain changes in terrestrial water storage reserves. Therefore, after subtracting the change in surface water storage, the change in groundwater storage can be deduced, which provides a new way to study the change in regional water resources [8]. Many scholars have conducted extensive studies using GRACE data and hydrological model data and have mainly focused on groundwater storage change monitoring [9,10,11,12], while some have carried out technical method improvements, accuracy testing, and spatial resolution downscaling [13,14,15]. Several studies have also shown that anthropogenic and natural factors play an important role in groundwater storage changes [16,17,18]. Guo et al. [19] found that anthropogenic factors had a greater impact on groundwater depletion in the Haihe basin than natural factors, with continued over-exploitation of groundwater, agricultural activities and irrigation, and urban expansion being the main factors in groundwater resource shortages. Li et al. [20] found that groundwater extraction had a greater impact on changes in groundwater storage than precipitation and evapotranspiration in the North China Plain. Zhang et al. [21] calculated the change in groundwater storage in China from 2003 to 2016 and concluded that the change in groundwater storage in China showed a clear interannual trend and seasonal variation, and the seasonal fluctuation in groundwater storage was consistent with precipitation. Zhou Miao et al. [22] studied groundwater storage in several typical regions of China using GRACE and meteorological data, which concluded that groundwater storage in most of Southeast China is increasing year by year, with the main recharge coming from precipitation. Despite previously studied groundwater storage changes in China, these studies did not consider gravity losses due to coal mining when estimating changes in terrestrial water storage reserves, resulting in poor accuracy of groundwater storage change estimates. At the same time, no comprehensive studies of human activities and climate change effects on groundwater storage change at longer time scales and larger spatial scales have been conducted.

Therefore, the main contents of this paper are: (1) To study the spatial and temporal variations in groundwater storage anomalies (GWSAs) in different regions of China (Xinjiang, Tibet, Inner Mongolia, North China Plain, South China, and Northeast China) based on the data from GRACE and Global Land Data Assimilation System (GLDAS), from 2003 to 2016. (2) Groundwater-level monitoring data were used to verify the accuracy of the GWSA estimates. At the same time, the gravity loss caused by the quality of raw coal to the land water storage was deducted when estimating the annual-scale GWSA variation, which further improved the accuracy of the GWSA estimation. (3) We analyzed the effects of different natural factors (precipitation and air temperature) and human activities (coal mining, industry, and agriculture) on average annual changes in groundwater storage.

2. Materials and Methods

2.1. Study Area

China has a wide geographical area with diverse climates and topography. Generally, Northeast China has a monsoon climate, and Northwest China has a continental climate [23]. In contrast, Southern China has a subtropical monsoon climate [24]. Additionally, precipitation’s spatial and temporal distribution is uneven, with more precipitation in summer and autumn and less precipitation in winter and spring. Further, the spatial distribution of precipitation gradually decreases from the southeast coast to the northwest inland [25]. For example, July and August receive 60–70% of the annual precipitation in North China and Northeast China. In the coastal areas of South China, there is a heavy rain period in April and May and a relatively dry period in June and July. In the drought regions of Northwest China, precipitation is scarce throughout the year [26]. Due to the influence of complex factors, the variation in groundwater storage in different regions is different. Due to sufficient precipitation, Southern China generally has more groundwater recharge than water consumption. Therefore, groundwater storage has shown a positive trend. However, in Northern China, groundwater consumption by industrial and agricultural activities is much greater than the amount of groundwater recharge, causing a groundwater deficit [27]. Therefore, GWSA plays a crucial role in water resource management in these areas. In order to analyze the regional features of GWSA change, different subareas of China were selected as research areas: Xinjiang, Tibet, Inner Mongolia, North China Plain, South China, and Northeast China (Figure 1).

2.2. Data

2.2.1. GRACE Data

The GRACE project launched two satellites about 220 km apart. As the mass migrates on Earth’s surface, it causes a change in gravity. Satellites are subjected to different gravity, and their state of motion will change, resulting in a slight change in the distance between the two satellites. By measuring this slight variation, we can infer the gravity in different regions [28,29,30]. In this paper, we use the terrestrial water storage anomalies (TWSAs) provided by the newly released GRACE dataset (the GRACE CSR RL06 Mascon V2) (http://www2.csr.utexas.edu/grace, accessed on 15 February 2022) [31] from the University of Texas Center for Space Research (CSR). Compared to the traditional spherical harmonic coefficient (SH) method, the Mascon method reduces leakage errors and eliminates the north–south stripe error, with a higher spatial resolution and more effective signals [32]. The TWSA dataset was relative to the baseline average from January 2004 to December 2009 and expressed in the equivalent height of water (EWH). The time was January 2003 to December 2016, and the spatial and temporal resolution of the CSR mascon solutions was 0.25° and one month, respectively. For individual missing months during the study period, the average value of adjacent months is used for interpolation [33,34].

2.2.2. GLDAS Data

GLDAS is a collaborative project developed by NASA and the National Oceanic and Atmospheric Administration (NOAA) [35]. GLDAS uses advanced data assimilation techniques to integrate satellite and ground observations into a unified model, including four land surface models. In this study, the output of the Noah model was selected, which contains all kinds of surface information in the grid-data form, such as snow water equivalent, canopy water storage, soil moisture, etc. Due to the high-resolution simulation performed by GLDAS-Noah, the output has been used in many studies addressing water storage changes [36,37]. In this study, we use the output of Noah 2.7.1 (https://podaac-tools.jpl.nasa.gov/, accessed on 20 February 2022) with a spatial resolution of 1° × 1°, resampled to 0.25° × 0.25°, and a temporal resolution of monthly scale, with a time length of January 2003 to December 2016.

2.2.3. Monitoring Groundwater Levels in Wells

Groundwater level (GWL) monitoring data from observation wells in different areas of China were compiled for the period of 2005~2013 in the “Chinese Geological, and Environmental Monitoring Groundwater Level yearbooks” [38] and were published by China Land Press. These data were used to verify the accuracy of the GWSA data. Due to the lack of data and uneven spatial distribution in the original statistics of GWL, the annual average GWL can be obtained by screening some wells with relatively complete data and supplementing the few missing data with a linear interpolation of values from adjacent months. Finally, we selected data from 123, 22, and 73 wells in the North China Plain, South China, and Northeast China, respectively, to verify the accuracy of the estimated GWSA (the method of estimation is presented in the following section) in each of the three regions. The locations of observation wells are shown in Figure 1.

2.2.4. Meteorological Data: Rainfall/Temperature

Precipitation and temperature are the meteorological factors affecting the variability in GWSA. The changes in air temperature can give rise to variations in evapotranspiration from the ground, which can cause a loss of gravity/groundwater. Precipitation data estimated from TRMM3B43 satellite remote sensing were obtained. The time was January 2003 to December 2016, and rainfall data’s spatial and temporal resolution was 0.25° and one month, respectively. Air temperature data were obtained from the National Earth System Science Data Center. The time is from January 2003 to December 2016, and the temperature data’s spatial and temporal resolution was 1 km and one month, respectively. Areal precipitation and temperature data in each region can be obtained by averaging the values at all regional grid boxes. For consistency with GRACE data, meteorological data were selected relative to the baseline average from January 2004 to December 2009.

2.2.5. Statistics of Human Water Consumption

Human water consumption mainly includes industrial, agricultural, and coal-mining water, etc. Except for water used for coal mining, the data for other types of water use can be found in the China Statistical Yearbook [39]. Coal-mining water consumption can be obtained from the annual production of raw coal in the China Statistical Yearbook, and the conversion method is as follows [40,41].

$$\begin{array}{c}{M}_w=\mu M_c\end{array}$$

(1)

where M_c denotes the annual production of raw coal (t/year). μ refers to the water consumption coefficient per ton of coal (m³/t) in which μ = 0.87 [42], and M_w is the water consumption associated with mine drainage (m³/year).

2.3. Methods

2.3.1. Calculation of Groundwater Storage Anomalies

The terrestrial water storage anomalies (TWSAs) mainly include groundwater storage anomalies (GWSAs), soil moisture anomalies (SMSAs), surface water anomalies (SWAs), snow water equivalent anomalies (SWEAs), canopy water storage anomalies (CWSAs), and biomass anomalies (BMAs). In the long-term aspect, the impacts of surface runoff and biomass changes on terrestrial water storage are negligible [43]. For regions rich of coal resources, such as Xinjiang, Inner Mongolia, Northeast China, and North China Plain, the annual production of raw coal can reach hundreds of millions of tons each year. Since the TWSA includes the gravity loss caused by raw coal (GLCRA), the gravity of raw coal needs to be deducted from the TWSA, and the actual TWSA changes in the study area are as follows:

$$\begin{array}{c}TWS{A}^{\prime }=TWSA-\Delta coal\end{array}$$

(2)

where Δcoal refers to the GLCRA; only the annual mass of raw coal production (AMRCP) can be obtained. Thus, only the annual GLCRA was deducted from the annual TWSA series. The mass of mined raw coal needs to be converted to water height. The specific conversion method is as follows:

$$\begin{array}{c}{h}_m=\frac{m_i}{s_i\ast \rho }\end{array}$$

(3)

where h_m is the equivalent height of water converted from the AMRCP, m_i is the coal-mining output of each region, S_i is the area of each region, and ρ is the density of water 1000 kg/m³.

Therefore, using datasets of GRACE and GLDAS, GWSA can be separated from TWSA according to Equation (4):

$$\begin{array}{c}GWSA=TWS{A}^{\prime }-SMSA-CWSA-SWEA\end{array}$$

(4)

where SMSA, CWSA, and SWEA are provided in the GLDAS dataset, which was produced by the Noah2.7.1 Land hydrology model. TWSA is provided by GRACE CSR RL06.

First, this paper obtains the TWSA from the GRACE data. Since the TWSA includes the gravity losses due to coal mining, the GLCRA must be deducted from the TWSA before calculating the change in GWSA. The raw coal mass is converted to an equivalent water height according to Equation (3), which obtains the GLCRA. The GLCRA is subtracted from Equation (2) to obtain a more accurate TWSA. Finally, TWSA subtracts SMSA, CWSA, and SWEA from GLDAS to obtain GWSA, according to Equation (4).

2.3.2. Correlation Analysis Method

In this paper, Pearson correlation analysis (R) was used to study the relationship between GWSA and GWL and that between the GWSA and the AMRCP, respectively, as well as that between GWSA and natural factors (rainfall, temperature), and human factors (the water consumption of coal mining, industry, and agriculture) were calculated. For any pair of samples X_i (i = 1, 2, …, n) and Y_i (i = 1, 2, …, n), the Pearson correlation coefficient R can be calculated by the following Equation (5):

$$\begin{array}{c}R=\frac{{{\displaystyle \sum }}_{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(X_i-\overline{X}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}}\end{array}$$

(5)

where: n is the number of samples, X_i and Y_i are the observed values of two time series, and $\overline{X}$ and Y are the average values of X_i and Y_i, respectively.

2.3.3. Mann–Kendall Trend Test

The Mann–Kendall [44,45,46] test is non-parametric and is also known as a no-distribution test. Its advantage is that the sample is not required to follow a certain distribution and is not disturbed by a few outliers. This study tested the significance of the long time series of GWSA. For any data to be tested X_i, the statistic is defined according to the following equation:

$$\begin{array}{c}S={\displaystyle {\displaystyle \sum }_{i=1}^{n-1}}{\displaystyle {\displaystyle \sum }_{j=i+1}^n}Sgn\left(X_j-X_i\right)\end{array}$$

(6)

$$\begin{array}{c}Sgn=\{\begin{array}{ll}1& (X_j-X_i)>0\\ 0& (X_j-X_i)=0\\ -1& (X_j-X_i)<0\end{array}\end{array}$$

(7)

$$\begin{array}{c}Var=\frac{n\left(n-1\right)\left(2n+5\right)}{18}\end{array}$$

(8)

$$Z=\{\begin{array}{ll}\frac{S-1}{\sqrt{Var}}& S>0\\ 0& S=0\\ \frac{S+1}{\sqrt{Var}}& S<0\end{array}$$

(9)

where X_j and X_i are the observation data at different times, i and j refer to the time step, n represents the length of the time series, and Sgn is the sign function. In determining the level of significance of α, if statistic |Z| > U_1−α/2, this indicates that there is a significant trend of change in the sequence; otherwise, there is no significant trend of change. In general, the value of α is 0.01 or 0.05, which indicates whether the trend of change in the time series of the elements is significant, at a confidence level of 0.01 or 0.05. p < = 0.05 is “significant” and p < = 0.01 is “very significant”.

3. Results

3.1. Trends of Groundwater Storage Changes

In this paper, six typical regions in China were studied: Xinjiang (73°40′~96°23′ E, 34°22′~49°10′ N), Tibet (26°50′~36°53′ E, 78°25′~99°06′ N), Inner Mongolia (40°20′~50°50′ E, 106°~121°40′ N), North China Plain (110°23′~122°71′ E, 29°68′~42°67′ N), South China region (104°29′~117°19′ E, 25°31′~29°05′ N), and Northeast China (118°53′~135°05′ E, 38°43′~53°33′ N). The GWSA for each region is obtained by calculating the average over the grid points in each region. The results were expressed in the form of equivalent water height. Figure 2a–f show the trend of groundwater storage in the study area. During 2003–2016, only GWSA in South China shows a significant trend of increase (p < 0.01), with mean rates of 4.79 mm/a. The other five regions showed a significant trend of decrease (p < 0.01), among which the GWSA decreased relatively fast in the North China Plain and Xinjiang region, approximately to −7.35 mm/a, −7.35 mm/a and −6.24 mm/a, respectively, and the GWSA decreased more slowly in Northeast China, about −0.75 mm/a. GWSA decreased at −3.33 mm/a and −3.17 mm/a in Tibet and Inner Mongolia, respectively.

3.2. Spatial Variation in Trends of GRACE-Based Groundwater Storage Change

The mean spatial variations in GRACE-based GWSA in study regions from 2003 to 2016 were calculated and analyzed. The spatial trends of GRACE-based GWSA in each region are shown in Figure 3a–f. The fastest declining region of groundwater storage in Xinjiang was mainly in Aksu and Yili, with a maximum rate of decline of about −50.11 mm/a. This might be due to low annual rainfall, limited groundwater recharge, and the extraction of large amounts of groundwater (Figure 3a). In Southeastern Tibet, the decreasing trend of GWSA was obvious, while it showed an increasing trend in the northeast, and Inner Mongolia also showed the same spatial variation (Figure 3b,c). The fastest declining areas of GWSA in Tibet were mainly concentrated in the Linzhi region, with a maximum decline rate of about −56.24 mm/a. The study in reference [47] mentioned that water consumption in the Linzhi region was extensive, and the serious depletion was mainly used for agricultural irrigation, accounting for 86.68% of the total water consumption. Therefore, anthropogenic factors might be the main reason for the decrease in groundwater in the region.

The fastest decline in GWSA in the North China Plain mainly occurred in the central area of the region, including the north of Henan Province, the southeast of Shanxi Province, and the southwest of Hebei Province, with a maximum decline rate of about −31.46 mm/a (Figure 3d). The variation trend of groundwater drought in the North China Plain in reference [48] was also similar to spatial trends of GWSA in this paper. The severe decline in GWSA might be due to the over-exploitation of groundwater. Shanxi and Henan are rich in coal-mining resources, and annual coal mining also consumes a large amount of water resources. In addition, with the rapid development of urbanization, the consumption of groundwater by industry is also likely to increase gradually. The decreasing trend of GWSA from the center of the North China Plain to the northern part of Beijing–Tianjin–Hebei has slowed down, which might be related to the South-to-North Water Transfers Project, which alleviated the groundwater shortage to some extent. In contrast, the GWSA in the Southwest of the North China Plain showed an upward trend with a maximum rate of increase of about 14.07 mm/year. Groundwater storage in South China generally showed a spatial decline from the center to the north–south direction (Figure 3e). In the northeast of Guangxi, the north of Guangdong, Hunan, and the south of Hubei were areas with increasing GWSA, with a maximum increasing rate of about 12.15 mm/a. The southwest of Guangxi, the southwest of Guangdong, and the north of Hubei were areas where GWSA decreased, with a maximum rate of decrease of about −1.27 mm/a. The GWSA of Northeast China showed a decreasing trend from north to south (Figure 3f). The maximum rate of decrease was about −11.50 mm/a, and the decreasing areas were mainly in the north of Liaoning, the south of Jilin, and the east of Heilongjiang. The maximum rate of increase was about 9.88 mm/a, mainly in the northwest of Heilongjiang Province. It is worth noting that this area in Northeast China has well-developed agriculture. From 2003 to 2016, agricultural water consumption in this area increased from 31.9 billion m³ to 48.9 billion m³, and the total groundwater supply increased from 21.13 billion m³ to 27.37 billion m³. Therefore, agricultural water consumption might have some impact on groundwater reduction in Northeast China.

To further analyze the declining trend of groundwater in the study area, we divided the change rate of groundwater storage according to the equidistant division method. The change rate of GWSA was divided into six levels, which are: dramatic decrease (<−30), rapid decrease (−30 to −15), a slight decrease (−15 to 0), a slight increase (0 to 15), rapid increase (15 to 30), and dramatic increase (>30) (Figure 4). As shown in Figure 4, GWSA decreased more severely in the northwest of Xinjiang, southeast of Tibet, and the central part of the North China Plain. The GWSA mostly shows a slight upward trend at the border between Xinjiang and Tibet, the border between Heilongjiang and Inner Mongolia, South China, and the south of the North China Plain. The percentage of each change level accounting for the area of each region is shown in

Table 1. Percentage of area for different grades of change (%).

Region	Dramatic Increase	Rapid Increase	Slight Increase	Slight Decrease	Rapid Decrease	Dramatic Decrease
Xinjiang	0.00%	0.00%	12.30%	77.30%	6.81%	2.78%
Tibet	0.00%	2.41%	40.61%	44.77%	8.43%	3.78%
Inner Mongolia	0.00%	0.00%	12.92%	87.08%	0.00%	0.00%
North China Plain	0.00%	0.00%	40.29%	44.83%	14.28%	0.60%
South China	0.00%	0.00%	99.32	0.68%	0.00%	0.00%
Northeast China	0.00%	0.00%	46.22%	53.78%	0.00%	0.00%

. From the data in Table 1, it could be concluded that the areas of groundwater reduction in Xinjiang, Inner Mongolia, the North China Plain, Northeast China, and Tibet were all more than 50%. Only groundwater was increased in most areas of South China. It reflected the significant decreasing trend of groundwater resources in China. Therefore, it was necessary to understand the factors that affect the decline in groundwater resources and to take appropriate measures to improve the condition of groundwater scarcity.

3.3. Verification in GWSA Results

Annual-scale GWL data were used to validate the GWSA results. Figure 5a shows the time-series trends of GWSA and GWL in the North China Plain, where GWSA is not deducting GLCRA on an annual scale. In contrast, GWSA in Figure 5c deducts GLCRA. Figure 5b,d show the correlations between GWSA and GWL, before and after deducting GLCRA. The results show that the correlation between GWSA and GWL becomes higher after deducting GLCRA, with the correlation coefficient increasing from 0.8793 to 0.8804.

Similarly, we analyzed the trends of the annual-scale GWSA and GWL time series in South China and Northeast China regions, comparing the correlation between GWSA and GWL before and after subtracting GLCRA, and calculating the RMSE of both. As the GWSA and GWL units are different, the GWSA and GWL were first normalized before calculating the RMSE. In South China, the correlation between GWSA and GWL became higher after deducting GLCRA, with the correlation coefficient increasing from 0.8398 to 0.8410 (Figure 6). In Northeast China, the correlation between GWSA and GWL also increased after deducting GLCRA, with the correlation coefficient increasing from 0.7885 to 0.7898 (Figure 7). Meanwhile, as shown in

Table 2. Root mean square error of annual average GWSA and GWL before and after deduction of gravity losses from raw coal.

Region	RMSE (Before Deducting GLCRA)	RMSE (After Deducting GLCRA)
South China	0.0765	0.0761
Northeast China	0.1301	0.1298

, the RMSE of both GWSA and GWL decreases after TWSA minus GLCRA in the South China and Northeast China regions.

In the Inner Mongolia and Xinjiang regions, we also validated GWSA using GWL. In the Inner Mongolia region, the correlation between GWSA and GWL became higher after deducting GLCRA, with the correlation coefficient increasing from 0.7536 to 0.7531 (Figure 8). In the Xinjiang region, the correlation between GWSA and GWL also increased after deducting GLCRA, with the correlation coefficient increasing from 0.6455 to 0.6458 (Figure 9). In the Xinjiang region, the correlation between GWSA and GWL is relatively small, which may be due to the small number of groundwater-level monitoring wells in Xinjiang and their uneven spatial distribution.

The correlation between GWSA and GWL is greater than 0.75 in North China, South China, Northeast China, and Inner Mongolia, where the number of groundwater-level monitoring wells is high, and the spatial distribution is relatively homogeneous. The correlation between GWSA and GWL is smaller in Xinjiang due to the small number of groundwater-level monitoring wells and the uneven spatial distribution. It is concluded that the GWSA has better consistency with the groundwater-level observation data. Therefore, estimating GWSA using GRACE and GLDAS data is reliable and feasible. The accuracy of the GWSA estimates could be further improved after deducting GLCRA.

In areas where more coal is mined, more gravity falls. Therefore, AMRCP and GWSA were used for comparison in areas with high coal mining volumes (Xinjiang, Inner Mongolia, South China, and Northeast China). Figure 10 shows the annual change trend of GWSA and AMRCP in Xinjiang, Inner Mongolia, South China, and Northeast China, where GWSA removed the influence of the GLCRA. The result showed that AMRCP was significantly and negatively correlated with GWSA in Xinjiang and Inner Mongolia, with correlation coefficients of −0.94 and −0.98, respectively, indicating that the more coal was mined, the lesser the groundwater storage was. In South China and Northeast China, it was more difficult to determine the relationship between the AMRCP and GWSA.

Table 3. Correlation between GWSA and the AMRCP.

Region	Before Deducting GLCRA	After Deducting GLCRA
Xinjiang	−0.9385	−0.9386
Inner Mongolia	−0.9772	−0.9781
South China	−0.1630	−0.1646
Northeast China	−0.1979	−0.2011

shows the correlation between GWSA and the AMRCP before and after removing GLCRA from GWSA in each region. The results shows that the correlation between GWSA and the AMRCP improved in all regions after removing GLCRA, and it indicated that the GWSA could better reflect the changes in groundwater storage after considering the influence of the raw-coal-mining activity.

3.4. Factors Influencing Changes in Groundwater Storage

In recent years, complex climate change and the frequent occurrence of extreme weather should have led to changes in groundwater storage. In addition, groundwater storage was also very sensitive to many human activities, such as water consumption from coal mining, industry, and agriculture. In this paper, we choose natural factors (rainfall and temperature) and human factors (annual water consumption by coal mining, agriculture, and industry) to analyze their effects on the changes in groundwater storage.

3.4.1. Relationship between Natural Factors and Groundwater Storage

Figure 11a–f show the monthly variation in GRACE-based GWSA with rainfall and average air temperature in the six areas. It could be seen from the figure that precipitation and temperature showed a clear seasonal cycle and a clear relationship between changes in groundwater storage and rainfall. During periods of heavy precipitation, the GWSA increased accordingly. The precipitation gradually decreased from November to March of the following year, and the GWSA also showed a decreasing trend. It was calculated that GWSA in Tibet, Inner Mongolia, and South China lagged behind rainfall by one month, three months, and three months, respectively (Figure 11b,c,e). This was probably due to the infiltration of rainfall into the ground, eventually forming groundwater. The GWSA in Figure 11a shows a clear downward trend in 2008 and 2014, but there was no clear downward trend in rainfall in the same year, indicating that there might be other anthropogenic factors, such as water consumption by coal mining, industry, and agriculture, which led to a continuous decline in GWSA. The GWSA in the North China Plain showed a general upward trend from 2003 to 2005, while rainfall showed a decreasing trend (Figure 11d). It was found that although rainfall was low and temperatures were comparatively high in 2004, relatively heavy rainfall and low air temperature occurred in 2003 and 2005. This might be the reason for the rising GWSA in the North China Plain from 2003 to 2005. From 2008 to 2016, the GWSA in the North China Plain continuously decreased. Annual precipitation showed a decreasing trend from 2008 to 2014 and an increasing trend from 2014 to 2016. The decrease in annual precipitation will lead to a decrease in groundwater recharge and result in a decrease in GWSA. However, as annual rainfall increased after 2014, GWSA continued to decrease, indicating that anthropogenic factors such as water consumption by agriculture and coal mines might be the main reason for the decline in GWSA. Monthly rainfall in Southern China showed an increasing trend with a rate of 0.016 mm per month. In addition, the annual rainfall increased at a rate of about 25.21 mm/a. Air temperatures showed an increasing trend with an increased rate of about 0.0036 °C/month (Figure 11e). Both evaporation and rainfall increased, indicating that rainfall may be the main reason for the increase in GWSA in South China.

3.4.2. Relationship between Anthropogenic Factors and Changes in Groundwater Storage

We analyzed the influence of anthropogenic factors, such as annual water consumption by coal mining, industry, and agriculture, on the changes in groundwater storage in six regions of China. Figure 12, Figure 13 and Figure 14 show the relationship between average annual GWSA and the annual water consumption of coal mining, industry, and agriculture in each study area from the perspective of time series, respectively, where the GWSA deducted the effect of the GLCRA.

Figure 12a–e show the relationship between GWSA and annual water consumption of coal mining in the study area. Due to the small production of raw coal in Tibet, the impact on GWSA was small, so it was not analyzed. The results showed that the changing trends of GWSA in Xinjiang, Inner Mongolia, and North China Plain were similar, showing a decreasing trend of fluctuation from 2003 to 2016, and the annual water consumption of coal mines showed an increasing trend in fluctuation (Figure 12a–c). The two were significantly negatively correlated, and the correlation coefficients were −0.94, −0.98, and −0.71, respectively. It was concluded that the annual water consumption of coal mining might be directly related to the reduction in groundwater reserves in Xinjiang, Inner Mongolia, and the North China Plain, which might be an important factor leading to the continuous decrease in GWS in these three regions.

However, there was no significant correlation between GWSA and coal-mining water consumption in South China and Northeast China (Figure 12d,e). The correlation coefficients only reached −0.16 and −0.2, respectively. Hence, it indicated that the effect of coal-mining water consumption on GWSA in South China and Northeast China might be minimal.

Figure 13a–e show the variation in the trend of average annual GWSA and annual water consumption by industry in each study area, respectively. The annual water consumption of industry showed an overall trend of fluctuation increasing in Xinjiang, Inner Mongolia, and North China Plain, which was opposite to the trend of GWSA (Figure 13a,c,d). They also showed a significant negative correlation with correlation coefficients of −0.87, −0.74, and −0.72, respectively. Hence, this indicated that the annual water consumption of industry might be the main factor in the reduction of GWSA in these areas.

However, there was no significant correlation between GWSA and annual industrial water consumption in Tibet, South China, and Northeast China (Figure 13b,e,f). Among them, the correlation coefficient was −0.47 in Tibet, indicating that the annual industrial water consumption might have minimal impact on the change in GWSA in Tibet. However, the correlation coefficients in South China and Northeast China were 0.11 and 0.03, respectively, indicating that annual industrial water consumption might not be the main driver of GWSA changes in these two regions.

Figure 14a–e show the trend of variation between average annual GWSA and annual water consumption of agriculture in each study area, respectively. The annual water consumption of agriculture showed an overall trend of fluctuation increasing in Xinjiang, which was opposite to the trend of GWSA (Figure 14a). They also showed a significant negative correlation with correlation coefficients of −0.81. Therefore, this indicated that the annual water consumption of agriculture might be the main driver of the decrease in GWSA in Xinjiang. However, there was no significant correlation between GWSA and annual agricultural water consumption in North China Plain, South China, and Northeast China (Figure 14d–f). The correlations of these three regions were −0.16, −0.47, and −0.38, respectively. Hence, this indicated that the effect of agricultural water consumption on GWSA in the North China Plain, South China, and Northeast China might be minimal.

Agricultural water consumption formed a significant positive correlation with GWS change in Inner Mongolia, whereas no remarkable correlation was observed in Tibet (Figure 14b,c). The correlation reached 0.72 and 0.27 in Inner Mongolia and Tibet, respectively. However, this association was contrary to the theoretical relationship between these variables. Therefore, agricultural water consumption might not be the main factor affecting the reduction in GWSA in Tibet and Inner Mongolia.

4. Conclusions

Based on data from the GRACE gravity satellite and the data in GLDAS, this study applied the water storage balance equation to calculate and analyze the characteristics of spatial and temporal variations in groundwater in six typical regions of China from 2003 to 2016. We analyzed the correlation between GWSA and groundwater-level data before and after deducting GLCRA to verify the accuracy of GRACE-based GWSA. Based on the meteorological data and the annual statistical data on human water consumption in different regions, we analyzed the effects of climate change and human activities on groundwater storage. The following conclusions were drawn:

• Temporally, the trend of GWSA was significantly decreasing in Xinjiang and the North China Plain, with decrease rates of, respectively, −6.24 mm/a and −7.35 mm/a. The decreasing trend of GWSA was slight in Tibet, Inner Mongolia, and Northeast China, and the rate of decrease was, respectively, −3.33 mm/a, −3.17 mm/a, and −0.75 mm/a. However, South China showed an upward trend with an increase rate of 4.79 mm/a. Spatially, it was obvious that northwestern Xinjiang, Southeastern Tibet, and the central part of the North China Plain have declining trends, and their maximum rates of change, respectively, reached −50.11 mm/a, −56.24 mm/a, and −31.46 mm/a. The border between Xinjiang and Tibet, Northern Tibet, Southern North China Plain, South China, and the north of Northeast China was in a region of increasing trend, and their maximum rates of change, respectively, reached 17.90 mm/a, 21.32 mm/a, 14.07 mm/a, 12.15 mm/a, and 9.88 mm/a.
• In the North China Plain, South China, Northeast China, Inner Mongolia, and Xinjiang regions, the correlation between GWSA and GWL data showed that the correlation between annual mean GWSA and annual mean GWL increased after deducting GLCRA, and all were greater than 0.65. Inner Mongolia, South China, and Northeast China revealed that the correlations between GWSA and ACMP all increased after subtracting GLCRA from GWSA. It was concluded that the estimation of GWSA in this study could reflect the actual changes in GWSA in China, and the accuracy will be improved after deducting GLCRA.
• The annual-scale GWSA deducted GLCRA in each region and was compared with the data of various factors, and it was concluded that rainfall and temperature were the reasons for the periodical fluctuation in GWSA. Among them, in Xinjiang, the annual water consumption by coal mines, industry, and agriculture might be the main drivers of the continued decline in GWSA. In Inner Mongolia and the North China Plain, the annual water consumption by coal mining and industry might be the main drivers for the continued decline in GWSA. However, in South China, rainfall might be the main reason for the continuous increase in GWSA, and rainfall recharge was larger than groundwater consumption.