Simulating the Shallow Groundwater Level Response to Artiﬁcial Recharge and Storage in the Plain Area of the Daqing River Basin, China

: Water shortage and overexploitation of groundwater (GW) have become the key factors restricting the development of the Xiongan New Area (XNA), the environmental management of Baiyangdian Lake, and the social and economic development of surrounding areas. This study used a numerical GW ﬂow model to quantitatively analyze the changes to the shallow GW level and GW reserves of the plain area of the Daqing River Basin over the next 15 years (2021–2035) under different artiﬁcial recharge schemes with the south to north water diversion project (SNWDP) acting as the GW recharge source. The results showed increasing GW storage and GW levels and that the regional GW resources are in a positive equilibrium state. The rates of change of the well irrigation supply scheme and the joint river-well irrigation supply scheme in the XNA will increase by 14.56% and 11.04% by 2035 as compared with the current situation. The well irrigation recharge scheme for the XNA was determined to be the most effective when comparing with the effects of the different artiﬁcial recharge schemes on the GW levels and recharge. This study provides a reference for the management and protection of aquifers in other areas suffering serious GW overexploitation. The results showed that all the recharge schemes resulted in increases to the shallow GW level in the depression cone area of Xiongxian, among which well GW recharge had the largest effect followed by combined river–well GW recharge. This result showed that all the recharge schemes are helpful for the recovery of the GW level in the Xiongxian depression cone area.


Introduction
Groundwater (GW) is an essential source of water for the environment, domestic water supply, agriculture, and industry due to its generally good quality and widespread occurrence [1]. While GW has supported rapid socioeconomic development in China and globally, overexploitation of groundwater has caused many environmental and geological problems. The negative environmental impacts of GW depletion from overexploitation are well known, and include deteriorating GW quality, ecosystem degradation, and in some cases, land subsidence and/or seawater intrusion [2][3][4][5]. The excessive depletion of GW over the last 50 years has become a global problem affecting major regions of North Africa, the Middle East, South and Central Asia, North China, North America, and Australia [6][7][8].
The North China Plain (NCP) is recognized to be one of the areas experiencing serious GW depletion [9,10]. The Chinese government established the Xiongan New Area (XNA) on the 1st April, 2017. The establishment of the XNA is of great practical and far-reaching historical significance [11]. However, achieving sustainable development in the XNA is currently being hindered by water deficits, severe water pollution, and overexploitation of GW. Although the GW system has a certain natural regulation ability, this has been increasingly degraded under the strong influence of human activities. Therefore, artificial GW storage has become an important measure implemented to prevent the continuous shortage in GW resources and to maintain the vitality of the GW system. Artificial recharge can be defined as the artificial regulation and storage of water resources in space and time to ensure their full and reasonable utilization. An important prerequisite for artificial GW replenishment is the availability of a rechargeable source of abundant and good water quality. Water for GW replenishment in the plain area of the Daqing River Basin mainly includes regulated water from the middle route of the south to north water diversion project (SNWDP), reservoir-regulated water, and rainfall runoff. The middle route of the SNWDP has provided ecological water replenishment to Hebei Province through the Fu, Bao, Beiyishui, Nanjuma, Hutuo, and other rivers since September, 2018. To date, the total volume of replenished water has exceeded 1.2 billion m 3 [12]. The Wangkuai, Xidayang, Angozhuang, and Longmen reservoirs have been built in the plain area of the Daqing River Basin since the 1950s. In recent years, each reservoir has continuously replenished water through each river to the plain area of the Daqing River Basin to meet ecological requirements. Approximately 420 million m 3 of regulated water from each reservoir enters the plain area of the Daqing River Basin every year [12]. Therefore, the study area is characterized by conditions conducive to artificial GW recharge. The Chinese government has issued a series of policies to reduce GW exploitation and increase recharge to solve the problem of continuous GW overexploitation. These developments are likely to lead to complex changes in the GW level in this area.
The numerical simulation of GW is the main technical method used to study the response of GW processes system to GW regulation and storage measures. GW flow models can assist in the formulation of GW regulation and storage management policies [13]. Numerical GW models are clearly appropriate for the assessment of sustainable GW use given the complexity and heterogeneity of aquifer systems and the ability of numerical modeling to integrate different data and to evaluate regional GW dynamics [14][15][16]. Several numerical models have been constructed to evaluate the GW resources and flow dynamics in the NCP. Zhang et al. [17] simulated and analyzed the initial state of GW exploitation in plain area of the Haihe River Basin, with their simulation results providing a basis for water resource management and restoration of the water environment. Wang et al. [18] evaluated the GW depletion in the NCP from January 2002 to December 2003 with their results estimating a total inflow and outflow of 49.4 km 3 and 56.5 km 3 , respectively, for a budget deficit of 7.1 km 3 . Cui et al. [19] used GW modeling to evaluate the effects of reductions in GW extraction over a ten-year period in response to the additional water provided by the SNWDP. Their study concluded that the SNWDP plays an important role in the recovery of GW in the NCP. Xue et al. [20] use numerical simulation to comprehensively analyze the change in the GW system and the impact of long-term GW exploitation of the shallow aquifer in the plain area of the Haihe River Basin. Liu et al. [21] simulated GW flow dynamics from 1960 to 2004 to evaluate the effects of urbanization of rural areas in the vicinity of Shijiazhuang City and the SNWDP. Hu et al. [22] integrated a crop-growth model and a GW model to evaluate the effects of a reduction in irrigation in Shijiazhuang, with the results of their study showing that a reduction in irrigation of 140 mm/yr could reverse GW depletion in this area.
The implementation of policies mandating GW artificial recharge and GW compression mining will result in changes to the GW-dependent water supply of the NCP. These changes will relieve water deficits to some extent, thereby allowing an opportunity for the recovery of the GW system. The present study established a GW flow numerical model to evaluate the recovery of the GW level and the change in GW reserves. Artificial recharge of GW can directly and effectively increase GW resources, regulate and store surface water, prevent or control the decline in GW levels, control the formation and further expansion of depression cones, and prevent or control land subsidence and other environmental geological problems. It is hoped that the results of the present study can guide a suitable GW recharge and storage scheme to alleviate the environmental risk resulting from the continuous decline of the GW level.

