Quantification of seepage in a multi-layered disconnected river-aquifer system

Quantification of seepage in disconnected river-aquifer systems is significant for local water management and groundwater pollution control, especially in areas with water shortage or contamination. The vadose zone under riverbeds usually exhibits a multi-layered structure, particularly when paved with low permeability liners. To evaluate the impact of engineering solutions to seepage under such conditions, we proposed an approach by combining GIS and the minimum flux in saturation layer (MFSL) method. MFSL can calculate the stable seepage rate by assessing the dominant low permeability layers (including but not limited to the liners) in multilayered disconnected river-aquifer systems. We used MFSL to calculate local seepage rate, and used GIS to extend the results to a regional scale. The reliability of MFSL is discussed by comparing the results with the double ring infiltration test, the numerical simulation by HYDRUS, and the methods of stream package in MODFLOW, including its improved form. A case study was conducted in the Yongding River with river-aquifer seepage calculated under various conditions, including different river water levels (i.e., under the designated water level, drought stage level, flood stage level and flood inundation level) and with/without low permeability liners (i.e., ecological membranes or geomembrane). Results showed that low permeability liners could greatly reduce the seepage rate. However, if an unlined inundation area exists, the seepage rate may increase greatly. The results indicated that the proposed method was reliable and convenient for calculating long-term, large area seepage in disconnected river-aquifer systems especially those paved with liners.


Introduction
Infiltration occurs when vertical seepage of surface water flows through riverbeds to the groundwater aquifer [1] . This seepage can result in the loss of valuable surface water resources and harm the riverine ecology, if short of water recharge [2] . Also, if the river water is contaminated, the polluted water may seep into the aquifer and put groundwater at risk [3] . Therefore, quantification of seepage in a river-aquifer system is crucial for both the protection of aquatic ecosystems and control of groundwater contamination.
Connected, transitional and disconnected are terms used to identify the interaction between surface water and groundwater [4] . A disconnected system refers to a surface water and groundwater system where the infiltration rate is independent of the water table position or identified by an unsaturated zone under the river [4] . When clogging layers exist (hereafter we refer to the dominant low permeability layers as liners) in the riverbed and the groundwater table is far below the riverbed, a disconnected system will usually develop.
The seepage rate in river-aquifer systems can be calculated theoretically according to Darcy's law, as has been done in the MODFLOW [5] stream package, whereas for a multi-layered system, the saturated hydraulic conductivity of the riverbed cannot be obtained directly similar to the homogeneous system. Fok [6] used the effective hydraulic conductivity, i.e., the thickness weighted harmonic mean hydraulic conductivity to account for the multi-layered impacts. Note that this approach assumes all soil layers between river and aquifer are saturated, which does not apply to disconnected riveraquifer systems where an unsaturated zone, usually induced by the upper clogging layer, exists. Bouwer [7] included a critical pressure head at the base of the clogging layer in Darcy's law, as did Rovey [8] . The main difficulty is how to determine the value of that critical pressure head. Osman and Bruen [9] suggested that the maximum suction head h mis at the base of the clogging layer should satisfy the condition that the calculated seepage rate according to the Darcy's law from the above clogging layer should equal the unsaturated hydraulic conductivity of the lower layer. Based on this assumption, they proposed an improved form for calculating stream seepage in MODFLOW. Using the above approach, Crosbie et al. [10] combined the simplified infiltration model and the field survey for estimating infiltration from Billabong Creek in south-eastern Australia. Wang et al. [11] also used similar assumption to modify the Green-Ampt model for estimating infiltration in fine soil with coarse sand interlayer.
The above assumption has a solid theoretical basis, which can improve the accuracy of estimation and avoid the complicated numerical simulations. However, it is important to note that it requires the quantitative description of the unsaturated hydraulic conductivity for the soil under the clogging layer, which may not be convenient for practical use.
Numerical models based on Richards' equation, describing variably saturated flow, have developed rapidly in recent years, For example, the HYDRUS series [12,13] and FEFLOW [14] are able to simulate saturated-unsaturated flow with water interaction between river and aquifer. Using HYDRUS, Wang et al. [15] investigated ponded infiltration into fine-textured soils with coarse interlayer, which is similar to the seepage in a disconnected groundwater-aquifer system. However, this numerical approach requires even more detailed hydraulic parameters, which are difficult to obtain in practice. Besides, numerical simulations are time consuming and require large computer memory when applied to large areas.
Multi-layered structures commonly exists in/under the riverbed due to both natural and artificial processes. Mostly, studies have not considered the impact of multilayered structure on seepage, including the possibility that other layers can be dominant in controlling seepage in addition to the surface layers in riverbeds. This raises the question of how to determine seepage when there are several clogging layers. For this purpose, Mao and Shang [16] proposed the method of minimum flux in saturation layer (MFSL) to assess clogging layers and to calculate the seepage in a multi-layered soil system. The main aim of this paper is to couple MFSL with the powerful spatial data analysis tool, GIS, to evaluate the impact of engineering solutions to seepage through a multi-layered vadose zone in the regional disconnected river-aquifer system. A case study was conducted on the Yongding River under various scenarios and the accuracy of the results is discussed by comparing with the other commonly used methods.

