Quantitative Assessment of Shallow Groundwater Sustainability in North China Plain

Abstract

The depletion of shallow groundwater has seriously affected the sustainable development of water resources in the North China Plain (NCP). Based on 556 well monitoring observations over a period of 13 years, we quantitatively evaluated the shallow groundwater sustainability in the NCP via various indices (e.g., the reliability, resilience, vulnerability, and sustainability indices), and further discussed the contribution of different drivers (including climatic and non-climatic factors). The main conclusions are summarized as follows: (1) the yearly trend of shallow groundwater shows a serious long-term deficit in the Piedmont Plain but is not significant in the East-Central Plain. (2) As for the sustainability of shallow groundwater in the NCP, the reliability is below the medium level (reliability < 0.5) in most areas and the ability of shallow aquifers to restore groundwater is very weak (resilience < 0.2), while the lack of groundwater storage in most shallow aquifers is not serious (vulnerability < 0.4). The final sustainability index (<0.1) shows the poor sustainability of most shallow aquifers in the NCP. (3) The non-climatic factor is the dominant driver of shallow groundwater depletion in the NCP when compared to the climatic factor. This result is helpful to formulate the water management policies for sustainable shallow groundwater storage in the NCP.

1. Introduction

The North China Plain (NCP) has become one of the most water-stressed regions in China due to population explosion, climate change, agricultural activities, and urbanization [1,2]. In recent years, groundwater resources have accounted for more than 70% of water resource utility in the most densely populated regions of China, especially in the NCP [3,4,5]. The contradiction between the supply and demand of water resources is very prominent. Unsustainable high-intensity groundwater extraction has resulted in a significant decrease in groundwater of 0.5–2.0 m/year [6,7]. As a result, a series of geological and environmental problems such as land subsidence, soil salinization, and aquifer depletion have become increasingly serious, greatly restricting the social and economic development of the NCP [6,8]. To better optimize the allocation of water resources in the NCP and alleviate the “groundwater crisis” [9], many scholars devote themselves to studying the main factors leading to groundwater depletion to further explore the sustainable utilization of groundwater resources in the NCP.

At present, most studies on groundwater storage deficiency in the NCP are based on GRACE observation data to assess the trend of the yearly decline of groundwater storage [10,11,12,13]. Moreover, the studies in which groundwater decline was typically attributed to continuous groundwater extraction [14,15,16], changes in precipitation, and hydrological drought [17,18,19] did not separate groundwater storage changes caused by groundwater extraction and precipitation. However, it has been observed that the water balance in Asia showed a stable trend during the moist climate regime from 2010 to 2013 [20]. This illustrates the need to isolate and assess the individual contributions of factors that influence groundwater losses. In other words, previous analyses of groundwater trends have not presented the actual pressure on shallow groundwater caused by climate change and increasing human demand. Thus, to better differentiate and quantify the variation in precipitation and the contribution of human-induced influences to shallow groundwater, we used the in situ monitoring well observations and the Global Land Data Assimilation System (GLDAS) model from 2005–2018 to assess the impact of human demands and changing precipitation patterns on shallow groundwater storage variations in the NCP to complement the existing work.

In addition, the current research on the sustainable development of groundwater in NCP. On the one hand, the sustainable utilization and evaluation of groundwater resources in some areas of the NCP (such as Beijing, Hebei, Shandong, Henan, etc.) lack overall and systematic analysis [21]. On the other hand, the sustainable analysis of the overall groundwater resources in the NCP based on a single index only considers the impact of groundwater extraction on the sustainable utilization of groundwater [22], but does not comprehensively consider the impact of aquifer structure and climate change on groundwater change; it is easy to be restricted by a single index and the evaluation is limited. Therefore, some studies choose to use the random performance index of water resource systems to evaluate the sustainability of groundwater. After considering the water storage capacity of the aquifer and the serious situation of groundwater deficiency, the case of aquifer resilience achieves an easily quantified multi-angle sustainability assessment [23,24]. Subsequent studies have added the groundwater deficiency index (GDI index) based on this method to obtain more standard quantitative results of water storage [25] and to conduct a more accurate sustainability evaluation, which is of great significance for the sustainable utilization of subsequent groundwater in the study area.

In conclusion, under the comprehensive consideration of aquifer structure, human demand, and the influence of climate change, we chose to use the random performance index of the water resource system and groundwater deficiency index (GDI) to analyze the reliability, elasticity, and vulnerability index (RRV) of shallow wells in the NCP to evaluate shallow groundwater depletion from multiple perspectives. On this basis, we quantitatively evaluate the sustainable development of shallow groundwater from the perspective of space.

