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

: Balancing groundwater supply and food production is challenging, especially in large regions where there is often insu ﬃ cient 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 e ﬃ cient 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 6 ha and 3.21 × 10 6 ha, respectively, or by reducing both by 0.93 × 10 6 ha. In the NECP, the groundwater can sustain an additional area of 0.62 × 10 6 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 e ﬀ ectiveness of using GRACE to understanding the nexus between groundwater supply and food production at large scales.


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], 3 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 GRACEobserved 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)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(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.

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]: Remote Sens. 2020, 12, 2151 4 of 21 This study uses in situ GWL anomalies (GWLAs) to validate GRACE-derived GWSA.
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. 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].

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].

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.

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.

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

In Situ GWL Observations
To verify the accuracy of GRACE-derived GWSA, monthly in situ GWL data ( Figure 1

In Situ GWL Observations
To verify the accuracy of GRACE-derived GWSA, monthly in situ GWL data (  Figure 5. In situ GWSA is calculated by multiplying in situ-measured GWLA with a specific yield ( ) of each well: In this study, the spatial distribution map of the in the NCP was provided by the China Geological Survey Bureau, and the ranges from 0.025 to 0.29 [26]. In the NECP, a mean value of 0.095, determined by pumping tests, was used to convert GWLA to GWSA [27,55].

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 In situ GWSA is calculated by multiplying in situ-measured GWLA with a specific yield (S y ) of each well: GWSA = S y × 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].

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

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).

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 and Figure 6): I, reduce the sown area of winter wheat and summer maize within the range of −3.00 × 10 6~0 ha at an interval of 1.0 × 10 4 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

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).

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 and Figure 6): I, reduce the sown area of winter wheat and summer maize within the range of −3.00 × 10 6~0 ha at an interval of 1.0 × 10 4 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 6 to 2.00 × 10 6 ha with an interval of 1.0 × 10 4 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).
where GWSA 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]. α 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]. β 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]. ∆A crop is the variation in the sown area of grain crops, and A region is the total area of the study region.  -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 --

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. 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.

Groundwater Storage Changes Monitored by GRACE
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), 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 Figure 9 and Table 3 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 6 ha, while the total grain yield gradually increased from 42.0 × 10 6 t to 64.2 × 10 6 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 6 ha in 2003 to 24.91 × 10 6 ha in 2016 and the grain yield increased from 76.3 × 10 6 t to 146.6 × 10 6 t, while the annual groundwater abstraction increased from 72.2 mm to 91.1 mm (Figure 9b).  Figure 9 and Table 3 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 6 ha, while the total grain yield gradually increased from 42.0 × 10 6 t to 64.2 × 10 6 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 6 ha in 2003 to 24.91 × 10 6 ha in 2016 and the grain yield increased from 76.3 × 10 6 t to 146.6 × 10 6 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.

Long-Term Groundwater Abstraction Variations and Grain Production
where ̅̅̅̅̅̅ , ̅̅̅̅̅̅̅ and ̅ 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   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.
where GSA, GWA and 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. 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 6 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 6 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 6 ha, by reducing the area of maize by 3.21 × 10 6 ha or by reducing both by 0.93 × 10 6 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 6 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 6 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 6 ha and 1.45 × 10 6 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 6 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 6 ha for rice while decreasing the same amount of the sown area of maize, could be implemented to maintain the GWS balance. Remote Sens. 2020, 12, x FOR PEER REVIEW

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 6 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 6 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 3 of groundwater, respectively. Based on the FYAP, we

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 6 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 6 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 3 of groundwater, respectively. Based on the FYAP, we designed the following scenario simulation: the wheat fallow area is 0.157 × 10 6 ha, the efficient water-saving irrigation area is 0.378 × 10 6 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). 17 designed the following scenario simulation: the wheat fallow area is 0.157 × 10 6 ha, the efficient watersaving irrigation area is 0.378 × 10 6 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 6 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 6 ha of the maize sown area, change the irrigation area of rice (−2 × 10 6 −2 × 10 6 ha). The simulated results (Figure 12b) show that the irrigation area of rice can be increased by 0.67 × 10 6 ha ( Figure 12d) on the premise of groundwater balance.

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 Table 2). The final trend in error for in situ GWS in both NCP and NECP is 1.0 mm/yr. According to the current agricultural Supply-Side Structural Reforms (referred to as SSSR hereafter) in NECP, 0.813 × 10 6 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 6 ha of the maize sown area, change the irrigation area of rice (−2 × 10 6 −2 × 10 6 ha). The simulated results (Figure 12b) show that the irrigation area of rice can be increased by 0.67 × 10 6 ha (Figure 12d) on the premise of groundwater balance.

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 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 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.

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 6 ha, by reducing the area of maize by 3.21 × 10 6 ha or by reducing both by 0.93 × 10 6 ha. For the NECP, our results suggest that the GWS balance can be maintained given a decrease of 0.62 × 10 6 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 6 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.