Groundwater Recharge Modeling under Water Diversion Engineering: A Case Study in Beijing

: The of surface water resource exploitation and utilization projects on groundwater has been widely studied. Surface water diversion projects lead to a reduction in river discharge, which changes the recharge of groundwater systems. In this study, the numerical simulation method is used to predict the variation in groundwater level under different diversion scale scenarios. The Zhangfang water diversion project in Beijing, China, was chosen for the case study. The downstream plain area of the Zhangfang water diversion project is modeled by MODFLOW to predict the inﬂuence of reducing water diversion on the dynamic change in the groundwater level in the downstream plain area. The model results show that the difference in groundwater recharge and discharge on the downstream plain of Zhangfang is 9,991,900 m 3 /a, which is in a negative water balance state, and the groundwater level continues to decrease. Reducing the amount of water diverted by the Zhangfang water diversion project to replenish groundwater is beneﬁcial to the rise of the groundwater level in the downstream plain area. The results indicate that the groundwater ﬂow model in the downstream plain area of Zhangfang performed well in the inﬂuence assessment of surface water resource exploitation and utilization projects on groundwater. This study also provides a good example of how to coordinate the relationship between surface water resources and groundwater resources.


Introduction
Groundwater plays an important role in social development and in maintaining the ecological environment [1]. River leakage recharge is an important part of groundwater recharge, so the interaction between the river and aquifer has become a research hotspot [2][3][4]. The interaction between surface water and groundwater is constantly changing due to the different spatial and temporal distribution of the exchange between surface water and groundwater [5,6]. A water diversion project will change the interaction between the surface water and groundwater, and directly reduce the discharge of the river channel, which is one of the man-made factors causing the change of regional groundwater level and quality [7,8]. Among them, there are many studies on the impact of external water transfer in the water receiving area, including the impact of the South-to-North water transfer project on the groundwater systems of the North China plain [9] and Beijing [10], the impact of watershed ecological water transfer on the hydrological situation in the lower reaches of the Heihe River basin [11], and the impact of the ecological water supplement on the section of Yongding River, Beijing plain [12]. Research shows that water transfer is conducive to the recovery of the groundwater level in the water receiving area. However, there have been few studies on the impact in the water transfer area. The water transfer project directly leads to the reduction of groundwater discharge in the water transfer area [13], resulting in a decline of the downstream groundwater level [14].
The Zhangfang water diversion project is located in the Fangshan district of Beijing, supplying water from the Juma River to Beijing through the Shengtian Canal. In September 1981, Beijing decided to divert water from the Juma River to relieve the pressure of a water shortage. In 2004, the Zhangfang water diversion project was completed, and supplied water to Beijing. Since the construction of the project, the Zhangfang water diversion project has diverted approximately 100 million m 3 of water from the Juma River every year. The large amount of water diversion has reduced the surface runoff of the Nanjuma River and the Beijuma River downstream of the Zhangfang water diversion project, which has caused a serious impact on downstream water resources and the water environment. Therefore, research on the influence of the Zhangfang water diversion project on downstream water resources is conducive to rational water diversion projects and can enrich the research on the influence of water diversion on the water diversion area.
The water diversion project inevitably leads to changes in hydrogeological conditions along the river and changes the supplement and discharge relationship between downstream groundwater and the surface. A variety of methods can be used to quantitatively describe the process: the water balance method [15], Darcy law [16], the groundwater level dynamic method [17], the isotope method [18] and the numerical simulation method [19]. Numerical simulation achieves the purpose of simulating the actual groundwater system by solving the mathematical model. It is one of the most popular methods for solving groundwater-related problems and has been widely used in groundwater resource evaluation [20], the evaluation of the impact of groundwater development [21], the impact of water conservancy projects on groundwater [22], the mutual conversion of groundwater and surface water, etc. Popular groundwater flow numerical simulation codes include the finite element groundwater flow system FEFLOW [23] and modular 3D groundwater flow model MODFLOW based on the finite difference method [24]. MODFLOW-2005 [25] is a recent popular core version of MODFLOW, which uses a block-centered finite difference method to simulate confined and unconfined groundwater flow, as well as external stresses that include regional recharge, pumping wells, evapotranspiration, rivers, etc.
In this study, we use MODFLOW-2005 as a solver to construct a transient threedimensional groundwater flow numerical model in the downstream plain of Zhangfang and study the current groundwater flow of this area. Based on this model, different water diversion amounts from the Zhangfang water diversion project are designed to evaluate the impact of the Zhangfang water diversion project on the groundwater level in the downstream plain area. Additionally, this study reveals the importance of river infiltration for groundwater recharge and proposes suggestions for the rational development and utilization of downstream water resources.

