Comparing Groundwater Storage Changes in Two Main Grain Producing Areas in China: Implications for Sustainable Agricultural Water Resources Management

Abstract

Balancing groundwater supply and food production is challenging, especially in large regions where there is often insufficient information on the groundwater budget, such as in the North China Plain (NCP) and the Northeast China Plain (NECP), which are major food producing areas in China. This study aimed to understand this process in a simple but efficient way by using Gravity Recovery and Climate Experiment (GRACE) data, and it focused on historical and projected groundwater storage (GWS) changes in response to changes in grain-sown areas. The results showed that during 2003–2016, the GWS was depleted in the NCP at a rate of −17.2 ± 0.8 mm/yr despite a decrease in groundwater abstraction along with an increase in food production and a stable sown area, while in the NECP, the GWS increased by 2.3 ± 0.7 mm/yr and the groundwater abstraction, food production and the sown area also increased. The scenario simulation using GRACE-derived GWS anomalies during 2003–2016 as the baseline showed that the GWS changes in the NCP can be balanced (i.e., no decreasing trend in storage) by reducing the area of winter wheat and maize by 1.31 × 10⁶ ha and 3.21 × 10⁶ ha, respectively, or by reducing both by 0.93 × 10⁶ ha. In the NECP, the groundwater can sustain an additional area of 0.62 × 10⁶ ha of maize without a decrease in storage. The results also revealed that the current groundwater management policies cannot facilitate the recovery of the GWS in the NCP unless the sown ratio of drought-resistance wheat is increased from 90% to 95%. This study highlights the effectiveness of using GRACE to understanding the nexus between groundwater supply and food production at large scales.

1. Introduction

Water shortages and food security are important issues in the world. As the largest source of freshwater, groundwater is critically important for irrigated agriculture, and hence, for global food security [1]. As irrigated agriculture constitutes approximately 40% of global food production [2], increased groundwater abstraction for irrigation has contributed substantially to the expansion in the global food supply [3]. However, this process has also led to groundwater depletion in many parts of the world, such as northwestern India, the North China Plain (NCP), and the central United States and California [4].

The NCP is one of the most important grain production regions in China [5]. However, it is also a region that suffers from very serious groundwater overexploitation, which is mainly due to agricultural irrigation [6,7,8]. The groundwater in NCP has been overexploited since the 1970s [9]. Previous studies have shown that one of the most effective ways to reduce the amount of groundwater exploitation in the NCP is by implementing seasonal fallow and appropriately reducing the sown area of winter wheat [5,6,10]. It was reported that the net groundwater consumption in the Hebei Plain (part of the NCP) for grain production (mainly winter wheat) was 113 km³ during the past 25 years [11]. A previous study showed that the optimal production of winter wheat could be guaranteed by 300 mm of irrigation [6]. Therefore, there is an urgent need to develop optimal irrigation water management to avoid further overexploitation of groundwater and to maintain sustainable crop production. Wu et al. found that the decline rate of the groundwater table had a decreasing trend before and after 2001, as did the planting area of winter wheat in the Northern Hebei Plain and Southern Hebei Plain [12]. From 1998 to 2010, the sown area of winter wheat in the Hebei Plain decreased by 49.62 × 10⁴ ha, and 15.96 × 10⁸ m³ of irrigation water was saved [13]. In order to protect groundwater resources and reduce groundwater exploitation, the sown area of winter wheat in the High Plains of the United States decreased by 36% from 1980 to 2010 [14]. In addition, several studies on supplemental irrigation and limited or deficit irrigation have been conducted to improve grain yield and increase irrigation water use efficiency (WUE) [5,6,15]. Additionally, Zhang et al. found that significant improvements in WUE have increased the crop yield in the NCP by 50% in the past 20 years [16].