Site description
The Yongding River is the largest river in the Beijing municipality, with a length of 170 km and a drainage area of 3168 km 2 for the river reach in Beijing. Overexploitation of water resources upstream led to the reach in Beijing drying out in the 1980s. To restore the ecology and environment, in 2009 the Beijing government launched the program called, Green Ecological Corridor in Yongding River [17] . Given the shortage of natural upstream inflow [18] , reclaimed water is to be used as the water supply to theYongding River. The program plans to fill six lakes with an area of 6.8 km 2 along the Yongding River, with a combined channel length of 50 km and area of 2.7 km 2 . To evaluate seepage loss and control groundwater contamination, quantification of river-aquifer seepage is required under a range of possible circumstances. Here we focus on the river reach passing the Beijing's western suburbs, with the length of 37 km from Sanjiadian Village in Mentougou District to Efang Village in Daxing District (Fig. 1).
The study area includes four lakes, Mencheng Lake, Lianshi Lake, Xiaoyue Lake, Wanping Lake, and two channels, Channel Shijingshan and Channel Daxing, and the riverbank area where inundation is possible in a flooding period (generally along Mencheng Lake and Lianshi Lake). Three boreholes were drilled in the study area in order to understand the geological conditions in the vadose zone [19] . For seepage calculation in each subregion, geological information was obtained from the nearest borehole (Table 1). It is recognized that the location and limited number of boreholes is not ideal because geological characteristics in each sub-region are not evenly distributed. Therefore, we expect more detailed information will improve the estimation accuracy in the future.

Principle of MFSL
To calculate the seepage rate in a multi-layered system that may have more than one clogging layer, Mao and Shang [16] proposed the method of minimum flux in saturation layer (MFSL). The principle of this approach is briefly illustrated below.
Assuming a soil column has n soil layers, each layer has a thickness of L i (m) and saturated hydraulic conductivity of K si (m$d -1 ), i = 1,2,...,n. The bottom of the soil column is at a fixed groundwater table and initially it is at the hydrostatic condition, i.e., with the same soil water potential and no water flow along the column. Then a fixed water pressure head H w (m) is applied to the top of the soil (Fig. 2).
For fine-textured soil overlying a relatively coarse layer, previous study has shown that the wetting front tends to pause temporarily at their interface, then the infiltration rate slows down to a constant value [21] , and the lower coarse layer remains unsaturated with stable moisture content [11] . However, for a soil column with several layers, it is difficult to tell at which interface this phenomenon would occur and if stable infiltration could be permanently achieved. The following method was proposed to assess this for soil column with multiple interfaces. Assuming the first k layers are saturated by seepage and the sublayer, k1, remains unsaturated, when neglecting the water suction at the interface of saturated and unsaturated soil layers, the infiltration rate through the saturated zone, q sk (m$d -1 ) can be calculated according to Darcy's law, i.e., Fig. 1 The location of the Yongding River reach passing the Beijing plain area and schematic of river subdivisions Note: a , The designated water level is estimated by investigation and discussion with the local engineer, both to satisfy the ecological water demand inside the river and to ensure there is enough water for scenery purpose; b , obtained from the borehole data monitored by Beijing Institute of Geological and Prospecting Engineering; c , obtained from the empirical value [20] .
where K k is the effective hydraulic conductivity of the layers 1-k. K k can be calculated from [6] There exists an m that is less than or equal to n and meets Then it can be concluded that the saturated wetting front of infiltration will actually temporarily pause at layer m, and the seepage reaches a stable state with the rate of q sm (m$d -1 ). The detailed proof of this conclusion can be found in Mao and Shang [16] .
MFSL is used here for calculating stable seepage rate in each subregion of the Yongding River, with the information of liners mentioned above and the hydrogeological information in the vadose zone shown in Table 1. Groundwater level is assumed to be at the bottom of the boreholes. This is in accord with data presented in the Beijing Water Resources Bulletin that averaged groundwater depth in the plain area of Beijing at over 20 m in recently years. Also, when low permeability upper layers exist, the wetting front would generally stay in the upper layers and reach a stable seepage state. According to Eqs.
(1) -(3), this would not be influenced by groundwater fluctuations in the lower layer. This also matches the definition of a disconnected system, in which the infiltration rate is independent of the groundwater table position, as indicated by Brunner et al. [4] .
Also, it is worth noting that MFSL is used for calculating the stable seepage rate. It does not address the unstable seepage process in the early period. The error for long-term total seepage caused by this method will be discussed later.

