China is located east of Eurasia and west of the Pacific Ocean, at latitudes and longitudes ranging from 73°22′E to 135°03′E and 3°51′N to 53°31′N. It covers a total land area of about 9.6 million km² and comprises 34 provincial administrative units (Figure 1). As the largest developing country in the world, mainland China had a total population of 1.395 billion at the end of 2018, which accounted for about 1/5th of the global population. Its GDP, at about $13.6 trillion, was ranked second in the world in 2018 (NBSC, 2019). However, China has been identified as “water short” by the United Nations, and the average water resource per capita, at 2085.72 m³ between 1988 and 2022, is 1/4 ^th of the global average and is ranked 128 out of 183 countries worldwide (FAO, AQUASTAT Database, https://www.fao.org/aquastat/en/databases/). Large amounts of groundwater are extracted to satisfy the water demand and this has resulted in serious GWS decline. For example, the GWS in China decreased from 850.186 billion m³ in 2000 to 824.65 billion m³ in 2018 (MWR, 2018). The number of groundwater extraction wells increased from 4 million in 1988 (Chen and Ma, 2017) to 97.49 million in 2011 (MWR, 2011). The climate in China varies significantly between different regions, and can be divided roughly into the eastern monsoon climate zone, northwest temperate continental climate zone, and the Qinghai-Tibet alpine climate zone. Similarly, socioeconomic development is much higher in the east and south than in the north and west.

The effects of both natural and human activities on GWS at national, regional, and provincial scales were analyzed to ensure the results were complete and represented the climatic, geographical, and administrative variations. Mainland China is generally divided into seven geographical regions, namely, Northern China, Northeast China, Eastern China, Central China, Southern China, Southwest China, and Northwest China (Figure 1) (Tian and Zhang, 2020; Chen et al., 2021). Moreover, because of the variability in socioeconomic activities across the nation, decision-making about socioeconomic development is commonly targeted at provincial and regional administrative units. Therefore, to account for the national, regional, and provincial scales, the temporal changes in GWS for the entire mainland of China was quantified, for the seven geographical regions, and for 31 mainland provinces, respectively. Hong Kong, Macao, and Taiwan were not included because no data were available.

Northern China covers four provincial-level administrative units (Beijing, Tianjin, Hebei, and Shanxi) and eight prefecture-level cities in the Inner Mongolia Autonomous Region (Hohhot, Baotou, Hulunbuir, Xingan League, Tongliao, Chifeng, Xilin Gol League, and Ulan Qab League); Northeast China covers three provincial-level administrative units (Liaoning, Heilongjiang, and Jilin); Eastern China covers seven provincial-level administrative units (Shanghai, Jiangsu, Zhejiang, Anhui, Jiangxi, Fujian, and Shandong); Central China covers three provincial-level administrative units (Henan, Hubei, and Hunan); Southern China covers three provincial-level administrative units (Guangdong, Guangxi Zhuang Autonomous Region, and Hainan); Southwest China covers five provincial-level administrative units (Chongqing, Sichuan, Guizhou, Yunnan, and Tibet Autonomous Region), and Northwest China covers five provincial-level administrative units (Shaanxi, Gansu, Qinghai, Ningxia Hui Autonomous Region, and Xinjiang Uygur Autonomous Region) and four prefecture-level cities in the Inner Mongolia Autonomous Region (Alxa League, Bayannaoer, Wuhai, and Ordos).

The framework for investigating the spatial-temporal variation of GWS change and driving forces across China was developed as shown in Figure 2.

The GRACE mission, launched in March 2002 by NASA and Das Deutsche Zentrum für Luft-und Raumfahrt (DLR), monitors the changes in the earth’s time-varying gravity field, which reflected variations in the terrestrial water storage (TWS) (Long et al., 2020). Therefore, GRACE data does not provide information about the total (absolute) amount of water stored at a location, but represent anomalies with respect to the mean state [i.e., the average surface mass during the baseline time (January 2004–December 2009)], known as Terrestrial Water Storage Anomalies (TWSA), which characterize the relative TWS. In this study, the monthly TWSA products from April 2002 to December 2016 at a resolution of 0.25 ° × 0.25 ° from the Center for Space Research’s TELLUS website (http://www2.csr.utexas.edu/grace) were obtained. There was a 19-month gap in the data of the study period because of issues with the GRACE satellite sensor. The missing values for this period were interpolated by averaging the values of the 2 months before and 1 month after the month without data (Long et al., 2015; Yin et al., 2018).

In addition, other properties of terrestrial water from GLDAS data and WGHM (WaterGAP Global Hydrology Model) data were obtained, including soil moisture (SM), snow water (SN), and surface water (SW). SM and SN data were obtained at a spatial resolution of 0.25 ° × 0.25 ° from the GLDAS-2/NOAH model (monthly product) (https://search.earthdata.nasa.gov/search?q=GLDAS). The SW data from information about the amount of reservoir, lake, river, and wetland storage from the monthly product of the WGHM 2.2 days model were obtained at a spatial resolution of 0.5 ° × 0.5 ° (Huang et al., 2019; Müller Schmied et al., 2021) (https://doi.pangaea.de/10.1594/PANGAEA.918447). The WGHM 2.2 days data were re-gridded to a horizontal grid resolution of 0.25 ° × 0.25 ° to match with the resolution of the GRACE and GLDAS data using ArcGIS 10.6 software (ESRI, Environmental Systems Research Institute). As it was to ensure comparisons against the TWSA are consistent, the temporal mean of each dataset were computed over a common time period and subtracted from the respective time series to obtain soil moisture anomalies (SMA), snow water anomalies (SNA), and surface water anomalies (SWA) (Huang et al., 2015; Chen et al., 2019).

According to the terrestrial water balance equation, TWS is composed of SW, GWS, SN, SM, and biomass water (which is negligible in most cases) (Liesch and Ohmer, 2016). Therefore, the TWSA based on GRACE data can be represented by the anomalies of the different terrestrial water components, as follows: T W S A = S W A + S M A + S N A + G W S A (1)

Of these, TWSA and SWA were derived from the GRACE and WGHM data, respectively, and SNA and SMA were derived from GLDAS data. Eq. 1 therefore becomes: G W S A = T W S A − S W A − S M A − S N A (2)

The GWSA represents the anomalies in GWS relative to the baseline temporal average (i.e., the relative GWS), usually expressed by the equivalent water thickness (cm). For different spatial scales, the GWSA was calculated as the mean value of all the grids summarized by the spatial domain of the mainland, geographic regions, or administrative provinces accordingly.

The GWS change (∆GWS) is the time derivative of the GWSA, which is referred specifically to as dGWSA/dt at the monthly scale, using the simple derivative in Eq. 3: ∆ G W S t ≈ d G W S A d t ≈ G W S A t − G W S A t − 1 ∆ t (3)

The linear change trend analysis was used to identify the GWS change during the period of 2002–2016. Moreover, we highlighted the difference in the GWS change between the different three time periods (2002–2006, 2007–2011, and 2012–2016), which approximately aligned with the Chinese national 5-year development programs. These national development plans have deeply affected the socioeconomic development of China, and the related development policies were often substantially adjusted every 5 years.

The accuracy of the GWS estimates was verified against in-situ groundwater water level observations, which were obtained from 289 official national-level monitoring stations (Figure 1). Because of issues with data availability, data were not available for the entire study period. The annual average groundwater level data were obtained from 2005 to 2016 in the China Geological Environment Monitoring Groundwater Level Yearbook. To match the annual time scale of this study, the monthly groundwater level data as annual data was averaged and processed into anomalies, with 2005–2009 as the baseline. The estimated GWS anomalies derived from GRACE were verified by comparing with the observation data using the Pearson correlation coefficient (Nanteza et al., 2016; Gao et al., 2020).

Referring to previous studies (Chen et al., 2018; Chen et al., 2020; Zhao et al., 2020), 10 natural and anthropogenic factors that were probably directly or indirectly related to GWS change were selected, namely, precipitation, evapotranspiration, total population, agricultural water consumption, industrial water consumption, domestic water consumption, GDP (Gross Domestic Product) of primary industries, GDP of secondary industries, GDP of tertiary industries, and total GDP. All the data were extracted from Statistical Yearbooks and Water Resources Bulletins from 2002 to 2016.

The Grey Relational Analysis (GRA) to reflect the relationships between GWS change and driving factors was used (Chen et al., 2020). GRA is based on geometrical mathematics and complies with the principles of normality, symmetry, entirety, and proximity, and has been useful for resolving questions about the complex interrelationships between multiple factors and research objects. Thus, GRA can show the similarity in the changes and trends of different subsystems, and is more suitable for systems with a small number of samples than traditional regression analysis (Yi et al., 2021).

In GRA, the closeness of the relationship between the changes in GW storage and multiple influences are judged by determining the geometrical similarity between the reference sequences Y = {Y(k)|k = 1, 2, … , n} and the different comparison series Xi = {Xi(k)|k = 1, 2, … , n}, i = 1, 2, … , m. In this study, the change in the GWS was determined as the reference sequence, while the influential factors were determined as the comparison sequences. Each sequence was initialized separately to unify the dimension (Deng, 1990). The grey relational coefficient was calculated as below. The grey relational grade represents the level of correlation between the reference and comparability sequences and, as the grey relational grade increases, the correlation between the comparability and reference sequences strengthens. r i = 1 n ∑ k = 1 n ξ i k (4) Where r represents the grey relational grade, and i represents the spatio-temporal sequence of the i-th driver, and i = 1, 2,..., 10 represent the 10 drivers. The time length and the time node of the sequence are represented by n and k, respectively. The study period was from 2002 to 2016, so n = 15, and k = 1, 2, … , 15. ξi(k) is the correlation coefficient, that is, the degree of correlation between the reference sequence and the comparison sequence at each time node: ξ i k = min i min k y k − x i k + ρ max i max k y k − x i k y k − x i k + ρ max i max k y k − x i k (5) Where, ξ_i(k) is the grey relational coefficient, ρ is the resolution coefficient, and the strength of the resolution increases as the value decreases, usually ρ = 0.5 (Mo, 2021; Reddy et al., 2021).

The change in GWS inverted by GRACE fitted well with the in-situ observation data, with a high significant correlation (value of 0.87, p < 0.01) between the GRACE-based GWS change and the in-situ observations. The results at the grid scale showed obvious spatial and temporal changes in GWS from 2002 to 2016 (Figure 3).

At the national scale (Figure 4), the mean annual GWS change (i.e., relative GWS) declined by 4.13 cm from 2002 to 2016 at a rate of −0.27 cm/a (p < 0.01). The GWS change fell into three different periods, namely 2002–2006, 2007–2011, and 2012–2016. The rate of decline in the GWS across the entire mainland was highest for 2007–2011 [−0.42 cm/a (p < 0.05)], followed by 2012–2016 (−0.22 cm/a), and was lowest for 2002–2006 (−0.15 cm/a).

Our analysis across the seven geographic regions showed that the GWS changes varied considerably across the seven regions (Figure 5). From 2002 to 2016, the GWS decreased significantly in Northern, Northwest, Eastern, Northeast, and Southwest China while it increased significantly in Southern China (0.32 cm/a, p < 0.01), but not in Central China. The rate of decline was most rapid in Northern China (−0.64 cm/a, p < 0.01) and was slowest in Eastern China (−0.11 cm/a, p < 0.05).

The regional GWS change also varied between different periods (Figure 6). Over the three time periods, the GWS decreased continuously in Northern and Northwest China, increased continuously in Central China, and fluctuated in Northeast, Southwest, Eastern, and Southern China. The rate of decrease in the GWS during the middle period (2007–2011) was most rapid in Northern China (−1.19 cm/a, p < 0.01), while the decrease in the GWS in the later period (2012–2016) was most rapid in Northwest China (−0.71 cm/a, p < 0.01). In contrast, the region with the most rapid GWS increase in the early period (2002–2006) was Central China (0.41 cm/a, p < 0.1). In the regions where the GWS changes fluctuated, the GWS first decreased and then increased in Northeast, Eastern, and Southern China, but first increased and then decreased in Southwest China. The acceleration or deceleration of the GWS in five of the seven regions show that the regional-scale GWS improved significantly in the later period.

The GWS changes varied considerably among the different provinces within each region, especially in Southwest, Central, and Eastern China (Figure 7). For example, the GWS decreased in Southwest China at an overall average rate of −0.12 cm/a (p < 0.05) for 2002–2016. Out of the five provinces in the region, only the GWS for Tibet decreased (−0.53 cm/a, p < 0.01), but the GWS for the other four provinces increased significantly, with rates of increase ranging from 0.21 to 1.03 cm/a. In Central China, the change in GWS was not significant at the regional scale, but was significant for some of the individual provinces within it. The GWS of Hunan Province increased at a rate of 0.48 cm/a (p < 0.01) and the GWS of Henan province decreased at a rate of −0.66 cm/a (p < 0.01). There was an overall regional decrease in Eastern China (−0.12 cm/a, p < 0.05), but the GWS increased significantly in two provinces, Jiangxi and Shanghai, at rates of 0.33 (p < 0.05) and 0.17 cm/a (p < 0.05), respectively. The trends in the GWS changes in the provinces in the other four regions (Northern, Northeast, Northwest and Southern China) were consistent with the overall change for the region.

The GWS changes in the provinces varied in the different time periods (Figures 7B–D). The GWS change at the provincial-scale for 2002–2006, 2007–2011, and 2012–2016 either increased continuously, decreased continuously, or fluctuated. The GWS increased continuously in three of the 32 provinces, namely Chongqing (Southwest China), Guangdong (Southern China), and Anhui (Eastern China). The GWS decreased continuously from 2002 to 2016 in eight provinces, mainly in Northern (e.g., Beijing) and Northwest (e.g., Xinjiang) China. The GWS fluctuated and increased in 2012–2016 in 13 of the 21 provinces in Northeast China (e.g., Liaoning), Eastern China, Central China, and Southern China (e.g., Guangxi), but showed decreasing trends in the remaining eight of these 21 provinces, mainly distributed in Southwest China (e.g., Sichuan) and Northwest China (e.g., Qinghai). Overall, the increase in the GWS for 2012–2016 was greatest in Liaoning and Guangxi, and the decrease was greatest in Beijing, followed by Shanghai, Jiangsu, Xinjiang, and Tibet. Also, fewer provinces had declines in GWS for 2012–2016 than for 2002–2006, which implies that the stress on groundwater exploitation had been relieved in some provinces. The provincial- and regional-scale results were consistent.

The GRA showed that the most important drivers of GWS change at all three spatial scales were precipitation, agricultural water consumption, and the total population, and that some other factors were important in some regions and provinces (Figure 8; Table 1). At the national scale, the grey relational coefficient (GRC) was highest between the GWS change and precipitation (0.87), followed by the agricultural water consumption (0.86) and the total population (0.86). The drivers showed variation at the regional scale. For example, precipitation had the highest GRC in Eastern China (0.83), Central China (0.84), and Southern China (0.84), agricultural water consumption had the highest GRC in Northern China (0.86), and industrial water consumption had the highest GRC in Northeast China (0.77). In contrast, it was difficult to assess the drivers of GWS change in Southwest and Northwest China. For example, precipitation, evapotranspiration, total population, and different types of water consumption had GRC values of 0.76 in Southwest China, while precipitation, industrial water consumption, and domestic water consumption had GRC values of 0.86 in Northwest China. Additionally, there were some differences between the main drivers at the provincial and regional levels. The GRC of the highest ranked drivers at the provincial and regional levels were different, and the primary drivers of 12 provinces were inconsistent with those of the region to which they belonged, mainly in Eastern, Southern, Northern, and Northwest China. For example, industrial water consumption had the highest GRC in Shandong province (Eastern China), while precipitation had the highest GRC at the regional scale. The reason was there existed obvious variability in natural, social, and economic characteristics between different provinces, which was inconsistent with those of the upper regional scale. This resulted in the variety of dominant driving factors across different spatial scales.

China’s GWS declined continuously over the study period, despite the great awareness of the need to protect groundwater. At the national scale, the trend in the GWS was generally inversely related to the changes in the groundwater usage, and a decrease in the GWS was usually accompanied by an obvious increase in the groundwater usage. For example, the national-scale GWS showed an overall decrease (−0.27 cm/a, p < 0.01) from 2002 to 2016, but the groundwater usage increased at a rate of 0.47 billion m³/a (p < 0.05). The decline in the GWS was fastest from 2007 to 2011, mainly reflecting the 3.74% increase in the groundwater usage. The slowdown in the GWS depletion after 2012 shows the results of China’s efforts to control groundwater extraction. For example, China’s groundwater withdrawal decreased by 5.73% from 2012 to 2017, which reflects the effectiveness of the strict groundwater protection measures adopted by the Chinese government. However, relative to the gains reported in other countries, there is still much room for improvement in groundwater conservation in China. For instance, the amounts of groundwater used in Egypt and Ukraine decreased by 13.33% and 47.33% from 2012 to 2017, respectively (FAO, AQUASTAT Database, https://www.fao.org/aquastat/en/). The groundwater will continue to be an important source of water for China into the future. Thus, groundwater in China will face ongoing pressure.

The differences in the GWS changes between the north and the south, and the east and the west, increased. For example, the GWS continued to increase in Southern China from 2002 to 2016 (0.32 cm/a, p < 0.01), but continued to decrease in Northern China (−0.64 cm/a, p < 0.01). Even so, the gap in the GWS change between Northern and Southern China narrowed, because the GWS increased at a faster rate in Southern China and decreased at a much slower rate in Northern China from 2012 to 2016 than in 2007–2011. The accelerated increase in the GWS in Southern China resulted from a 16.21% increase in the average precipitation and 8.31% decrease in the groundwater usage. In contrast, the ongoing decline in the rate of decrease in the GWS in Northern China reflects the combination of an increase of 15.21% in the average precipitation and a decrease of 13.73% in the groundwater usage. So, the net effect of the increase in precipitation and the decrease in the groundwater usage caused the slowdown in the GWS depletion in Northern China. Similarly, the imbalance between groundwater exploitation and replenishment in the west of China intensified from 2002 to 2016, causing an increase in the rate of GWS depletion, and a widening of the gap in the GWS decline between the east and the west. For example, in the absence of significant precipitation changes from 2002 to 2016, groundwater exploitation increased by 44.76% in Northwest China, but decreased by 27.79% in Eastern China (MWR, 2018). The groundwater usage, and so the reliance on groundwater, only increased in Northwest China. Therefore, measures to protect the groundwater usage in the Northwest of China should be urgently minimized.

There was significantly more variation in the GWS changes at the provincial scale than at the regional scale. For example, the GWS in Southwest China decreased (−0.12 cm/a, p < 0.05) from 2002−2016, but increased significantly in four of the five provinces in this area, with a maximum rate of increase of 1.03 cm/a (p < 0.01) in Chongqing. Similarly, the GWS showed a decrease in Eastern China (−0.12 cm/a, p < 0.05), but increased significantly in Jiangxi (0.33 cm/a, p < 0.05) and Shanghai (−0.17 cm/a, p < 0.05). In Central China, there was no significant change in the GWS, but there were significant changes in the GWS in some provinces in this region. These results remind us that some of the variability will be hidden by examining data at a single scale, and that the variability in the GWS changes should be examined at a range of scales.

Some of the drivers of GWS changes were common across different scales. At the national scale, the main drivers of the GWS change were precipitation, agricultural water consumption, and the total population, while other factors, such as industrial and domestic water consumption and evapotranspiration, were also significant at regional and provincial scales. This shows that changes in these factors can significantly affect the GWS changes at different scales. For example, the strong correlation between precipitation and GWS change was evidenced by the fact that the average precipitation across China decreased from 2002 to 2011 (by 11.94%) and then increased from 2011 to 2016 (by 25.98%), and the increase in the later period greatly slowed the decline in the GWS. Additionally, China is a major agricultural country and so has to address the challenge of high-water consumption due to inefficient irrigation in agriculture. The agricultural sector in China accounts for more than 60% of the water withdrawal, but the effective use coefficient of farmland irrigation water (that is, the ratio of the actual effective use of water irrigated into the field, excluding deep seepage and field loss to the total water introduced from the canal) in China was only 0.55 in 2019, compared with 0.89 in Israel (Martinho, 2020). These figures indicate that the agricultural water demand in China could be reduced by developing high-efficiency water-saving irrigation techniques and adjusting agricultural planting structures (Qiu, 2010).

At regional and provincial scales, other anthropogenic and natural factors, such as water consumption by different sectors and evapotranspiration, noticeably affected the changes in the GWS. For example, in Hebei (Northern China), agricultural water consumption had the highest GRC, followed by industrial water consumption and domestic water consumption. From 2002 to 2016, Hebei’s surface water resource accounted for only 0.23% of the whole nation’s water resource, so about 78.58%, 78.59%, and 83.65% of Hebei’s agricultural water, industrial water, and domestic water came from groundwater, respectively (DWRHP, 2016). In Gansu (Northwest China), the evapotranspiration increased by 20.44% from 2002 to 2016, which potentially reduced the recharge from surface water to groundwater. The differences in the drivers of the GWS change suggest that specific polices are needed for different scales.

The variations in the GWS change and the stresses at different spatial scales across China also reflected the effect of policies. The slowdown in the GWS decline in China from 2012 to 2016 resulted from the implementation of a series of policies, such as water saving, water diversion, and water replenishment. For example, at the national level, policies such as “The National 12th Five-Year Plan” (https://www.gov.cn/) and “The National Land Subsidence Prevention and Control Plan (2011–2020)” focused on promoting a water-saving society, endorsing efficient irrigation technology, and curbing the exploitation of groundwater. With these plans, China’s agricultural water consumption and groundwater extraction decreased by 3.45% and 6.81%, respectively, from 2012 to 2016. Accordingly, specific solutions were adopted by regions and provinces to achieve the national targets. For example, transregional water diversion programs were proposed to alleviate the pressure of groundwater depletion in the north, such as the South-to-North Water Diversion that spans Eastern, Central, and Northern China. By diverting water to Beijing via the South-to-North Water Diversion, the groundwater depletion has decreased by about 3.7 km³, accounting for 40% of the total GWS recovery from 2006 to 2018 (Long et al., 2020). At the provincial level, some provinces have adopted their own groundwater protection policies, such as the “Administrative Measures for the Water Intake Permit System of Hebei Province” (http://slt.hebei.gov.cn/) and the “Groundwater Regulations of Shaanxi Province” (http://www.shaanxi.gov.cn/). In Hebei (North China), 1196 urban wells were closed in 2016, and the ecological water supplement in Shaanxi (Northwest China) increased by 68.28% from 2016 to 2020. China has implemented stricter measures to protect groundwater, and groundwater pumping during the period from 2016 to 2020 was 50% lower than that for 2012 to 2016; also, from 2016 to 2020, the ecological water supplement increased to be five times greater than it was in 2012–2016. The effects of these groundwater protection policies in more detail need to be evaluated in the future.

There were some limitations to this study. For example, 2002–2016 was chosen as the study period because WHGM data were only available up to 2016. As shown in Figure 1, the distribution of national-level stations is not uniform, which would have impacted the accuracy of the results, for example, the accuracy verification of the estimates of GWS derived from remote sensing data. As the groundwater monitoring system is extended, and more in-situ monitoring data become available, the verification accuracy can improve. Additionally, also because of data availability, the hydrological flux variables associated with the physical processes of groundwater changes were not considered, such as groundwater recharge from precipitation and surface water, withdrawals and pumping, capillary rise, and groundwater discharge to streams and lakes. The investigation was focused on the most relevant drivers of GWS change rather than quantifying the precise contribution of each factor. Therefore, when exploring the contribution of climatic and anthropogenic drivers, the collinearity of the variables should be considered. Moreover, as the main aim of this study was to identify the different impacts of human activities on GWS change, the study was based on administrative boundaries rather than natural aquifers.