Meteorology of the Study Area
The study area is located between 39 • 35 57" and 39 • 4 50" east longitude and 115 • 30 40" and 116 • 15 26" north latitude. The total area is approximately 2603 km 2 . The research scope is shown in Figure 1. The climate in this area is a temperate continental monsoon climate. The perennial mean temperature is 13.4°C, the winter is severely cold from December to January, and the summer is hot from July to August. The mean annual precipitation is 489.9 mm.

Hydrogeological Conditions in the Study Area
The study area is located on multistage geological faults in front of the alluvialproluvial fan of the Juma River. The main fault structures are the Huairou-Laishui deep fault and Donglutou-Qiaoliufan fault. West of the Huairou-Laishui deep fault, quaternary deposits directly cover Proterozoic carbonate rock, and the buried depth of the bedrock is generally 10-40 m. East of the fault to the line of the Donglutou-Qiaoliufan fault, underlying the quaternary strata are tertiary cemented conglomerate and glutenite, and the buried depth of the bedrock is generally 20-60 m. To the east of the Donglutou-Qiaoliufan fault line, quaternary deposits suddenly thicken to 170-550 m. The hydrogeological profile in the study area is shown in Figure 2b.

Hydrogeological Conditions in the Study Area
The study area is located on multistage geological faults in front of the alluvial-proluvial fan of the Juma River. The main fault structures are the Huairou-Laishui deep fault and Donglutou-Qiaoliufan fault. West of the Huairou-Laishui deep fault, quaternary deposits directly cover Proterozoic carbonate rock, and the buried depth of the bedrock is generally 10-40 m. East of the fault to the line of the Donglutou-Qiaoliufan fault, underlying the quaternary strata are tertiary cemented conglomerate and glutenite, and the buried depth of the bedrock is generally 20-60 m. To the east of the Donglutou-Qiaoliufan fault line, quaternary deposits suddenly thicken to 170-550 m. The hydrogeological profile in the study area is shown in Figure 2b.
The aquifer is mainly recharged by means of lateral runoff in the piedmont zone, and precipitation infiltration and water infiltration of the Nanjuma River and Beijuma River; it is discharged by means of artificial mining and lateral outflow of the aquifer. Affected by topography, groundwater flow in the study area mainly flows from northwest to southeast. In recent years, as a result of the large scale of water through the Zhangfang water diversion project and human exploitation of groundwater, the aquifer has almost drained.

The Water Diversion Information of the Zhangfang Water Diversion Project
The Zhangfang water diversion project supplies water to Beijing through the Shengtian Canal, with an annual flow of approximately 100 million m 3 from the Juma River. The water diversion amounts can be estimated through the observation data of hydrologic stations. The location of three hydrological stations on the Juma River is shown in Figure  1. Zhangfang hydrological station is located downstream of the Duya hydrological sta-  The aquifer is mainly recharged by means of lateral runoff in the piedmont zone, and precipitation infiltration and water infiltration of the Nanjuma River and Beijuma River; it is discharged by means of artificial mining and lateral outflow of the aquifer. Affected by topography, groundwater flow in the study area mainly flows from northwest to southeast. In recent years, as a result of the large scale of water through the Zhangfang water diversion project and human exploitation of groundwater, the aquifer has almost drained.

