Yan ZHU, Yuan-yuan ZHA, Ju-xiu TONG, Jin-zhong YANG*

State Key Laboratory of Water Resources and Hydropower Engineering Science,Wuhan University, Wuhan 430072, P. R. China

1 Introduction

Understanding the interaction between soil water and groundwater dynamics is of great importance in water resources planning and management in an extensively distributed unsaturated-saturated zone. Many models have been developed to simulate the complex processes of water flux regulation between unsaturated and saturated zones. To estimate the recharge flux from the unsaturated zone, some groundwater models consider the recharge flux to be the model input, while groundwater table dynamics are generally considered in a very simplified form, or not considered at all (Carrera and Medina 1999). However, water movement in the unsaturated zone is sensitive to climate change, soil utilization, and human activities. Thus, it is necessary to calculate the hydraulic change in the unsaturated zone and estimate the transient recharge flux to the saturated zone.

Nowadays, there are two kinds of coupled unsaturated-saturated water flow models commonly used in simulation of subsurface systems. One type is fully three-dimensional (3-D)unsaturated-saturated water flow models, such as FEFLOW (Diersch and Kolditz 1998), which uses a 3-D equation to represent the water flow in the subsurface system. However, it is a little difficult to apply a fully 3-D model to a large area with an unsaturated zone because it spends too much time disposing the hydraulic parameters and solving the 3-D unsaturated water flow equation due to its high nonlinearity. On the basin scale, the direction of unsaturated flow averaged over large grid cells is usually vertical (Mantoglou 1992; Chen et al. 1994), which makes it possible to simplify the 3-D soil water flow model into a one-dimensional (1-D)vertical flow model and couple it with a fully 3-D groundwater model to develop a more efficient unsaturated-saturated water flow model. The earlier coupled models were developed on the field scale. The model developed by Pikul et al. (1974)is a typical one. This model combined Richards’ 1-D equation for the unsaturated zone with the two-dimensional (2-D)Boussinesq’s equation for the saturated zone, and used the storage coefficient as the linkage between the two systems. But subsequent research shows that it is very difficult to determine the storage coefficient and not convenient to use this coupled model (Vachaud and Vauclin 1975). Another problem is that the earlier coupled model used the 2-D or quasi 3-D sub-model to simulate the groundwater flow, which cannot be done for a complex groundwater system.Considering these factors, most of the coupled unsaturated-saturated water flow models recently developed adopts 3-D groundwater models. MODFLOW is one of the most popular 3-D groundwater models, with high efficiency and accuracy (McDonald and Harbaugh 1988).Researchers have developed different soil water packages for simplifying soil water simulation and combined them with MODFLOW to form coupled models for the unsaturated-saturated zone (Havard et al. 1995; Facchi et al. 2004; Niswonger et al. 2006; Lin et al. 2010). These models usually have some constraints in their application. For example, the UZF1 package established by Niswonger et al. (2006)was restricted to a homogeneous medium, and the simplified coupled model developed by Lin et al. (2010)could only be applied to an aquifer with a gentle head gradient where the groundwater flow could be simplified into 2-D flow. On the other hand, although MODFLOW is widely used all over the world, some limitations exist.For instance, significant deformation of meshes may lead to inaccurate simulation results or even induce some errors (McDonald and Harbaugh 1988).

The purpose of our work is to establish an efficient coupled unsaturated-saturated water flow model that can be applied to large irregular areas with complex boundary conditions. A 1-D flow module for simulating the main infiltration processes in the unsaturated zone was developed and coupled with a simplified 3-D groundwater flow model developed by Zhu et al.(2010)and could be applied to an irregular aquifer system. The source code for the flow package in the unsaturated zone was extracted and modified from Hydrus 5.0 (Vogel et al.1996). For the unsaturated-saturated system, most attention is paid to the recharge flux occurring at the phreatic surface. The recharge flux is seen as the connection of the two separate systems, which is expressed by the hydraulic head difference between the adjacent nodes in the unsaturated and saturated zones. To couple the two separate systems, the source codes of the two sub-models are modified according to the expression of the recharge flux.When the program begins to run, the hydraulic head of the nodes in the saturated zone is assumed first, and the depth of the 1-D soil column and the groundwater table are determined.The recharge flux is calculated and the soil water flow equation for each 1-D soil column is solved by taking the recharge flux as the bottom boundary condition. Similarly, the nodal hydraulic heads in the groundwater system are calculated by considering the recharge flux as the upper boundary condition. Iteration continues until the discrepancy between the calculated and assumed hydraulic heads is less than the accuracy control value. It should be noted that the soil water model and groundwater model were coupled as a united program. Some special sub-modules were established to simulate other important flow processes, e.g., crop cover,rainfall/irrigation, and pumping wells, in subsurface systems.

