North China Plain (NCP) in China is an important agricultural region increasingly dependent on groundwater to meet the demands of water for irrigation which consequently has resulted in groundwater depletion. Quantifying spatio-temporal variations of groundwater storage (GWS) is important in NCP for monitoring groundwater depletion. Gravity Recovery and Climate Experiment (GRACE) satellite data provide the potential for quantifying regional GWS changes. However, its coarse spatial resolution and errors in disaggregation have limited the application of GRACE for localized groundwater studies, which are essential for effective groundwater management. We, therefore, implement a Random Forest (RF) Machine Learning (ML) model to establish an empirical relationship between GRACE-derived Terrestrial Water Storage variations (TWS), hydro-meteorological variables, and available in situ groundwater level data for shallow and deep aquifers. In-situ and RF modeled groundwater level variations show a high correlation during training and validation. Therefore, the modeled empirical relationship was extended to the whole of NCP to produce monthly GWS variations at 5 km resolution. This resolution is similar to previous downscaling studies. Deep aquifers show rapid GWS losses compared to shallow aquifers suggesting a relatively slow recharge process of the deep layers of the groundwater reservoirs. The methodology presented in this paper shows an effective downscaling of GRACE mass change observations for localized GWS assessment which can also be replicated in other regions.
Research Article
Machine Learning application for modeling high-resolution groundwater storage variations in North China Plain
https://doi.org/10.21203/rs.3.rs-2062965/v1
This work is licensed under a CC BY 4.0 License
Version 1
posted
You are reading this latest preprint version
Groundwater, an important component of the global water cycle, is also a vital freshwater resource for many arid and semi-arid regions around the world. Approximately half of the drinking water and 40% of the industrial water needs are met through the groundwater (Wada et al., 2014; Zekster and Everett, 2004). Groundwater constitutes 43% of consumptive irrigation water usage (Siebert et al., 2010), making the availability of groundwater essential for food sustenance in many regions.
Approximately 1.5 billion people use groundwater for drinking, therefore population growth can have a direct correlation with groundwater usage (Famiglietti et al., 2014). Further, rapid industrial and agricultural growth, and climate-induced droughts, among others have led to increased groundwater dependence in many aquifers around the world. As a result, groundwater recharge rates from precipitation and inflows have been unable to match the discharge rates from withdrawals and outflows. This condition has resulted in the loss of groundwater storage (GWS) over time (Wada et al., 2010). The impacts of such groundwater depletion are already visible in the form of degradation in groundwater quality and quantity such as saltwater intrusion, soil salinization, and land subsidence. These groundwater-related issues could also have huge socio-politico-economic impacts on a local, regional, and transboundary scale (Famiglietti et al., 2014).
Studies in the past have used diverse methods to reliably estimate GWS and monitor the health of aquifers (Bierkens and Wada, 2019). Monitoring groundwater level in wells provides critical information on the aquifers and can help quantify GWS. Groundwater levels can be multiplied by storage coefficients or specific yield to obtain storage changes. Such field data is also important for the numerical modeling of aquifers. However, groundwater levels are point observations having limited spatio-temporal coverage, making their applicability for regional studies difficult. Further storage coefficients also have uncertainties which might propagate into uncertainties in GWS.
GRACE is a unique satellite mission that has enabled a continuous measurement of Terrestrial Water Storage (TWS) from 2003 to 2018 at a spatial resolution of 333 km half-wavelength (Frappart et al., 2018). TWS integrates water storage components related to groundwater, surface, soil, snow, and canopy and is an important hydrologic quantity. GWS variations can be obtained by disaggregating other components using external data sources from GRACE-derived TWS (Famiglietti et al., 2011; Feng et al., 2013). However, numerous error sources related to disaggregation and coarse resolution of GRACE data limit its applicability for local groundwater studies (Scanlon et al., 2012). For instance, soil moisture storage is modeled from Land Surface Models (LSMs) which is subtracted from GRACE-derived TWS to obtain GWS variations. However, irrigation water usage and recharge are not simulated by the LSMs which causes errors in derived soil moisture.
In addition, Machine Learning has been proposed as an alternative approach to downscale GRACE derived TWS and make it more relevant for local-scale GWS estimation. This modeling approach relies on establishing a good empirical relationship between the input and output variables (Dramsch et al. 2020). Therefore, detailed knowledge of physical properties of aquifer is not required for accurate modeling of the GWS variations (Adamowski and Chan, 2011; Miro and Famiglietti, 2018; Rahaman et al., 2020). Several ML models, ranging from Multiple Linear Regression (MLR), Artificial Neural Network (ANN), Random Forest, Boosted Regression Tree, and Deep Learning have been used in different groundwater-related modeling applications and achieved good accuracy (Bai et al., 2022; Chen et al., 2019; Rasool et al., 2022).
Quantifying GWS is crucial for North China Plain (NCP). This region is an important grain base for China where increasing urban and industrial water usage along with a focus on growing water-intensive crops has overstressed the aquifer and depleted the groundwater resources (Liu et al., 2008; Zheng et al., 2010). Several negative effects such as severe falling groundwater levels, land subsidence, drying of rivers, seawater intrusions, shrinking of wetlands have been reported (Changming et al., 2010; Chen and Li, 2004; Zhang et al, 2009). In response, several studies have assessed regionwide groundwater depletion in NCP using groundwater monitoring wells (Cao et al., 2013), numerical modeling (Cao & Zheng, 2016; Li et al., 2017) and GRACE (Feng et al., 2013, 2018), or a combination of the above techniques (Gong et al., 2018). Besides, some studies have used geodetic information from Interferometric Synthetic Aperture Radar (InSAR) (Zhou et al., 2016) and GPS (Wang et al., 2017) to quantify GWS.
Several approaches have been implemented to improve the accuracy of GWS in NCP by mathematical and statistical techniques. These techniques seek to integrate ground-truth data with GRACE and modeling data in NCP to improve accuracy as well as improve the resolution of GWS. Huang et al., (2015) used the forward modeling approach using a priori information derived from in situ groundwater level measurements and other field data to estimate GWS variations in NCP. Yin et al., (2018) implemented a statistical downscaling approach based on the simple linear regression between coarse resolution GRACE-derived GWS variations and high-resolution evapotranspiration. However, this approach is only applicable where a strong correlation between evapotranspiration and GWS variations exists. The results were therefore not satisfactory in all regions within NCP. Yin et al., (2020) assimilated GRACE-derived TWS in a hydrological model to improve the spatial resolution and accuracy of GWS estimates for NCP. However, their GWS estimates at 0.5° resolution could be improved in relation to GWS from in-situ wells.
Chen et al. (2019) used Random Forest (RF) Machine Learning (ML) model to establish an empirical relationship between hydrological variables and GRACE-derived TWS variations as the output variable in NCP. They disaggregated various hydrological components from the modeled TWS to derive the GWS at 0.25° resolution in NCP from 2003 to 2016. However, they only validated the GWS variations in Handan prefecture of Hebei province and encountered significant challenges in modeling groundwater variations. Similarly, Zhang et al., (2021a) used different ML methods to downscale the monthly GWS variations to 0.25° resolution in NCP. However, they used only 12 wells from shallow aquifers for training and 6 wells for validating the results. Specifically, they constructed models for individual groundwater wells limiting the applicability of their method to cases with few wells. Their results could only capture the long-term trends and not model groundwater levels numerically, reflecting poor accuracy. Besides, these studies used the GRACE data without consideration of leakage and other correction, which are quite important. Zhang et al., (2021b) downscaled annual GWS over the whole of China, using a variety of datasets at 1 km resolution using GRACE and hydro-meteorological data.
The above studies prove that empirical machine learning methods to downscale the GRACE-derived TWS variations are effective. However, no study has been performed to integrate in-situ groundwater level data ML to model localized spatio-temporal GWS variations in NCP. More specifically we aim to fill the following science gaps: (i) Investigate the usefulness of the RF ML model for establishing an empirical relationship between the GRACE-derived TWS variations and other hydrological/meteorological variables as model input while in-situ groundwater level variations as target variables, (ii) obtain regional estimates of GWS from ML modeling and test its consistency with past studies. (iii) obtain the GWS trend with a resolution as fine as 0.05° (~ 5 km) over the whole of NCP.
North China Plain (NCP), with an area of approximately 140,000 km2, is the largest alluvial plain and one of the most important grain bases in China (Fig. 1). It includes plain areas of Beijing and Tianjin municipality, most of the Hebei provinces, and parts of the Henan and Shandong provinces to the north of the Yellow River. Based on the hydrogeological properties, NCP is usually divided into the Piedmont plain (PP) and the East Central Plain (ECP).
NCP consists of a shallow unconfined aquifer and three deeper confined aquifers with varying depths (40–60 m, 120–170 m, 250–350 m, and 400–600 m (Sakura et al., 2003). Since precipitation is lower in NCP compared to regions in South China, irrigation is important from spring to early summer which are important growing seasons. The aquifers in NCP have been primarily tapped to increase the maize and wheat production in NCP.
In Piedmont Plains, upper unconfined and sometimes confined aquifers are significant sources of groundwater for irrigation intensively used in this region (Shu et al., 2012). In East Central Plains, deeper confined aquifers are primarily tapped for the domestic and industrial water needs as shallow aquifers mostly contain brackish water owing to the seawater intrusion from the Bohai Gulf (Wang et al., 2007). Most of the recharge in NCP occurs in the shallow unconfined aquifers of Piedmont Plains from the precipitation or the lateral recharge from the neighboring Taihang Mountains (Cao et al., 2013). Other sources of groundwater recharge in the NCP include irrigation return flow and the leakage from the Yellow River and irrigation canals.
3 Data and Methods
The methodology adopted in this paper is shown in Fig. 2.
2.1 Random Forest Model
Random Forest (RF) model has been used in various geoscience applications such as global soil property mapping (Hengl et al., 2017), downscaling GRACE data (Milewski et al., 2019, Rahaman et al., 2019), or assessing the groundwater overextraction related contamination (Smith et al., 2018). The algorithm is well-described in Breiman, (2001) and we discuss it concisely here.
It is a robust ML model consisting of an ensemble of decision trees. Each decision tree consists of several decision-making elements called nodes architectured in the form of a tree. The nodes use different criteria to split their data into two child nodes. The node at the top in a decision tree is called the root node, while the child node at which splitting stops is known as the leaf node. Splitting continues until the following user-defined condition is met, such as the maximum number of splits, the depth of decision tree, or the minimum number of samples in the leaf node. These user-defined values are called hyperparameters. For a good split at a node, the input variable should be chosen to ensure child nodes are “purer” (i.e., more homogeneous) in terms of the outcome variable than the parent nodes.
While a single decision tree can have low bias, it can also have high variance, a condition known as overfitting. Ensembling several trees together is useful to reduce the variance compared to a single decision tree. RF algorithm is built on the concept of ensembling (Breiman, 2001). It consists of the following two steps:
1. Select a random sample of observations (with replacement) from the entire dataset. This process is known as bagging. For each decision tree, samples selected for model building are known as “in bag” samples, while the “out of bag” (OOB) samples can be used as independent data set for model validation.
2. Randomly select a few features for splitting at each decision tree to reduce the correlation between the decision trees and prevent overfitting.
Different random splitting criterion applied at each decision tree results in a different model being built for each decision tree. Steps 1 and 2 are repeated/iterated until the stopping criterion, similar to a single decision tree, is achieved.
Finally, “predictions” from each of the “in bag” samples of different decision trees are aggregated for the final output, and errors on “out of bag” samples can additionally be used as an accuracy measure. This process increases the robustness of modeled results and protects against overfitting.
2.2 Input Data to the ML model
2.2.1. Precipitation and Temperature
We obtained monthly precipitation (PPT) and temperature (TEMP) from China Meteorological Data Service Center (CMDSC) at a grid resolution of 0.25° (~ 25 km).
2.2.2. Groundwater Storage
We used monthly GRACE Level 2 (L2) global gravity field data version 06 in the form of spherical harmonic coefficients (SHCs) available at the Center for Space Research (CSR) at the University of Texas at Austin. The SHCs are truncated to the degree and order(d/o) of 60 to suppress the noise of higher d/o SHCs (representing shorter wavelengths). This truncation represents low pass filtering in the spatial domain, causing a limited spatial resolution of GRACE data due to signal dampening. We further replace the degree 1 coefficients and the zonal degree 2 coefficients following the GRACE science data product recommendations (Swenson et al., (2008); Cheng & Ries, (2018)). The effect of Glacial Isostatic Adjustment (GIA) geophysical process is corrected using a GIA forward model (A et al., 2013). We have used a non-isotropic DDK4 filter to reduce the correlated errors in the original spherical harmonics (Kusche, 2007), which cause North-South stripping errors in the subsequent TWS changes.
The mean SH field was first computed between 2005–2013 and it was removed from each monthly field to obtain monthly anomalous SHCs. Anomalous SHCs were then converted to gridded TWS anomaly in Equivalent Water Height (EWH) (Wahr et al., 1998) at an oversampled grid resolution of 0.25°. The 0.250 interpolated EWH resolution is justified by the published data grid size of the CSR masons (Save et al., 2016) which with empirical constraints, extended the GRACE EWH solutions to a grid size of 25 km. To account for leakage-out (amplitude damping caused because of signal leakage from the target region to the surrounding areas) and leakage-in (signal leaking in the target region), we apply the signal leakage correction based on the iterative forward modeling without a priori information (Chen et al., 2014). The procedure is critical when the basin size approaches the GRACE spatial resolution size of 200,000 km2 and the process for accounting for attenuation is mentioned in Agarwal, (2021) and Agarwal et al., (under review). We used the root zone soil moisture content for the depth of 0-200 cm at 0.25° resolution from Global Land Data Assimilation System (GLDAS) Noah Land Surface Model L4 V2.1 (GLDAS_NOAH025_M) [accessed October 2020] for NCP (Rodell et al., 2004). Long-term mean of the soil moisture content during the study period of 2005–2013 is removed to obtain the soil moisture anomaly to be consistent with the GRACE TWS anomaly. We used the difference, TWSA- SMA, a rough approximator of GWS anomaly for NCP as one of the input variables for the RF ML model.
Here we have used a lag of 0, 1, 2, and 3 months for GWS, precipitation and temperature as input variables in ML modeling. The corresponding input variables have been coded with acronym name of the variable and lag month, for instance, precipitation has been coded as PPT0, PPT1, PPT2, PPT3 for lags of 0, 1, 2, and 3 months.
3. Saturated hydraulic conductivity (K)
K is estimated globally using data-driven empirical functions by Zhang et al., (2018). In their approach, easily obtainable information about soil such as global soil texture and bulk density is used to derive K at 1 × 1 km2 resolution globally. Spatially K decreases from the west to the east (Fig. 1), a patch of high K can be found in the central plains in and around Hengshui and Xingtai prefectures which is also the region where agriculture is prominent.
4. Topographic Slope
We used Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) with 1° (~ 30-m) ground resolution (Rabus et al., 2003). The landform is primarily flat, with an elevation of less than 100 m. Topography inclines from north, west, and southwest towards the Bohai Sea.
2.3 Target (output) variable to the ML model: Groundwater well data
Groundwater well observations in NCP were obtained from the yearly handbook compiled by the Institute of China Geological Environment Monitoring. In -situ groundwater level data was available from 2005–2013. We selected the groundwater wells based on the spatial distribution and continuity of the observations. Accordingly, we only used observation from 74 wells for building our model, 38 of which were shallow wells. According to our segregation, the shallow and deeper wells have depth < = 25 m and > 25 m, respectively.
2.4 ML model building and feature importance
We use a lag term for each input variable except the topographic slope and saturated hydraulic conductivity (K), which are static variables. For instance, PPT0, PPT1, PPT2 refer to the total precipitation at the current month, previous month, and the month before, respectively. Similar naming conventions are applied to other variables and a lag of 3 months has been used for each dynamic variable based on our empirical tests.
To prepare the data for ML modeling, we separated 15% of the original data from the model building process. This test data was separated randomly such that it includes different months throughout our study period 2005–2013. Rest 85% of the data is. We used 5-fold cross-validation that uses 5 different combinations of training and validation data to ensure the robustness of the final model (Agarwal, 2021). We implement the Bayesian Hyperparameter Optimization technique to identify the most optimal RF model hyperparameters using the cross-validated training and validation dataset (Snoek et al., 2012). The accuracy of the trained model is assessed on the independent test data, which is not used for model-building process.
Further, to identify the most important input variable, we use the drop column method where an input variable is dropped and the percentage increase/decrease in error on validation data relative to the baseline model is computed (Parr et al., 2020). Baseline model is the “best” possible RF model after Bayesian Hyperparameter Optimization using all the input variables. Other models using dropped input variables are also trained from scratch using Bayesian Hyperparameter Optimization. The relative increase in RMSE error is used to assess the importance of input variables. The input variables which cause any reduction in RMSE error are dropped from the RF model.
To test the accuracy of the modeled results against the actual in-situ values, we have used the following measures in this study: Root Mean Squared Error (RMSE), scaled RMSE (RMSE*), Correlation coefficient (r), Nash-Sutcliffe Efficiency (NSE). For a detailed description of these measures, the reader is referred to the Supplementary document.
2.5 Conversion of groundwater level to groundwater storage anomaly (GWSA)
The groundwater level GWL changes were converted to GWSA by first removing the long-term mean of the GWL (Eq. 1) and then multiplying it with the specific yield or the storage coefficient of the aquifers (Eq. 2).
The spatial distribution map of the specific yield and storage coefficient in the NCP is provided by the China Geological Survey Bureau (Zhang et al., 2009). The distribution of the specific yield and storage coefficient in NCP is such that the lowest values are found near the coastal plain around the margin of the Bohai Sea and are greatest in the piedmont plain region of the Taihang Mountains.
3.1 Overall results
The overall results for shallow and deep wells suggest that RF is suitable for modeling groundwater variations in NCP (Fig. 3a and b). For shallow wells, correlation, RMSE, NSE and RMSE* for training (test) data are 0.99(0.97), 0.59 (1.35), 0.97(0.91), and 0.10 (0.25) respectively. For deep wells, correlation, RMSE, NSE and RMSE* for training (test) data are 0.95(0.90), 1.43 (2.89), 0.94 (0.88), and 0.15 (0.29). The accuracy of the training and test data area are comparable and indicate that overfitting has been avoided during the training phase. For shallow wells, the accuracy seems to be good throughout, while for deep wells, accuracy is worsened slightly for deeper wells (Fig. 3a). A similar analysis is also done for out-of-bag modeled results, showing excellent agreement with the in-situ groundwater level (Figure S1).
3.2 Feature importance
For shallow wells, precipitation is the most important input variable followed by GWSA with the significance of 33% and 31%, respectively. Hydraulic conductivity and topographic slope show relative importance of 18% and 13%, respectively, while temperature is the least significant (5%). For deep wells, GWSA is the most important input variable followed by hydraulic conductivity with a significance of 47% and 21%, respectively. Temperature and slope show lesser importance of 14.5% and 12%, respectively. Precipitation is the least significant (~ 5.5%).
Saturated hydraulic conductivity (K) is an important hydrogeological parameter that determines water flow in soils and sediments and helps estimate groundwater transport and recharge. It holds a high importance for NCP, though its importance varies slightly across shallow and deep wells. The slope determines the direction of groundwater flow and affects the groundwater recharge behavior (Bai et al., 2022). It has similar importance for both shallow and deep wells. These two variables are also an important component of integrated hydrological models (Faunt, 2009).
Temperature is an important variable for describing the climatic conditions which might affect the groundwater usage patterns. For instance, groundwater provides drinking water year-round, while also providing water for irrigation and heating during winter months. Temperature holds high importance for both Precipitation is also a critical input variable as it determines incoming flux to the groundwater system. The lesser influence of precipitation in the deeper wells than the shallow wells is reasonable and reflects that the surface water due to rainfalls cannot reach the deep compartments of the aquifer systems in the region.
GWSA from GRACE-derived TWSA and soil moisture anomaly has high importance both shallow and deep wells. GRACE-derived TWSA incorporates storage variations due to human activities such as irrigation which became prominent during droughts of 2003–2009 in NCP (Yin et al., 2020). Soil moisture from the GLDAS does not model the irrigation water usage. Groundwater discharge for irrigation is higher for deeper wells in East Central Plains which might explain the higher importance of GWSA for deeper wells.
The fact that feature importance within the RF models for shallow and deeper wells are different implies that different processes might be determining the groundwater level patterns. It justifies our usage of more complicated input patterns and possibly additional geophysical/hydrogeological information for better modeling of GWS variations in deeper aquifer systems. However, till additional information is available, the results presented here provide a reliable estimate of groundwater level variations.
3.3 Groundwater level variations for different regions in NCP
Visual analysis of Fig. 4 and Table 1 proves that RF can model the groundwater level variations at the scale of individual wells for both shallow and deep wells (See Supplementary section for additional examples). The seasonality of groundwater level variations is captured adequately by the model (Fig. 4), while the trend is underrepresented (Table 1). Here we applied the unweighted 6-parameter least-squares fitting with annual and semi-annual terms included in GRAMAT software to estimate the trend (Feng et al., 2018).
Groundwater level for shallow wells shows primarily an increasing trend in Beijing municipality and decreasing trend in Tianjin municipality. Prefectures in Hebei province show diverse groundwater level variations (Fig. 4 and S2). Baoding and Cangzhou show a large negative trend, Shijiazhuang and Xingtai show a moderate negative trend, while Hengshui and Qinhuangdao show less negative or positive trends. In Henan province, wells in Xinxiang and Jiaozuo show significant declines while the other cities show increasing trends in groundwater level. Groundwater level for deep wells shows a significant decreasing trend in Beijing and Tianjin municipality. Similarly, in Hebei province, almost all wells show significant groundwater level declines (except the wells in Fig S3f, k), the highest decline noted in Hengshui, related to intensive irrigation withdrawal. Groundwater level in shallow wells started to rise after 2010–2011 (shallow wells in Fig. 4a, b, e), whereas it continued to decline for deeper wells throughout the study period. However, a few of them, such as those at Hengshui and Jiaozuo, have started to recharge after 2012. This is most likely due to higher precipitation being recorded for the years 2010 and 2012–2013 (Fig. 2).
| S. No. | REGION | Correlation | RMSE (m) | RMSE* | NSE | GWS trend (cm yr− 1) | |
|---|---|---|---|---|---|---|---|
| In-Situ | Modeled | ||||||
| Shallow wells | |||||||
| 1 | Beijing | 0.95 | 0.26 | 0.37 | 0.87 | -15.3 ± 3.7 | -13.4 ± 2.4 |
| 2 | Xingtai (Hebei) | 0.92 | 0.55 | 0.48 | 0.77 | -20.2 ± 7.2 | -19.4 ± 4.1 |
| 3 | Hengshui (Hebei) | 0.90 | 0.36 | 0.52 | 0.73 | -0.1 ± 4.6 | -1.2 ± 2.1 |
| 4 | Tangshan (Hebei) | 0.91 | 0.44 | 0.54 | 0.70 | 10.8 ± 5.7 | 3.8 ± 2.7 |
| 5 | Jiaozuo (Henan) | 0.96 | 0.53 | 0.30 | 0.90 | -59.6 ± 5.5 | -52.6 ± 4.8 |
| Deep wells | |||||||
| 1 | Beijing | 0.95 | 1.89 | 0.29 | 0.90 | -229.3 ± 6.2 | -187.8 ± 11.4 |
| 2 | Tianjin | 0.95 | 1.53 | 0.35 | 0.88 | -165.5 ± 4.9 | -128.7 ± 6.4 |
| 3 | Baoding (Hebei) | 0.97 | 0.58 | 0.25 | 0.93 | -79.0 ± 4.9 | -69.86 ± 4.5 |
| 4 | Hengshui (Hebei) | 0.96 | 1.39 | 0.30 | 0.90 | -80.6 ± 17.4 | -72.0 ± 9.4 |
| 5 | Cangzhou (Hebei) | 0.96 | 0.81 | 0.29 | 0.91 | -44.5 ± 11.5 | -39.±7.3 |
3.4 Spatial variations in GWS
GWS trend map for shallow and deep wells, respectively at grid resolution of 0.05° (~ 5 km) is shown in Fig. 5a and b. 4–5 km is also the modeled resolution by previous studies such as Agarwal, (2021) and Miro et al. (2018). Handan city in Hebei province and Xinxiang city in Henan are the key regions where shallow wells show high groundwater declines in Piedmont plains. In the central flood plain region lies the main irrigated production zone and slight groundwater depletion can be seen in Hengshui city for shallow wells, which has been studied before (Tamanyu et al., 2009; Hu et al., 2010). Groundwater from shallow wells might also be used for small cities in other parts of Hebei province and Beijing. The coastal areas do not show any significant depletion for shallow wells primarily because groundwater in such wells might be saline and not useful for industries or drinking (Wang et al., 2009) and more water is available in the confined aquifer in these areas (Jingli et al., 2013). Further, irrigation makes a significant contribution (as high as 80%) to agriculture production in most of the piedmont plain due to the less precipitation and more intensive agriculture, while the proportion of groundwater consumed for irrigation in coastal areas is low (Fig. 5c) which might result in lesser depletion for shallow wells (Zhang et al., 2012). However, depression cones related to groundwater over withdrawal have been observed in piedmont plain and coastal areas. The appearance of these cones might also be related to the long-term history of groundwater withdrawal and local geology of the area that the field studies can only confirm.
Groundwater depletion in deeper wells is more widespread throughout the piedmont and East-Central Plain, and can be seen in the Baoding, Shijiazhuang, Cangzhou, Hengshui, Beijing, and in coastal areas of Tianjin (Fig. 5b). These are also the locations where the depression cones have been reported (MLR, 2009, Hu et al., 2004), validating the accuracy of the modeled results. Typically, deep, confined aquifers supply drinking water and industrial activities for big cities, while shallow, unconfined aquifers supply water for small cities and agriculture, leading to more subsidence in the middle east section of NCP (Changming et al., 2010).
3.5 Temporal variations in GWS
As seen in Fig. 2, precipitation is typically highest in the summer months of May-July, while TWSA shows seasonal variations and has the minimum value in June and reaches a peak in August-September. TWSA shows a declining trend during 2005–2009 and stable-to-increasing character thereafter. Such variations are dictated by the balance between evapotranspiration and precipitation (Cao et al., 2013). Further, soil moisture anomaly is the lowest in May-June, while the wettest surface soil is observed in August-September as the soil gets recharged with the precipitation. The soil moisture anomaly shows a declining trend during 2005–2008, stable character during 2008–2011, and an increasing trend after that, which might be closely linked with the variations in precipitation.
Figure 6a shows the in-situ and modeled GWS variations considering only the observations for the in-situ wells (point measurements) shown in Fig. 1b and c. The well observations are the groundwater level records at specific days in each of the months, whereas the input data of the ML models are either monthly (such as GRACE TWSA) or aggregated to monthly values using proper tools. Therefore, differences in results might be due to temporal aliasing as expected, modeling errors as well as issues related to the spatio-temporal coverage of the groundwater wells.
GWS derived from the model reaches a peak in January (shallow wells)- February (deep wells), which is similar to the behavior of the in-situ groundwater wells. Since soil moisture storage reaches a peak during August- September in NCP, the lag between soil moisture and groundwater peaks is 4–5 months and precipitation is most effective for recharging root-zone soil moisture in the region (Han et al., 2008). GWS shows a rapid decline from March-June as intensive groundwater pumping is required during March-May for wheat growth when little precipitation is available. When precipitation arrives during June- September, it sustains the water demand for maize (Pan et al., 2017) and groundwater starts to rebound slowly through the autumn.
During our study period, in-situ (modeled) GWS trend from the 38 shallow wells is -11.5 ± 4.9 (-17.3 ± 4.6) and − 12.7 ± 4.1 (-10.4 ± 1.3) mm EWH yr− 1 respectively over 2005–2013 and 2007–2009 (Fig. 6a). Here 2007–2009 is the period with low precipitation which shows higher decline in GWS. In-situ (modeled) GWS trend from the 36 deep wells is -56.8 ± 6.9 (-49.9 + 3.6), -36.9 + 8.6 (-51.9 + 10.7) mm EWH yr− 1 respectively for 2005–2013 and 2007–2009 respectively.
Figure 6b shows the modeled GWS variations for the grid at 0.05° resolution as shown in Fig. 5. The modeled GWS shows a decline of -15.1 ± 1.8 mm EWH yr− 1 during 2005–2013 while the highest decline of -30.1 ± 8.3 mm EWH yr− 1 is observed during 2007–2009.
For deep groundwater wells, the situation is rather different; the groundwater level continued declining from 2005 to 2013, showing that the GWS variability in deeper wells in NCP is less influenced by variations in precipitation. The GWS trend is -36.9 ± 4.0 mm EWH yr− 1 during 2005–2013, while the highest rate is observed between 2007–2009 as -59.7 ± 12.1 mm EWH yr− 1. The rate of decline in GWS is much higher in deeper wells compared to the one observed for the shallower wells, which agrees with reported results of Hu et al., (2016). GWS trend considering both shallow and deep wells is -24.0 ± 2.3 and − 44.2 ± 10.3 mm EWH yr− 1 over 2005–2013 and 2007–2009, respectively.
Huang et al., (2015) estimated a GWS trend of -28.1 ± 8 mm EWH yr− 1 during 2005–2010 for the whole of the NCP based on the downscaled GWS derived from GRACE, hydrological and in-situ groundwater level data. During the above period, the GWS trend is -22.1 ± 2.8 mm EWH yr− 1 and − 42.9 ± 4.2 mm EWH yr− 1 for shallow and deeper wells, respectively, based on our modeled results yielding an average GWS rate of -31.9 ± 5.5 mm EWH yr− 1 for the whole of the NCP. Thus, our results are comparable to that of Huang et al., (2015).
For the computation of GWS rate in terms of km3 yr− 1, we considered the whole of the NCP covering 144,000 km2 for shallow well and only the 100,000 km2 area for the extent of the deeper wells (Fig. 5). The GWS rate between 2005–2013 is -2.1 ± 0.3 km3 yr− 1 for shallow aquifer and − 3.7 ± 0.5 km3 yr− 1 for deeper wells. The total rate of GWS loss is thus − 5.9 ± 0.8 km3 yr− 1 over 2005–2013. Other estimates such as the one by Feng et al., (2018) is -7.2 ± 1.1 km3 yr− 1 for a different period 2002–2014 and from GWSA derived from GRACE and other hydrological data. However, our estimates are still within the same order of magnitude as their estimates and within the uncertainty range suggesting that our model can estimate the correct magnitude of GWS changes in the region.
It is important to note that the implementation of the South- to- North Water Diversion (SNWD) project since 2003 and the regulation imposed by the government for controlling groundwater withdrawal in NCP since 2002 has also been suggested to play a part in the lowering of the groundwater depletion in recent times (Long et al., 2020). Interactions of diverse factors like these might play a role in deviation of results from one approach to another.
This study uses a ML-based approach to analyze groundwater depletion during 2005–2013 in the NCP region, China’s most important grain base. We used GRACE-derived TWSA and other hydro-meteorological data as input variables and in situ groundwater level measurements as the output variable for the RF model. We created two separate models for shallow wells and deeper wells because of the different groundwater fluctuating patterns observed. Both models show good accuracy and feature importance plot shows that RF-based ML approach is suitable for downscaling GRACE-derived TWSA for modeling local groundwater changes.
We validated our modeled results with in-situ groundwater level variations in different provinces and municipalities within NCP. We obtained good accuracy and the modeled trend of groundwater level for the shallow wells agreed with the actual trends; however, for the deeper wells, the modeled trend was underestimated compared to actual observations. Further improvements in designing models for deeper wells can be undertaken in later studies by using alternative time lags and/or incorporating more geophysical information about the subsurface variations, if available.
We also used the modeled results from the RF model to generate a GWS trend map at 5 km resolution for both the shallow and deep aquifer systems over the entire NCP. The map captured the localized regions undergoing groundwater depletion and can be a useful resource for management authorities to identify local changes in groundwater resources. For instance, the map shows that depletion is more prominent in shallow aquifers in the Piedmont plain region in the south, while both Piedmont plain and East-Central Plain display prominent groundwater depletion for deeper confined aquifers which are more resistant to recharge.
Our results are comparable within the uncertainty range with the past studies, most of which were based on the GRACE-derived GWSA obtained after removing the surface water mass changes from recovered TWSA. We also estimated that NCP lost − 5.7 ± 0.8 km3 yr− 1 of GWS from 2005 to 2013. Our model, therefore, can model local and regional variations in GWS over NCP and none of the past studies incorporated such approaches in the region.
Uncertainties still exist in our model, however with improvements in data sources, including GRACE-FO, and the development of efficient ML models such as ensemble models for correcting the residuals, we can expect better performances. Uncertainties in the quantification of GWS changes due to lack of sufficient information about storage coefficient can be offset by using GPS or InSAR data from which we can derive the long-term vertical deformation and combine with the long-term groundwater level changes from this study to obtain reliable estimates of storage coefficient. Notwithstanding these suggested improvements, our ML model undeniably can create robust results and presents a significant advancement for studying the local and regional groundwater variations in NCP which can also be replicated for other stresses aquifers.
Acknowledgement
Wei Feng is supported by the National Nature Science Foundation of China (41874095). We declare no other conflict of interest.
A, G., Wahr, J., & Zhong, S. (2013). Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada. Geophysical Journal International, 192(2), 557–572. https://doi.org/10.1093/gji/ggs030
Agarwal, V, Akyilmaz, O., Shum, C. K., Feng, W., Yang, T.-Y., Forootan, E., Syed, Tajdarul H., & Uz, M. (n.d.). Effective Machine Learning-Aided Downscaling of Satellite Gravimetry Estimated Groundwater Level in Central Valley, California. Journal of Hydrology.
Agarwal, Vibhor. (2021). Machine Learning Applications for Downscaling Groundwater Storage Changes Integrating Satellite Gravimetry And Other Observations. The Ohip State University.
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32. https://doi.org/10.1023/A:1010933404324
Cao, G., & Zheng, C. (2016). Signals of short-term climatic periodicities detected in the groundwater of North China Plain. Hydrological Processes, 30(4), 515–533. https://doi.org/10.1002/HYP.10631
Cao, G., Zheng, C., Scanlon, B. R., Liu, J., & Li, W. (2013). Use of flow modeling to assess sustainability of groundwater resources in the North China Plain. Water Resources Research, 49(1), 159–175. https://doi.org/10.1029/2012WR011899
Changming, L., Jingjie, Y., & Kendy, E. (2010). Groundwater Exploitation and Its Impact on the Environment in the North China Plain. Water International, 26(2), 265–272. https://doi.org/10.1080/02508060108686913
Chen, J., Li, J., Zhang, Z., & Ni, S. (2014). Long-term groundwater variations in Northwest India from satellite gravity measurements. Global and Planetary Change, 116, 130–138. https://doi.org/10.1016/j.gloplacha.2014.02.007
Cheng, M. K., & Ries, J. R. (2018). Monthly estimates of C20 from 5 SLR satellites based on GRACE RL06 models. https://podaactools.jpl.nasa.gov/drive/files/allData/grace/docs/TN-11_C20_SLR.txt
Dramsch, J. S. (2020). 70 years of machine learning in geoscience in review. Advances in Geophysics, 61, 1-55. https://doi.org/10.1016/bs.agph.2020.08.002
Famiglietti, J.S., 2014. The global groundwater crisis. Nature Climate Change 2014 4:11 4, 945–948. https://doi.org/10.1038/nclimate2425
Faunt, C. C. (2009). USGS Professional Paper 1766: Groundwater Availability of the Central Valley Aquifer, California. https://pubs.usgs.gov/pp/1766/
Feng, W., Shum, C., Zhong, M., & Pan, Y. (2018). Groundwater Storage Changes in China from Satellite Gravity: An Overview. Remote Sensing, 10(5), 674. https://doi.org/10.3390/rs10050674
Feng, W., Zhong, M., Lemoine, J.-M., Biancale, R., Hsu, H.-T., & Xia, J. (2013). Evaluation of groundwater depletion in North China using the Gravity Recovery and Climate Experiment (GRACE) data and ground-based measurements. Water Resources Research, 49(4), 2110–2118. https://doi.org/10.1002/WRCR.20192
Han, S., Yang, Y., Lei, Y., Tang, C., & Moiwo, J. P. (2008). Seasonal groundwater storage anomaly and vadose zone soil moisture as indicators of precipitation recharge in the piedmont region of Taihang Mountain, North China Plain. Hydrology Research, 39(5–6), 479–495. https://doi.org/10.2166/NH.2008.117
Hengl, T., Jesus, J. M. de, Heuvelink, G. B. M., Gonzalez, M. R., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., & Kempen, B. (2017). SoilGrids250m: Global gridded soil information based on machine learning. PLOS ONE, 12(2), e0169748. https://doi.org/10.1371/JOURNAL.PONE.0169748
Hu, R. L., Yue, Z. Q., Wang, L. C., & Wang, S. J. (2004). Review on current status and challenging issues of land subsidence in China. Engineering Geology, 76(1–2), 65–77. https://doi.org/10.1016/J.ENGGEO.2004.06.006
Hu, Y., Moiwo, J. P., Yang, Y., Han, S., & Yang, Y. (2010). Agricultural water-saving and sustainable groundwater management in Shijiazhuang Irrigation District, North China Plain. Journal of Hydrology, 393(3–4), 219–232. https://doi.org/10.1016/J.JHYDROL.2010.08.017
Huang, Z., Pan, Y., Gong, H., Yeh, P. J.-F., Li, X., Zhou, D., & Zhao, W. (2015). Subregional-scale groundwater depletion detected by GRACE for both shallow and deep aquifers in North China Plain. Geophysical Research Letters, 42(6), 1791–1799. https://doi.org/10.1002/2014GL062498
Jingli, S., Ling, L., Yali, C., & Zhaoji, Z. (2013). Groundwater Flow Simulation and its Application in Groundwater Resource Evaluation in the North China Plain, China. Acta Geologica Sinica - English Edition, 87(1), 243–253. https://doi.org/10.1111/1755-6724.12045
Kusche, J. (2007). Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models. Journal of Geodesy 2007 81:11, 81(11), 733–749. https://doi.org/10.1007/S00190-007-0143-3
Li, X., Li, G., & Zhang, Y. (2014). Identifying Major Factors Affecting Groundwater Change in the North China Plain with Grey Relational Analysis. Water 2014, Vol. 6, Pages 1581-1600, 6(6), 1581–1600. https://doi.org/10.3390/W6061581
Li, X., Ye, S.-Y., Wei, A.-H., Zhou, P.-P., & Wang, L.-H. (2017). Modelling the response of shallow groundwater levels to combined climate and water-diversion scenarios in Beijing-Tianjin-Hebei Plain, China. Hydrogeology Journal 2017 25:6, 25(6), 1733–1744. https://doi.org/10.1007/S10040-017-1574-4
Long, D., Yang, W., Scanlon, B. R., Zhao, J., Liu, D., Burek, P., Pan, Y., You, L., & Wada, Y. (2020). South-to-North Water Diversion stabilizing Beijing’s groundwater levels. Nature Communications 2020 11:1, 11(1), 1–10. https://doi.org/10.1038/s41467-020-17428-6
Milewski, Thomas, Seyoum, & Rasmussen. (2019). Spatial Downscaling of GRACE TWSA Data to Identify Spatiotemporal Groundwater Level Trends in the Upper Floridan Aquifer, Georgia, USA. Remote Sensing, 11(23), 2756. https://doi.org/10.3390/rs11232756
MLR. (2009). Groundwater regime in major cities of China. http://www.mlr.gov.cn/
Moiwo, J. P., Tao, F., & Lu, W. (2013). Analysis of satellite-based and in situ hydro-climatic data depicts water storage depletion in North China Region. Hydrological Processes, 27(7), 1011–1020. https://doi.org/10.1002/hyp.9276
Pan, Y., Zhang, C., Gong, H., Yeh, P. J.-F., Shen, Y., Guo, Y., Huang, Z., & Li, X. (2017). Detection of human-induced evapotranspiration using GRACE satellite observations in the Haihe River basin of China. Geophysical Research Letters, 44(1), 190–199. https://doi.org/10.1002/2016GL071287
Parr, T., Wilson, J. D., & Hamrick, J. (2020). Nonparametric Feature Impact and Importance. ArXiv, 1–17. http://arxiv.org/abs/2006.04750
Rabus, B., Eineder, M., Roth, A., & Bamler, R. (2003). The shuttle radar topography mission—a new class of digital elevation models acquired by spaceborne radar. ISPRS Journal of Photogrammetry and Remote Sensing, 57(4), 241–262. https://doi.org/10.1016/S0924-2716(02)00124-7
Rahaman, M., Thakur, B., Kalra, A., Li, R., & Maheshwari, P. (2019). Estimating High-Resolution Groundwater Storage from GRACE: A Random Forest Approach. Environments, 6(6), 63. https://doi.org/10.3390/environments6060063
Rodell, M., Famiglietti, J. S., Chen, J., Seneviratne, S. I., Viterbo, P., Holl, S., & Wilson, C. R. (2004). Basin scale estimates of evapotranspiration using GRACE and other observations. Geophysical Research Letters, 31(20), 20504. https://doi.org/10.1029/2004GL020873
Sakura, Y., C, T., Yoshioka, R., & Ishibashi, H. (2003). Intensive use of groundwater in some areas of China and Japan. In M. Llamas & E. Custodio (Eds.), Intensive Use of Groundwater: Challenges and Opportunities.
Save, H. (2019). CSR Grace RL06 Mascon Solutions. Texas Data Repository Dataverse, 1.
Save, H., Bettadpur, S., Tapley, B.D., 2016. High-resolution CSR GRACE RL05 mascons. J. Geophys. Res.-Sol. Earth 121 (10), 7547–7569.
Shu, Y., Villholth, K. G., Jensen, K. H., Stisen, S., & Lei, Y. (2012). Integrated hydrological modeling of the North China Plain: Options for sustainable groundwater use in the alluvial plain of Mt. Taihang. Journal of Hydrology, 464–465, 79–93. https://doi.org/10.1016/J.JHYDROL.2012.06.048
Siebert, S., Henrich, V., Frenken, K., & Burke, J. (2010). Update of the Global Map of Irrigation Areas to version 5. http://www.fao.org/3/I9261EN/i9261en.pdf
Smith, R., Knight, R., & Fendorf, S. (2018). Overpumping leads to California groundwater arsenic threat. Nature Communications, 9(1), 1–6. https://doi.org/10.1038/s41467-018-04475-3
Snoek, J., Larochelle, H., & Adams, R. P. (2012). Practical Bayesian optimization of machine learning algorithms. Advances in Neural Information Processing Systems, 4, 2951–2959. https://arxiv.org/abs/1206.2944v2
Swenson, S., Chambers, D., & Wahr, J. (2008). Estimating geocenter variations from a combination of GRACE and ocean model output. Journal of Geophysical Research: Solid Earth, 113(B8), 8410. https://doi.org/10.1029/2007JB005338
Tamanyu, S., Muraoka, H., & Ishii, T. (2009). Geological interpretation of groundwater level lowering in the North China Plain. Bulletin of the Geological Survey of Japan, 60(1–2), 105–115. https://doi.org/10.9795/BULLGSJ.60.105
Wada, Y., van Beek, L.P.H., van Kempen, C.M., Reckman, J.W.T.M., Vasak, S., Bierkens, M.F.P., 2010. Global depletion of groundwater resources. Geophysical Research Letters 911 37. https://doi.org/10.1029/2010GL044571
Wada, Y., Wisser, D., Bierkens, M.F.P., 2014. Global modeling of withdrawal, allocation and consumptive use of surface water and groundwater resources. Earth System Dynamics 5, 914 15–40. https://doi.org/10.5194/ESD-5-15-2014
Wang, J., Huang, J., Blanke, A., Huang, Q. Q., & Rozelle, S. (2007). The agricultural groundwater revolution: opportunities and threats to development. In M. Giordano & K. G. Villholth (Eds.), The agricultural groundwater revolution: opportunities and threats to development.
Wang, L., Chen, C., Du, J., & Wang, T. (2017). Detecting seasonal and long-term vertical displacement in the North China Plain using GRACE and GPS. Hydrology and Earth System Sciences, 21(6), 2905–2922. https://doi.org/10.5194/HESS-21-2905-2017
Wang, S., Song, X., Wang, Q., Xiao, G., Liu, C., & Liu, J. (2009). Shallow groundwater dynamics in North China Plain. Journal of Geographical Sciences, 19(2), 175–188. https://doi.org/10.1007/S11442-009-0175-0
Yin, W., Han, S., Zheng, W., Yeo, I., Hu, L., Tangdamrongsub, N., & Ghobadi-Far, K. (2020). Improved water storage estimates within the North China Plain by assimilating GRACE data into the CABLE model. Journal of Hydrology, 590, 125348. https://doi.org/10.1016/j.jhydrol.2020.125348
Zhang, G., Fei, Y., Wang, J., Yan, M., Yu, C., Wang, G., Wang, H., Feng, H., Liu, Z., & Liu, K. (2012). Adaptation between Irrigation Agriculture and Groundwater in North China Plain.
Zhang, G., Zheng, W., Yin, W., & Lei, W. (2021). Improving the Resolution and Accuracy of Groundwater Level Anomalies Using the Machine Learning-Based Fusion Model in the North China Plain. Sensors, 21(46). https://doi.org/10.3390/s21
Zhang, J., Liu, K., & Wang, M. (2021). Downscaling Groundwater Storage Data in China to a 1-km Resolution Using Machine Learning Methods. Remote Sensing 2021, Vol. 13, Page 523, 13(3), 523. https://doi.org/10.3390/RS13030523
Zhang, Y., Schaap, M. G., & Zha, Y. (2018). A High‐Resolution Global Map of Soil Hydraulic Properties Produced by a Hierarchical Parameterization of a Physically Based Water Retention Model. Water Resources Research, 54(12), 9774–9790. https://doi.org/10.1029/2018WR023539
Zhang, Z., Fei, Y., Chen, Z., Zhao, Z., Xie, Z., Wang, Y., Miao, J., Yang, L., Shao, J., & Jin, M. (2009). Investigation and Assessment of Sustainable Utilization of Groundwater Resources in the North China Plain. Geological Publishing House: Beijing, China.
Zheng, C., Liu, J., Cao, G., Kendy, E., Wang, H., & Jia, Y. (2010). Can China Cope with Its Water Crisis?—Perspectives from the North China Plain. Groundwater, 48(3), 350–354. https://doi.org/10.1111/J.1745-6584.2010.00695_3.X
Zhou, C., Gong, H., Chen, B., Zhu, F., Duan, G., Gao, M., & Lu, W. (2016). Land subsidence under different land use in the eastern Beijing plain, China 2005-2013 revealed by InSAR timeseries analysis. GIScience & Remote Sensing, 53(6), 671–688. https://doi.org/10.1080/15481603.2016.1227297