Hydrogeological Settings
The plain area of the Daqing River Basin is in the piedmont plain of the Taihang Mountains and covers an area of 13,338 km 2 ( Figure 1). The elevation increases slightly from northwest to southeast, the terrain is flat, and the surface slope is 1-4‰. The Daqinghe River Basin falls within a temperate semihumid and semiarid continental monsoon climate zone with four distinct seasons. The annual average temperature of the region is 12.3 • C (2015-2018a) [12]. Precipitation in the plain area shows both intra-and interannual heterogeneity and a regional distribution with an annual average of 526 mm (2015-2018a) [12]. Annual average evaporation in the plain area is 1309.7 mm (2015-2018a) [12]. The stratigraphy of the study area consists mainly of unconsolidated quaternary sediments. These sediments contain most of the GW [23]. Previous studies [24][25][26] have shown that the aquifer structure varies from a single aquifer composed of gravels and pebbles in the upper parts of the piedmont fan in the west to multiple aquifers composed of sand, silt, and clay in the east. The GW system can be divided into four aquifer groups based on this stratigraphy (I-IV). Aquifer groups I and II represent shallow aquifers, whereas groups III and IV represent deep aquifers [27] (Figure 2). Aquifer groups I (10-50 m depth) and II (≤210 m) represent unconfined and semiconfined aquifers, respectively, and are exploited for irrigation and ecosystem services. Aquifer groups III (≤310 m) and IV (>400 m) represent confined aquifers. The shallow GW system is characterized by a faster water cycle as compared to the deep GW system, and mainly includes the aquifer groups I and II. GW exploitation in the study area is mainly of the shallow GW system. Overexploitation of the shallow GW system has resulted in water shortages, GW fall cones, and a series of environmental and geological problems. Artificial GW recharge is urgently required to assist in the recovery of this aquifer. Therefore, the present study concentrated on the shallow aquifer system under aquifer groups I and II. The deep GW system (aquifer groups III and IV) mainly consists of confined GW. This aquifer group primarily receives recharge from lateral runoff and shallow water. Artificial exploitation is currently the main discharge route for the deep GW system. GW recharge processes are mainly precipitation infiltration, lateral flow in the mountains, and irrigation return flow. GW discharge processes are mainly phreatic evaporation and human GW exploitation. However, GW exploitation in the central and eastern parts of the plain results in flow of GW into the center of depression cones and in the flow of shallow GW into the deep aquifers.
The study area has a per-capita water resources availability of 282 m 3 /a, which is only 7% of the average per-capita water resources availability in China. There is a high reliance on GW in the study area to support agriculture, industry, and domestic water use. Shallow GW is mainly exploited, accounting for 86.33% of total GW exploitation. Although a series of measures have been introduced in Hebei Province since 2015 to reduce GW exploitation, agricultural production activities in the study area remain reliant on GW.
The middle route of the SNWDP has been an important measure to alleviate water deficits and the deterioration of the ecological environment in the plains of the Daqing River Basin. The middle route of the SNWDP was officially opened in December, 2014 [28]. Under current water transfer levels, the SNWDP can be used to artificially replenish 160 million m 3 of GW per year in the plains of the Daqing River Basin.

