An improved method for the calculation of unsaturated–saturated water flow by coupling the FEM and FDM

Numerical modeling of water movement in both unsaturated soils and saturated groundwater aquifers is important for water resource management simulations. The development of efficient numerical algorithms for coupling unsaturated and saturated flow has been a long-lasting challenge in hydrologic modeling, especially for regional-scale modeling. In this study, a new method coupling the Finite Element Method (FEM) and Finite Difference Method (FDM), FE-FDM, is developed to solve Richards equation for simulating unsaturated–saturated water flow. The FEM is adopted to discretize the governing equation in the horizontal direction, while the FDM is used in the vertical direction. This method combines the advantages of FEM in domain discretization, especially for complex computational domain, and the advantages of FDM in modeling simplicity and efficiency. The validity of the new method is demonstrated with three test cases that show that the FE-FDM solutions are accurate and are applicable for regional scale problems. In the test cases, the FE-FDM solutions are compared with experimental data and numerical results calculated with common software packages including FEFLOW and COMSOL. This study verified that the FE-FDM is applicable for simulating water flow in the unsaturated–saturated zone.

Understanding the hydrologic cycle is critical for proper management of groundwater resources. In the context of surface-water and groundwater interaction, the hydrologic interactions among soil, vegetation, atmospheric processes, and groundwater dynamics should be considered, which depends to a large extent on the characteristics of the unsaturated zone 1,2 . The commonly used equation for flow in the unsaturated zone, the Richards equation 3,4 , can be troublesome to solve because it is highly nonlinear often causing convergence issues with solution schemes, especially with scale differences between horizontal and vertical dimensions in regional models 1,5 . Because of these difficulties, there are limitations to the coupling of unsaturated and saturated flow in numerical simulations 6 . Therefore, it is essential to develop a mathematical/computational method that can efficiently simulate the coupled unsaturated-saturated flow to manage groundwater resources, especially at the regional scale 2,7 .
The solutions methods for simulating coupled saturated and unsaturated water flow have been studied by many researchers. Rubin 8 developed a transient numerical model integrating the saturated and unsaturated zones, which solved Richards equation for two-dimensional, transient groundwater flow in a rectangular saturated-unsaturated soil domain. Freeze 9 further developed a three-dimensional, transient, saturated-unsaturated flow model to solve the saturated-unsaturated flow equation in the unsaturated zone and the saturated flow equation in underlying unconfined and confined aquifers. These early solution methods were prone to mass-balance errors. This issue was investigated by Mitly 10 and he suggested lumping procedures and methods of time averaging of the storage term to reduce mass-balance errors. Of course, a number of researchers simulate the saturated-unsaturated flow to solve elliptic problems that are robust for high contrasts in material properties and heterogeneity while providing locally conservative velocities for transport over the years [11][12][13][14][15][16][17] . However, variably saturated flow simulation has proven challenging 6 .
Since the 1970s, a number of commonly used groundwater flow and transport models have been developed that simulate saturated and unsaturated flow conditions. FEFLOW was first introduced in 1979 by the Institute for Water Resources Planning and Systems Research Inc. (WASY GmbH) of Berlin, Germany a part of The Danish Hydraulic Institute (DHI) group 18 . A 2D finite-element simulation model called SUTRA for saturated-unsaturated, fluid density-dependent groundwater flow with energy transport was developed by Voss 19 . Yeh and Ward 20 further developed a 3D finite-element model for variably saturated flow, 3DFEMWATER, in  1987. The DHI 21,22 developed a more comprehensive saturated-unsaturated flow model, based on the FDM, MIKE  SHE (Système Hydrologique Européen), where the unsaturated flow is simplified to 1-D in Richards equation while saturated flow is controlled by 3-D Boussinesq equation. COMSOL Multiphysics 23 is a numerical simulation software based on finite element method which builds models on the basis of general partial differential equations or partial differential equations. Richards equation is built into subsurface flow module of COMSOL to simulate unsaturated-saturated flow. FDM and FEM are recognized as the most popular numerical methods for common groundwater flow models 24 . Although FDM is simpler and more efficient, it is not flexible and efficient for representing irregular boundaries and using refined local grid spacing. FEM requires more computer memory than FDM, though FEM has good mesh adaptability for representing complex boundaries and mesh refinement.
In saturated and unsaturated flow models, the horizontal dimension is generally much larger than the vertical dimension, especially in regional models. In addition, the magnitude of the hydraulic conductivity varies greatly in the vertical direction of the unsaturated zone in the process of rainfall recharge. Thus, a much finer mesh is required in the vertical direction than in the horizontal. Using Different Methods to discrete Equations in horizontal and vertical directions respectively is one of ways to resolve this difference.
This study developed a new numerical method, FE-FDM by coupling FEM and FDM, to solve the Richards equation for simulating unsaturated-saturated water flow. The FEM is adopted to discretize the governing equation in horizontal direction, while the FDM is used in the vertical direction. This method combines the advantages of FEM in domain discretization especially for the irregular boundaries and the advantages of FDM in modeling simplicity and efficiency.

