A saturated–unsaturated coupling model for groundwater flowing into seepage wells: a modeling study for groundwater development in river basins

As a major and popular groundwater extraction structure, seepage wells are often used to transfer river water into aquifers for harvesting water resources. It can help ameliorate the imbalance between supply and demand, in particular, in areas of water scarcity. Large drawdowns due to pumping may cause the river to disconnect from the groundwater and form an unsaturated zone, which seriously affects the efficiency of seepage wells. However, most of the current models of extraction structures of non-tube wells account for only saturated flow and do not consider unsaturated conditions. To address this limitation, a saturated–unsaturated coupling model was developed using the exchange flow rate between the well pipe and the aquifer as the coupling point. Moreover, the model was evaluated with physical simulation test data. The statistical results indicated that the model could estimate the drawdown and pumping rate well with root-mean-square deviations of 0.0114 m and 0.0079 L s−1, respectively, for a river with strong leakage capacity; and 0.0129 m and 0.0099 L s−1, respectively, for a river with weak leakage capacity. The critical drawdown, where the river disconnects from the aquifer, as well as variations of the unsaturated zone, was also discussed. The present study provides important information for the design of seepage wells with reasonable drawdown, while being able to predict the potential water yield, and at the same time help protect the groundwater environment.


Introduction
Groundwater is one of the most important available water resources (Qiao et al. 2020;Wang et al. 2020) and is often harvested from aquifers in river basins where a huge population often resides (Li and Qian 2013;Li et al. 2014). However, some river basins cannot supply a large quantity of water using traditional wells because of low-yield capacity or thin aquifers. Studies have shown that seepage wells constructed near rivers can effectively solve this problem (Fu 2012;Wang 2013). Seepage wells are complex groundwater collecting structures, consisting of a shaft, galleries, chambers, and radiate pipes (Wang and Zhang 2008), which can purify the river and drive river water into aquifers via the sandy layer (as shown in Fig. 1). Compared to a traditional well, seepage wells and other similar water extraction structures have many advantages, such as high production, low consumption, convenient management, and relatively low cost of water supply. Therefore, they are widely used as a means to increase water resources (Ameli and Craig 2018;Maroney and Rehmann 2017;Wu et al. 2019). Often, drawdown of the well is increased in the effort to obtain large water yields. However, it may lead to a disconnection between the groundwater aquifer and an associated river, thereby forming an unsaturated zone (Brunner et al. 2011;Cognac and Ronayne 2020;Glennon 2007;Sun 2017;Wang et al. 2009;Xie et al. 2015). In this way, the water yield would not increase significantly, but the water extraction 1 3 711 Page 2 of 13 efficiency would decrease with the appearance of the unsaturated zone (Chang et al. 2019).
Previous studies mainly focused on the problem of water intake from non-tube wells under saturated flow conditions. For these studies, models have been developed based on various assumptions. For example, non-pipe well water intake structures were characterized as an equal-flow or equal-head boundary to calculate the water yield (Hantush 1961;Park and Zhan 2003;Zhan et al. 2001). However, Chen et al. (2003) noted that the head in the well was not equally distributed due to head loss and established a coupled seepagepipe flow model based on the equivalent hydraulic conductivity. The coupled seepage-pipe flow model considered the head loss and since then it has been widely applied (Lee et al. 2010;Sun 2015;Zhang et al. 2008). Wang and Zhang (2008) established a calculation model of seepage wells based on the theory of equivalent hydraulic conductivity and a coupled seepage-pipe flow model. However, this model could only be used to calculate a single seepage well due to the need for fine subdivisions. Therefore, Wang et al. (2013) improved the model by generalizing seepage wells as a series of nodes and pipes, which made the calculation of multiple seepage wells working simultaneously possible and reflected the coexistence of multiple flow regimes scientifically. Unfortunately, this improved model cannot correctly simulate how the infiltrated recharge from the river becomes a leaching recharge when the unsaturated zone occurs.
Studies about water intake in non-tube wells under unsaturated flow conditions have rarely been reported. Relatively few studies on unsaturated seepage of groundwater have provided an important guidance for the development of the present research. Richards (1931) applied Darcy's law to unsaturated flow and established the unsaturated seepage equation, later coined Richards' equation. Neumen (1973) considered the state of saturated-unsaturated seepage flow and applied a finite element method to solve the numerical problem. Twarakavi et al. (2008) used MODFLOW (USGS, USA) and HYDRUS (Simunek et al. 1998) to simulate an infiltration experiment and a hypothetical regional-scale groundwater flow problem. The results showed that the coupling model could well solve the unsaturated-saturated flow problem. Zhu et al. (2012) established a saturated-unsaturated coupling numerical model based on the vertical flow between the saturated and the unsaturated zone. Subsequently, the saturated-unsaturated seepage flow theory has been widely used (Chen et al. 2020;Diouf et al. 2020).
The overall objective of the present study was to develop and verify a saturated-unsaturated coupling model of groundwater flowing into seepage wells. Specific objectives include: (1) to carry out a physical simulation experiment, using a sand tank, and obtain the water yield, as well as the aquifer and shaft head information during pumping of seepage wells; (2) to establish a saturated-unsaturated coupling model and simulate the water intake of seepage wells; (3) to explore the critical drawdown of seepage wells when the river disconnects from the aquifer; and (4) to reveal a variation tendency of the unsaturated zone between river and groundwater with the drawdown of the shaft. The saturated-unsaturated coupling model can simulate water intake from seepage wells adequately both in saturated and unsaturated conditions, which can help predict water yield and guide the design of seepage wells.