In this paper, the details of developing the model are described. Case studies of model application to different conditions of water flow, as well as the simulated results, are presented and discussed. In addition, the benefits and limitations of the coupled model in practical application are also demonstrated.

2 Development of coupled model system

2.1 Overview of Hydrus 5.0 and groundwater models

Hydrus 5.0 (Vogel et al. 1996)is a numerical model for simulating water, heat, and solute movements in 1-D variably saturated media. We chose the module for water flow in Hydrus 5.0 and modified it to couple it with the groundwater model. It is assumed that soil water flows in the vertical direction and the model can deal with the prescribed head boundary, flux boundary, and atmospheric boundary conditions. The governing flow equation is solved numerically using Galerkin-type linear finite element schemes. As a 1-D unsaturated-saturated model, in most conditions, Hydrus 5.0 is ineffective in dealing with groundwater flow, because groundwater flow is in three directions and mainly in the lateral direction. To simulate groundwater flow, a 3-D groundwater model has been established based on the results of water balance analysis in the finite volume, which is suitable for groundwater simulation in large areas, especially in irregular aquifer-aquitard systems (Zhu et al. 2010). It has been validated as being more efficient than the fully 3-D finite element model FEFLOW and can obtain higher simulation accuracy than the finite difference model Visual MODFLOW in irregular aquifers.Although the model can be applied to simulation of groundwater flow under different boundary conditions, such as the specified head boundary condition, flux boundary condition, and sink terms, it has difficulty computing the distributed groundwater recharge. Therefore, if the recharge flux obtained based on the hydraulic heads of the saturated and unsaturated zones is taken as the input of the groundwater model, the spatial-temporal variation of hydraulic characteristics in the unsaturated and saturated zones will be properly depicted.

A schematic diagram of the coupled system is shown in Fig. 1. The unsaturated zone is divided into different sub-areas in the horizontal direction according to the soil material characteristics, atmospheric boundary conditions, land use conditions, and crop types, each of which is represented by a 1-D soil column. The 1-D soil water equation is discretized by the finite element method, described in Vogel et al. (1996). The saturated zone is divided into different layers in the vertical direction, and the total matrix for the groundwater equation is developed by the method presented by Zhu et al. (2010). Each sub-area has its own recharge flux to the saturated zone, which has no relationship with other sub-areas.

Fig. 1 Schematic diagram of coupled soil water and groundwater model

2.2 Structure of integrated unsaturated-saturated water flow model

The coupled model begins with reading the input data which are required to initiate the soil water sub-model and the groundwater sub-model. First of all, water flow in the unsaturated zone is calculated with the assumed groundwater table. Then, the water flow in the saturated zone is calculated based on the results of the recharge flux from the unsaturated zone. Some modules are programmed to link the two water flow systems. The iteration continues if differences between the assumed and simulated hydraulic heads in the saturated zone are beyond the accuracy control value. The simulation processes of the coupled model are shown in Fig. 2, where t is the simulation time, Δt is the time step, and tmaxis the maximum simulation time. The average unsaturated water flow in a sub-area is represented by a 1-D soil column. When there are more than one 1-D soil columns, we can form and solve the matrices for 1-D soil water columns one by one.

Fig. 2 Flow chart of computation of coupled model

2.3 Vertical flow in unsaturated zone

Modules related with the soil water flow processes in Hydrus 5.0 (Vogel et al. 1996)are extracted and utilized as the unsaturated zone flow package in this study. This soil water flow package can be applied to different 1-D soil columns with variable height under different conditions, extending from the soil surface to the groundwater table. Some changes to the items in the matrix for the nodes adjacent to the groundwater table should be made, so that the soil water flow package can be coupled with the 3-D groundwater sub-model. The details for the changes are explained in the subsequent part of this section.

The equation describing soil water flow in the unsaturated zone can be written as

where θ is the volumetric water content, K is the unsaturated hydraulic conductivity as a function of the pressure head h, z is the coordinate in the vertical direction, and S is the sink or source term.