Numerical form of Richards Equation Discretized by FE-FDM
Governing equation. The Richards equation, which is commonly used to describe saturated-unsaturated groundwater flow, is expressed as 2,3 : where Kx(h), Ky(h), and Kz(h) are x-, y-, and z-directions hydraulic conductivity as a function of pressure head; h is the pressure head; w is a source/sink term; z is the elevation head; C(h) is the specific moisture capacity as a function of pressure head; Se is the water saturation; S is the specific storage; t is time; and H is hydraulic head calculated by: The specific moisture capacity C(h) and the hydraulic conductivity K(h) both are functions of the pressure head h in the governing Eq. (1), which leads equation to be nonlinear. In order to solve the nonlinear governing equation, functional relations must be obtained between two model parameters, C(h) and K(h), with unknown variable h. Van Genuchten 25 developed a model for soil water retention curve and hydraulic conductivity function based on earlier work from Mualem 26 .
The saturation-pressure relation is expressed by: where θ is the water content; θs the saturated water content; θr the residual water content; α, m and n are coefficients specifying a particular medium type, and m is restricted as m = 1 − 1/n. The hydraulic conductivity function in both saturated and unsaturated soils is written as: where Ks is the saturated hydraulic conductivity, Se is degree of saturation which could be expressed by: The specific moisture capacity can be given by: Finite element method in the horizontal direction. The governing equation for horizontal flow is discretized by the FEM. According to the Galerkin method, the discrete form in horizontal direction is constructed as follows: where N p is node number in a horizontal layer; N i , N j is basis function which is a function of x and y; k stands the kth horizontal layer; q is specified flux in the boundary. Symbols were introduced as follows: Then the Eq. (7) is expressed as: Finite difference method in the vertical direction. In this section, Eq. (12) is discretized by FDM in terms of the space (z direction) and time. A central difference method and a fully implicit difference method are developed for the z-direction and time, respectively. Equation (12) discretized in the z direction is expressed as: where n is nth time step; − K z k 1 2 is harmonic mean of hydraulic conductivity in (k-1)th layer and kth layer; + K z k 1 2 is harmonic mean of hydraulic conductivity in kth layer and (k + 1)th layer.
The method of full implicit discretization is adopted in time discretization. Then the expression is obtained as: According to Eq. (14), the final global equation of the system is formed assembling all the horizontal layers.
= DH F (15) where D is global matrix consisting of the coefficients in the Eq. (14); H is the unknown hydraulic head values in the nodes for the current time step and current iteration level; F is the known values in the nodes from the last time step and last iteration level. Equation (15) is a non-linear equation which is solved by Picard iteration.

Verification Examples
In this study, three test problems are used to verify the accuracy and reliability of the method. The results calculated by FE-FDM are compared with popular models and empirical data.