Laboratory experiment
A rectangular sand tank (with a length of 3.63 m, a width of 0.5 m, and a height of 1.13 m) was designed to simulate the groundwater flow into the seepage well, as shown in Fig. 2. The sand tank was made of 10 mm thick steel and 5 mm thick angle steel (reinforced on the front and back sides). Sand with particle sizes in the range from 0.25 to 1.0 mm was used to simulate aquifer media. One hundred and eighty observation holes on the front and back sides of the sand tank were designed to observe the water level of the aquifer. Four water intakes were set up at the bottom of the sand tank. Before the test, the sand tank was filled with water from the bottom to the top to remove the air in the sand tank. During the test, four water intakes were closed. A hole was made on the back of the sand tank to link the shaft. A water tank on the top of the sand tank was used to continuously recharge the sand tank during the test. The setup also included an overflow tank, to keep the water level at the top of the sand tank constant. When the water level exceeded the overflow height, the water flowed out through the overflow hole. On top of the sand, two types of permeable materials, namely, color banded cloth and gauze, were applied to simulate riverbeds with different leakage capacities. A color banded cloth represented a river with weak leakage capacity and gauze represented a river with good leakage capacity. A layer of pebbles was laid on top of the permeable materials to simulate the riverbed media. The bottom of the sand tank was paved with a layer of pebbles for better water filling. Tensiometers were placed at the bottom of the river to record the occurrence of the unsaturated zone. Distilled water was used in the experiment.
Seepage wells were designed with five galleries, a single shaft, and five chambers, and five radiant pipes were installed in each chamber (Fig. 2). The length of each radiate pipe was 25 cm, four of which had an elevation angle of 45°, while the angle of the fifth was 90°. The length of each gallery between the two adjacent chambers (with a side length of 10 cm) was 50 cm. The experiments involved three different drawdowns of the shaft for each type of river with different leakage capacities, giving a total of six test schemes. The flow rate was recorded by a measuring cylinder and a triangular weir (Table 1).
After installing all the required instruments, the sand was saturated slowly from the bottom of the sand tank to the top by opening water intakes at the bottom. The saturation process was carried out over 12 h to remove air in the sand. After saturation, the water level in the tank was 8 cm higher than the top surface of the sand. Before the start of the experiment, the water intake at the sand tank bottom was closed. The water used in the experiment was supplied from the water tank. Moreover, the initial head of the observation holes was recorded. Then, seepage wells began pumping. When the pumping was stable, the flow rate of the seepage wells and the head of the observation holes were recorded. Seepage wells continued pumping until the final drawdown test was completed.
After finishing the experiment, the hydraulic conductivity of sand was measured by the standpipe method (Chen 2000), and moisture characteristic parameters of sand were determined by the sandy funnel method (Li et al. 2011). These measured parameters can provide input information for the numerical model.

