In-situ pumping tests and numerical simulations of seepage wells in the Yellow River valley, China

Seepage wells that can convert surface water into groundwater are often constructed near river valleys to obtain more water and lead to smaller drawdown compared with traditional wells. Seepage wells have been widely used, whilst the groundwater and river-level variations caused by seepage wells are still unclear, and numerical models are rarely verified due to the lack of in-situ observational data, which may lead to results that are quite different from the actual conditions. To address those limitations, a large-scale pumping test was carried out near the Yellow River valley in China and a coupled seepage–pipe flow model was established using the exchange yield between the aquifer and pipe as the coupling key in this research. The coupled model was evaluated with in-situ measurement. The field observation showed that both the river and groundwater had a positive response to the pumping of the seepage wells. The simulation results indicated that our model can well estimate the pumping rate and drawdown with root-mean-square deviations of 158.235 m/d and 0.766 m, respectively. Further, it is also found that the groundwater showed the obvious characteristics of three-dimensional flow under the influence of seepage wells and the maximum drawdown should be less than 15 m to ensure exploitation efficiency. These findings provide important information that can guide the design and construction of seepage wells to improve the rational exploitation of groundwater.


GRAPHICAL ABSTRACT INTRODUCTION
Understanding the relationship between groundwater and surface water is crucial for the management and utilization of water resources in arid regions . Previous studies have observed and recorded the variation of water level by using the traditional vertical wells pumping test and analyzed the relationship between the river and the aquifer Huang et al. 2012). There are few studies on water-level observation of rivers and aquifers during the pumping process of seepage wells. Some physical simulation experiments have explored the mechanisms responsible for changes in the groundwater flow (Zheng 2013;Liu 2018;Gao et al. 2020). However, the simulation tests failed to accurately describe the actual situation due to the limitations of both test equipment and conditions. Further in-situ observations are therefore necessary to reveal the mechanisms responsible for the water intake of seepage wells.
Previous studies have focused on pumping engineering, including horizontal wells, radiating, and seepage wells, for which models have been developed. For example, Hantush & Papadopulos (1962) characterized horizontal wells using the line sink method with equal strength distribution. Zhan et al. (2001) and Park & Zhan (2003) then used this method to generalize the non-tube-well water-collecting structures.  noted that the water head (H) in the well was not distributed in equal strength due to head loss in the well.  proposed equivalent hydraulic conductivity and coupled seepage-pipe flow models that considered a head loss in the well, which were widely accepted. Wang & Zhang (2007) established a calculation model of seepage wells based on the theory of equivalent hydraulic conductivity and coupled seepage-pipe flow model. Wang et al. (2013) improved the model by taking the exchange yield between the radiant pipe and the aquifer as the coupling key, which realized the calculation of multiple seepage wells working simultaneously. Zhang et al. (2019) proposed that the flow rate of injection wells could be calculated by the pumping theory of radiant wells and established a mathematical model of the rising curve of water level around a radial well. Maroney & Rehmann (2017) developed a mathematical method to calculate the stream depletion rate for a radial collector well in an unconfined aquifer near a fully penetrating river. Yang (2020) used GMS software to simulate the water yield of seepage wells, and proposed engineering optimization suggestions for seepage wells. Mitrinovićet al. (2021) calculated the drawdown in radial wells including the effects and number of regenerations and rate of increase of the local hydraulic losses for constant flow or constant level using mathematical and software models. Due to the lack of in-situ observation data, most of the models lack verification, which results in a gap between the predicted results and the actual conditions.
The transition relationship between river and groundwater from the pumping of traditional vertical wells has been studied extensively. However, the research on seepage wells mainly focuses on simulation tests and model development. Most of these studies are only ideal cases, and the models are hard to verify. This research was undertaken at YR valley as the considerable water resources, which are of great importance to arid and semi-arid regions, and the hydrogeological conditions are suitable for large-scale in-situ observation. The objectives of this study include: (1) observe and record water levels of the aquifer and the YR during pumping; (2) establish a three-dimensional numerical model for seepage wells and verify it; (3) simulate the distribution and variation of the groundwater flow field; (4) discuss the applicability of the model under different drawdowns and obtain the maximum pumping drawdown. Our study can guide the future design and construction of seepage wells that can improve the rational exploitation of groundwater, particularly in arid and semi-arid regions.