Changes in the GW Level
Interannual fluctuations in the shallow GW level are directly affected by atmospheric precipitation and artificial exploitation. Changes to GW levels in the study area over time can be roughly categorized into five periods: (1) the natural flow stage prior to 1964; (2) a stage of a localized declines in GW levels (1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972); (3) a stage of rapid decline in regional GW levels (1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985); (4) a stage of fluctuating declines in GW level  and; (5) a stage of gradual GW level rise post 2015 [1]. Shallow GW in the study area showed a general decline from 1982 to 2015 ( Figure 3). While the rate of decline was small in the piedmont area and the area surrounding Baiyangdian Lake, larger declines in GW were evident in the Wangdu Baoding Gaoyang and Boye Lixian areas due to the large-scale exploitation of GW. GW levels in areas of large-scale exploitation dropped dramatically, resulting in the formation of depression cones. The regional GW level has gradually increased since 2015 due to the implementation of a policy restricting GW exploitation and promoting artificial recharge ( Figure 4).   Figure 1 shows a map of the plain area of the Daqing River Basin. A model was constructed to represent the complex multilayer aquifer system, including two aquifer groups and one aquitard. Aquifer groups I and II were combined into one aquifer group representing the shallow aquifer system due to the presence of several incomplete aquitards. These incomplete aquitards showed enhanced hydraulic connections due to the presence of many pumping wells. Aquifer group III was taken to represent the deep aquifer system. Aquitards were located between the aquifer groups. These aquifers were either characterized by strong horizontal flow due to their wide distribution and thickness or slight flow in the vertical direction. The recharge and discharge of this area changes over time along with the GW level, showing features of transient flow. At the same time, GW parameters show spatial variation, representing heterogeneous flow. As a result, the GW system of the Daqing River Basin could be described using a conceptual hydrologic model comprising three tiers and a heterogeneous, three-dimensional transient flow system with horizontal isotropy.