2. Data and Methods

2.1. Data

2.1.1. Groundwater Storage Estimation from Well Observation

Water resources in China are not evenly distributed from north to south. About half of the population and two-thirds of the farmland are located in northern China, while the water resources there only account for one-fifth of the national total [26,27]. The NCP (Figure 1) covers an area of 140,000 sq km, has a population of 350 million, and is responsible for 10% of China’s food production; however, it is short on water resources per capita and relies heavily on groundwater for human water supply [28,29].

The NCP is mostly distributed with loose rock pore groundwater. The lithology of the aquifer changes from gravel, medium-coarse sand, and medium-fine sand in the Piedmont Plain to medium-fine sand and fine sand in the central plain, and gradually to fine sand and silty sand in the Marine plain. The shallow aquifer thickness transients from 120~210 m in the Piedmont Plain to 50~120 m in the East-Central Plain [23]. The thickness of the deep aquifer gradually increases from 100 m in the Piedmont Plain to about 600 m in the East-Central plain, but the average storage coefficient is 0.00125, about 2% of the area-weighted ratio yield of the shallow aquifer (0.075) [30], so the release of deep groundwater is much lower than that of the shallow aquifer. This is neglected in the calculation of water storage change [31]. The high heterogeneity of the shallow aquifer system indicates the need for a more detailed assessment of the response of different shallow aquifers to individual components of groundwater loss. At the same time, the high heterogeneity of aquifer systems suggests a need for a more detailed assessment of the response of different aquifers to individual components of groundwater depletion. Therefore, shallow aquifers with more wells and wider coverage were selected for detailed analysis.

The monthly groundwater level measurements from 556 unconfined wells during 2005–2018 were obtained from the Chinese groundwater management authority and the groundwater yearbook. The locations of 556 monitoring wells are presented in Figure 1. Due to the uneven spatial distribution of monitoring wells, the regional mean groundwater level is estimated by the area-weighted Thiessen polygon method. The Thiessen polygons for these unconfined wells in the NCP are presented in Figure 2.

The monthly groundwater storage is calculated as the multiplication of the groundwater level measurement, the area of each Thiessen polygon, and the corresponding specific yield [10]. The range of specific yield in the NCP is 0.025~0.29, and the spatial data of comprehensive specific yield were obtained from Cao et al. (2013) [32].

2.1.2. Precipitation Data

Monthly precipitation data from 2005 to 2018 were derived from the catchment land surface model (CLSM) product of the Global Land Data Assimilation System (GLDAS) released jointly by NASA and the US National Oceanic and Atmospheric Administration. The product provides estimates of shallow groundwater storage. However, the results only represent natural variations in groundwater storage [33,34]. In this study, the CLSM model is also exploited to estimate the non-climatic portion of shallow groundwater storage in the NCP.

2.2. Groundwater Performance Indicators

RRV indicators (reliability, resilience, and vulnerability) were used to assess the overall situation of shallow groundwater storage in the NCP. The RRV index was initially used to evaluate the random performance of water resource systems [34], which was convenient for the comprehensive analysis of water resource systems from multiple perspectives. In order to achieve RRV analysis of shallow groundwater systems in the NCP, we used the GDI index to calculate the groundwater performance indicators [25]. The first step of GDI index calculation is to eliminate the seasonality of groundwater storage anomaly (${\mathrm{GWSA}}_{\mathrm{i}}$) by deducting the seasonal mean (${\mathrm{CM}}_{\mathrm{i}}$) [35].

$$C{M}_i=\frac{{{\displaystyle \sum }}_{j=1}^{n_i}GWS{A}_j}{n_i}$$

(1)

where i represents the month of observation and ${\mathrm{n}}_{\mathrm{i}}$ represents the number of observations. The abnormal groundwater storage after removing ${ \mathrm{CM}}_{\mathrm{i}}$ represents the net deviation of groundwater storage change, namely groundwater storage deviation (GSD) [36]. GDI is normalized by using the standard scoring method with mean (${\mathsf{\mu }}_{\mathrm{GSD}}$) and variance (${\mathsf{\sigma }}_{\mathrm{GSD}}$), as shown in the following equation:

$$GDI=\frac{GS{D}_t-{\mu }_{GSD}}{{\sigma }_{GSD}}$$

(2)