A 1-D soil column is first discretized into n-1 adjoining elements, with the last element located at the soil surface. There are n nodes for the 1-D unsaturated column. It is assumed that the vertical coordinate axis is directed downward. We give immediately the following matrix equation for the discretized 1-D soil water equation (Vogel et al. 1996):

where A is the coefficient matrix in the global matrix equation for water flow, B is the vector in the global matrix equation for water flow, and h is the vector in the global matrix equation of pressure head.

Since it is assumed that there are n nodes for one 1-D soil column, the coefficient matrix is an n×n dimensional matrix, and Eq. (2)has the following form:

Eq. (4)represents water balance at node i, which is adjacent to the groundwater table in a 1-D soil column as shown in Fig. 3:

where zNis the elevation of node N in the phreatic surface adjacent to node i, zi+12is the upper surface elevation of the control element of node i, qi+12andare the vertical flow fluxes through the upper and lower surfaces of the control element, and θiis the volumetric water content of node i.

Two items, (zN- zi)ΔθiΔt andqN, can be given as follows:

Fig. 3 Water balance at node i in unsaturated zone

where Ciis the soil water capacity at node i;andare the pressure head of node i at the (k+1)th and kth time steps, respectively;is the hydraulic head of node N at the(k+1)th time step; andis the mean unsaturated hydraulic conductivity of node i and node N.

Substituting Eq. (5)and Eq. (6)into Eq. (4), we have

The matrixes A and B can be modified according to Eq. (7), and the final matrix is obtained in Eq. (8):

It should be noted that the number of the nodes in the 1-D soil column is fixed, but not all the nodes are included in the calculation. Only those with their elevations over the groundwater table can be considered active nodes, and those under the groundwater table are taken as inactive nodes and not used for calculation. As calculation runs to the next time step,if the groundwater table decreases, some of the nodes under the groundwater table in the previous time step will be over the groundwater table, and they will be calculated in the present time step. In contrast, when the groundwater table rises, some of the nodes over the groundwater table in the previous time step will be under the groundwater table, and they will not be calculated in the present time step. In a word, the number of nodes in the calculation in the 1-D soil column changes with the groundwater table.

2.4 Water flow in groundwater system

The groundwater model developed by Zhu et al. (2010)was used in this study. It is a simplified 3-D groundwater model based on the theory of the quasi 3-D model. The model runs according to the following steps: (1)the aquifer-aquitard system is stratified into some layers in the vertical direction and some triangular prism elements (Fig. 4)in the horizontal direction within each horizontal layer based on the aquifer extension and distribution; (2)an average head gradient surface of the triangular prism element is determined; (3)the element stiffness matrix is established from the results of water balance calculation in the control volume based on the average head gradient surface; (4)the flow flux between upper and lower layers is computed with the finite difference method; and (5)the global stiffness matrix is formed by combining the results of the finite element analysis in each layer and the finite difference analysis between the layers. This model’s suitability for groundwater modeling of large aquifer-aquitard systems has been validated by comparison of its results with those of other common software, such as FEFLOW and Visual MODFLOW. The matrix formation of the nodes in the phreatic surface needs to be modified when it is combined with the equation for the unsaturated zone. Details about coupling the groundwater model with the unsaturated zone equation are provided below:

In the groundwater model, the water balance equations for the nodes under the water table are different from those in the phreatic surface. There is no upper flow flux to the nodes in the phreatic surface. They are analyzed respectively.

The control volume for any node under the groundwater table is shown in Fig. 5.

Fig. 4 Triangular prism element ijki′ j′ k′

Fig. 5 Control volume of node M

The water balance equation for node M is expressed as

where μMis the elastic storage coefficient at node M; VMis the control volume of node M;QMis the net lateral flux into the control volume of node M; qM+12and qM-12are the vertical flow fluxes of the upper and the bottom surfaces, respectively; and SMis the surface area of the control volume.

Eq. (11)represents the water balance equation of the nodes in the 3-D groundwater system except for the nodes in the phreatic surface, for which there is no upper flow flux. The recharge flux to the groundwater system is considered a source term in the groundwater model.

For the groundwater system, if there are three layers and the number of nodes in each layer is n, the matrix of the saturated model has the following form:

However, when the unsaturated zone and saturated zone are coupled as a united system, the recharge flux should be considered the upper flux for the groundwater system. The nodes in the phreatic surface have a direct connection with nodes in the unsaturated zone. The water balance analysis for nodes in the phreatic surface can help to form the coupling matrix for the groundwater system.