Saturated-unsaturated coupling model
The coupled seepage-pipe flow model was developed by (Chen et al. 2003) to calculate the groundwater flowing to horizontal well. However, the differential equation of saturated seepage no longer applies when the unsaturated zone appears. Therefore, to describe the unsaturated groundwater flow caused by pumping, the governing equation for the saturated groundwater flow should be modified into the "mixed form" of Richards' equation, as follows (Thoms et al. 2006): where H is the groundwater head; K(ψ) is the unsaturated hydraulic conductivity which is a function of the pressure head ψ; Θ(ψ) is the effective saturation as a (1) The description of the experimental equipment function of the pressure head; C(ψ) is the specific moisture capacity; W denotes the volumetric flux per unit volume representing source(s) and/or sink(s); S s is the specific storage, and t is time.
As the aquifer becomes fully saturated, the Richards' equation converges to the general flow equation used by MODFLOW (Harbaugh 2005) for saturated subsurface flow: where K x, sat is the saturated hydraulic conductivity (m d −1 ).
The soil characteristic functions K(ψ), Θ(ψ), and C(ψ) have been represented in the literature studies by several different empirical and theoretical methods (Haverkamp et al. 1977;Millington and Quirk 1961;Mualem 1976;van Genuchten 1980;White and Broadbridge 1988). In this study, van Genuchten soil characteristic functions were used. The functions are represented as follows: where Θ r is the residual soil saturation; n v , m are the soil parameters (m = 1 − 1/n v ) representing the degree of poresize uniformity; and α is the parameter representing the inverse characteristic length of the soil pores. A detailed description is described in a literature study (van Genuchten 1980). Equation 1 combined with other definite conditions can be solved by the finite difference method. The saturated-unsaturated coupling process in MODFLOW (developed by USGS, USA) provided details on how to incorporate the pressure-head dependency within the hydraulic conductivity, as well as storage terms for the unsaturated cells, and the capability and availability of incorporation (Thoms et al. 2006). Munson et al. (2002) described the flow in the pipe using the Darcy-Weisbach equation (Eq. 6) which is suitable for laminar and turbulent flow: Then, the flow velocity and volumetric flow rate are expressed as: where l is the length of the flow path; g is the gravity acceleration (m s −2 ); u is the average flow velocity in the pipe (m s −1 ); A is the cross-sectional area of the pipe (m 2 ); d is the diameter of the pipe (m); and f is the friction factor determined by different equations according to different flow regimes (Colebrook and White 1937;Munson et al. 2002).
Therefore, the exchange flow rate between the aquifer and pipes and the flow rate in the pipes is expressed as follows  where Q e is the exchange flow rate between the aquifer and pipe; Q p is the flow rate in seepage wells (m 3 d −1 ); C is the conductance of seepage wells (m 2 d −1 ); H p is the head of pipe flow (m); π is the mathematical constant pi (unitless); ΔH is the head loss (m); v is the average flow velocity in the well pipe (m s −1 ); e is the coarseness degree of the inner wall of the well pipe, and υ is the kinematic viscosity of water (m 2 d −1 ).
Based on the discretized grid of the seepage well, the water balance is satisfied at the radiate pipe node i when the flow reaches the steady state, as follows: where j is the number of well nodes connected to well node i, n p is the total number of well nodes are connected to well node i, and Q p,i.j is the flow rate from the well node i to the well node j. Then, the seepage flow in the aquifer is coupled with the pipe flow in the radiate bore by the exchange flow between the radiate bore and the aquifer ( Table 2).
Three variables (the head difference in the pipe, the average flow velocity, and the Reynolds number) are dependent on each other. Therefore, an iterative scheme must be used to solve these equations. Further, these equations are coupled into Eq. 2 by treating the exchange flow rate as the sinks of finite difference cells.
This mathematical model can be solved numerically. At the start of computation, the elevation of the top of the well nodes is set to the initial head of the well nodes, and the initial flow regime is assumed to be laminar. Then, exchange flow rates are computed using the initial head of the finite difference cells. Furthermore, the distribution of the flow rates through the seepage well is determined. Based on the distribution of flow rates and the diameter of the well pipe, the Reynolds numbers can be obtained, which are used to select the appropriate equation to calculate the new heads in the seepage well. The new heads of the seepage well can be calculated, and are used to repeat the calculation using the above mentioned Equation. The exchange flow rates after these iterations are added to finite difference equations as a sink to solve the new heads of the aquifer. Then, the new exchange flow rates are determined according to the new heads of both the aquifer and the well nodes, which are used to repeat the above mentioned calculation. The iteration for solving groundwater heads stops when the maximum absolute head differences of two successive iterations meet the precision criteria. Finally, the saturated-unsaturated coupling model is used to solve the problem of the seepage well with a large drawdown in which the head loss inside the well pipe and unsaturated flow beneath the streambed layer are considered. A detailed description can be found elsewhere (Wang et al. 2017;Wang and Zhang 2008).