Method for calculating regional river-aquifer seepage
To obtain the regional river-aquifer seepage, the seepage rate at specific locations was first calculated by MFSL, then the local seepage calculations were extended to the regional scale by GIS. The procedure was as follows: (1) Input the data of study area to GIS, including the outline, area and water head distribution of surface water, the soil layered structure and hydraulic characteristic of each drilling point, the river bed lining plans.
(2) Divide the study area into subregions, i.e., lakes, channels and related river banks according to the hydrogeology information and surface water head. Then calculate the area of each subregion using GIS.
(3) Calculate the stable seepage rate of each subregion by MFSL and then multiply the results by subregion area to obtain the regional river-aquifer seepage.
(4) If applicable, calculate the seepage under different water depths (e.g., different seasons or hydrological years) and under different low permeability liners.

Other methods for river-aquifer seepage calculation for comparison with MFSL
The in situ double ring infiltration test, numerical simulation by HYDRUS, the stream package in MOD-FLOW and its improved form were applied to test the reliability of MFSL results. A brief introduction of these follows.
(1) The double ring is a commonly used method for measuring the field infiltration rate. It includes two nested steel rings set into the soil layer, with the outer ring only serves as a screen to minimize lateral outflow and the inner ring to measure the seepage.
(2) The HYDRUS series [12] , is numerical software for simulating water, heat and solute movement in variably saturated porous media. It uses Richards' equation (Eq. (4)) for simulating water flow, which was solved numerically using Galerkin type linear finite element schemes.
where h is the soil water matrix potential (negative) in the unsaturated zone or pressure head (positive) in the saturated zone; t is time; z is the vertical coordinate, positive downward; q is the volumetric soil water content; and K is hydraulic conductivity. HYDRUS provides several functions to depict the relationships among q, K and h. Here the commonly used Van Genuchten [22] model for soil water retention curve and the associated model for soil hydraulic conductivity [23] was adopted.
(3) The stream package in MODFLOW calculates seepage in a river-aquifer [5] as where q s is the seepage rate, K b is the saturated hydraulic conductivity of riverbed, M is the thickness of riverbed, H r is the elevation of the river water surface, and h 0 is the aquifer head.
(4) Osman and Bruen [9] proposed an improved form for calculating river seepage in MODFLOW, where Y bot is the elevation of the riverbed (i.e., the elevation of the base of the clogging layer). Thus h mis can be calculated by iteratively solving Eq. (7), where K 2 (h mis ) is the unsaturated hydraulic conductivity of soil under the clogging layer at suction head h mis , which can be calculated from a soil hydraulic property model, e.g., the van Genuchten-Mualem type model [24] .

Results of river-aquifer seepage calculated by MFSL and GIS
As an important influencing factor for river-aquifer seepage, surface water pressure head (i.e., water depth) was estimated for each season according to the annual rainfall, the artificial recharge plan and the field survey ( Table 2). Table 3 shows the calculated results for three cases, without liners, and with EM or G3, on the riverbed. The total annual seepage without a liner was estimated to reach 1.84 Â 10 10 m 3 , which is far above the annual runoff of the Yongding River. Actually this value could not be reached because the river would dry out as soon as there is not enough inflow. In other words, the designated water level cannot be maintained in natural conditions without a liner. For the two cases with liners, the seepage can be several orders of magnitude smaller than without a liner. Also, seepage decreases proportionally with the decrease of the hydraulic conductivity of the liner (data not shown). When the hydraulic conductivity of the geomembrane reaches 10 -11 m$d -1 , the total annual seepage is only 1.10 Â 10 2 m 3 , over 7 orders of magnitude smaller than the case without a liner. For G1 and EM, the total annual seepage are  1.10 Â 10 4 Flood stage water level 3.78 Â 10 -6 6.00 Â 10 -6 6.00 Â 10 -6 6.00 Â 10 -6 3.78 Â 10 -6 3.78 Â 10 -6 Drought stage water level 2.11 Â 10 -6 4.33 Â 10 -6 4.33 Â 10 -6 4.33 Â 10 -6 2.11 Â 10 -6 2.11 Â 10 -6 1.1 Â 10 8 m 3 and 7.9 Â 10 8 m 3 , respectively. These values are of the same magnitude as the results of Hao et al. [25] , showing that the total annual seepage in the Yongding River could be 3.06 Â 10 8 m 3 when the Middle Route Project for South-to-North Water Transfer is in service.