The Water Diversion Information of the Zhangfang Water Diversion Project
The Zhangfang water diversion project supplies water to Beijing through the Shengtian Canal, with an annual flow of approximately 100 million m 3 from the Juma River. The water diversion amounts can be estimated through the observation data of hydrologic stations. The location of three hydrological stations on the Juma River is shown in Figure 1. Zhangfang hydrological station is located downstream of the Duya hydrological station, with Shengtian Canal between them. Figure 3 shows the monthly average flow data of each hydrological station. According to Figure 3, the discharge data of the Duya Hydrological station is significantly higher than that of the Zhangfang hydrological station, which is caused by the interception of a large amount of water from Juma River by the Shengtian Canal. The Luobaotan hydrological station is located on the Juma River and downstream of the Zhangfang hydrological station. Except for flood season, the flow of Luobaotan hydrological station is almost 0, indicating that the Nanjuma River is in a state of perennial flow failure. The observation data of the hydrological stations show that the Zhangfang diversion project has caused a serious impact on water resources and aquatic ecology in the downstream.

Groundwater Flow Numerical Model
The western boundaries of the study area are the junction of a mountainous area and plain area, which mainly receives lateral recharge in front of the mountain and is generalized as prescribed flow boundaries. The northern boundaries of the study area are perpendicular to the multiyear average groundwater level contour, which can be generalized as no-flow boundaries. The eastern and southern boundaries are mainly the administrative boundaries of cities and counties. Observation wells are distributed near the boundaries, and the variation range of the groundwater level is small; therefore, the boundary can be generalized as prescribed head boundaries. The generalization of boundary conditions is shown in Figure 2a. Based on the study objectives and the distribution characteristics of the aquifer, we generalized the model as a single-layer model. The upper boundary is the phreatic surface, which accepts external replenishment. Due to the continuous increase in groundwater exploitation, the depth of the groundwater table is greater than

Groundwater Flow Numerical Model
The western boundaries of the study area are the junction of a mountainous area and plain area, which mainly receives lateral recharge in front of the mountain and is generalized as prescribed flow boundaries. The northern boundaries of the study area are perpendicular to the multiyear average groundwater level contour, which can be generalized as no-flow boundaries. The eastern and southern boundaries are mainly the administrative boundaries of cities and counties. Observation wells are distributed near the boundaries, and the variation range of the groundwater level is small; therefore, the boundary can be generalized as prescribed head boundaries. The generalization of boundary conditions is shown in Figure 2a. Based on the study objectives and the distribution characteristics of the aquifer, we generalized the model as a single-layer model. The upper boundary is the phreatic surface, which accepts external replenishment. Due to the continuous increase in groundwater exploitation, the depth of the groundwater table is greater than 10 m; therefore, the evaporative water loss of the diving surface can be ignored. The bottom boundary of the phreatic aquifer is selected as the model boundary and regarded as a nonflow boundary.
In this study, the flow field on 1 January 2017 was taken as the initial condition of the model. Due to topography, groundwater mainly flows from northwest to southeast. Figure 4a shows the initial flow field in the study area. The study area is located on the alluvial fan of the Juma River, and the groundwater type is quaternary loose rock pore water. According to the survey data, the northwest of the study area is mainly gravel, medium sand, and coarse sand, while the southeast is interbedded with medium sand, fine sand, and silty clay. According to formation lithology and the geological structure, the study area is divided into five regions of hydrogeological parameters. Figure 4b shows the five regions of the hydrogeological parameters. Table 1 shows the value range of hydrogeological parameters in different partitions.   Groundwater recharge in the study area includes precipitation infiltration, agricultural irrigation regression, lateral inflow, and river leakage. The precipitation infiltration coefficient in the study area is about 0.11-0.22 and the agricultural irrigation coefficient is about 0.1-0.15. According to the data of the Water Resources Bulletin, the average annual agricultural irrigation water consumption in the study area is about 3.74 million m 3 /a, and the agricultural irrigation regression water amount is about 56 million m 3 /a. The lateral inflow is calculated by Darcy′s law according to the contour map of the groundwater level. The lateral inflow of piedmont and other inflow boundaries is calculated to be 85 million  Groundwater recharge in the study area includes precipitation infiltration, agricultural irrigation regression, lateral inflow, and river leakage. The precipitation infiltration coefficient in the study area is about 0.11-0.22 and the agricultural irrigation coefficient is about 0.1-0.15. According to the data of the Water Resources Bulletin, the average annual agricultural irrigation water consumption in the study area is about 3.74 million m 3 /a, and the agricultural irrigation regression water amount is about 56 million m 3 /a. The lateral inflow is calculated by Darcy s law according to the contour map of the groundwater level. The lateral inflow of piedmont and other inflow boundaries is calculated to be 85 million m 3 /a. According to the investigation data, the river leakage coefficient in the study area is 0.45. Groundwater discharge in the study area includes exploitation and lateral outflow. According to the Water Resources Bulletin, the average groundwater exploitation in the study area is 461.56 million m 3 /a. The lateral outflow is calculated by Darcy's law. The distribution of precipitation infiltration, agricultural irrigation and groundwater exploitation is divided according to the administrative boundaries of each county. This is mainly determined by the data of the Water Resources Bulletin that we collected. Figure 1 shows the administrative divisions of Zhuozhou, Yuyuan, Gaobeidian and Dingxing.
According to the variation in the groundwater flow field, the model is extended to a three-dimensional transient groundwater flow model with heterogeneous anisotropy. To simulate the change in groundwater level in the study area, MODFLOW-2005 is used to solve the three-dimensional numerical groundwater flow model in the study area. MODFLOW-2005 uses the finite difference method to calculate the hydraulic conduction between cells. First, the asymmetric conduction matrix is generated, and then the conjugate gradient method is used to solve the problem. MODFLOW-2005 is realized through the GMS platform developed by the Environmental Modeling Research Laboratory of Brigham Young University in the United States. In this study, we set the basic grid size of the model to 500 m×500 m. In the vertical direction, the region is discretized into one layer by using a deformed mesh. A total of 2592 grids exists in the model.
The simulation period is from January 2017 to December 2023, which is divided into 29 stress periods. Each stress period from 2017 to 2018 is one month, with 24 stress periods in total and a time step of 10 days. This time period is used as the identification and verification period of the model. From 2019 to 2023, each stress period is one year, and the time step is 30 days.