GDI indices calculated from monitoring wells form the basis of various statistical indicators [37,38], and a positive GDI demonstrates a satisfactory state, while a negative GDI denotes an unsatisfactory state [25].

2.2.1. Reliability

Reliability (Rel) represents the probability that the water resource system is in a satisfactory state [24], while in the groundwater resource system, it represents the possibility that aquifer storage is higher than a certain threshold (the water storage does not increase or decrease corresponding to zero GDI index: greater than zero indicates abundant groundwater; a value less than zero indicates a lack of groundwater). The Rel value is calculated as the quotient of the number of satisfactory scenarios divided by the total number of scenarios [37,38]. $\mathrm{Rel}$ closer to 1 represents more satisfactory events, indicating stable and reliable groundwater storage in the region. In contrast, $\mathrm{Rel}$ closer to 0 represents more unsatisfactory events, indicating that the reliability of groundwater is poor in this region.

2.2.2. Resilience

Resilience (Res) represents the ability of a water resource system to overcome an unsatisfactory state [24], while in the groundwater resource system, it represents the likelihood that aquifers will return to normal values after overuse of groundwater or prolonged groundwater drought. Resilience is quantified as the quotient of the number of satisfactory scenarios following unsatisfactory scenarios divided by the number of unsatisfactory scenarios [37,38]. The value ranges from 0 to 1, and closer the Res is to 1, the better is the ability of the aquifer system to recover groundwater storage.

2.2.3. Vulnerability

Vulnerability (Vul) represents the occurrence possibility of water resource system failure (unsatisfactory state) events [24]. In this paper, Vul is the sum of the GDI value corresponding to the unsatisfactory events multiplied by the probability (${\mathrm{p}}_{\mathrm{j}}$) of the unsatisfactory events occurring during the study period:

$$Vul=\sum {GDI}_{j\times p_j}$$

(3)

The closer the Vul is to 1, the more serious the loss of groundwater reserves.

2.2.4. Sustainability Index

The sustainability coefficient (Si) describes each region as a function of Rel, Res, and Vul. In this study, the product of three performance indicators was used to form the sustainable development coefficient. Compared with the equal-weight summative coefficient of the three performance indicators [25], it has the advantage of increasing the weight of the statistical measure with small values and producing a higher index value only when all statistical measurements are large.

$$Si=Rel\times Res\times \left(1-Vul\right)$$

(4)

2.3. Multiple Linear Regression

2.3.1. Functional Model

Since groundwater is an important component of global fresh water demand, it is critical to distinguish the drivers that cause significant changes in groundwater storage. In this study, the least squares multiple regression (MLR) method with multiple datasets was used to quantify the impacts of climatic and non-climatic factors on shallow groundwater storage variations. Meanwhile, the MLR-based approach paves the way for dominance analysis, which helps to distinguish the dominant drivers of groundwater storage variations in NCP. The proposed groundwater functional model has a general form:

$$SG{W}_{At}={\beta }_0+{\beta }_1{P}_{1t}+{\beta }_2{P}_{4t}+{\beta }_3{P}_{12t}+{\beta }_{4}N{P}_t$$

(5)

where ${\mathrm{SGW}}_{\mathrm{At}}$ is the measured monthly change in shallow groundwater storage, ${\mathsf{\beta }}_0$ to ${\mathsf{\beta }}_4$ is the model coefficient, ${\mathrm{P}}_{1\mathrm{t}}$ is the precipitation of the previous month at observation time t, ${\mathrm{P}}_{4\mathrm{t}}$ is the total precipitation in 4 months before observation time t, ${ \mathrm{P}}_{12\mathrm{t}}$ is the total precipitation in 12 months before observation time t, and ${\mathrm{NP}}_{\mathrm{t}}$ is the influence of non-climatic factors.

Since the dominant drivers of groundwater variations in NCP are human activities and precipitation change [8], we divided the simulated shallow groundwater storage change into climatic (${\mathrm{P}}_{1\mathrm{t}}$, ${\mathrm{P} }_{4\mathrm{t}}$, ${\mathrm{P}}_{12\mathrm{t}}$) and non-climatic (${\mathrm{NP}}_{\mathrm{t}}$) in Equation (5). In the regression model, precipitation (climatic) components are considered on different time scales (one month, four months, and twelve months) to reflect the time delay of natural recharge in the precipitation process.