The NCP and the Northeast China Plain (NECP) (Figure 1) are major food production regions in China [7,8,17]. The NCP (140,000 km²) has a typical temperate and monsoonal climate with annual precipitation of 500–1000 mm. Groundwater accounts for ~75% of the total water consumption [18]. The cropping system features two harvests per year, mostly with a wheat–corn rotation [17]. According to the Haihe River Basin Water Resources Bulletin and the statistical yearbooks of Beijing, Tianjin, Hebei, Henan and Shandong, the average groundwater abstraction in the NCP during 2003–2016 was 131 mm/yr, of which more than 70% was used for agricultural (mainly winter wheat) irrigation [5,6,19,20,21]. However, the annual average groundwater recharge is only about 100–120 mm/yr [19,22,23]. Since the amount of groundwater that is exploited is larger than the recharge, the GWS is continuously depleted. The NECP (350,000 km²) encompasses four provinces (Liaoning, Jilin, Heilongjiang and Inner Mongolia) and is characterized by a continental and monsoon climate with annual precipitation of 350–700 mm. The planting system is one crop-harvest per year, with wheat, corn, and rice being the major crops. In the NECP, groundwater irrigation accounts for ~47% of the total water consumption [24]. According to the Songliao River Basin Water Resources Bulletin, the average groundwater exploitation from 2003 to 2016 was 83 mm/yr, with ~40% used for agricultural irrigation, while the average annual groundwater recharge is ~120 mm/yr [24].

In recent decades, the conflict between agricultural water demand and groundwater sustainability have been widely observed. In addition to the widely reported groundwater depletion in the NCP [7,8,25,26], some regions of the NECP had also suffered from groundwater depletion [27,28]. With regard to the current situation in the NCP and NECP, it is important that we balance the GWS in the NCP and avoid serious GWS depletion in the NECP.

The monitoring and analysis of groundwater storage (GWS) changes are of great importance for sustainable use and effective management of water resources. In situ well measurement is a traditional method for monitoring GWS. However, estimating the GWS is a great challenge because of the complexity of the aquifer system and the uneven distribution of in situ monitoring wells, especially in some areas where there are insufficient monitoring wells. Moreover, water withdrawals from pumping wells are often unrestricted and unmonitored, which makes it even more difficult to estimate rates of groundwater consumption [29]. The Gravity Recovery and Climate Experiment (GRACE) satellite mission, launched in 2002, provides an excellent opportunity to monitor total water storage anomalies (TWSA) [30,31]. Wahr et al. provided the basic theory and method for estimating TWSA using synthetic GRACE data [32]. Swenson and Wahr improved the methods for extracting regional TWSA information from GRACE gravity coefficients [33]. By subtracting independent estimates of the water storage anomalies in the soil, snow and surface reservoirs from the GRACE-observed TWSA, the GWS anomaly (GWSA) can be estimated [34,35,36,37]. GRACE data has been successfully used to quantify GWS depletion in different irrigated regions worldwide, such as in northwestern India [36,38,39], the High Plains aquifer [37] and Central Valley [29,37,40,41] in the United States, and the NCP [7,8,25,26,42].

In this study, we used the GRACE data and groundwater level (GWL) measurements to estimate GWSA (2003–2016) in the NCP and NECP, and to investigate the relationship between agricultural water demand (i.e., groundwater supply) and food production with a focus on the appropriate sown areas. Based on the GRACE-derived GWS changes, the GWS response to changes in sown areas of grain crops was simulated using a linear relationship between GWS changes and grain-sown areas. The historical and projected GWS changes in the NCP and NECP were compared and discussed with a focus on how to balance agriculture water demand and groundwater sustainability at large scales with the assistance of GRACE data.

2. Data and Methods

The GRACE-derived TWSA includes several components of mass change, including snow water equivalent storage anomalies (SWESA), surface water storage anomalies (SWSA), soil moisture storage anomalies (SMSA), and GWSA. Using in situ or modeled estimates of SMSA, SWESA and SWSA, the GWSA can be calculated as a residual of the disaggregation equation [34,38,43]:

$$\mathrm{GWSA}=\mathrm{TWSA}-\mathrm{SMSA}-\mathrm{SWESA}-\mathrm{SWSA}$$

(1)

This study uses in situ GWL anomalies (GWLAs) to validate GRACE-derived GWSA.

2.1. GRACE-Derived TWSA

