Simulation-optimization model for groundwater replenishment from the river: a case study in the Hutuo River alluvial fan, China

Groundwater is one of the most valuable natural resources for industry, agriculture, and human livelihoods (Khosravi et al. 2018). However, as the principal source of water supply in the North China Plain (NCP), groundwater is inevitably overused, which has led to a series of geological and environmental problems (e.g., groundwater levels drawdown and land subsidence) (Kinzelbach et al. 2022). On the other hand, long-term overexploitation of groundwater caused a constantly declining in groundwater levels and the formation of extremely thick vadose zone, which provides storage space for groundwater replenishment from the river (Fei et al. 2011). Moreover, the operation of the middle route of the South-to-North Water Diversion Project (m-SNWDP) has provided opportunities to restore groundwater storage in the NCP by increasing water supply and reducing groundwater exploitation since December 2014 (Zhang & Li 2014; Long et al. 2020).

Under the circumstances, ecological water replenishment (EWR) has been carried out in the NCP since September 2018. Over 22 rivers along the water transfer channel receive water from the m-SNWDP to replenish groundwater and improve water quality in the receiving areas. (Li et al. 2018; Ding et al. 2019; Sun et al. 2021). Hutuo River, as one of the first pilot rivers for EWR, is referred to as ‘the Mother River’ by the natives of Shijiazhuang, a provincial city of Hebei. However, due to anthropogenic activity and climate change, the river channel flow has been continuously interrupted since the 1970s. In addition, the groundwater quality of the Hutuo River alluvial fan continues to deteriorate as a result of the leaching and nitrification reactions promoted by the overexploitation of groundwater (Omer et al. 2016; Zhang et al. 2019). Previous research indicated that the Hutuo River has the top priority for groundwater replenishment (Zhu et al. 2020). However, there is little research on the optimal groundwater recharge from different reaches of the Hutuo River, especially in the case of large-scale EWR.

Although river infiltration replenishing aquifer is effective for restoring groundwater resources, excessively high groundwater levels will cause corresponding environmental and engineering safety issues, such as soil salinization (Shokri-Kuehni et al. 2020) and interference with underground structures (Gattinoni & Scesi 2017). Therefore, the selection of a reasonable recharge strategy is crucial in the process of groundwater replenishment. However, designing the feasible and efficient operation of groundwater replenishment is far from simple, because the aquifers’ characteristics are usually complex, and the response of water levels to recharge is nonlinear (Hao et al. 2018). Therefore, to find a satisfactory groundwater recharge scheme, it is significant to figure out the states and characteristics of the groundwater flow and various land use constraints.

Coupling analytical or numerical groundwater simulation models with optimization algorithms is a common approach to solving such problems. The simulation-optimization (S/O) method is a powerful tool employed to solve groundwater resources management problems in many fields (Singh 2014). Previous studies suggest that the S/O model was used in groundwater remediation strategy (Mategaonkar & Eldho 2012), seawater intrusion management for coastal aquifers (Dentoni et al. 2014), and optimal management of the artificial recharge-pumping system (Kawo et al. 2018). However, coupled S/O models may present particular challenges due to the long-time scales and model runtimes that are typically involved. Furthermore, the technical complexity introduced by using linked software packages and data-exchange libraries is a common issue in conventional approaches for the coupled model of such groundwater systems (Castilla-Rho et al. 2015). To address these limitations, integration of pillars in this study is performed via Python programming language, with extensive use of the Python FloPy library for writing inputs, and post-processing the majority of outputs. Compared with a low-level language (Fortran or C, for example), complex tasks can be achieved with a few lines of readable code in Python (Bakker et al. 2016).

Traditional optimization techniques tend to prematurely converge to a local optimum in solving complex non-linear non-convex problems, such as linear programming (Veintimilla-Reyes et al. 2018), non-linear programming (Bajany et al. 2021), and dynamic programming (Yu et al. 2017). As a stochastic global search optimization method, a genetic algorithm (GA) shows advantages in handling discontinuous and non-differentiable problems (Goldberg 1989; Yeh 2015). GA is widely used in the research of optimization problems because of its advantages, such as easy convergence and high accuracy (Ayaz 2017). Rahimi et al. (2014) combined the analytic hierarchy process and GA to identify the most suitable areas for the artificial groundwater recharge by ranking all the alternatives and weighing the criteria involved. Pramada et al. (2018) coupled SEAWAT with GA to determine the optimum pumping strategy, reducing the risk of pumping contaminated water. Hassan & Khalaf (2020) integrated the numerical groundwater model and GA to obtain the optimal use of groundwater resources in the Karbala Desert, Iraq. All the works mentioned above show the advantages of GA in the procedure of groundwater management, especially for the optimization of the facilities’ locations and pumping strategy, but the optimization of groundwater replenishment from the river has rarely been attempted in previous studies.