Conceptual Model
The mountains form the western boundary of the plain study area and also act to provide lateral GW recharge; therefore, the mountains are treated as a flow boundary in the model. The study area falls within the Daqing River Basin, with the eastern part of the study area bounded by the Daqing and Zhongting rivers and the northeast by the Mangniu River. The southern, northern, and eastern parts of the model were based on the surface watershed obtained from hydrological analysis conducted in the ArcGIS software. Since there are two possibilities for inflow and outflow of GW at the boundary, these boundaries were treated as general head boundaries in the model. The GW table was regarded as the top boundary in the model within which exchange between the GW and the other systems occurred through infiltration of precipitation and evaporation of phreatic water. The bottom of the system was regarded as impermeable due to the presence of loam and clay beneath the quaternary strata.
Since infiltration of precipitation is the dominant source of GW recharge in the study area, GW recharge is affected by the amount of precipitation and the infiltration coefficient. Within the model, the monthly precipitation was evenly distributed across the study area, whereas the infiltration coefficient differed spatially according to six zones, with the coefficient ranging from 0.12 to 0.3 [29]. Although lateral flow from the west is a significant source of recharge in the study area, it is poorly quantified [30]. The recharge rates originating from lateral flow used in the model were based on previous studies in which they were estimated by the profile method [31] and hydrological analysis [32,33]. GW is the primary source of water for irrigation in the study area, accounting for 90% of the total irrigation water [33,34]. The infiltration coefficient of irrigation used in the model was less than that of precipitation considering land use since evaporation is higher during the irrigation season. The irrigation infiltration used in the present study was distributed in 42 zones based on the monthly irrigation and infiltration coefficients. Groundwater recharge through surface water infiltration decreased after the 1960s in the study area, and by the 1980s, the rivers had nearly dried up [20,35]. Before the establishment of the XNA, river GW recharge in the basin mainly occurred through infiltration recharge from the Nanjuma River. Ephemeral streams only influence GW recharge rates in the study area during the flood season. The recharge from ephemeral streams during the flood season was incorporated into total recharge in the model. Most of the rivers in the study area have dried up or have become ephemeral streams during the flood season due to decreasing precipitation and upstream reservoirs retaining runoff [35]. Under the current situation, seepage from Baiyangdian Lake recharges the surrounding GW, and the average lake level in 2017 used in the model was 6.5 m.
The current study generated the evapotranspiration rate by Thiessen polygons using daily potential evaporation data measured by four meteorological stations. These data were provided by the China Meteorological Data Sharing Service system and were based on the phreatic evaporation coefficient, which is decided by the characteristics of the vadose zone. The extinction depth of the entire modelled study area was set to 4 m [35]. Artificial discharge from the aquifers in the study area occurs through well pumping [35,36]. Data for GW exploitation were obtained from the hydrological bureau of each municipality in the study area and previous field research [35][36][37][38][39]. Although detailed information on the locations and depths of the wells was not available, the layers of exploitation were known.

Numerical Model
The GW flow model was developed using MODFLOW [40]. MODFLOW includes a set of stress packages that allows the simulation of external flow stresses, such as wells, areal recharge, evapotranspiration, drains, and rivers [41,42].
The simulation was divided horizontally into 100 rows and 100 columns with a cell size of 1.5 km × 1.5 km. The model was divided into three vertical layers based on the hydrogeological conditions and data on the porous aquifers at the study site. Layer 1 included the aquifer groups I and II, layer 3 represented aquifer group III, and layer 2 represented the aquitards. The simulation period extended from December 2016 to December 2017 due to data limitations. Each calendar month represented a stress period with constant source and sink terms. Contour maps for December 2016 were used as the initial head.
All sources and sinks were input using MODFLOW packages such as well (WEL), recharge (RCH), and evapotranspiration (EVT). Hydrogeological parameters such as the hydraulic conductivity, specific yield, and storage coefficient were imported by subarea [43].

Model Calibration
The model calibration period extended from December 2016 to June 2017, whereas model validation extended from July 2017 to December 2017. The calibration process was divided into two parts: (1) fitting of the observed and simulated GW flow fields and (2) comparison between the observed and simulated hydrographs at typical observation wells [31]. The simulated flow fields were similar to the observed water level contours in 2017 for the shallow and deep aquifers (layers 1 and 3; Figure 1). A total of 45 observation wells were selected for comparison of the GW hydrographs. Figure 5 shows the goodnessof-fit between simulated and observed water levels at all monitoring wells. Figure 6 shows a comparison between the simulated and observed water levels for five observation. These results showed that the model simulations provided a good representation of the actual GW flow field and the features of the GW regime.   Table 1).   An equilibrium analysis of the transient GW flow model was used to compute the GW equilibrium in the Daqing River Basin from December 2016 to December 2017. The recharge value for the entire GW system was 1.858 billion m 3 and the recharge-discharge shortfall was 1.008 billion m 3 . These results suggest that the GW discharge exceeded recharge for several years, resulting in the continuous reduction in GW storage.