This study used the monthly GRACE mascons (i.e., mass concentration) products (2003–2016) from the Center for Space Research (CSR) (http://www.csr.utexas.edu/grace/) [44,45] at the University of Texas, NASA’s Goddard Space Flight Center (GSFC) (https://earth.gsfc.nasa.gov/geo/) [46] and the Jet Propulsion Laboratory (JPL) (http://grace.jpl.nasa.gov) [47], United States of America. The data processing of all the above mascon products included the following: the degree-2 order-0 (C20) coefficients were replaced with the solutions from satellite laser ranging [48], the degree-1 coefficients (Geocentre) were estimated using the method from Swenson et al. [49], and a glacial isostatic adjustment (GIA) correction was applied based on the ICE6G model [50].

Although the mascon solutions greatly reduced leakage errors from land to ocean [44,45,46,47], the gain factors provided by JPL were applied to the JPL RL06M data to account for leakage errors. The uncertainties in the GRACE data were quantified by propagating the errors associated with each GRACE product. The errors in JPL and GFSC mascons were officially provided along with the TWSA data. The error in CSR mascon was assumed to be 2 cm as suggested by CSR. The missing values in the GRACE data were conservatively estimated through cubic-spline interpolation.

Figure 2 shows the time series of the TWSA from 2003 to 2016 in the NCP (Figure 2a) and NECP (Figure 2b). A comparison of the monthly time-series of TWSA derived from the three mascon solutions indicates high consistency, while notable deviations still exist, especially for the period of 2014–2015 in the NCP. This might have resulted from a combination of (1) poor data quality (i.e., an increase in missing data) since 2011 that led to increased errors, and (2) extreme drought/wet events (e.g., 2014–2015 drought) that exaggerate the deviations. The average of the three products is used in this study to reduce the uncertainty [51].

2.2. SMSA and SWESA

Soil moisture and snow cover are important components of the TWSA and can be obtained through in situ measurements or model simulation. Previous studies have shown that the models (e.g., Community Land Model (CLM) and Mosaic, Noah and Variable Infiltration Capacity (VIC) models) from the Global Land Data Assimilation System (GLDAS) have relatively accurate simulations of SMSA and SWESA [52,53]. GLDAS was developed by NASA’s GSFC and National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Prediction (NCEP). Its goal is to generate optimal fields of land surface states and fluxes through ingesting satellite and ground-based observational data products using advanced land surface modeling and data assimilation techniques [54].

This study used the average SMSA and SWESA simulated by the four land surface models (CLM, Noah, Mosaic and VIC) from GLDAS. Furthermore, the modeled SMSA data from GLDAS were validated against in situ observations provided by the China Meteorological Administration, including 97 monitoring sites in the NCP and 45 monitoring sites in the NECP (Figure 1). The soil depth in the CLM, Noah, Mosaic and VIC is 3.43 m, 2.0 m, 3.5 m, and 1.9 m, respectively, while the soil depth of in situ observation is 0.75 m. Figure 3 shows that the average values of the four GLDAS models are basically consistent with the in situ data. The correlation coefficients of SMSA from in situ and GLDAS-AVE in the NCP and NECP were 0.59 and 0.72, and the RMSE values were 28.43 mm and 19.82 mm, respectively. For the SWESA, we also use the average values of the four GLDAS models (Figure 4), which can reduce uncertainty.

2.3. In Situ SWSA

Surface water mainly includes water in rivers, lakes, reservoirs and other water bodies. In this study, the SWSA refers to the changes in water stored in large and medium-sized reservoirs. We collected the reservoir storage data from the Haihe River Water Resources Bulletin, the Songliao River Water Resources Bulletin and the Hydrological Information Annual Report of China.

2.4. In Situ GWL Observations

To verify the accuracy of GRACE-derived GWSA, monthly in situ GWL data (Figure 1) from 844 wells in the NCP and 651 wells in the NECP during 2003–2016 were obtained from the Ministry of Water Resources of China. Also, monthly in situ GWL data from 81 wells in the NCP and 130 wells in the NECP during 2005–2016 were obtained from the GWL Almanac compiled by the China Institute of Geological Environment Monitoring. The GWLA and the mean series in the NCP and NECP are shown in Figure 5.

In situ GWSA is calculated by multiplying in situ-measured GWLA with a specific yield ($S_y$) of each well:

$$\mathrm{GWSA}=S_y\times \mathrm{GWLA}$$

(2)

In this study, the spatial distribution map of the $S_y$ in the NCP was provided by the China Geological Survey Bureau, and the $S_y$ ranges from 0.025 to 0.29 [26]. In the NECP, a mean $S_y$ value of 0.095, determined by pumping tests, was used to convert GWLA to GWSA [27,55].

2.5. Sown Area of Main Grain Crops

The main grain crops in the NCP and NECP are maize, wheat and rice. The sown area data for the main grain crops in the NCP were collected from the statistical yearbooks (2003–2016) for Beijing, Tianjin, Hebei, Henan and Shandong provinces. The sown area data for the main grain crops in the NECP were collected from the statistical yearbooks (2003–2016) for Heilongjiang, Jilin, Liaoning and Inner Mongolia provinces. Figure 6a shows that the main grain crops in the NCP are wheat and maize, and the total area of grain crops has remained basically unchanged since 2003. The main grain crops in the NECP are rice and maize, and the sown area of maize has increased significantly since 2003 (Figure 6b).

2.6. Auxiliary Data

Groundwater abstraction and precipitation are closely related to GWSA and food production, and were used for joint analysis. Annual groundwater abstraction data (2003–2016) were collected from the Water Resources Yearbook of the Haihe River basin and Songliao basin. Monthly gridded (0.5° × 0.5°) precipitation data (version 2.0) were obtained from the China Meteorological Administration (http://data.cma.cn).

2.7. Irrigation Area Scenario Simulation Settings

We designed four sets (1250 kinds) of scenarios accounting for the sown areas of both winter wheat (high irrigation demand) and summer maize (low irrigation demand) (

Table 1. The irrigation water requirements under the different precipitation probabilities based on previous studies and irrigation water quotas from the Beijing, Tianjin, Hebei, Henan, Shandong, Heilongjiang, Jilin, Liaoning and Inner Mongolia provinces.

Region	Irrigated Crops	Probability of Precipitation	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov
NCP	wheat (mm)	75%	75	75	95	-	-	-	-	-	75
		25%	75	75	95	-	-	-	-	-	-
	corn (mm)	75%	-	-	-	65	65	-	-	-	-
		25%	-	-	-	65	-	-	-	-	-
NECP	rice (mm)	75%	-	-	190	185	115	45	40	-	-
		25%	-	-	163	160	100	38	34	-	-
	corn (mm)	75%	-	-	12	50	50	30	30	-	-
		25%	-	-	8	37	37	22	22	-	-

and Figure 6): I, reduce the sown area of winter wheat and summer maize within the range of −3.00 × 10⁶~0 ha at an interval of 1.0 × 10⁴ ha; II, reduce only the sown area of winter wheat; And III, reduce only the sown area of summer maize. The reduction range and interval for Scenario II and III is the same as that of Scenario I. Scenario IV is reduce the sown area of winter wheat and increase the sown area of summer maize.

In order to optimize the agricultural planting structure and the sustainable utilization of groundwater, we set up the following four groups (1600 kinds) of scenarios (Table 1 and Figure 6): V, change the sown area of rice and maize; VI, change only the sown area of rice; VII, change only the sown area of maize; and VIII, reduce the sown area of rice and increase the sown area of maize by the same amount, or increase the sown area of rice and reduce the sown area of maize by the same amount. We only included negative changes in the sown area in the NCP, because the GWS balance can only be achieved by reducing the sown area in the overexploited NCP. The sown area of rice and maize was changed within the approximate range of −2.00 × 10⁶ to 2.00 × 10⁶ ha with an interval of 1.0 × 10⁴ ha.

Changing the sown area of grain crops will lead to changes in irrigation water demand and thus it will also change the groundwater abstraction. In this study, a linear response by GWSA to changes in sown areas, as shown in Equation (3), was developed to estimate the changes in GWSA, as shown in Equation (4).

$$\Delta GWS{A}_{t+1}=\Delta GWS{A}_t+\left(\alpha \times I_{t+1}\right)\frac{\Delta {\mathrm{A}}_{crop}}{{\mathrm{A}}_{region}}-\left(\left(1-\alpha \right)\times \beta \times I_{t+1}\right)\frac{\Delta {\mathrm{A}}_{crop}}{{\mathrm{A}}_{region}}$$

(3)

$$GWSA{A^{\prime }}_{\mathrm{t}}=GWS{A}_t+\Delta GWS{A}_t$$

(4)

where $GWS{A}^{\prime }$ is the updated GWSA following a change in sown area, and I is the irrigation (mm). Generally, there are 4 or 5 irrigations applied during the winter wheat season and 0 to 2 irrigations during the maize season in the NCP. There are ~5 irrigations applied during both the winter wheat season and the maize season in the NECP. The water requirements for each irrigation event were referenced from previous research results and the irrigation water quota (Table 1) [5,6,15,16,56,57,58,59]. $\mathsf{\alpha }$ is the recharge coefficient. The value of α was taken to be 0.1 for irrigations less than 90 mm and 0.15 for irrigations between 90 and 250 mm [5]. $\mathsf{\beta }$ is the percentage of groundwater irrigation in the total irrigation. In the NCP and NECP, the percentage of groundwater irrigation accounts for 75% and 47%, respectively [18,24]. $\Delta {\mathrm{A}}_{crop}$ is the variation in the sown area of grain crops, and ${\mathrm{A}}_{region}$ is the total area of the study region.

2.8. Uncertainty Estimation

The following describes the uncertainty estimation for each water storage component. The error for GRACE-derived GWSA comes from both TWSA and the mass change compartments of SWSA, SMSA, and SWESA. We estimated the TWSA uncertainty by propagating the errors associated with the three GRACE data products provided by CSR, JPL and GFSC. The standard deviation of the SMSA from four LSMs in GLDAS were used as the uncertainties in SMSA and SWESA. A value of 20% was assumed for the uncertainty in the SWSA estimates. All the errors were considered to estimate the uncertainties of GWSA using the error propagation method [60]. The uncertainty of $S_y$ is estimated to be 20%, which is used to estimate the uncertainty of GWS based on the in situ groundwater level data.

3. Results

3.1. Groundwater Storage Changes Monitored by GRACE

Figure 7 shows the monthly GWSA in the NCP and NECP detected by GRACE from 2003 to 2016. The reliability of GRACE-derived GWSA was verified by using in situ GWSA (Figure 7a,c). The results show that (1) the long-term trend in the GWSA estimated based on GRACE is similar to that of the in situ GWSA, (2) in the NCP, the correlation between the GRACE-derived and in situ GWSA is relatively good (r = 0.83), which indicates that the GRACE inversion results are reliable, and (3) in the NECP, the annual phases of the GRACE-derived and in situ GWSA are inconsistent and the correlation is poor (r = 0.22), which may be because the surface water storage changes in the rivers, lakes and wetlands are not well considered when subtracting the SWSA from GRACE TWSA. This seasonal inconsistency has also been revealed by previous studies [25,26]. Despite the low correlation, the seasonal inconsistency has little impact on the estimated long-term trend in GWSA, which is our major focus. Although the low correlation should be further investigated, it is considered acceptable in this study.

The long-term trends in GWSA from GRACE and in situ observations in the two areas were obtained by least squares linear fitting. From 2003 to 2016, the depletion rate of GWS detected by GRACE in the NCP is −17.2 ± 0.8 mm/yr (Figure 7a,

Table 2. The annual amplitude and phase, semiannual amplitude and phase, and the trend in equivalent water height (EWH) anomalies in GRACE-derived GWS (2003–2016).

Region	EWH Anomaly	Annual Amplitude (mm)	Annual Phase (°)	Semi-Annual Amplitude (mm)	Semi-Annual Phase (°)	Trend (mm/yr)
NCP	GWS (GRACE)	23.20	275.67	14.21	72.31	−17.2 ± 0.8
	GWS (in situ)	40.00	339.81	20.69	99.15	−16.8 ± 1.0
NECP	GWS (GRACE)	25.55	102.14	1.25	348.61	2.3 ± 0.7
	GWS (in situ)	39.73	253.74	2.58	44.13	2.1 ± 1.0

), and the in situ GWS drawdown rate is −16.8 ± 1.0 mm/yr, which is well comparable to the depletion rate (−18.2 ± 0.2 mm/yr) of GWS in 2003–2015 estimated by Gong et al. [24]. During the study period, the GWS in the NECP increased slightly, with a rate of 2.3 ± 0.7 mm/yr, and the in situ GWS had an increase rate of 2.1 ± 1.0 mm/yr. The GWSA shows the expected responses to the variations in annual precipitation (Figure 7a,b). For example, during 2005–2008 and 2014–2015, the GWS declined quicker due to precipitation deficits. In the NECP, the GWS shows more rapid responses to precipitation changes. Rapid decline in the GWS can be found in dry years such as 2004, 2007, 2011 and 2014–2015, with rapid recovery in the GWS found in wet years such as 2005, and 2012–2013.

Figure 8 shows the spatial distribution of the GWS trend in the NCP and NECP from 2003 to 2016. The GWS depletion occurs widely across the whole NCP, while the NECP generally represents stable or positive trend in GWS. The difference is a result of the combined impacts of climate and human activities.

3.2. Long-Term Groundwater Abstraction Variations and Grain Production

Figure 9 and

Table 3. Annual grain yield (GY), grain sown area (GSA), groundwater abstraction (GWA) and precipitation (P) in NCP and NECP from 2003 to 2016.

Year	NCP				NECP
	GY (106 t)	GSA (104 km2)	GWA (mm)	P (mm)	GY (106 t)	GSA (104 km2)	GWA (mm)	P (mm)
2003	42.00	15.05	144.64	660.45	76.31	16.75	72.23	494.53
2004	44.58	15.10	141.43	588.90	88.70	17.77	74.13	413.83
2005	47.87	15.20	136.93	529.17	95.89	19.40	75.52	533.85
2006	51.33	15.09	139.07	443.03	101.04	19.96	62.25	454.85
2007	52.22	15.06	139.81	514.17	100.66	20.28	80.94	385.61
2008	54.10	15.19	138.63	564.67	110.57	20.28	81.89	463.09
2009	54.65	15.15	132.40	597.85	103.86	22.61	85.11	483.00
2010	55.56	15.19	130.32	548.05	117.79	23.23	88.66	547.98
2011	58.46	15.27	130.79	562.45	131.65	23.74	88.57	425.36
2012	59.74	15.27	129.65	661.29	137.03	24.12	90.94	621.45
2013	61.30	15.26	124.19	562.17	145.24	24.47	91.89	604.29
2014	61.59	15.23	122.12	414.69	142.82	24.96	92.07	492.90
2015	62.59	15.25	119.20	546.76	148.01	25.23	88.43	514.53
2016	64.19	15.20	108.42	615.88	146.57	24.91	91.06	589.48

shows the evolution of groundwater abstraction, grain sown area, and grain yield in the NCP and NECP. From 2003 to 2016, the total grain sown area remained stable at around 15.2 × 10⁶ ha, while the total grain yield gradually increased from 42.0 × 10⁶ t to 64.2 × 10⁶ t, which represents an increase of 53% in the NCP. However, annual groundwater abstraction gradually decreased from 144 mm to 108 mm (Figure 9a). During the same period in the NECP, the sown area of grain increased from 16.75 × 10⁶ ha in 2003 to 24.91 × 10⁶ ha in 2016 and the grain yield increased from 76.3 × 10⁶ t to 146.6 × 10⁶ t, while the annual groundwater abstraction increased from 72.2 mm to 91.1 mm (Figure 9b).

We set GY as dependent variable, GSA, GWA and P as independent variables, and carried out multivariate linear regression analysis. The regression equations for grain yield in the NCP and NECP are as shown in Equation (5) and Equation (6), with r-square equal 0.88 and 0.92, respectively.

$$GY=0.349\times \overline{GSA}-0.661\times \overline{GWA}-0.130\times \overline{P}$$

(5)

$$GY=0.920\times \overline{GSA}+0.004\times \overline{GWA}+0.065\times \overline{P}$$

(6)

where $\overline{GSA}$, $\overline{GWA}$ and $\overline{P}$ are standardized values (using zero mean normalization method) of GSA, GWA and P, respectively. It can be seen from the regression equation that GY is negatively correlated with GWA in NCP, while in NECP, GY is positively correlated with GSA. This shows that the increase in GY in NCP is mainly due to the improvement in irrigation water use efficiency, while in NECP, the increase in GY is mainly due to the increase in GSA.

4. Discussion

4.1. The Projected GWSA under Different Scenarios

Figure 10 shows the simulated monthly GWSA based on Scenarios I–IV for the NCP from 2003 to 2016 and the GWS trend fitted from monthly time series. It was found that reducing the sown area of grain crops, especially winter wheat, can significantly recover GWS. Given a reduction of per 1 × 10⁶ ha in the sown area of both winter wheat and summer maize (Scenario I), only winter wheat (Scenario II) and only summer maize (Scenario III), there is an increase in GWS at the rate of 18.5 mm/yr, 13.1 mm/yr and 5.4 mm/yr, respectively. For Scenario IV, i.e., reducing the sown area of winter wheat while increasing the sown area of summer maize by the same amount, a reduction of per 1 × 10⁶ ha in the sown area of winter wheat leads to an increase in GWS at the rate of 7.8 mm/yr. The scenario simulation using GRACE-derived GWS anomalies during 2003–2016 as the baseline shows that a balance (i.e., no decreasing trend in storage) in GWS changes in the NCP can be achieved by reducing the area of winter wheat by 1.31 × 10⁶ ha, by reducing the area of maize by 3.21 × 10⁶ ha or by reducing both by 0.93 × 10⁶ ha. This finding is comparable to that from previous studies that use other methods. For example, Xu et al. found that the sustainable utilization of groundwater could be achieved by reducing the sown area of winter wheat by 48.3% in the Hebei Plain [10]. Hu et al. used the decision support system for agro-technology transfer (DSSAT) method to simulate groundwater irrigation in the Shijiazhuang irrigation district, and found that a 29.2% reduction in irrigation could stop groundwater drawdown and lead to a decrease of 8.42% and 2.11% in the yield of wheat and maize, respectively [61]. As the NCP is the main grain base in China, the grain yield cannot be greatly reduced. Therefore, it may be better to reduce the sown area of winter wheat while increasing the sown area of maize or other drought-resistant grain crops.

Figure 11 shows the simulated monthly GWSA based on Scenarios V–VIII in the NECP from 2003 to 2016 and the trend in the GWS. Given an increase of per 1 × 10⁶ ha in the sown area of both rice and maize (Scenario V), only rice (Scenario VI), and only maize (Scenario VII), there is a decrease in GWS at the rate of −6.9 mm/yr, −5.3 mm/yr, and −1.6 mm/yr, respectively. For Scenario VIII, i.e., increasing the sown area of rice while reducing the sown area of maize by the same amount, an increase of per 1 × 10⁶ ha in the sown area of rice leads to a decrease in GWS at the rate of −3.7 mm/yr. It was found that the GWS can sustain an additional increase in rice or maize in the sown area of 0.43 × 10⁶ ha and 1.45 × 10⁶ ha, respectively, without depleting the aquifers in the NECP. The GWS balance can also be maintained given a simultaneous increase in the sown area of 0.33 × 10⁶ ha for rice and maize. Since the maize production in the NECP is required to be decreased according to the Supply-Side Structural Reforms, a more rational change, i.e., increasing the sown area of 0.62 × 10⁶ ha for rice while decreasing the same amount of the sown area of maize, could be implemented to maintain the GWS balance.

4.2. Evaluation of Policies on Recovering GWS in the NCP

According to the “Five Year (2018–2022) Action Plan for Integrated Governing of Groundwater Overexploitation in the Hebei Province” (referred to as FYAP hereafter), three strategies will be implemented by 2022: (1) 0.157 × 10⁶ ha of winter wheat will be fallowed in the whole province; (2) the development of efficient water-saving irrigation area will be 0.378 × 10⁶ ha; and (3) the sowing ratio of drought resistant and water-saving wheat varieties will be over 90%. These three strategies will save 0.387, 0.311 and 1.003 billion m³ of groundwater, respectively. Based on the FYAP, we designed the following scenario simulation: the wheat fallow area is 0.157 × 10⁶ ha, the efficient water-saving irrigation area is 0.378 × 10⁶ ha, and change the sown area proportion of drought resistant and water-saving wheat varieties (0~100%). The simulated results (Figure 12) show that the groundwater will be balanced in the NCP if the sown ratio of drought-resistant and water-saving wheat is increased (as estimated by the FYAP) to 95% (Figure 12c).

According to the current agricultural Supply-Side Structural Reforms (referred to as SSSR hereafter) in NECP, 0.813 × 10⁶ ha maize sown area was reduced, and the irrigation area of rice was appropriately increased in some areas. So, we designed the following scenario simulation: reduce 0.813 × 10⁶ ha of the maize sown area, change the irrigation area of rice (−2 × 10⁶−2 × 10⁶ ha). The simulated results (Figure 12b) show that the irrigation area of rice can be increased by 0.67 × 10⁶ ha (Figure 12d) on the premise of groundwater balance.

4.3. Uncertainties

The uncertainty in GRACE-derived GWSA was estimated by propagating the error components in TWSA, SMSA, SWESA, and SWSA. The estimated uncertainties are showed in

Table 4. The monthly errors for GWSA,TWSA,SMSA and SWSA derived from GRACE, GLDAS and in situ observation.

Region	TWSA (mm)	SMSA (mm)	SWESA (mm)	SWSA (mm)	GWSA (mm)
NCP	50.50	12.41	0.07	0.73	52.01
NECP	32.71	10.55	1.70	2.11	34.48

. The monthly error for GRACE-derived GWSA is 52.01 mm and 34.48 mm in NCP and NECP, respectively, from 2003 to 2016. The monthly error for GRACE-measured TWSA is the largest component (50.50 mm/month for NCP and 32.71 mm/month for NECP from 2003 to 2016). Another major source of error is SMSA from the hydrologic models used (12.41 mm/month for NCP and 10.55 mm/month for NECP from 2003 to 2016). The total trend in the error for GWS from 2003 to 2016 is calculated as 0.8 mm/yr in NCP and 0.7 mm/yr in NECP (Table 2). The final trend in error for in situ GWS in both NCP and NECP is 1.0 mm/yr.

Based on the law of error propagation and the error of GRACE-derived GWSA, the monthly error of GWSA in the simulated scenario is calculated as 52.45–53.31 mm in NCP and 34.94–35.04 mm in NECP (

Table 5. The monthly errors of GWSA from simulated scenarios (2003–2016). FYAP and SSSR stand for Five Year (2018–2022) Action Plan for Integrated Governing of Groundwater Overexploitation in the Hebei Province and Supply-Side Structural Reforms, respectively.

NCP	GWSA (mm)	Scenario Ⅰ	Scenario Ⅱ	Scenario Ⅲ	Scenario Ⅳ	Scenario FYAP
		53.31	52.97	52.59	53.31	52.45
NECP	GWSA (mm)	Scenario Ⅴ	Scenario Ⅵ	Scenario Ⅶ	Scenario Ⅷ	Scenario SSSR
		35.00	34.97	34.94	34.95	35.04

). The error of GWSA in the simulated scenario and GRACE-derived is almost the same. This shows that the error in simulated scenarios mainly comes from the error of GRACE-derived GWSA.

5. Conclusions

In this study, we provided a case study to understand the historical and projected GWS changes at the large scale by using GRACE data, with a focus on the cause-effect relationship between agricultural water demand and the trend (balance) of GWS in two major food producing areas of China, i.e., the NCP and NECP. A linear response equation was developed to simulate the GWS changes in response to varied sown areas of major grain crops based on GRACE-derived GWSA. The simulation focused on investigating how much agricultural water demand should be reduced to achieve a balance in the GWS in the NCP, as well as how much agricultural water demand can be increased while maintaining a balance in the GWS in the NECP.

The GRACE observations showed that the GWS was being depleted at a rate of −17.2 ± 0.8 mm/yr in the NCP, while it generally remained stable (2.3 ± 0.5 mm/yr) in the NECP, during the period of 2003–2016. The GWS trends estimated from GRACE were comparable to those from in situ well measurements, despite a seasonal mismatch between each of these in the NECP, which is probably due to biased input of surface water storage changes. The contrasting trends in the GWS changes in the NCP and NECP reflected the differences in the water demand and water supply. Increasing grain production along with decreasing groundwater abstraction in the NCP during 2003–2016 indicated the critical role of agricultural practices in this region, while groundwater mainly sustained the increase in grain production in the NECP.

The scenario simulation based on GRACE-derived GWSA showed that the declining trend in the GWS can be recovered by reducing the sown area of grain crops in the NCP, while additional grain crops can be sustained by extracting more groundwater without a declining trend in the NECP. It was found that a balance (i.e., no decreasing trend in storage) in GWS changes in the NCP can be achieved by reducing the area of winter wheat by 1.31 × 10⁶ ha, by reducing the area of maize by 3.21 × 10⁶ ha or by reducing both by 0.93 × 10⁶ ha. For the NECP, our results suggest that the GWS balance can be maintained given a decrease of 0.62 × 10⁶ ha in the sown area of maize while increasing the same sown area of rice.

The simulation with realistic scenarios showed that the GWS balance in the NCP could be achieved by implementing the Five Year (2018–2022) Action Plan for Integrated Governing of Groundwater Overexploitation in the Hebei Province, under the condition that the sown ratio of drought-resistance and water-saving wheat is increased from 90% to 95%. For the NECP, the GWS balance could be maintained with a maximum increase of 0.67 × 10⁶ ha in the sown area of rice under the framework of the Supply-Side Structural Reforms.

This study demonstrated the usefulness and potential of GRACE data in understanding the food–water nexus at large scales. The linear response equation jointly used with the GRACE-derived GWSA provides a simple but effective way to simulate or evaluate groundwater responses to changes in agriculture practices. It is of significant value when an overall understanding of agricultural water demand and groundwater sustainability is needed for regions where little information is available for complex model simulations. However, it should be noted that the uncertainties associated with the simplified parameterization and the spatial-temporal variabilities in the cause-effect relationship (including climate change), may lead to biased estimates of GWS responses.