Model Calibration and Validation
Model calibration is the process of adjusting and improving the model parameter structure and parameter values. The model in this study is corrected by the trial and error method to continuously adjust the hydraulic conductivity and recharge input parameters [26,27] to minimize the error between the simulated value and the observed value [28,29]. The measured groundwater flow field on 26 May 2017 is used as the identification and verification flow field. Figure 5a shows that the simulated flow field and the measured flow field have the same trend and flow patterns. Except on the southwestern side of the study area, which has a poor simulation effect, the other parts all reflect the actual groundwater flow trend.
A total of 27 monitoring wells are used in this calibration process, and their distribution locations are shown in Figure 2a. In this study, Yihezhuang, Xiyi an, Nangaoguanzhuang, Chenzhuang, Liyuzhuang and Chencunying are selected as six typical monitoring holes for display (the rest are not shown). Figure 5b shows the process fitting curve between the calculated water levels and the measured water levels of the selected monitoring wells. The simulated groundwater level process line of the typical monitoring wells is consistent with the measured groundwater level, which accurately reflects the process of groundwater before and after water replenishment.
The parameter partition after identification and correction is shown in Table 2. The hydraulic conductivity changes show a clear trend of gradually decreasing from the upper part of the alluvial fan to the downstream plain region, and the corrected parameter is within its value range. Therefore, the established numerical model can reflect the variation

Groundwater Budget Calculation
The groundwater budget calculation in the study area from January 2017 to December 2018 is obtained according to the operation results of the groundwater flow numerical model, as shown in Table 3. During the simulation period, the total recharge of groundwater in the study area is 466.8 million m 3 , the total discharge is 476.79 million m 3 , and the supplementary discharge difference is −9.99 million m 3 . The groundwater system is in a negative water balance state. Among the recharge items, the precipitation infiltration is 257.35 million m 3 , followed by river leakage of 68.77 million m 3 . They account for 69.86% of the total groundwater recharge in the study area. The discharge is mainly assembled by artificial mining, which is 461.56 million m 3 , which accounts for 96.81% of the total discharge.