The intensity of precipitation plays an important role in maintaining groundwater balance. Due to the high heterogeneity of the aquifer system in the NCP, the recharge effect of precipitation also varies with different regions. From the Piedmont Alluvial Plain, Central Alluvial Lacustrine Plain, and Coastal Alluvial Plain, the thickness of the aquifer becomes thinner, and the permeability becomes lower, and the effect of the aquifer receiving precipitation is different. In the NCP, high-intensity precipitation mainly occurs in summer (June to September), which is conducive to rapid recharge of the Piedmont Aquifer with high permeability. The low-intensity precipitation is favorable for the slow recharge of the low permeability aquifer in the middle and east. Therefore, considering the change in recharge time required by different aquifers, ${\mathrm{P}}_{1\mathrm{t}}$, ${\mathrm{P}}_{4\mathrm{t}}$ and ${\mathrm{P}}_{12\mathrm{t}}$ are considered in Equation (5).

The first part (${\mathrm{P}}_{1\mathrm{t}}$) represents the shallow groundwater storage variations caused by precipitation in the previous 1 month of observation time T. The second part (${\mathrm{P}}_{4\mathrm{t}}$) represents the shallow groundwater storage variations caused by precipitation in the first 4 months of observation time T. As mentioned above, precipitation in North China is mainly seasonal precipitation from June to September. In order to capture the changes in shallow groundwater storage caused by seasonal high-intensity precipitation, we carefully selected the total precipitation from June to September each year. This part is especially important for aquifers with high permeability and fast recharge rate. Selecting other months will not capture the effect of seasonal precipitation on changes in shallow groundwater storage. The third part (${\mathrm{P}}_{12\mathrm{t}}$) represents the change in shallow groundwater storage caused by precipitation in the previous 12 months of observation time T. This part is especially important for aquifers with low permeability and slow recharge rate.

The non-climatic component (${\mathrm{NP}}_{\mathrm{t}}$) was estimated by subtracting the shallow groundwater storage from the CLSM surface model. Since the CLSM estimation only reflects the changes in shallow groundwater storage caused by climate change, the non-climatic component can be obtained by deducting this estimate, which has also been applied in previous studies such as [35].

2.3.2. Dominance Analysis Using Proportional Reduction of Error

To further identify the dominant factors leading to changes in shallow groundwater storage, we used proportional reduction of error (ProRE) to perform dominance analysis of the predictors. Dominance analysis is a tool for calculating the relative contribution of each predictor in multiple regression models [39]. The dominance analysis was performed via different models with a single unit of measured observations, and the ProRE of each component was obtained by removing the corresponding components. For example, the ProRE of ${\mathrm{P}}_{1\mathrm{t}}$ is calculated as

$$ProR{E}_{P_{1t}}=\frac{\left(RSS of G{W}_{At2}-RSS of G{W}_{At1}\right) }{RSS of G{W}_{At2}}$$

(6)

where RSS denotes the residual sum of squares of a time series. ${\mathrm{GW}}_{\mathrm{At}1}$ is the model containing all components as

$$G{W}_{At1}={\beta }_0+{\beta }_1{P}_{1t}+{\beta }_2{P}_{4t}+{\beta }_3{P}_{12t}+{\beta }_{4}N{P}_t$$

(7)

In addition, ${\mathrm{GW}}_{\mathrm{At}2}$ is the model excluding ${\mathrm{P}}_{1\mathrm{t}}$:

$$G{W}_{At2}={\beta }_0+{\beta }_2{P}_{4t}+{\beta }_3{P}_{12t}+{\beta }_{4}N{P}_t$$

(8)

3. Results and Discussions

3.1. State of Shallow Groundwater Storage

To preliminarily access the state of shallow groundwater storage in the NCP, we calculated the yearly trends of well monitoring (from 2005 to 2018) observations. It should be noted that the positive trend in groundwater is shown in shades of red, while the negative trend is represented in shades of blue. As shown in Figure 3, the shallow groundwater in the Piedmont Plain region adjacent to the Taihang Mountains shows an obvious blue color, which means an obvious deficit, indicating that the Piedmont Plain region had a long-term deficit from 2005 to 2018. In the East-Central Plain, the trend is close to zero, and the deficit is not obvious, meaning that the groundwater deficit is not serious. This is consistent with the reality of shallow groundwater exploitation in the NCP: agricultural irrigation is intensive in the Piedmont Plain, and the mining conditions of shallow aquifers are excellent, so a large amount of groundwater is extracted for agricultural irrigation, which is the main exploitation area of shallow groundwater, and the shallow groundwater is severely deficient. In contrast, the shallow groundwater in the East-Central Plain is unusable salt water [40], thus the shallow groundwater depletion in the central and eastern regions is not as severe.

