Analysis of Spatiotemporal Groundwater-Storage Variations in China from GRACE

Abstract

In China, where some regions are over-reliant on groundwater, groundwater consumption is faster than replenishment, which results in a continuous decrease in the groundwater level. Here, we applied spatial and temporal methods to analyze the spatiotemporal variations in groundwater in China from GRACE, GRACE-FO, and GLDAS data. From a national perspective, groundwater storage showed a decreasing trend in northern China and an increasing trend in southern China. The results showed that the rates of groundwater depletion in North China, the Loess Plateau, and Northwest China were −10.09 ± 0.94, −10.05 ± 1.05, and –4.91 ± 0.28 mm y⁻¹ equivalent height of water from 2003 to 2019, respectively. Furthermore, the groundwater in South China, the middle-lower Yangtze River, and the Ch-Yu region had a positive trend, with rates of 7.26 ± 1.51, 7.73 ± 1.35, and 3.61 ± 0.53 mm y⁻¹ equivalent height of water, respectively. We also found that groundwater storage fluctuated slightly before 2016 on the Qinhai-Tibet Plateau and in Northeast China and decreased significantly after 2016. The Yun-Gui Plateau had a fluctuating trend. Investigating the spatiotemporal variation in groundwater storage in China can provide data for initiating regional ecological and environmental protection.

1. Introduction

Groundwater is a key resource supplying water to billions of people and sustaining agricultural, industrial, and domestic activities [1]. It is often the last freshwater resource available for supplying water for domestic use and irrigation after the depletion of surface water in semiarid areas and densely populated countries [2]. More than 1.5 billion people worldwide rely on groundwater as their primary source of drinking water [3]. Due to extreme climatic episodes, population growth, and extensive groundwater reduction and overexploitation, long-term groundwater depletion has occurred in many regions of the world, exacerbated by various adverse environmental impacts, including groundwater-quality deterioration, land subsidence, seawater intrusion, and soil salinization [4,5]. Groundwater depletion has been considered a serious threat to the sustainability of water supplies [6].

However, monitoring groundwater storage (GWS) is limited in many regions [7]. The absence of in situ monitoring of well networks, low-quality observation data, and uncertainties in storage coefficients for translating groundwater-level changes to groundwater storage variations restrict our research on spatiotemporal variations in GWS [8,9]. In this context, long-term reliable and frequent observations of aquifer levels are necessary to monitor the variations in groundwater storage. The launch of the Gravity Recovery and Climate Experiment (GRACE) mission in March 2002 offered new perspectives for the monitoring of terrestrial parts of the hydrological cycle, including groundwater resources [10].

Feng et al. found that the rate of groundwater depletion based on GRACE was 2.2 ± 0.3 cm/a from 2003 to 2010 [11]. Rodell et al. used terrestrial water storage change observations from GRACE data and simulated soil-water variations from a hydrological modeling system to show that, from August 2002 to October 2008, groundwater depletion in northwestern India was equivalent to a net loss of 109 km³ of water, which is double the capacity of India’s largest reservoir [12]. Scanlon et al. compared updated estimates of groundwater storage changes in the Central Valley of California using GRACE satellites with storage changes from groundwater-level data. They found that groundwater depletion totaled 31.0 ± 3.0 km³ in this area [13]. Yin et al. used GRACE, hydrological data, and Global Hydrological Model (WGHM) to analyze the groundwater storage changes in Tasmania (Australia), and their results showed that this area experienced a pronounced GWS decline during 2003–2010 at a depletion rate of −2.57 mm/y, mainly due to decreasing precipitation [14]. Alexandra et al. demonstrated that a groundwater depletion rate in the Northwest Sahara Aquifer System of 2.69 ± 0.8 km³/y would result in the aquifer being depleted to 90% of its total storage in as few as 50 years [15]. Gholamreza et al. found that across the Middle East, the largest groundwater depletion occurred in Iran based on GRACE data, land surface models, and well observations, with a mean loss rate of 25 ± 3 Gt/y from February 2003 to December 2012 [16]. However, the lack of temporal and spatial continuity of data, such as monitoring well data and hydrological models, poses a challenge to using GRACE to retrieve groundwater storage changes. The application of GLDAS data, to a large extent, solves the issue of spatiotemporal discontinuity, algorithm inconsistency, and availability of data, which makes it possible to isolate other component variables from terrestrial water storage (TWS).