Groundwater Budget Calculation
The groundwater budget calculation in the study area from January 2017 to December 2018 is obtained according to the operation results of the groundwater flow numerical model, as shown in Table 3. During the simulation period, the total recharge of groundwater in the study area is 466.8 million m 3 , the total discharge is 476.79 million m 3 , and the supplementary discharge difference is −9.99 million m 3 . The groundwater system is in a negative water balance state. Among the recharge items, the precipitation infiltration is 257.35 million m 3 , followed by river leakage of 68.77 million m 3 . They account for 69.86% of the total groundwater recharge in the study area. The discharge is mainly assembled by artificial mining, which is 461.56 million m 3 , which accounts for 96.81% of the total discharge. Through model correction and validation, the groundwater flow model can be used to predict the dynamic changes of groundwater level in the downstream plain of Zhangfang in the future. According to the monitoring data of the hydrology station, the Zhangfang water diversion project diverted a large amount of water, which caused a serious impact to the downstream water resources and river ecology. Water diversion projects directly lead to a significant reduction in river discharge, coupled with a large degree of groundwater exploitation; surface water and groundwater interaction gradually developed into a single form of river recharge groundwater. Therefore, the Zhangfang water diversion project should consider the shortage of water resources in the downstream as a whole and make a reasonable allocation of water resources.
In this study, we considered three schemes to predict groundwater level changes in the study area from 2019 to 2023, we set the stress period as one year and the time step as one month in the prediction period.
Scheme 1 maintains the current water diversion volume of the Zhangfang water diversion project at 100 million m 3 . This scheme can be used to evaluate the change of groundwater resources of the downstream plain of Zhangfang. Scheme 2 considers reducing the amount of water diverted from the Zhangfang water diversion project by 50%, that is, to discharge 50 million m 3 /a of water into the South Juma River and North Juma River. The discharge is divided between the South Juma River and North Juma River in a natural water ratio of 7:3. The direct utilization coefficient of surface water is 0.5, and the infiltration coefficient of the river is 0.45. Therefore, the river infiltration replenishment in the downstream plain area increases by 11.25 million m 3 /a. Scheme 3 assumes that the water diversion volume of the Zhangfang water diversion project is 0 and analyzes the dynamic influence of groundwater recharge on the groundwater level of the downstream plain after the water volume of the Beijuma River and Nanjuma River increases. In Scheme 3, groundwater infiltration recharge is increased to 22.5 million m 3 /a.