The node N illustrated in Fig. 6 is taken as an example for water balance analysis in the phreatic surface. The water balance equation can be written as follows:

Fig. 6 [0]Control volume of node N

The recharge flux from the unsaturated zone qNcan be expressed as Eq. (6). Then, the water balance equation of the node N in the groundwater table is determined by Eqs. (6)and (13):

The modified matrix can be expressed as follows:

3 Simulation performance

Different examples were used to test the reliability of the coupled model by comparing the obtained results with those from common software, such as Hydrus-1D, SWMS-2D, and FEFLOW. Most boundary conditions and source or sink items listed in Fig. 1 were tested in the case studies.

3.1 Column infiltration test

The soil profile was 5 m in depth, and the soil was assumed to be homogenous and isotropic with a saturated hydraulic conductivity of 0.6 m/d. The residual water content θrwas 0.057, and the saturated water content θswas 0.35. The coefficient in the soil water retention function α was 4.1 m-1, and the exponent in the soil water retention function β was 2.28. It was assumed that the initial pressure head was in a steady state, and the initial phreatic surface was at a depth of 3.3 m. The infiltration rate was 0.05 m/d at the soil surface, and there was no flow around it, a situation which could be considered a 1-D vertical water flow. The simulation was carried out for 10 d, as shown in Fig. 7 and Fig. 8.

Figs. 7 and 8 show the pressure head and the volumetric water content of the soil profile at different simulation times calculated by the coupled model and Hydrus-1D. The results produced by the coupled model agreed quite well with those obtained by Hydrus-1D.Therefore, the coupled model is effective.

Fig. 7 Soil profile pressure head comparison of coupled model and Hydrus-1D at different simulation time

Fig. 8 Soil profile volumetric water content comparison of coupled model and Hydrus-1D at different simulation time

3.2 Water flow with crop growth

This example considered 1-D water flow in a grassland profile, and the precipitation,evaporation, and transpiration were also considered in the modeling. The simulation results obtained with the coupled model were compared with those from Hydrus-1D. The soil profile was 3 m in depth and the depth of the root zone was within 0.3 m from the ground surface.The soil water parameters were the same as those in Section 3.1. The surface fluxes were calculated by daily precipitation rates (shown in Fig. 9)and the potential transpiration rate(assumed to be 0.002 m/d). A no-flow boundary condition was assumed at the bottom boundary of the soil column. The initial pressure head was from -1.3 m at the ground surface to 1.7 m at the bottom. The pressure head and volumetric water content of the soil profile obtained by the coupled model and Hydrus-1D are shown in Figs. 10 and 11.

Fig. 9 Daily precipitation rate curve

Fig. 10 Soil profile pressure head comparison of coupled model and Hydrus-1D at different simulation time

Fig. 11 Soil profile volumetric water content comparison of coupled model and Hydrus-1D at different simulation time

The pressure head in the root zone changed dramatically because it was influenced by the upper boundary condition in this example. The coupled model was sensitive to the changing boundary condition as shown in Fig. 10 and Fig. 11 and could obtain reasonable results.

3.3 2-D water flow

The coupled model was used to simulate the water flow within two rivers while there was precipitation infiltration from the upper boundary. The distance between the rivers was 40 m.The infiltration rate was 0.004 m/d. The soil profile was 3 m, and it was assumed to be homogenous with a saturated hydraulic conductivity of 0.5 m/d. The residual water content θrwas 0.02, and the saturated water content θswas 0.3. The coefficient in the soil water retention function α was 4.1 m-1, and the exponent in the soil water retention function β was 1.964. The flow system was plotted in Fig. 12. The analytical solution of the hydraulic head when the water flow was steady is given in Eq. (18):

where H is the hydraulic head at the location with a distance x from the left river; H1andare the initial hydraulic heads in the left and right rivers, respectively; l is the distance between the two rivers; W is the precipitation rate; and K′ is the saturated hydraulic conductivity.

Fig. 12 Profile of flow system

Because the flow flux from the unsaturated zone to the saturated zone equals the precipitation rate when the water flow is steady, we can use one 1-D soil column to represent the whole unsaturated zone. For the 2-D unsaturated-saturated water flow model, it is difficult to output the recharge rate between the unsaturated and saturated zones. We just compared the recharge rate obtained from the coupled model with that from SWMS-2D in the steady state(Fig. 13). The groundwater table calculated from the coupled model was compared with the analytical solution and that from SWMS-2D (Fig. 14). The comparison of the pressure head in the soil profile obtained from the coupled model and SWMS-2D at different locations is shown in Fig. 15.