Artificial Recharge Scenarios and the Response of the GW Level
The calibrated and verified model was used to investigate scenarios to predict and evaluate GW levels for the next 15 years (2021-2035). The present study only analyzed the shallow GW levels as this GW layer was the focus of the current study. The study area was in a significantly negative equilibrium state under the current conditions, which could be attributed to the local overexploitation of GW, resulting in the formation of several GW depression cones. The regional GW main recharge and discharge processes were precipitation recharge and artificial exploitation, accounting for 70% and 90% of total recharge and discharge, respectively. This result indicated the need to conduct artificial GW recharge through the use of river channels and recharge wells in the study area while also regulating GW exploitation. Under this approach, water consumed by exploitation can be counterbalanced by recharge to the aquifer. Four scenarios were developed based on these assumptions.

Scenario 1: River Recharge across the Entire Study Area
The action plan for comprehensive treatment of GW overexploitation in North China issued by the Ministry of Water Resources in January 2019 sets a target plan for future GW replenishment by major rivers. Therefore, the Zhulong, Tang, Baohe, and Nanjuma rivers will be recharged in the future. These rivers will become a major supply channel for GW resources and for surface water resources in Baiyangdian Lake. The comprehensive river recharge scheme under the pressure mining conditions in the study area is designed to provide each of the Nanjuma, Bao, Tang, and Zhulong rivers with 100 million m 3 /a of recharge water.

Scenario 2: River Recharge in the XNA
The average GW recharge for each river was set with reference to ecological water recharge and GW recharge in the implementation plan of the Ministry of Water Resources. The focus of the model prediction was on the GW recharge schemes for the Nanjuma and Baigou rivers flowing through the main GW level cone within the XNA. Since the water source of the middle route of the SNWDP mainly recharges the Nanjuma river, the Nanjuma River and Baigou Diversion River were designed to facilitate GW recharge of 40 million m 3 /a and 10 million m 3 /a, respectively on the basis of natural river infiltration.

Scenario 3: Well Recharge in the XNA
A number of recharge wells were set within 1 km of the Baigou Diversion River according to the regional GW flow field and GW level cone distribution characteristics. The water source of the middle route of the SNWDP was used to recharge the regional GW. The focus of the scenario was on the recovery of the GW table in the district depression cone in Xiongxian and Rongcheng. The evaluation of the regional suitability showed that natural infiltration near the Baigou Diversion River is poor. Therefore, the pressurized method was selected for injection of recharge water within the reinjection wells. Recharge water was injected into the phreatic aquifer. Artificial recharge of 50 million m 3 /a was implemented and single-well recharge of 5000 m 3 /d was required for a minimum of 27 recharge wells.

Scenario 4: Combined Supply of River and Well GW Recharge in the XNA
A combined supply of river and well recharge (river-well GW recharge) was established in the XNA according to the regional GW flow field and GW level cone distribution characteristics. A river recharge volume of 10 million m 3 /a and well recharge volume of 40 million m 3 /a were implemented. The pressurized method was used for artificial recharge at the recharge wells due to the poor natural infiltration near the Baigou Diversion River. The recharge water was injected into the phreatic aquifer. A single well recharge rate of 5000 m 3 /d was implemented and a minimum of 22 recharge wells were required.