In this study, we present a S/O-GA model for groundwater replenishment from the Hutuo River. The calibrated simulation model was performed as a script-based model using the FloPy library (Bakker et al. 2022), which allows the MODFLOW-USG (Panday et al. 2013) model to be executed in a Python programming environment. Specifically, we optimized the volume of groundwater recharge from each river reach during the management period to maximize the recovery of the groundwater storage capacity under the upper limits for groundwater levels.

The principal aquifers in the study area are the Quaternary porous aquifers (Figure 1(b)). Based on the geological stratification of the Quaternary, the four strata Q_h-Q_p1 from new to old correspond to three aquifers. The uppermost two layers, Holocene and upper Pleistocene have coalesced into aquifer I, which is a phreatic-micro confined aquifer. The middle and lower Pleistocene correspond to aquifer II and aquifer III, respectively, and both are deep confined aquifers. From west to east, the thickness of the Quaternary alluvial deposits increases, and the grain size decreases. These Quaternary aquifers consist of gravels, coarse sand, locally fine sand, and discontinuous clay characterized by a stratified distribution. Groundwater flows from the northwest to the southeast in this area. The main recharge of groundwater is precipitation, river infiltration and lateral inflow, while the groundwater discharge is mainly artificial exploitation. The western aquifers boundaries are composed of the division lines between mountain and plain, processed as flux boundaries. Parallel to the streamline, the north-western boundary is a zero-flux boundary. The eastern and southern boundaries are administrative boundaries, which are processed as specified head boundaries, controlled by observed data (Figure 1(a)).

It is an essential objective of EWR to allocate the abandoned water before the flood period in the south, and to restore groundwater resources of the receiving areas in the north by using the surplus water supply capacity of the m-SNWDP. From September 13 to December 27, 2018, the cumulative water of 420 × 10⁶ m³ flowed into the Hutuo River, including 365 × 10⁶ m³ from m-SNWDP, 34 × 10⁶ m³ from the reservoir, and 21 × 10⁶ m³ from reclaimed water.

The numerical simulation model was based on MODFLOW-USG (Panday et al. 2013), which supports unstructured grid (U-grid) types, including nested grids. U-grid is more flexible than traditional structural grid generation method and is used to subdivide the grids along the main channel of the Hutuo River. With initial cell size of 2,000 × 2,000 m, the cells near the river are refined to 500 × 500 m, and the total cell number is 15,312. Based on the hydrology and geology conditions of the study area, the model consists of three layers, and the boundaries are shown in Figure 1(a). Particularly, specified flux conditions are used to simulate natural and artificial recharge from each river reach to groundwater. The distribution and values of hydraulic conductivity, regional topographic and geological conditions are supported by previous research (Zhang et al. 2009). Model simulations were performed for the period from January 2000 to December 2019. The first 15 years were used for model calibration, and the last 5 years were used for validation. The calibration of the groundwater flow model was manually adjusted by the trial-and-error method.

After the calibration and validation of the groundwater model, the simulation model for EWR is set up for the period from January 2020 to December 2029. The management period is 10 years divided into 40 stress periods, further subdivided into time steps of 90 days. The end of 2019 hydraulic heads are used to set the initial conditions for the S/O model. The same recharge and discharge over the previous decade (2010–2019) are considered.

The Hutuo River between Huangbizhuang reservoir and Xian county is divided into seven reaches (R1-R7) as shown in Figure 1(a). Table 1 lists the basic characteristics and infiltration capacity of each reach. The maximum recharge rates are given according to previous research (Du et al. 2013) and the current condition of the river channel. Furthermore, the minimum recharge rate of each reach is set to one-tenth of the maximum recharge rate.

Considering the safety of underground infrastructures, the controlled limit value of groundwater table depth in urban areas is restricted to 15–20 m. In addition, to prevent the accumulation of salts in the soil, the limit value in the agricultural area is set to 2–5 m (Mu et al. 2020).

A GA realizes the heuristic search of complex space that finds the global optimal solution with a high probability by imitating the recombination and mutation of genes, which make them constantly have new traits to be retained in the process of biological evolution in nature. In this study, an improved GA with elite reservation is adopted that replaces the moderate of non-optimal individuals of the parent generation with the same number of individuals of the offspring population to get a new generation population (Deb et al. 2002). Hence, the convergence speed of the algorithm can be effectively improved.