Construction of the numerical model
The finite difference method from MODFLOW was employed to solve the problem of water intake in seepage wells. The sand tank was divided into 91 rows (along the X axis) and 15 rows (along the Y axis) on the plane and 25 layers vertically. The river at the top and the sandy pebble at the bottom were set as two separate layers (the thickness was 8 cm), and the middle sand layer was evenly divided into 23 layers. The galleries and chambers were set on the 20th layer. The bottom and sidewalls of the sand tank were set as zero flow boundaries as there was no water exchange with the outside. Rivers with strong and weak leakage capacity were set as constant head and river boundaries, respectively. The head of the shaft was set as the constant head boundary during pumping (as shown in Fig. 3). From this, a hydrogeological conceptual model was established.
Before the simulation, seepage wells must be treated reasonably to complete the coupling of seepage and pipe flow. First, the column and row numbers of each node in the MODFLOW grid were counted based on the top view of the seepage wells. Second, nodes and pipes of seepage wells were numbered using the specific numbering rules to ensure the connectivity of water flow in seepage wells. Then, the elevation of each node was calculated according to the position of chambers and elevation angle of radiant pipes. Noteworthy, nodes of one radiant pipe must be continuous in layers, otherwise, the elevation of model layers or nodes of radiant pipe should be adjusted appropriately. Finally, a.cfp file that can be recognized by MODFLOW must be created, and the format can be found in Shoemaker et al. (2008). After finishing all of the above work, MODFLOW can now be used to simulate the water intake of seepage wells.
To evaluate the performance of the model, statistical parameters, such as relative errors (RE), correlation coefficient (R 2 ), and root-mean-square deviation (RMSD), were used (Zheng et al. 2020). Higher values of R 2 and smaller values of RE, RMSE show a better performance of the model. The equations are shown as follows: where N is the total number of the simulation or observation data; i is the ith simulation or observation data; sim denotes the simulation data; obs is the observation data; obs is the average of the observation data; sim is the average of the simulation data. Finally, the contents of this study and the relationship among different parts are shown in Fig. 4.