Double ring infiltration test
According to the in situ double ring infiltration test data obtained near Lianshi Lake and Efang Village, the stable infiltration (seepage) rates were 2.15 and 1.50 m$d -1 , respectively [19] . MFSL was applied to calculate the stable seepage rate with the same water head (0.1 m) for double ring test, and with the hydrogeological data of the nearest borehole (Table 1). MFSL calculated results were 1.53 and 2.6 m$d -1 . The difference between observed and estimated value was 0.62 m$d -1 near Lianshi Lake and 1.10 m$d -1 in Efang Village. It should be noted that the hydrogeological information at the specific experimental points is not available, therefore this comparison can only be used as a general comparison of the seepage monitored and calculated in that region. Considering the spatial heterogeneity, the calculated data are in reasonable agreement with the monitored data because they are of the same magnitude and did not deviate greatly.

HYDRUS
The seepage at Mencheng Lake with an ecological membrane at flood stage (i.e., water depth is 3 m, Table 2) was simulated by HYDRUS. The associated soil hydraulic parameters in the numerical simulation are shown in Table 4.

MODFLOW package and its improved form
In the MODFLOW stream package, when the aquifer head drops below the riverbed, i.e., the disconnected river-aquifer system discussed here, the head of aquifer, h 0 in Eq. (5), cannot be used directly. Instead it is assumed to be equal to the elevation of the bottom of the riverbed (i.e., the bottom of the clogging layer). While for the improved form MODFLOW developed by Osman and Bruen [9] , h mis , the suction under the clogging layer is considered as shown in Eq. (6) and obtained by solving Eq. (7).
The difficulty with these two methods is how to determine the clogging layer(s). With the same case as used for HYDRUS simulation, there are large increases in hydraulic conductivity from layer 2 to layer 3 and from layer 3 to layer 4. It is most likely that the saturation layers are layers 1 and 2, or layers 1-3. Given that the clogging layers may include several layers, the saturated hydraulic conductivity K b in Eqs. (5) and (6) is calculated according to Fok [6] .
The calculation results of the left and right hand side terms of Eq. (7) are shown in Fig. 3. Obviously the intercept points of right hand term and the Y axis indicate the solution to Eq. (5), i.e., the solution by MODFLOW stream package. The intercept points between the left and right hand terms indicate the solution h mis to Eq. (7), i.e., the improved form by Osman and Bruen [9] .