Fig. 13 Comparison of recharge rate for groundwater system in steady state obtained from coupled model and SWMS-2D

Fig. 14 Comparison of groundwater table calculated from coupled model with analytical solution and that from SWMS-2D

Fig. 15 Comparisons of pressure heads calculated from coupled model and SWMS-2D at different distances from left river

In this example, the groundwater table increased from 2.0 m to 2.68 m with the highest value occurring at the middle location. Flow exchange between soil water and groundwater was complex and drastic, changing from 0.0 m/d at the start to 0.004 m/d in the end. As shown in Figs. 14 and 15, good results were obtained by the coupled model as compared with those from SWMS-2D. The largest absolute error and relative error for the hydraulic head in the simulation were 0.02 cm and 0.65%, respectively. The simulation results were reasonable.However, it was shown that the coupled model had no obvious advantage over the SWMS-2D model considering the simulation cost.

3.4 3-D water flow

The simulation domain has a size of 5 000 m × 5 000 m × 30 m. The groundwater depth was 10 m from the soil surface. There were two precipitation sub-areas divided by the midline. The precipitation rates of the two sub-areas were 0.001 m/d (sub-area #1)and 0.001 5 m/d (sub-area#2). A no-flow boundary condition was assumed around the simulation zone. The soil was assumed to be homogenous, and the saturated hydraulic conductivity was 1.08 m/d. The residual water content θrwas 0.065, and the saturated water content θswas 0.41. The coefficient in the soil water retention function α was 7.5 m-1, and the exponent in the soil water retention function β was 1.89. The simulation was carried out for 200 d. The unsaturated zone and saturated zones were divided into 50 layers and eight layers, respectively. We adopted 2 000 nodes in each layer, so there were 118 000 nodes in FEFLOW and 18 102 nodes in the coupled model.

The hydraulic head change processes at different depths obtained from the coupled model and FEFLOW were shown in Fig. 16. This indicates that the coupled model almost had the same simulation accuracy as the 3-D model. However, it reduced the CPU run time significantly. The simulation time was 0.33 h for the coupled model and 4.5 h for FEFLOW.

Fig. 16 Comparison of hydraulic heads at different depths obtained from coupled model and FEFLOW

4 Comparison and discussion

Resultsfrom the coupled model were compared with those from the Hydrus-1D,SWMS-2D, and FEFLOW. Case studies demonstrate the capability of the coupled model in simulating water flow in the unsaturated and saturated zones on a large scale. The simulation results are in good agreement with those from Hydrus-1D, SWMS-2D, and FEFLOW. As shown in case 1 and case 2, the coupled model runs well when linking water flow between the unsaturated and saturated zones and determines the groundwater table accurately. Besides that,the model can also capture the large head change in the root zone. As we know, it is difficult to obtain the recharge flux from the unsaturated zone to the saturated zone when using the fully 1-D or 2-D unsaturated-saturated water flow model such as Hydrus-1D and SWMS-2D.However, in our coupled system, the unsaturated zone and saturated zone are strictly separated,and it is easier to locate the unsaturated and saturated zones accurately and calculate the recharge flux for the groundwater system. It is of importance for groundwater estimation on a large scale.

Although the fully 3-D unsaturated-saturated water flow models such as FEFLOW can be used to model large areas, the computational cost is too expensive in some cases. Although the thickness of the unsaturated zone is much less than that of the groundwater zone, most of the computational cost is spent on the unsaturated zone. Simplification of the unsaturated zone is necessary to decrease the computational cost. According to the results of case studies, the technology used to couple soil water with groundwater as described in this paper is feasible and can satisfy the accuracy requirement. The more important aspect is that it reduces the computational cost significantly, especially for a simulation domain with a large number of nodes.

5 Conclusions

A coupled unsaturated-saturated water flow model was established using the 1-D vertical Richards’ equation to represent the soil water flow and a 3-D numerical model to represent the groundwater flow. The model was applied to case studies of water flow simulation with different boundary and water flow conditions in the unsaturated-saturated zone. The simulation results were compared with those of other models. The main conclusions are as follows:

(1)The exchange flux between the unsaturated zone and saturated zones is always involved in the evaluation of the coupling technique. Case studies demonstrate the capability of the coupling method in calculating the recharge rate. The simulation results are in good agreement with those from Hydrus-1D, SWMS-2D, and FEFLOW.

(2)Compared with other fully 1-D, 2-D, or 3-D saturated-unsaturated water flow models(e.g., Hydrus-1D, SWMS-2D, and FEFLOW), it is much easier for the coupled model to obtain the recharge rate by strictly separating the unsaturated and saturated zones.

(3)It is reasonable to simplify soil water flow into 1-D vertical flow to reflect the hydraulic characteristics of the unsaturated zone in a large area. The technology is useful for decreasing the computational cost for the unsaturated-saturated water flow model.

(4)The coupled model cannot be used for the condition when the lateral flow plays the same important role as the vertical flow or is even greater than the vertical flow in the unsaturated zone.

(5)As the flow flux between soil water and groundwater is calculated by the head difference of the two systems, which is disposed of with an explicit scheme, the solution may be divergent when the water table changes frequently and strongly.

Carrera, J., and Medina, A. 1999. A discussion on the calibration of regional groundwater models. Feyen, J.,and Wiyo, K., eds., Proceedings of International Workshop of EurAgEng’s Field of Interest on Soil and Water, Modelling of Transport Processes in Soils at Various Scales in Time and Space, 629-640.Wageningen: Wageningen Pers.

Chen, Z. Q., Govindaraju, R. S., and Kavvas, M. L. 1994. Spatial averaging of unsaturated flow equations under infiltration conditions over areally heterogeneous fields: 1. Development of models. Water Resources Research, 30(2), 523-534. [doi:10.1029/93WR02885]

Diersch, H. J. G., and Kolditz, O. 1998. Coupled groundwater flow and transport: 2. Thermohaline and 3D convection systems. Advances in Water Resources, 21(5), 401-425. [doi:10.1016/S0309-1708(97)00003-1]

Facchi, A., Ortuani, B., Maggi, D., and Gandolfi, C. 2004. Coupled SVAT: Groundwater model for water resources simulation in irrigated alluvial plains. Environmental Modelling and Software, 19(11),1053-1063. [doi:10.1016/j.envsoft.2003.11.008]

Havard, P. L., Prasher, S. O., Bonnell, R. B., and Madani, A. 1995. LINKFLOW, a water flow computer model for water table management, Part I: Model development. Transactions of the ASABE, 38(2),481-488.

Lin, L., Yang, J. Z., Zhang, B., and Zhu, Y. 2010. A simplified numerical model of 3D groundwater and solute transport at large scale area. Journal of Hydrodynamics, 22(3), 319-328. [doi:10.1016/S1001-6058(09)60061-5]

Mantoglou, A. 1992. A theoretical approach for modeling unsaturated flow in spatially variable soils: Effective flow models in finite domains and nonstationarity. Water Resources Research, 28(1), 251-267. [doi:10.1029/91WR02232]

McDonald, M. G., and Harbaugh, A. W. 1988. A modular three-dimensional finite-difference groundwater flow model. Techniques of Water-resources Investigations. Denver: U.S. Geological Survey.

Niswonger, R. G., Prudic, D. E., and Regan, R. S. 2006. Documentation of the unsaturated-zone flow (UZF1)package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005. Techniques and Methods, 6-A19. Denver: U.S. Geological Survey.

Pikul, M. F., Street, R. L., and Remson, I. 1974. A numerical model based on coupled one-dimensional Richards and Boussinesq equations. Water Resources Research, 10(2), 295-302. [doi:10.1029/WR010i002p00295]

Vachaud, G., and Vauclin, M. 1975. Comments on ‘A numerical model based on coupled one dimensional Richards and Boussinesq equations’ by Mary F. Pikul, Robert L. Street, and Irwin Remson. Water Resources Research, 11(3), 506-509. [doi:10.1029/WR011i003p00510]

Vogel, T., Huang, K., Zhang, R., and van Genuchten, M. Th. 1996. The Hydrus Code for Simulating One-Dimensional Water Flow, Solute Transport, and Heat Movement in Variably-satruated Media.Riverside: U.S. Salinity Laboratory, Agriculture Research Service, U.S. Department of Agriculture.

Zhu, Y., Yang, J. Z., and Tong, J. X. 2010. A simplified 3D numerical model for groundwater flow in declining aquifer-aquitard system. Journal of Sichuan University (Engineering Science Edition), 42(6), 43-50. (in Chinese)