In situ data and other data, as mentioned above, are also lacking in some parts of China. Groundwater is an important resource for production and life in most areas of China. Since the reform and opening up, with the development of the economy and a population explosion, groundwater has diminished in many areas of China. These changes have also attracted the attention of Chinese scholars. The research areas based on GRACE data to monitor the changes in groundwater storage are mainly distributed in the Haihe River Basin [11,17], Yellow River Basin [18,19], Hexi Corridor-Alxa Inner Flow Area, northwestern China [20], Liaohe River Watershed [21], and Songhua River area [22] in the northeastern and southwestern Karst areas [23]. Studies have focused on small-scale areas or zones to conduct research on groundwater in China. Therefore, it is necessary to grasp the changes in groundwater storage in China over a long-time scale and at a larger spatial scale.

In this study, we used GRACE and GRACE-FO terrestrial water storage anomalies and Global Land Data Assimilation System (GLDAS) data to analyze the variations in groundwater storage in China during 2003–2019 based on spatial and time-series analytical methods. We hope that our results can provide data support for others to grasp the temporal and spatial changes in groundwater storage in China and to carry out ecological environmental restoration. In Section 2, we introduce the basis of the study area and the data and methods. In Section 3, we analyze all of China utilizing hotspot analysis, regional analysis, and seasonal analysis of changes in groundwater storage, and a summary is provided in Section 4.

2. Materials and Methods

2.1. Study Area

China is a vast territory with different terrain topography, human activities, agricultural production, etc., that cause changes in groundwater regions. China has increasingly experienced severe water scarcity with rapid population growth and socioeconomic development. Especially in North China, the Loess Plateau, and Northwest China, due to the drying climate and intensive human activities in recent decades, groundwater has been used and even overexploited to meet ever-expanding water demand [19,20,24]. In southern areas of China, including southern China and the middle–lower Yangtze River, although there is a phenomenon of excessive reliance on groundwater, due to the impact of precipitation, these areas generally have more groundwater recharge than the water used; therefore, groundwater storage has shown a positive trend.

In this study, we created a zonation for China using Hu’s plan [25] based on 9 regions: North China (NC), Loess Plateau (LP), Middle-Lower Yangtze River (MLY), South China (SC), Northeast China (NEC), Northwest China (NWC), Qinhai-Tibet Plateau (TP), Ch-Yu Region (Ch-Yu) and Yun-Gui Plateau (Yun-Gui) (Figure 1). Due to the wide area and complex terrain in China, the distribution of groundwater in each region has certain characteristics. In NWC, NEC, and NC, the confined aquifers are mainly the synclinal or large fault zones of Paleozoic limestone, intermountain active basins, and piedmont alluvial fan zones, with the unconfined aquifers being mainly distributed in mountain bedrock cracks, partly karst landforms, deserts, permafrost, alluvium, and diluvium. In the LP region, the groundwater table is at a depth of more than 50 m owing to the thick loess strata, where the unconfined aquifer is the alluvium in the valley basin or the pebble bed at the bottom of the loess terrace. In Yun-Gui, MLY, Ch-Yu and SC, the confined aquifers are mainly synclinal structures of Paleozoic and Mesozoic strata and strata belts of intermountain basins, and the unconfined aquifers in these areas are alluvium in the valley, plains, and rock weathering zones in mountain land. In the TP region, the unconfined aquifers are concentrated in swamps, glaciers, and frozen soil landforms, with groundwater in some valleys and intermountain basins [26].

2.2. Data

2.2.1. GRACE Data