GA begins by initializing a population of N individuals, each representing a possible recharge scheme in the feasible solution space. Then the phenotypes of individuals are written into the MODFLOW-USG input file. Run the groundwater flow model and obtain the groundwater levels for each stress period. The groundwater level constraints were handled using a penalty function. Individuals that violate the constraints are assigned a penalty factor that reduces their fitness. The fitness of individuals is evaluated by a fitness function that synthetically considers the objective function and the penalty function. Elite with the highest fitness of each generation will be reserved. Moreover, N-1 individuals are selected by the tournament method (Chu et al. 2018). The optimal individual of the parent generation and the offspring constitute a new population generation for crossover and mutation. Likewise, the new population is subjected to the same evolutionary process until a convergence criterion is satisfied or the maximum evolution algebra is reached. All the methods and calculations were carried out in a Python programming environment.

The GA is applied to determine the optimal solution for the groundwater recharged by the Hutuo river. The population size is set to 500, and the maximum number of iterations is set to 150. Randomness plays a substantial role in keeping the algorithm searching in the solution space. In this study, the probabilities of crossover and mutation are set to 0.70 and 0.03, respectively. The summary of the GA results and the optimized groundwater recharge scheme is as follows.

The total amount of water recharged into the aquifers is 5.36 × 10⁹ m³ during the management period, and the percentages of recharge in each reach are listed in Table 1. The upper river reaches are the main recharge area, accounting for 67.58% of the total amount. This is because that the infiltration capacity of the upper reaches is better than that of the downstream channel. Moreover, due to its closer proximity to the cone of groundwater depression in Shijiazhuang, R1 has greater potential for recharge than R2 and R3. The annual average recharge from R1 is 1.91 × 10⁸ m³/a. According to the previous study, if the volume of artificial recharge water in R1 is less than 7.86 × 10⁸ m³/a, groundwater reservoir is a better storage method than the local reservoir (Du et al. 2013).

The inherent stochastic nature of genetic algorithms is an unneglectable factor that influences the fluctuation of the curves in the plateau period. For instance, the maximum recharge rate can be maintained without exceeding the upper limit water level due to the relatively low infiltration capacity at R4, R5, and R6. In this case, the fluctuation of the recharge rate comes from the randomness in the process of GA's selection, crossover, and mutation.

As presented in GW9 and GW15, the rises of groundwater levels between regions mutually affect each other. When the groundwater level in GW9 reaches the upper limit, the optimal recharge rate of R1 declines rapidly, and the rate of groundwater level rise in GW15 decreases from 3.4 to 2.2 cm/d. Compared with the constant-rate recharge, the optimal recharge scheme dynamically adjusts the recharge rate according to the storage space and groundwater regime, so as to ensure that the groundwater level is maintained below the upper limit. However, in the area close to the R4, R5, and R6, the upper limit of groundwater level does not affect the optimal recharge rate of reach. Until the end of the management period, the groundwater level in GW27 is 26.1 m below the upper limit. Therefore, the optimization result almost coincides with the groundwater level under the constant-rate recharge scenario, and the infiltration capacity of the river reach is the main constraint for recharging. The recharge rate of R7 drops significantly before the level in GW35 approaches the upper limit. So, the reason why the groundwater level in GW35 is not maintained close to the upper limit is that the recharge rate holding the water level in GW35 constant will cause the calculated water levels in other cells along the river to exceed the upper limits.

Based on the EWR project of the Hutuo River, which has been implemented since 2019, this study proposed an optimization framework that integrates a GA with the MODFLOW-USG model to facilitate decision-making for groundwater management. The groundwater level data collected by fieldwork were used to validate the groundwater model. With NSE of 0.90 and R of 0.97, the results suggest that the model can accurately reflect the variation of groundwater levels caused by changes in hydrological regimes and anthropogenic interferences. From the performance of the GA, the population's fitness values rise with iterations, until the convergence is reached after the 75th generation, where the global optimal solution is obtained. Finally, the spatio-temporal optimal groundwater recharge strategy is presented, taking into account the upper limits for groundwater levels in different areas. The main findings are as follows:

• (1)

  The optimization results show that the total amount of water recharged into the aquifers from the Hutuo river is 5.36 × 10⁹ m³ from 2020 to 2029, and the main replenishment is concentrated in the upper reaches R1, R2, and R3, totally accounting for 67.58%.

• (2)

  EWR significantly influences the groundwater flow field in the study area, and the groundwater recharge from the Hutuo River accounts for 26.20% of the total recharge. Nevertheless, the optimal recharge reduces as groundwater exploitation declines, which are 17.23% and 18.44%, respectively.

• (3)

  Compared with recharging at invariable rates, the optimal recharge scheme can efficiently augment the groundwater storage without exceeding the upper limits of the groundwater levels.

The optimization framework based on the 3D transient groundwater flow model and GA has been found to be very effective for groundwater optimization problems and can provide a meaningful reference for the EWR. However, it takes a considerable computational effort to solve problems with the large-scale model and massive decision variables. Therefore, parallel computing methods are required to improve computational performance in the future.

This work was financially supported by the National Natural Science Foundation of China (41702282) and China Geological Survey Project (DD20190303).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.