Results comparison
The numerical simulation results calculated using HYDRUS (Fig. 4) show that the seepage rate at the earlier stage is much higher than the stable seepage rate. However, it reaches the steady seepage stage quickly (only after about 0.2 d). Therefore, errors of total seepage caused by ignoring the initial unstable seepage rate, in the MFSL method are acceptable considering that the total seepage period is much longer than the unsteady seepage stage. HYDRUS results show that the seepage rate drops sharply and then fluctuates greatly over the first 0.05 d. The fluctuation is possibly caused by the numerical oscillation [26] , especially for our case when the wetting front passes through layers 1 and 2, and there are large difference in soil texture. The stable seepage rate calculated by HYDRUS is the lowest compared with all the other methods. Among results of other methods, the MFSL result and the MODFLOW package result assuming Note: q r and q s are the residual and saturated water content, respectively, a and n are empirical parameters, L is the pore-connectivity parameter, and K s is the saturated hydraulic conductivity. These parameters were obtained from representative values for specified soil texture and with model calibrations [12,20] . layers 1 and 2 to be clogging layers are equal, and both are the closest to the HYDRUS result, being only about 41% higher. This difference could be caused by the assumptions made in the MFSL and MODFLOW package, i.e., the layer where wetting front stops is fully saturated and the soil water potential is zero at the bottom of the clogging layer. Considering the results simulated by HYDRUS, the pressure head (i.e., the soil water potential when it is negative) at the base of the clogging layer is -0.41 m, which is lower than the assumed data (i.e., 0 m) by the MFSL and MODFLOW package, although it must be noted that a lower soil water potential could increase the pressure gradient and thus the seepage rate. However, in the meantime the hydraulic conductivity might not be fully saturated (especially in the lower part of the clogging layer). These two coupled effects result in a decrease of seepage rate in HYDRUS compared with that of MFSL and MODFLOW packages, which assumes the clogging layers are fully saturated. The value from MODFLOW package were 0.27 and 0.52 m$d -1 assuming layers 1-2 and 1-3 are the clogging layers, respectively. The result assuming layers 1-2 are the clogging layers is equal to the MFSL result. However the result assuming layers 1-3 are the clogging layers is 174% higher than the HYDRUS result, indicating the wrong identification of clogging layers may induce a large error in seepage estimation. The results from the improved MODFLOW package show that soil water potential at the base of the clogging layer are -1.96 and -0.52 m when assuming layers 1-2 and 1-3 are clogging layers, respectively (Fig. 3), while the corresponding seepage rates are 0.42 and 0.55 m$d -1 (Fig. 3). It is surprising that although the improved MODFLOW package considers the soil water suction at the base of the clogging layer, it does not improve the calculation accuracy compared with the original MOD-FLOW package. In our opinion, this is mainly because the improved MODFLOW package uses the saturated hydraulic conductivity of the clogging layer when calculating the seepage rate based on the Darcy's law when considering a  higher pressure gradient (because it considers the water suction at the base of the clogging layer), which potentially overestimates the seepage rate. This also indicates that caution is needed when calculating seepage in a multilayered river-aquifer system with a low permeability liner using the improved MODFLOW package.
It must be noted that both versions of the MODFLOW package do not provided advice on selecting clogging layers. Therefore, large hidden error exists in these methods when an inappropriate clogging layer is selected.
In summary, taking HYDRUS results as the benchmark, most of the calculation methods overestimated the infiltration rate. The results of MFSL and the original MODFLOW package are the closest to those of HYDRUS when the clogging layer is known. If the clogging layer is unknown, MFSL could be the most reliable compared with the results of the two versions of the MODFLOW. Also, MFSL only needs the saturated hydraulic conductivity of each layer, without a requirement for complex computations, and detailed input variables and parameters in the numerical method and the improved MODFLOW package. This indicates MFSL is most useful for calculating large area river-aquifer seepage.

River seepage under flood condition
To investigate the seepage quantity under flood condition (assuming it lasts for 2 months, i.e., July and August), we calculated the total seepage with G3 on the riverbed as mentioned in Section 2.3. According to field investigation, two lakes, Mencheng Lake and Lianshi Lake, will overflow their banks into surrounding areas where there are no liners. Figure 5 shows that the seepage in both riverbanks without liners is much higher than the rivers or channels with liners. If flood inundation occurs in summer, the annual seepage would increase to 7.46 Â 10 8 m 3 with liner (G3) on the riverbed. The value is still lower than the total seepage (1.84 Â 10 10 m 3 ) with no liners on the riverbed, yet four orders of magnitude higher than the total annual seepage (1.10 Â 10 4 m 3 ) for the case with the same liner at flood stage (i.e., without flood inundation). This is mainly caused by the high seepage intensity in the unlined riverbanks. Therefore, special attention needs be given to reducing water seepage under flood inundation. For example, some seepage proof liners could be installed in areas subject to submergence.

Conclusions
This paper proposes a method for calculating seepage through a multi-layered vadose zone in disconnected riveraquifer systems, especially for rivers paved with low permeability liners. The MFSL method was applied to calculate the seepage rate at each specific location according to the river level and the hydrogeological properties of the vadose zone. Then the total amount of seepage in the research area was obtained by GIS based on regional information.
A case study was conducted on the river-aquifer seepage of the Yongding River, under various low permeability liners and different water levels. The calculated seepage in the Yongding River would reduce significantly if low permeability liners were installed in the riverbed. If there was flood inundation, the seepage from the Yongding River would be extremely high especially along riverbanks without liner protection.
Comparison of results between MFSL and field tests, HYDRUS and two versions of the methods of stream package, MODFLOW, show that MFSL is reliable for predicting seepage under multi-layered conditions. Overall, the proposed approach of coupling GIS and MFSL is suggested to be reasonable and convenient for practical applications, such as for evaluating river-aquifer seepage loss under various possible conditions, and for the water resource management when interactions of surface water and groundwater are involved.