Experimental results
The drawdown of the shaft and the flow rate of the seepage wells were recorded throughout the experiment (Fig. 5). For the simulated river with a strong leakage capacity, the flow rate increased linearly with the increase of the drawdown in the shaft, which indicated that the water yield of the seepage wells could be guaranteed due to a sufficient recharge source. However, for the simulated river with a weak leakage capacity, the previously increasing rate of flow rate decreased during the third stage (B2-D3). This was attributed to the limited recharge capacity of the river preventing the aquifer from being recharged quickly enough, which Fig. 4 The framework illustrating the logic of the study led to the river being disconnected from the groundwater, thereby forming an unsaturated zone. The unsaturated zone was also observed using tensiometers during this period. Furthermore, for the same drawdown, the flow rate of the river with a strong leakage capacity was larger than that of the river with a weak leakage capacity. Figure 6 shows the drawdown contours of the aquifer in each scheme when the flow rate became stable. The drawdown of the aquifer increased with the increase of the shaft drawdown and was always smaller than that in the shaft. Comparative analysis of the schemes in B1 and B2 indicated that that the drawdowns in the B2 schemes were larger than those in B1. This suggested that the timely recharge of rivers mitigated the decrease of aquifer levels. It was also found that the drawdown of the aquifer increased with proximity to the seepage wells. These results indicate that the water yield of the seepage wells mainly came from direct recharge of the nearby aquifer. Figure 7 shows that the drawdown of the aquifer increased with the increment of drawdown for the shaft but it was smaller than that in the shaft. The drawdown of the aquifer first increased and then decreased from top to bottom and reach the maximum near the seepage wells. It indicates that the drawdown of the upper aquifer was mitigated by the recharge of the river. Water to the seepage wells was also supplied from aquifer under the seepage wells. By combining combined analysis of Figs. 6 and 7 deduced that groundwater was supplied from the aquifer surrounding the seepage wells, and the groundwater flowing to the seepage wells showed clear three-dimensional flow characteristics. Six samples were obtained from different positions in the sand tank and hydraulic conductivities of the sand layer were measured by the standpipe method (shown in Tab.1). The hydraulic conductivities ranged from 31.84 to 47.82 m/d, with an average of 40.76 m/d. As shown in Fig. 8, the desorption process and the moisture absorption process were fitted using the soil characteristic function. A good hysteresis was observed between the sorption and desorption curves, and soil characteristic function fitted the measured data well. Therefore, the moisture characteristic parameters of sand were determined (Table 2).

Model and simulation evaluations
Six schemes were simulated using the saturated-unsaturated coupling model, and the model was verified by experimental data using the sand tank setup. Table 3 presents the simulated and experimental results of the flow rate and their errors. Clearly, the calculated and observed flow rates were in good agreement and that flow rates increased with increasing drawdown of the shaft. The largest absolute error between the calculated and the experimental flow rates was 0.0148 L/s, and the relative error was less than 4%.
To further evaluate the performance of the model, the simulated results (including the water yield of seepage wells and the drawdown of the observation holes) and the experimental data were compared in Fig. 9a and b. The model effectively simulated the flow rate of the seepage wells and drawdown. For the river with strong leakage capacity, the flow rate and drawdown had R 2 values of 0.99 and 0.98, respectively. The values of RMSD were 0.0079 L s −1 and 0.0114 m, respectively, and values of RE for flow rate and drawdown were 1.96% and − 0.13%, respectively. For the river with weak leakage capacity, the R 2 , RMSD, and RE for flow rate were 0.99, 0.0099 L s −1 , 2.54%, respectively, and for drawdown were 0.98, 0.0129 m, and − 2.58%, respectively. These results show the feasibility of using this model to simulate water intake in seepage wells both in saturated and unsaturated conditions. The final model parameters are K v = 45 m d −1 , K h = 29.5 m d −1 , and S y = 0.2.

Determination of the critical drawdown
Owing to the limitations of the experimental conditions, only six representative schemes were designed to simulate the water intake of seepage wells. Disconnection between the river and groundwater aquifer appeared in the scheme B2-D3, forming an unsaturated zone. However, the existing coupled seepage-pipe flow model was only suitable for saturated seepage flow, and the unsaturated unit was treated as an inactive unit without water exchange, leading to calculation errors. Therefore, it was of great significance to determine the critical drawdown of the shaft when the river Fig. 5 The relationship between experimental flow rate and drawdown. B1 and B2 represent rivers with strong and weak leakage capacities, respectively, and D1 to D3 represent the first to third drawdown, respectively and groundwater disconnected. To determine the critical drawdown, the existing coupled seepage-pipe flow model and the saturated-unsaturated coupling model were used to calculate the water yield of seepage wells under different drawdowns of the shaft. Figure 10 illustrates that the calculation results of the two models were very close when the drawdown of the shaft was less than 33.588 cm, which indicated that both models can be used to calculate the water intake of seepage wells under saturated conditions. However, the water yield calculated using the former model rapidly decreased to 0 once the drawdown went beyond 33.588 cm. When the saturated-unsaturated coupling model was used to calculate the water yield, the water yield continued to increase at first (but the growth rate decreased) after the drawdown of the shaft was input to be larger than 33.588 cm, and finally became stable at drawdowns above 44 cm. Compared to the calculation results of the two models, the critical depth of disconnection between the river and groundwater was finally determined to be 33.588 cm. Noteworthy, the reasonable drawdown of the shaft is very important in practical engineering. If the drawdown of the shaft is too small, the water yield cannot meet the demand. However, if the drawdown of the shaft is too large, the water yield of the seepage well does not increase significantly, and the water intake efficiency decreases.