The GRACE RL06 level 3 filtered harmonic products from the terrestrial water storage data during January 2003–December 2019 were employed in this study from the CSR (Center for Space Research at University of Texas, Austin, TX, USA). The data are available at https://search.earthdata.nasa.gov/search/ (accessed on 5 October 2020). The positive anomaly shows mass gain, and the negative anomaly indicates mass loss. The monthly land-mass grids contain water-mass anomalies given as equivalent water thickness that represents the total terrestrial water storage anomalies from soil moisture, snow, and surface water (including rivers, lakes, reservoirs, etc.), as well as groundwater and aquifers. Here, 1° × 1° grids and monthly data are used, and we used the scale factor to restore the signal. To ensure the continuity of the data sequence, we used the mean value of two adjacent months to interpolate the data of the missing months of GRACE, and we used the mean values of the same month of the two adjacent years to fill the 11-month gap of CRS TWS product. Although the solution is simple, it is common and efficient for a time-series analysis of TWS and GWS. Firstly, the water mass anomalies product was calibrated relative to the 2004.0–2009.999 time-mean baseline. Here, the TWS data from pre- and post-gap were time consistent, which indicates the mean approach had little impact on the trend estimate. In addition, the gap accounted for less than 5% of the whole study period (204 months), which also has little impact on trend estimates in the long time-series analysis.

2.2.2. GLDAS Data