3.2. Groundwater Performance Indicators

To assess the status quo of shallow groundwater from multiple perspectives, we calculated the performance indicators (reliability (Rel), resilience (Res), vulnerability (Vul), and sustainability index (Si)) for individual monitoring of wells based on the GDI index. The values of these indices are normalized based on their minimum and maximum values to compare the results on the same scale from 0 to 1. Figure 4a–d show the results of Rel, Res, Vul, and Si obtained from in situ well observations. Areas representing high values (close to 1) of Rel indicate the large storage of groundwater. The high values (close to 1) of Res indicate the better ability of the aquifer system to recover groundwater storage. The high values (close to 1) of Vul indicate the more serious loss of groundwater storage. The high values (close to 1) of Si indicate greater sustainability.

Figure 4a shows that the reliability is below the medium level (<0.5) in most areas of the NCP, which indicates that the shallow groundwater storage in most areas of the NCP is relatively minor after 13 years of extraction. At the same time, resilience analysis shows that most shallow aquifers are severely damaged after long-term groundwater extraction [8,41], and the ability of shallow aquifers to restore groundwater storage is extremely weak (<0.2) (Figure 4b). In Figure 4c, similar Vul calculations (nearly between 0.3 and 0.4) suggest that the lack of groundwater storage in most shallow aquifers is not significant in the NCP (<0.4), which is attributed to the fact that after realizing the importance of sustainable development of groundwater, the government adopted water transfer and water-saving measures to supplement groundwater storage in the later period [42]. Three indicators all affect the sustainable development of shallow groundwater storage in the NCP. Finally, the value of Si is less than 0.1, indicating that most shallow aquifers in the current NCP have poor sustainability and need to reduce extraction and increase water sources to promote further groundwater recovery. It is worth noting that there is a strong correlation (0.86) between resilience (Res) and sustainability (Si), which indicates that the ability of shallow aquifers to restore groundwater is the main factor affecting the sustainable development of shallow groundwater in the NCP. This means that irreversible damage to the aquifer structure, resulting in a continued deterioration of the aquifer’s ability to restore groundwater storage, must be avoided to better promote sustainable groundwater development.

Compared with previous studies, we have longer and more densely distributed in situ well observations. This paper provides spatial quantification of sustainable groundwater development over a longer time, from more angles, and with denser coverage points.

3.3. Influence of Different Drivers on Shallow Groundwater Variability

In order to determine the influence of climatic and non-climatic factors on shallow groundwater storage variations, we separated the main components through multiple linear regression simulation via Equation (5).

To assess the reliability of our multiple linear regression process, the root mean square errors (RMSEs) between the simulated results and the measured data were calculated. As shown in Figure 5, the RMSEs are all smaller than 0.3 cm, indicating the efficiency of the regression model in the whole study region.

The model coefficient results obtained by multiple linear regression are presented in Figure 6. The coefficients are very similar over the East-Central Plain but distinguishable over the Piedmont Plain. According to the distribution of coefficients, the climatic and non-climatic impacts on groundwater storage for the East-Central Plain are not clearly distinguished. However, for the Piedmont Plain, the coefficient results show that there is a negative correlation between precipitation and groundwater change in 1 and 12 months, and a positive correlation between non-climatic factors.

After that, we performed ProRE analysis on the fitting results to obtain the dominant drivers of shallow groundwater storage variations (Figure 7). The closer the area color is to red, the larger the driver is, and the more likely it is to affect the changes in shallow groundwater. As shown in Figure 7a, the non-climatic factor accounts for the largest proportion and plays a dominant role in the change in shallow groundwater storage, that is, due to the long-term and intense artificial exploitation of shallow groundwater for agricultural irrigation and daily water supply in the NCP. In contrast, the shallow groundwater is not sensitive to precipitation accumulated with different time delays (Figure 7b–d). Among the climate (precipitation) factors, shallow groundwater in the NCP is more sensitive to accumulated precipitation in the first four months.

In order to evaluate the variation in precipitation, a trend analysis was also conducted for precipitation at different time accumulative scales. Figure 8a–c show the accumulative precipitation trend one month, four months and twelve months before the observation time, respectively. The results show that precipitation intensity in the NCP presents a generally positive trend, which is not consistent with the groundwater storage variations derived from well observations (Figure 3). It also indicates that precipitation (climatic factor) is not the main driver for the long-term decline in shallow groundwater storage in the NCP.