Variation law of the unsaturated zone
To establish the variation tendency of the unsaturated zone, the saturated-unsaturated coupling model was used to simulate the water intake of seepage wells under different drawdowns of the shaft. Figure 11 shows the variation of the unsaturated zone in the plane and profile. The critical drawdown was determined to be 33.588 cm. Therefore, the drawdown of the shaft was set to 34, 36, 38, 40, 42, 44, 46, and 48 cm, respectively (shown in Fig. 10). As illustrated in Fig. 11, the drawdown of the aquifer increased with increasing drawdown of the shaft. Furthermore, it can be seen that the unsaturated zone area at first was very small and appeared in the center of the first layer of the model.  The area of the unsaturated zone increased and expanded to the lower aquifer with the increase of shaft drawdown. Subsequently, the unsaturated zone took up all the space in the first layer and reached the fourth layer. The unsaturated zone expanded from the center of every layer to the surrounding area and from top to bottom of the model, showing a coneshaped spatial distribution.

Conclusion
Seepage wells are major and complex groundwater extraction structures. A fundamental understanding of unsaturation caused by pumping in seepage wells is vital for guiding the design of seepage wells. In the present study, a physical simulation test of a sand tank for water intake of seepage wells was designed. The results showed that the water yield of a river with a strong leakage capacity was larger than that of a river with a weak leakage capacity for the same drawdown of the shaft. There was no unsaturated phenomenon observed in the test scheme of the river with strong leakage capacity, and the relationship between water yield and drawdown of the shaft was linear. However, the growth rate of the water yield began to decrease after the appearance of an unsaturated zone, which was reflected in the river scheme with a weak leakage capacity. The experiment also revealed that the groundwater flowing into seepage wells showed clear characteristics of three-dimensional flow.
A saturated-unsaturated coupling model was developed to describe how groundwater flows into seepage wells both in saturated and unsaturated conditions. The model was verified by experimental data. The validation results for the drawdown and water yield of the river with strong leakage capacity showed root-mean-square deviations, RMSD, equal to 0.0114 m and 0.0079 L s −1 , and the correlation coefficients, R 2 , were 0.98 and 0.99, respectively. For the river with weak leakage capacity, the R 2 and RMSD of flow rate were 0.99 and 0.0099 L s −1 , and for drawdown were 0.98 and 0.0129 m, respectively. These results suggest that the model can appropriately simulate water intake in seepage wells both in saturated and unsaturated conditions, and indicated a good performance of the model. Furthermore, the critical drawdown of the shaft at the time of disconnection between the river and groundwater aquifer (forming an unsaturated zone) was determined to be 33.588 cm. The unsaturated zone expanded from the center of each layer to the surrounding area and from top to bottom of the model, showing a cone-shaped spatial distribution.
In summary, the saturated-unsaturated coupling model described herein showed a great potential to help predict  . 9 Comparison of simulated and observed results. Q obs and Q sim are the observed and simulated flow rates of seepage wells, respectively, and S obs and S sim are the experimental and modeled drawdown of observation holes, respectively. The black lines in the graph represent the function y = x. The blue and red circles represent rivers with weak and strong leakage capacity, respectively Fig. 10 The relationship between the simulated flow rate and drawdown. The red and blue lines are the results calculated by the former model and the saturated-unsaturated coupling model, respectively. The green square represents the water yield at the critical drawdown