Test problem 1: 2D unsaturated-saturated water flow.
In this test problem, the study domain is a 2D rectangular unconfined aquifer between two rivers 40 m apart with a 3 m thick aquifer. The soil parameters are presented in the Table 1. The conceptual model of test problem 1 is shown in Fig. 1. The initial hydraulic head is The pressure heads calcualted by FE-FDM and FEFLOW at t = 5d, 15d, 25d and 35d are shown on Fig. 3. FE-FEDM results are in good agreement with FEFLOW solutions at the different times. At early time (i.e. t = 5d), the pressure heads are similar to initial conditions. The smaller the contour of the pressure head, such as −0.7 m, has the greater the change, since the moisture content mainly in the unsaturated zone. At t = 35d, all the pressure heads have been raised obviously. Figure 4 shows this process more clearly. In Fig. 4, pressure head change process is shown at the location of (20,2) m which is at the center of water table in the initial time. As shown in the Fig. 4, the pressure head begins to grow after about 0.1d. There is a lag time when groundwater is recharged, since infiltrated water should pass through the vadose zone. The lag time depends on the characteristics of the vadose zone.
Test problem 2: Sand tank experiment. In this example, FE-FDM and FEFLOW are used to simulate a transient sand tank experiment which was conducted by Koichi et al. 27 . The sand tank is 315 cm in length, 23 cm in width and 33 cm in height. A two-dimensional vertical profile model was used to simulate the sand tank experiment. The soil parameters for the tank experiments, which are the same as those used in the model, are presented in Table 2. Figure 5 shows the conceptual model for test problem 2. The hydraulic head on the left boundary is increased from 10 cm to 30 cm at time equal zero, and the hydraulic head on the right boundary is 10 cm, both of which are constant head boundaries. The top and bottom boundaries are impermeable. The initial condition is a water table at 10 cm. The same rectangular mesh is used for the FEFLOW and FE-FDM models. The domain is discretized with 41 rows and 31 columns of nodes. The transient simulation was run for 4,800 seconds. Figure 6 shows the water table evaluations measured in the experiment and the water table elevations Table 1. Soil parameters of test problem 1.  www.nature.com/scientificreports www.nature.com/scientificreports/ FE-FDM are closer to the experimental data than those calculated with FEFLOW. As shown in Fig. 6, the water tables gradually rise during the experiment. Figure 7 shows the lateral velocity and the vertical velocity distributions in the unsaturated zone at z = 20 cm segment. As the Fig. 7 shows, both the lateral velocity and the vertical velocity decrease with distance from the infiltration line. The lateral velocity is much larger than the vertical velocity near the infiltration line. Thus the horizontal movement of water flow has an important role in the unsaturated zone where lateral flow is evident. Some numerical models ignore the horizontal velocity, because the horizontal hydraulic gradient is usually smaller   www.nature.com/scientificreports www.nature.com/scientificreports/ than the vertical gradient 24,28,29 . This case, demonstrates why a fully 3D coupling of the saturated-unsaturated model is necessary. The initial hydraulic head is 20 m in the model. Although this model is highly simplified, the model can be used to verify whether the FE-FDM can be applied to large-scale sites.
In this model, the scale in the horizontal direction is much larger than that in the vertical direction. Because hydraulic conductivity varies greatly in vertical direction in the unsaturated zone, a fine mesh spacing of 2 m is specified in the vertical direction. Course meshing is specified in the horizontal direction which has 256 triangle elements in every layer. The average side length of meshes is 954 m in the horizontal direction. The average side length of the meshes in the horizontal direction is 477 times the length in the vertical direction in this model. Therefore, the horizontal dimension is much larger than the vertical direction not only for the whole model but also for a single mesh.
FE-FDM and COMSOL are used to solve this test problem. COMSOL is a multiphysical field numerical simulation software which is a stable program that covers a wide range of applications for those interested in ground water modeling based on FEM 30 . The same mesh is used for both FE-FDM and COMSOL. Figure 8 shows the water table elevations in the central section at 1d, 3d and 5d. There are differences between the values calculated by FE-FDM and COMSOL. The differences mainly are lower by 0.035 m and their ratio to COMSOL solution are mainly less than 0.2%. The changes of the water table elevations with time as calculated by FE-FDM and COMSOL are consistent. The water table elevations fluctuate near the boundary of the model as the mesh is not fine enough to accurately calculate the water table position near the boundary. It is obvious that the results of FE-FDM are less affected by the boundary than that of COMSOL. The FE-FDM has better adaptability to coarse mesh than COMSOL. This proves that FE-FDM can be applied to large scale models.