Study area
The study was performed in the YR valley of MaZhen, Shaanxi province (latitude 38°37 0 22″-38°39 0 12″N, longitude 110°-52 0 04″-110°53 0 39″E). The YR valley has a temperate continental semi-arid monsoon climate, with mean annual precipitation, evaporation, and temperature of 422.63 mm, 1,788.4 mm, and 8.6°C, respectively. The main water system is the YR, which flows through the eastern part of the study area ( Figure 1), with an average annual runoff of 170.53 Â 108 m 3 /a. The Quaternary (Q) and Triassic (T) strata contribute to the main strata of this region according to geological borehole data. The thickness of the Quaternary strata is 15.0-20.0 m and the upper region is interbedded with mudstone, sandy mudstone, and arkose sandstone of unequal thickness. The lower region is dominated by sandstone with undeveloped fissures. The thickness of the Triassic strata is 160-200 m, including pebbles, fine sand, and medium/coarse sand. According to the pumping test results of exploration boreholes, the main hydrogeological parameters of the aquifer are shown in Table 1. The groundwater is recharged by agricultural irrigation, precipitation infiltration, and the YR during the flood season, and is discharged by surface runoff, evaporation, and artificial mining.  Pumping test design Pumping tests can effectively estimate the hydraulic parameters of aquifers, identify aquifer boundaries, and evaluate well performance (Sun 2018). Pumping tests of seepage wells have been performed to determine the hydraulic parameters of the aquifer and to reveal the variation rules of the groundwater field through observational data. The seepage well is composed of a shaft, several galleries and chambers, and radiant pipe groups. In this study, seepage wells were constructed with a single shaft, five galleries, and five chambers, with each chamber constructed from 20 radiant pipes ( Figure 2). The length of the radiant pipes that passed through the Quaternary or weathered zone was ∼3 m. The dimensions of each chamber were 3 Â 3 Â 4 m. The length of each gallery (height and width of the gallery: 1.8 and 1.6 m) between the two adjacent chambers was 60 m. The diameter of the radiant pipe and shaft were 89 mm and 3 m, respectively. Five drawdown pumping tests of the shaft, including 3.32, 5.97, 10.58, 23.94, and 31.50 m, were designed to analyze the relationship between the flow rate and drawdown. As the thickness of the Quaternary formation where the shaft is located was 17 m, along with the galleries located 44 m below the ground, the Levellogger automatic water-level gauge was set to a depth of ∼7, 11, and 15 m to observe the water levels of the Quaternary, and to 22 and 47 m to observe the water levels of the Triassic. The water levels of the three layers of the Quaternary were simultaneously observed in S1, S2, S3, S4, S5, S6, S7, S8, and S9, respectively. The bedrock water levels were observed in the observational wells S1, S2, S3, S4, S5, S6, and S8, respectively. S2 was also monitored to observe the bedrock water level at 47 m. S9 was set to observe the groundwater on the other side of the YR ( Figure 3). The observation time for all wells was 25 days with an interval of 1 min.