Prediction of Groundwater Level of the Downstream Plain of Zhangfang in 2023 in Scheme 1
A groundwater flow field in the study area at the end of 2018 is shown in Figure 6a. The prediction of the groundwater level in 2023 is shown in Figure 6b. Figure 6c shows the simulation of groundwater level variation from 2019 to 2023 in Scheme 1. According to Figure 6c, if the amount of diversion of the Zhangfang water diversion project remains unchanged, the groundwater level in the downstream plain of Zhangfang will continue to decline. The water diversion of the Zhangfang water diversion project is about 100 million m 3 /a, accounting for 1/5 of the total replenishment of the downstream plain area. We think that the amount of water diverted is catastrophic. The region with the largest drop in groundwater level is the northwest of the study area, where the aquifer is thin and the hydraulic conductivity is large. In the case of unfavorable upstream water inflow conditions, the groundwater in this region is very easy to lose. The maximum drop of water level in this area is more than 10 m. water diversion project reduces the recharge source of groundwater and intensifies the drop of groundwater level. At the same time, the Zhangfang water diversion project leads to the interruption of downstream river flow, which causes serious harm to water ecological security. Therefore, we should consider reducing the amount of water diverted by the Zhangfang water diversion project and rationally allocate the water resources.  The current development and utilization pattern of water resources in the downstream plain of Zhangfang is extremely unreasonable. The diversion of the Zhangfang water diversion project reduces the recharge source of groundwater and intensifies the drop of groundwater level. At the same time, the Zhangfang water diversion project leads to the interruption of downstream river flow, which causes serious harm to water ecological security. Therefore, we should consider reducing the amount of water diverted by the Zhangfang water diversion project and rationally allocate the water resources. to 2023., Figure 8a,b shows the groundwater level variation from 2018 to 2023 of Schemes 2 and 3.
According to Figure 8, compared with Scheme 1, reducing the water diversion amount of the Zhangfang water diversion project to recharge groundwater can alleviate the decline rate of the groundwater level in the whole study area. Although the water level in the eastern part of the research area still continues to decline, the maximum decline range is reduced to 2-4 m. The water level in the northwest of the study area increased greatly, and the groundwater level in the water receiving areas of the South Juma River and North Juma River also increased to a certain extent. Figure 8 shows that the unreasonable allocation of water resources in the Zhangfang water diversion project has caused serious impacts on the downstream water resources. The pressures resulting from the water shortage in the downstream can be alleviated by reasonably reducing the water diversion.  According to the operation results of Schemes 2 and 3, the calculation results of the groundwater budget in the study area are obtained. According to Table 4, the river leakage in Scheme 2 increases by 11.25 m 3 /a, and the supplementary discharge difference is 0.13 million m 3 /a. The groundwater system is in a positive water balance state. In Scheme 3, the river leakage increases by 22.5 million m 3 /a, and the supplementary discharge difference is 8.31 million m 3 /a. The groundwater system is in a positive water balance state. According to Figure 8, compared with Scheme 1, reducing the water diversion amount of the Zhangfang water diversion project to recharge groundwater can alleviate the decline rate of the groundwater level in the whole study area. Although the water level in the eastern part of the research area still continues to decline, the maximum decline range is reduced to 2-4 m. The water level in the northwest of the study area increased greatly, and the groundwater level in the water receiving areas of the South Juma River and North Juma River also increased to a certain extent. Figure 8 shows that the unreasonable allocation of water resources in the Zhangfang water diversion project has caused serious impacts on the downstream water resources. The pressures resulting from the water shortage in the downstream can be alleviated by reasonably reducing the water diversion.
According to the operation results of Schemes 2 and 3, the calculation results of the groundwater budget in the study area are obtained. According to Table 4, the river leakage in Scheme 2 increases by 11.25 m 3 /a, and the supplementary discharge difference is 0.13 million m 3 /a. The groundwater system is in a positive water balance state. In Scheme 3, the river leakage increases by 22.5 million m 3 /a, and the supplementary discharge difference is 8.31 million m 3 /a. The groundwater system is in a positive water balance state. Therefore, in order to ensure the sustainable utilization of groundwater in the study area, we propose to reduce at the diversion water amount of the Zhangfang water diversion project by at least 50%to recharge the groundwater. According to Figure 8b, when the water diversion amount of Zhangfang water diversion project is 0, the groundwater level of the South Juma River and North Juma River receiving area rises significantly, but the water level near Baigou River still drops. The reason for this is that we only recharge groundwater though natural rivers. Therefore, when a large number of surface water sources can be used to recharge the groundwater, the replenishment amount of the surface water increased by connecting channels.

Conclusions
By means of groundwater monitoring and numerical simulation, this paper studied the impact of the Zhangfang water diversion project on the downstream groundwater level. Through research, the following conclusions are obtained.
The groundwater flow model in the downstream plain area of Zhangfang is established. With the flow model, it is concluded that the Zhangfang water diversion project reduces the river discharge and the river infiltration recharge in the study area, which places the supplement and discharge difference in a negative water balance state.
If the Zhangfang water diversion project maintains the current water diversion amount of 100 million m 3 /a, the water level of the downstream plain of Zhangfang will continue to decline. According to the prediction results, the groundwater level of the South Juma River and the North Juma River will rise to a certain extent by reducing the amount of the Zhangfang water diversion project to recharge the groundwater. The water level of the Baigou River will continue to decline, but the decreasing rate will be effectively alleviated.
For the sustainable development of the region, it is recommended to reduce the amount of water diverted by the Zhangfang water diversion project by at least 50%.
Through numerical simulation, the dynamic influence of the Zhangfang water diversion project on the change of groundwater levels in the downstream is discussed from the perspective of water resources allocation, which has practical significance. The prediction scheme of the model can provide suggestions for decision-making.
Author Contributions: M.Z., conceptualization, formal analysis, investigation, methodology, software, and writing-original draft; X.M., conceptualization, formal analysis, and supervision; B.W., data curation and methodology; D.Z., data curation and software; Y.Z., data curation and Investigation; R.L., software. All authors have read and agreed to the published version of the manuscript.