To further investigate this issue, we also assessed the trend of shallow groundwater storage in the CLSM. Figure 8d shows the trend of shallow groundwater storage estimated from the CLSM (climatic components). It was found that climate-driven shallow groundwater changes very little (note the different color bars in Figure 8d), and the Piedmont Plain shows an obvious positive trend, which also presents different features compared to the trends presented in Figure 3. This again confirms the minor contribution of the climatic component.

Instead, as compared to the climatic composition, these results present obvious consistency between trends of non-climatic components (Figure 8e) and those derived from 556 well observations (Figure 3). For instance, there are negative trends over the Piedmont Plain of the NCP, that is, there are long-term losses in shallow groundwater, while the losses in the central and eastern regions are not serious. The correlation coefficients between the trend of climatic components (Figure 8d) or non-climatic components (Figure 8e) and total groundwater storage (Figure 3) are presented in

Table 1. Correlation analysis for climatic component and non-climatic component.

Correlation Index	Climatic Component	Non-Climatic Component
Shallow groundwater storage	−0.36	0.99

. The quantitative results also support that the shallow groundwater storage deficit in the NCP is mainly derived from non-climatic factors.

For the study of groundwater driving factors in the NCP, the advantages of this paper are as follows: more comprehensive measured data, the realization of a long-term point-to-point analysis, and more accurate results. In addition, instead of attributing the cause of the groundwater deficit to the combined influence of man-made extraction and rainfall, this paper conducts an isolated analysis of the driving factors causing groundwater change, obtains independent non-climatic and climatic components of shallow groundwater, and further separates the climatic components of different time scales considering the delayed effect of precipitation on groundwater recharge. This was performed in more detail compared to other papers that directly deduct the climatic component estimated by the CLSM from the overall groundwater (deep + shallow) to obtain the non-climatic component of groundwater [25]. We consider that the CLSM estimates represent only the climatic component of shallow groundwater [33,34]. Therefore, the interference of deep groundwater is eliminated, and accurate results of shallow groundwater drivers are obtained in a more detailed stratification.

In addition, since the specific yield is needed in the calculation process of in situ wells, but the aquifer structure changes due to long-term groundwater extraction, its groundwater storage capacity changes correspondingly, but the specific yield data are not updated accordingly. Is there a higher real-time specific yield to obtain more accurate groundwater data? Secondly, the time range of the existing in situ data is only from 2005 to 2018, which brings the question of whether longer data can be obtained. Moreover, this paper only analyzes the changes in shallow groundwater, so whether it is possible to use additional data sources to conduct a detailed exploration of the overall groundwater in the NCP remains to be seen. Finally, the number of wells obtained in the present paper is the most intensive under the 13-year time scale, but that is only about 39 wells per 10,000 square km, which is still not enough for the NCP, and a more detailed analysis is flawed. We also hope to have more detailed data to support further research.

4. Conclusions

Based on the in situ well observations, we analyzed the spatial distribution, the sustainability, and dominant drivers of shallow groundwater variations in the NCP. The results demonstrate the following:

(1)

There is an obvious long-term shallow groundwater deficit in the Piedmont Plain, while the deficit is not significant in the East-Central Plain.

(2)

The reliability is below the medium level (reliability < 0.5) in most areas of the NCP, and there is very low ability (resilience < 0.2) of most shallow aquifers to recuperate from the decline in groundwater storage. Meanwhile, the deficit of groundwater storage in most shallow aquifers is not serious (vulnerability < 0.4) in the NCP. The low reliability and resilience values in most areas of the NCP result in the weak sustainable development ability (sustainable index < 0.1) of shallow groundwater storage.

(3)

The multiple linear regression model was used to separate the impacts of precipitation and non-precipitation components, revealing the sensitivity of climatic and non-climatic drivers of shallow groundwater in the NCP. The results demonstrate that the non-climatic factor is the dominant driver of shallow groundwater storage depletion in NCP.

Our study provides a quantitative analysis of shallow groundwater status and sustainability in the NCP, and comprehensively assesses the impacts of climatic and non-climatic factors on shallow groundwater decline in the NCP. It will be of high relevance to governing bodies to implement water management policies for sustainable shallow groundwater storage. Since we focus on shallow groundwater variations but not whole groundwater variations in the NCP, the GRACE and GRACE Follow-On observations are not considered in this work. Considering the good performance of satellite gravimetry in tracking groundwater in the NCP, we will further distinguish between the shallow groundwater and confined water (deep groundwater) variations in our upcoming work.