Numerical methods
To facilitate the generalization of seepage wells, studies were performed in the YR valley (with the converted geodetic coordinates X ¼ 19487488-19489294 m, Y ¼ 4276624-4280386 m). The YR has a total area of 2.48 km 2 . The simulation area was determined according to the actual hydrogeological conditions of the study area. The hydrogeological conceptual model was established by reasonably generalizing the boundary conditions, aquifer structure, and groundwater flow characteristics of the study area. As shown in Figure 4, YR as the eastern boundary of the simulation area was set as the river boundary. The western boundary formed the demarcation between the low hills and the YR valley, due to the poor permeability of the bedrock (set as a zero flow boundary). Because of the broad water surface on the YR and its longitudinal extension (set as a river boundary), the water intake project was distanced from the upper and lower reaches of the YR. As such, the upper and lower reaches of the YR were generalized as a constant head boundary. The top surface of the simulation area was set as the water table, where vertical water exchange occurred, inclusive of precipitation infiltration and evaporation. The bottom was generalized as a zero flow boundary due to the complete Triassic bedrock.  Under natural conditions, the main recharge sources of the groundwater included precipitation, agricultural irrigation, and the YR. The discharge included surface water discharge, artificial mining, evaporation, and transpiration. During the pumping period of seepage wells, the YR recharge formed the main recharge source of the groundwater.
In recent years, several models have been developed to evaluate water intake engineering in seepage wells and other similar groundwater collecting structures (Mohamed & Rushton 2006;Wang & Zhang 2007;Lee et al. 2012;Wang et al. 2013). Our study presented herein employed the method of Wang et al. (2013) to calculate water intake from seepage wells. Radiant pipes and galleries were conceptualized as a series of interconnected nodes and pipes, which were superposed on the finite difference grid of MODFLOW, a modeling environment for practical applications in three-dimensional groundwater flow, which has been used widely to overcome issues related to groundwater (Al-Hashmi et al. 2020;Vellando et al. 2020). According to the actual position of each node, its plane and vertical positions were determined to simulate the actual radiant pipes and galleries. Using the exchange flow rate between seepage wells and the aquifer as the coupling key, seepage wells could be described as the head boundary. Using this generalization, seepage flow in the aquifer and pipe flow in seepage wells were well coupled. This seepage-pipe coupling model could be calculated via the numerical method for a discrete pipe and discrete grid system. The mathematical model describing the groundwater flow in the aquifer can be written as: where K xx , K yy , and K zz are the hydraulic conductivity (m/d); H is the water head (m); x, y, and z are coordinate variables (m). Equation (2) (3)) was used to calculate the conductance of the seepage wells: where Q e is the exchange flow rate between the aquifer and pipe; C j,i,k is the conductance of seepage wells (m 2 /d); H p is the head of pipe flow (m); K j,i,k is the conduit wall permeability term (m/d) in MODFLOW cell j, i, k; π is the mathematical constant pi (unitless); d ip is the diameter (m) of the pipe ip connected to the node in; Δl ip is the straight-line length between nodes (m) of the pipe ip connected to the node in; τ ip is the tortuosity of the pipe ip connected to the node in (unitless); r ip is the radius of the pipe (m). The flow rate in the pipes can be calculated by using Equation (4): where Q p is the flow rate in the seepage wells (m 3 /d); ΔH is the head loss (m); l is the length of the flow path (m); g is the gravity acceleration (m/s 2 ); d is the diameter of the well pipe (m); v is the average flow velocity in the well pipe (m/s); e is the coarseness degree of the inner wall of the well pipe; υ is the kinematic viscosity of water (m 2 /d).
Based on the above methods, the flow rates both in the aquifer and pipes are obtained. When the pipe is filled with water, the sum of the inflow and outflow should be 0 in any node, which can be written as: where i is the number of the node which is linked to the node j; n i is the total number of nodes which are linked to the node j; and Q p,i.j is the total flow rate from the node i to the node j. The above method is used to couple the seepage flow and pipe flow well. Three interdependent variables (the average flow velocity, Reynolds number, and the head difference in the pipe) exist in the equation of Q p . Therefore, an iterative method is needed to solve this equation. Firstly, the elevation was set as the initial head of the well nodes and the initial flow regime was assumed to be laminar. Then, the exchange flow of the well nodes was calculated by the initial water head, and the distribution of the flow rate of the seepage wells was obtained. The Reynolds numbers were calculated according to the flow rate and diameter of the seepage wells, and then the corresponding equation was re-selected to calculate the head of the seepage wells. The above calculation process was repeated by using the new heads of the seepage wells. The iteration stopped when the maximum absolute value of head differences between two contiguous iterations met the criteria of precision. Finally, the exchange flow rate was added to the MODFLOW equations to solve the heads of the aquifer. Through the above calculation process, the water intake problem of seepage wells can be simulated. A detailed description of the calculation method is described by Wang et al. (2013).
Based on the hydrogeology conceptual model and mathematical model, the numerical model was established to compute the exchange flow rate amongst the aquifer, radiant pipe, and pumping yield of the seepage wells. During the process of these simulations, a large pumping flow of the seepage wells could cause the water in the pipes to be immediately drained and the cells where the pipes were located to be dry during the calculation process. The coupled seepage-pipe flow model could not be used for iterative calculations during the transient flow simulation process. Furthermore, the fourth and fifth periods showed incomplete stability during pumping tests of the seepage wells. A vast number of nodes, cells, and pipes were also depicted in the model, which involved a series of complex calculations and longer computing times. Therefore, steadystate flow was simulated during the first three test periods in this study.
To validate the simulated drawdown of each observation well, statistics including correlation coefficient (R 2 ), root mean square deviation (RMSD), and relative errors (RE) were used . Higher values of R 2 and smaller values of RE and RMSE indicate a better performance of model and method. The equations are as follows: where obs is the observation data; sim is the simulation data; i is the ith simulation or observation data; N is the total number of the simulation or observation data; obs is the average of the observation; sim is the average of the simulation.