The total land-water storage anomalies were aggregated from the Global Land Data Assimilation System (GLDAS) NOAH model (https://search.earthdata.nasa.gov/ (accessed on 24 October 2020). GLDAS outputs land-water content by using numerous land surface models and data assimilation. The aggregated land-water anomalies (sum of soil moisture, snow, and canopy water) provided here can be used for comparison against and evaluations of the observations of the Gravity Recovery and Climate Experiment (GRACE) and GRACE-FO over land. The monthly anomalies are computed over the same days during each month as GRACE and GRACE-FO data and are provided on monthly 1 × 1 grids and monthly scales. We performed the same preprocessing as for the GRACE data.

2.3. Methods

2.3.1. Groundwater Storage Change Estimation from GRACE

GRACE-derived TWS changes represent an integrated change in various water storage components, i.e., soil moisture storage (SMS), snow water equivalent storage (SWES), canopy-water storage (CWS), surface runoff storage (SRS), and groundwater storage (GWS). To separate GWS changes from GRACE-derived TWS changes, other water storage components of TWS must be estimated from hydrological models (GLDAS). Thus, changes in GWS can be expressed as follows in Equation (1) [27]:

$$\Delta GWS=\Delta TWS-\left(\Delta SMS+\Delta SWES+\Delta CWS+\Delta SRS\right)$$

(1)

where $\Delta GWS$ is the average change in groundwater storage; $\Delta TWS$ is the average change in terrestrial water storage; and the sum of $\Delta SMS,\Delta SWES,\Delta CWS,\mathrm{and}\Delta SRS$ is the surface-water storage.

2.3.2. Getis-Ord Gi*

For spatial analysis, a map of changes in groundwater storage was generated based on groundwater storage data obtained from Getis-Ord Gi* statistics [28]. The Getis-Ord Gi statistics for each feature in the dataset are z-scores and p-values, which indicate where features with either high or low values cluster spatially. For statistically significant, positive Z scores, the larger the Z score is, the more intense the clustering of high values (hot spot). For statistically significant, negative Z scores, the smaller the Z score is, the more intense the clustering of low values (cold spot). Here, each province’s groundwater storage is a 1 × 1 grid using area-weighted and aggregation methods. The value of a grid in the map is calculated using Equation (2) as follows:

$$\begin{gathered} G_i^{*}=\frac{{{\displaystyle \sum }}_{j=1}^n{\omega }_{i,j}{x}_j-\overline{x}{{\displaystyle \sum }}_{j=1}^n{\omega }_{i,j}}{S\sqrt{\frac{n{{\displaystyle \sum }}_{j=1}^n{\omega }_{i,j}^2-{({{\displaystyle \sum }}_{j=1}^n{\omega }_{i.j})}^2}{n-1}}}\left(i,j=1,2,\cdots ,n\right), \\ where,\overline{x}=\frac{{{\displaystyle \sum }}_{j=1}^n{x}_j}{n},S=\sqrt{\frac{{{\displaystyle \sum }}_{j=1}^n{x}_j^2}{n}-{\overline{x}}^2} \end{gathered}$$

(2)

where $x_j$ is the value of equivalent water height in grid j, and ${\omega }_{i,j}$ is the weight of the i-th grid and j-th grid. We used a fixed threshold of 5 grids. When the distance between the i-th and j-th grids is less than the threshold, ${\omega }_{i,j}$ = 1; otherwise, ${\omega }_{i,j}$ = 0. $n$ is the sum of the grids.

2.3.3. Trend Analytical Methods

In this study, two nonparametric methods (Mann–Kendall and Sen slope) were used to detect the changing trend of groundwater storage.

We used the nonparametric [29] procedure for estimating the slope of the trend in the sample of n pairs of data:

$$\beta =Median\left(\frac{x_j-x_i}{j-i}\right)$$

(3)

where $i=1\dots \dots \dots n$,$j>i$, $x_j\mathrm{and}{x}_i$ are the data values of groundwater storage at times$j\mathrm{and}i\left(j>i\right),\beta$ is used to judge the time series trend, and the groundwater storage has increased if $\beta$ > 0; otherwise, the groundwater storage has decreased [30].

The Mann–Kendall test statistic S [31,32] is calculated as:

$$S={\displaystyle \sum }_{k-1}^{n-1}{\displaystyle \sum }_{j=k+1}^n{S}_{gn}(X_j-X_k)$$

(4)

where n is the number of data grids, and$X_j\mathrm{and}{X}_k$ are the data values in time series j and k (k > j), respectively. In cases where the sample size n > 10, the standard normal test statistic U is computed using Equation (5):

$$U=\left\{\begin{array}{cc}\frac{S-1}{\sqrt{V_{ar}}},& ifS>0\\ 0,& ifS=0\\ \frac{S+1}{\sqrt{V_{ar}}},& ifS<0\end{array}\right\}$$

(5)

In this study, we combined the Sen slope and Mann–Kendall tests to analyze the groundwater storage trend, and the confidence interval is computed at one significance level (α = 0.05), $U_{1-\alpha /2}$ = $U_{0.975}$ = 1.96. If β > 0 and |U| > 1.96, the groundwater storage has a significant increase; otherwise, the groundwater storage has a significant decrease.

3. Results

3.1. The Getis-Ord Gi* Analysis of Groundwater Storage Changes

Initially, we used Moran’s I to measure the spatial correlation of changes in groundwater storage. From 2003 to 2019, the p-values of all months were below 0.1, and the absolute values of Z were greater than 2.56, indicating that groundwater storage changes had significant spatial clustering in China.

In addition, the Getis-Ord Gi* method was used to verify the specific spatial location of high-value and low-value clusters (Figure 2). From 2003 to 2007, the area with an increasing trend in groundwater storage accounted for approximately 40.73%, and an area with a decreasing trend accounted for approximately 34.72%. Hot spots (high-value clustering areas) were mainly distributed in NWC, NC, and southern TP; that is, in these regions, changes in groundwater storage increased. However, during this period, cold spot areas (low-value clustering areas) were mainly distributed in NEC and MLY, indicating that in these two regions, the changes in groundwater storage mainly showed a decreasing trend from 2003 to 2007. During the period from 2008 to 2019, hot spots accounted for approximately 41.31%, and cold spots accounted for approximately 38.65%. The hot spots of groundwater storage changes were centered on the southern area of the Qinling-Huaihe line and the northern part of the TP, while the cold spots were mainly distributed in northern Xinjiang, LP, and NC. That is, during this period, the groundwater storage in the southern part of the Qinling–Huaihe line and northern areas of the TP showed an increasing trend, while the groundwater storage in the northern Xinjiang, LP, and NC regions decreased.

Overall, from 2003 to 2007, the groundwater storage in NWC, NC, and southern TP showed a positive trend, while that in NEC and MLY showed a negative trend. During 2008–2019, the groundwater storage in the southern part of the Qinling-Huaihe line and northern part of the TP showed a positive trend; the northern Xinjiang region, LP, and NC showed a negative trend.

3.2. Changes in Groundwater Storage in Nine Regions

We found that the groundwater storage changes in NC, LP, and NWC had a significant, negative trend (Figure 3a), which was the same as the research results in Feng [11], Xie [19], and Tangdamrongsub [17]. The variations in groundwater storage in NC were similar to that of LP, where anthropogenic factors were the main reason for the decrease in groundwater in these two regions [19,24], with mean rates of −10.09 ± 0.94 mm y⁻¹ and −10.05 ± 1.05 mm y⁻¹ equivalent height of water, respectively. The rate of depletion of groundwater was −4.91 ± 0.28 mm y⁻¹ equivalent height of water in NWC. When studying groundwater changes in the Hexi Corridor, Wang [20] pointed out that, compared with temperature and precipitation, human activity has a greater impact on the depletion of groundwater. The changing trend of groundwater storage on the TP and in NEC was smaller than that in the three regions mentioned above. Taking 2016 as the split point, during 2003–2015, there were fluctuating changes in groundwater in the two regions with mean rates of −0.2 ± 0.49 mm y⁻¹ and −3.26 ± 1.06 mm y⁻¹ equivalent height of water, respectively. After 2016, with mean rates of 3.28 ± 3.7 mm y⁻¹ and −21.9 ± 0.96 mm y⁻¹, respectively, the groundwater storage in the two regions decreased. Furthermore, in the SC, MLY, and Ch-Yu regions, the groundwater increased significantly, with mean rates of 7.26 ± 1.51, 7.73 ± 1.35, and 3.61 ± 0.53 mm y⁻¹ equivalent height of water, respectively. There was no obvious changing trend of groundwater storage in the Yun–Gui region where groundwater had an indistinctively positive or negative trend but fluctuated with a mean rate of 4.13 ± 0.68 mm y⁻¹ equivalent height of water. These results are consistent with the groundwater storage change trends calculated from China Water Resources Bulletin during 2003–2019, for example, which were −9.65 ± 1.94, −7.09 ± 0.8, 8.61 ± 2.1 and 7.53 ± 1.05 mm y⁻¹ in the NC, LP, SC, and MLY, respectively.

Using Sen slope combined with the Mann–Kendall nonparametric test, we can visually reflect the changes in groundwater in each area (Figure 3b). We found that groundwater on the LP, northern NWC, southern NEC, and southern TP had a negative trend in the range from −71.08 to −9.46 mm y⁻¹. Groundwater storage in most areas of Ch-Yu, SC, and Yun-Gui, and most areas of MLY and the northern area of TP, had a positive trend in the range from −0.87 to 18.9 mm y⁻¹, and other regions had no significant changes. From a national perspective, groundwater storage had a decreasing trend in northern China and an increasing trend in southern China, which were closely related to precipitation in the north and south and the development of industry and agriculture.

3.3. The Seasonal Changes in Groundwater Storage

We analyzed the changes in groundwater storage in all regions from a seasonal perspective (Figure 4 and

Table 1. Seasonal changing rate of groundwater storage in nine regions (mm y⁻¹).

	NC	SC	LP	MLY	NEC	NWC	TP	Ch Yu	Yun Gui
Spring	−11.34 ± 1.57	10.09 ± 1.86	−10.13 ± 1.39	9.21 ± 1.21	−8.05 ± 2.06	−5.23 ± 0.39	−7.62 ± 2.6	5.29 ± 1.87	8.31 ± 1.79
Summer	−9.71 ± 1.15	8.89 ± 1.32	−10.83 ± 1.07	12.77 ± 1.51	−7.95 ± 1.66	−5.56 ± 0.42	−5.34 ± 1.43	3.18 ± 0.85	8.5 ± 2.12
Autumn	−11.05 ± 1.04	8.99 ± 1.27	−11.41 ± 1.3	15.25 ± 1.74	−8.36 ± 2.33	−5.11 ± 0.36	−4.36 ± 1.04	6.95 ± 2.2	11.79 ± 2.69
Winter	−11.45 ± 0.96	10.22 ± 1.53	−10.11 ± 0.83	11.91 ± 1.42	−6.04 ± 1.35	−4.4 ± 0.47	−4.74 ± 2.03	5.22 ± 0.85	9.91 ± 1.99

). We found that in NC and LP, where the groundwater in the four seasons had a negative trend, NC had the maximum rates of groundwater depletion occurring in winter, which was −11.45 ± 0.96 mm y⁻¹ equivalent height of water, and LP had the maximum rate of decrease in autumn with a mean rate of −11.41 ± 1.3 mm y⁻¹. Seasonal changes in groundwater storage on the LP are attributed to changes in precipitation [19]. The seasonal changes in groundwater storage in the NEC and NWC regions had similar trends; in spring and summer, the extent of change in groundwater storage was greater than that in autumn and winter, which may have been caused by agricultural irrigation and changes in precipitation; groundwater storage in winter, spring, and summer remained basically unchanged, but there was a negative trend in autumn. NEC had the maximum rate of decrease in autumn, with a mean rate of −8.36 ± 2.33 mm y⁻¹ equivalent height of water, and the four seasons in NWC had little difference in value. Groundwater storage on the TP increased slightly in autumn and decreased in the other three seasons. The groundwater storage in MLY had a significant seasonal change, which was manifested as a significant increase in autumn with a mean rate of 15.25 ± 1.74 mm y⁻¹ equivalent height of water, a significant decrease in spring and winter, and a constant change in summer. In SC and Yun–Gui, the seasonal changes in groundwater storage had a significantly increasing trend and had the largest increase in winter and autumn with mean rates of 11.79 ± 2.69 mm y⁻¹ and 10.22 ± 1.53 mm y⁻¹, respectively. As the two regions are located in the monsoon climate zone, there are obvious seasonal changes in precipitation. Together, the seasonal changes in groundwater storage in the Ch-Yu region were represented by increases.

4. Conclusions

In this paper, we applied hot spot analysis, Sen slope, and Mann-Kendall nonparametric methods to analyze the spatiotemporal variation in groundwater in China based on GRACE, GRACE-FO, and GLDAS data. Our results showed that during 2003–2019 in NC, LP, and NWC, groundwater storage continued to decrease, with mean rates of −10.09 ± 0.94, −10.05 ± 1.05, and −4.91 ± 0.28 mm y⁻¹ equivalent height of water, respectively. In NC, maximum rates of groundwater depletion occurred in winter; on the LP, a significant decrease occurred in autumn, and there was little difference in the four seasons in NWC. Furthermore, the groundwater in SC, MLY, and Ch-Yu had a positive trend, with mean rates of 7.26 ± 1.51, 7.73 ± 1.35, and 3.61 ± 0.53 mm y⁻¹ equivalent height of water, respectively. There was a visible increase in groundwater in winter in the SC region and in autumn in MLY, then Ch-Yu had maximum rates of increase in autumn. Taking 2016 as the time demarcation point, the groundwater storage fluctuated slightly before 2016 on the TP and in NEC and decreased significantly after 2016. The groundwater storage in Yun-Gui showed a fluctuating trend. Overall, we found that the changes in groundwater storage in China were generally reduced in the north and increased in the south. This feature was closely related to the precipitation, industrial and agricultural development, and population distribution in the different regions. Our results can provide data support for others to grasp the temporal and spatial changes in groundwater storage in China and to carry out ecological environmental restoration.