Changes in the GW Level
A comparison of the predicted flow field in 2035 under various recharge schemes with that for the same year under current conditions showed an upward trend in the GW flow field in the study area ( Figure 8). The trend in the GW flow field was consistent across the Western Piedmont area, the flow field in the XNA showed obvious changes, the area of the depression cone decreased, and there was an obvious rise in the GW level of 20 m. The GW level in the Rongcheng depression cone area increased by 11 m whereas that in the Xiongxian depression cone area increased by 4 m. The downward trend in the GW flow field had a certain inhibitory effect. These results showed that the different recharge schemes adopted in the Baigou Diversion River and within 1 km of the river could effectively act to recharge the GW level of the depression cone area on the eastern and western sides of the river. The storage condition of aquifer is the key factor affecting the effect of groundwater recharge. The storage condition of the aquifer is characterized by two parts. One is that the target aquifer should have enough space to store recharge water, that is, the size of underground water storage space. The other is that the groundwater can be rapidly distributed in the whole water storage space without forming a groundwater mound after entering the aquifer, that is, the water transmitting ability of the aquifer.
The results showed that all the recharge schemes will result in a rise in the shallow GW level of the Rongcheng depression cone area, with the most obvious rise in the GW level under the pressure extraction of the GW scheme. Similar increases in the shallow groundwater levels were achieved under well GW recharge and river-well GW recharge schemes in the XNA. This result could be attributed to the ratio of the river-well GW recharge to well GW recharge being 1:4. Therefore, river GW recharge across the entire study area is helpful for facilitating the recovery of the GW level in the Rongcheng depression cone area.
The results showed that all the recharge schemes resulted in increases to the shallow GW level in the depression cone area of Xiongxian, among which well GW recharge had the largest effect followed by combined river-well GW recharge. This result showed that all the recharge schemes are helpful for the recovery of the GW level in the Xiongxian depression cone area.

Changes to Supply
Four artificial GW recharge schemes were implemented in the study area within the case of pressure extraction. The results showed that the regional GW resources entered a positive equilibrium state ( Table 2). The relative changes in well GW recharge in the XNA, combined river-well GW recharge in the XNA, river GW recharge across the entire area, and river recharge in the XNA were 14.56%, 11.04%, −0.12%, and −3.18%, respectively. A comprehensive comparison of various recharge schemes from the perspective of the overall effect of artificial recharge on GW level and the increase in recharge showed that the well irrigation recharge scheme for the XNA was relatively reasonable.
The results indicated that the influence of the SNWDP helped to control the decrease in the GW levels. The results of the present study are consistent with those of other studies [1,44,45]. The present study can improve the understanding of the impact of the SNWDP on GW levels in the plain area of the Daqing River Basin. Beijing has less than 100 m 3 per capita water resources within its administrative area, 1/20 of the national average, 1/80 of the global average, and far below the internationally recognized danger level of 1700 m 3 and absolute scarcity of less than 500 m, the degree of groundwater depletion is similar to that in plain area of Daqing River Basin. Long et al. used physical based hydrological models to quantify the impact of SNWD in the context of the effects of climate extremes and water policies (e.g., reduction in irrigation pumpage) using physically based hydrologic modeling, and to project how these factors may influence GWS in the future. The results show that water diverted to Beijing reduced cumulative GW depletion by 3.6 km 3 , accounting for 40% of total GW storage recovery during 2006-2018. Increased precipitation contributes similar volumes to GW storage recovery of 2.7 km 3 (30%) along with policies on reduced irrigation (2.8 km 3 , 30%) [14]. The conclusion is similar to that of this study.

Conclusions
A conceptual GW flow model was developed by analyzing the hydrogeological conditions in the Daqing River Basin. A three-dimensional transient numerical flow model was then developed using MODFLOW. The hydrogeologic parameters were identified and the model was calibrated by fitting the simulated GW flow field with the hydrographs of observation wells. The model was then used to simulate the response of the GW levels to different artificial recharge schemes. The GW levels were simulated and analyzed under four scenarios. When compared with the current situation, the GW level increased to a certain extent in 2035 after the implementation of four artificial recharge schemes. The order of the different recharge schemes according to the relative rate of change in GW was well GW recharge in the XNA (14.56%) > combined river and well GW recharge in the XNA (11.04%) > river GW recharge across the entire area (−0.12%) > river GW recharge in the XNA (−3.18%). A comprehensive comparison of various recharge schemes from the perspective of the overall effect of artificial recharge on the GW levels and the increase in recharge showed that the well GW recharge scheme for the XNA was relatively reasonable. Therefore, the implementation of artificial GW recharge can alleviate the environmental risk resulting from the continuous decline in the GW level in the plain area of the Daqing River Basin. This study provides a reference for the management and protection of aquifers in other areas suffering serious GW overexploitation.
The present study considered the impact of artificial recharge on GW under only one climate scenario. The interactions between climate change, multisource recharge, and socioeconomic activity for integration into a GW flow model requires further research.