RESULTS AND DISCUSSION
Pumping test analysis From Figure 5(a), the relationship between the drawdown of the shaft and the pumping flow can be defined. The flow rate of seepage wells increased concomitantly with the shaft drawdown. When the drawdown was excessive, partial disconnection of the groundwater recharge could occur, which affected the flow rate of the seepage wells (Cognac & Ronayne 2020). Figure 5(b) shows the high responsiveness between the pumping flow rate and the drawdown of the shaft. The groundwater levels decreased synchronously with the increase in pumping flow rate, with stages 1-3 reaching a stable state. However, the fourth stage did not achieve stability. The fifth test stage can be regarded as an extension of the fourth stage through the analysis of the curve of the pumping tests in the fourth and fifth stages. The final pumping flow of the two stages was similar, so the two stages could be regarded as a single stage that finally reached stability. The water levels of each observation well in the Quaternary failed to reach stability across any period ( Figure 6). This may be due to the large specific yield of the Quaternary strata that made it difficult to reach a stable state in a short time and the  recharge of the YR had a great influence on it. The groundwater level was found to vary with the pumping rate and water levels of the YR. When the river level rose (or fell) and the pumping rate decreased (or increased), the groundwater level in the observation well of the Quaternary immediately rose (or fell). Taking S9 as an example, the groundwater level was consistent with the variation trend of the water level of the YR. During the flood peak of the YR, the groundwater level quickly responded, but a time lag was observed. These results reveal that the YR exerts a strong recharge effect on the groundwater. Figure 6 shows the head differences between the water levels of the layers that increased with the continuous growth of the pumping rate. Taking S1 as an example, the head differences between the observation layers (7 and 10 m) were 0.07, 0.1, 0.12, and 0.19 m at stages 1 to 4, respectively. The head differences between the observation layers (10 and 13 m) were 0.12, 0.21, 0.26, and 0.27 m, respectively. The curves of other observation wells also reflect this phenomenon. These data provide evidence that the YR recharges groundwater through a three-dimensional flow, regardless of the exploitation or natural conditions.
Excluding the S2 observation well, the variation curve of the groundwater drawdown in the bedrock was consistent with that in the shaft (Figure 6), indicating that the groundwater responded to the pumping rate in time. The drawdown of the bedrock was larger than that in the Quaternary, as the chambers are located at the bedrock formation, and radiant pipes in the bedrock form open holes. From stages 1-3, due to the small storage capacity of the bedrock, the water level of the observation wells was stable. However, the drawdown curves of the fourth and fifth periods fluctuated to different degrees, following changes in the pumping rate.
The drawdown curve (Figure 7) of S2 differed from that of the other observation wells and shaft. Groundwater levels in S2 present the overall trend of fluctuation (initially increasing then decreasing). These differences can be explained by the permeability coefficient of the bedrock, which is small and leads to a lagged response of the bedrock. It was necessary to drain the aquifer during the early stages of seepage-well construction and the groundwater levels would recover at the end of construction. The groundwater level continued to recover during the initial two periods due to the lag response of the deep bedrock. In addition, the pumping rate was lower than the recharge rate of the aquifer during the first two periods, causing the groundwater level to continue to rise. The water level during the last three periods decreased due to the large pumping rate, but a significant lag-time was observed.
To analyze the influence of pumping and YR recharge on drawdown, the drawdown curves (the observation layer below the ground surface of ∼11 m, Figure 8) of S4, S5, S7, and S9 (the observation wells were approximately located on the same northwest-to-southeast profile, Figure 3) were compared. The distances from S4, S5, S7, and S9 (S4 is located beneath YR and S9 is located on the other side of the YR) to the D2 chamber are 52, 32, 22.2, and 121 m, respectively. Figure 8 shows that the drawdown of S4, S5, S7, and S9 during each test period successively decreased.
Although the distance of S7 to D2 was longer than the distance of S5 to D2, the drawdown of S7 exceeded that of S5. This can be rationalized as due to the distance of S7 to the top of the radiant pipe being closer than the distance of S5 to the top of the radiant pipe. The main water yield is from the top of the radiant pipe, which is greatly influenced by the pumping. The drawdown of S4 was less than S5 as S4 was located below the YR and could be quickly recharged. The drawdown of S9 was the smallest, as it was far away from D2, and was affected by changes in the water level of the YR. In general, groundwater in the main water intake layer of the radiant pipe was influenced by pumping, and the YR played a pronounced role in weakening the effects of seepage-well pumping.

Model and simulation evaluations
To evaluate the performance of the model, the simulated drawdown values (S sim ) and the observational data (S obs ) in each observation well and simulated pumping flow during each period were compared (Figure 9). For pumping flow, the correlation coefficient (R 2 ), root mean square deviation (RMSD), and relative errors (RE) are 0.9995, 158.235 m 3 /d, À0.0007, respectively. For drawdown, the R 2 of stages 1-3 were 0.938, 0.961, and 0.966, respectively. The RMSD for these simulations were 0.345, 0.689, and 1.044 m, respectively. The RE for these simulations were À0.105, À0.176, and À0.115 m, respectively.
These results suggest that the simulated results are close to the observational data. However, some errors were observed between the simulated and observational data, with the maximum error of the drawdown and pumping flow (S2) during the third period showing errors of 1.45 m and 201.97 m 3 /d, respectively. These may be due to the presence of dry cells in the model and pipes, which do not participate in the calculation due to the large drawdown of the shaft. Meanwhile, the observational data may also display errors during each period. In addition, the model inevitably harbored errors, as it was difficult to describe the actual situation in full. The model was however able to simulate the drawdown and pumping flow of the wells during each period. The final hydraulic parameters of the model for the Q are K h ¼ 29.5 m/d, K v ¼ 5.5 m/d, S y ¼ 0.2, and for the T aquifer included The model verified by actual data can be used to reveal the groundwater drawdown field under the pumping conditions. Figures 10(a), 10(b), and 11 indicate the drawdown of layers 15 (Quaternary), and 20 (Triassic), and the cross-sectional view of the shaft, respectively. Upon comparison of Figure 10(a) and (b), it can be concluded that although layer 15 forms the main water intake layer, its drawdown is lower than that of layer 20. The reason is that the YR had recharged the  upper aquifer in time, which reduced the impact of pumping. Figure 11 suggests that the drawdown of the radiant pipe and shaft was greater than that of the surrounding aquifer. These figures suggest that the groundwater forms a cone of depression centered on the radiant pipe and that the groundwater in the aquifer is collected through the radiant pipe to the chamber and shaft under the impact of the hydraulic gradient during the pumping process. Our results also reveal that the groundwater flowing to the seepage wells shows an obvious characteristic of three-dimensional flow.

Applicability of the model and determination of maximum pumping drawdown
This study provided evidence that the water yield of seepage wells increases with increments in the drawdown of the shaft, but the growth speed decreases ( Figure 5). The reason is that a disconnection zone (unsaturated zone) between the river and groundwater occurs when the drawdown of the shaft is increased to some extent and the pumping flow is too large, whilst the growth speed of the water yield decreases (Cognac & Ronayne 2020). Therefore, the water yield of seepage wells under different depth-reduction conditions was simulated. The results showed that the water yield initially increased then decreased, with the increase of the shaft drawdown (Figure 12), which suggests that an unsaturated zone formed between the river and groundwater.  On the one hand, when the drawdown of the shaft exceeded 15 m, the existing model was no longer applicable to these conditions. As such, the model must be improved to calculate the unsaturated flow formed from seepage wells. In addition, it is necessary to improve and simplify the calculation model so that the transient flow can be simulated when seepage wells are complex. The dynamic variation of the groundwater could also be obtained during the working process of the seepage wells. On the other hand, the maximum pumping drawdown should be less than 15 m, otherwise, it may increase the cost of water intake and cause negative environmental problems without a significant increase in water yield.

CONCLUSION
In-situ observation and the verification of the calculation model on water intake engineering in seepage wells have been poorly studied, despite seepage wells being widely used. In this study, field observation of a pumping test of seepage wells was conducted to explore the interaction among groundwater, river, and seepage wells and to evaluate the performance of the coupled seepage-pipe flow model. The pumping test results showed that the water yield of seepage wells increases, but the increased speed decreases after the drawdown of the shaft exceeds 15 m. Both the Yellow River and groundwater responded to the variation of shaft drawdown, and the timely recharge of the river weakened the impact on the upper aquifer. Our model showed a good performance in simulating water intake of seepage wells under saturated conditions with maximum root-mean-square deviations and correlation coefficients of 1.044 m and 0.966, respectively. The results of the model and experimental observations both indicated that the groundwater presented obvious three-dimensional flow characteristics during the process of seepage-well pumping. Further analysis indicated that an unsaturation zone was formed between the river and the groundwater and the maximum pumping drawdown should be less than 15 m.
There are two major limitations in this study. One is for the pumping test. Although the water level in different spatial wells was observed, it was difficult to add more observation wells due to the constraints of funding, so the actual spatial flow field under the pumping condition of the seepage well was not obtained. Another limitation is for the model. On the one hand, it is difficult for the model to simulate the transient flow when the structure of the seepage wells is complex. On the other hand, the model is not applicable to unsaturated conditions. The simplification and improvement of the calculation model will also form the focus of our further studies. Nevertheless, these results can guide the design and construction of seepage wells and can improve the rational development of groundwater, particularly in arid and semi-arid regions.