Regional Quasi-Three-Dimensional Unsaturated-Saturated Water Flow Model Based on a Vertical-Horizontal Splitting Concept

Yan Zhu 1, Liangsheng Shi 1, Jingwei Wu 1, Ming Ye 2, Lihong Cui 1 and Jinzhong Yang 1,* 1 State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China; 2011102080004@whu.edu.cn (Y.Z.); liangshs@whu.edu.cn (L.S.); jingwei.wu@whu.edu.cn (J.W.); 2015102060057@whu.edu.cn (L.C.) 2 Department of Scientific Computing, Florida State University, Tallahassee, FL 32306, USA; mye@fsu.edu * Correspondence: jzyang@whu.edu.cn; Tel.: +86-27-6877-5432; Fax: +86-27-6877-6001


Introduction
Accurate assessment of water flow in the unsaturated-saturated porous media is critical for groundwater exploitation, groundwater resources management, crop uptake estimation, and pollution remediation [1][2][3][4].Numerical models are effective tools to handle complicated sub-surface water flow problems of hydrologic systems, with complex and dynamic boundary conditions and heterogeneous characteristics.The complexity of subsurface systems entails not only models to reliably represent the system of interest, but also numerical solutions to efficiently solve the models for recurrent hydrological issues.A number of 3-D numerical models have been developed based on the basic idea of solving the 3-D unsaturated-saturated Richards' equation, which is used to describe the water flow in the subsurface system.While using the 3-D model is considered to be the most rigorous approach for solving real hydrological problems [5], transient 3-D modeling is computationally expensive, especially at the regional scale, because a large amount of computational nodes and time steps are necessary to simulate the systems with sufficient accuracy [6,7].Moreover, the numerical schemes used to solve the 3-D Richards' equation still face a number of computational challenges, which is another hamper to the practical applications of the 3-D unsaturated-saturated model [8].Specifically speaking, for the most popular schemes of the finite difference method (FDM) and the finite element method (FEM) adopted by MODFLOW [9] and FEFLOW [10], respectively, FDM is not flexible to handle complex domain boundaries and may yield inaccurate heads at the vicinity of the boundaries [11,12], and FEM does not maintain mass balance at the element level [13,14].Therefore, it is necessary to simplify the 3-D models conceptually and numerically.
Considering the complexity and nonlinear intrinsic characteristics of subsurface flow systems, a variety of concepts led to the development of simplified flow models.Quasi 3-D models such as the loosely or fully coupled unsaturated-saturated models have been developed.Differing from the fully 3-D unsaturated-saturated models, which consider the water flow in three directions, the quasi-3-D models only consider the 1-D vertical infiltration in the unsaturated zone, but consider the groundwater system as 3-D flow [5,15,16].A common feature of the quasi 3-D models is to loosely couple a 1-D vadose zone model with 3-D groundwater flow model such as the MODFLOW-type models [9].For example, the 1-D models of SVAT, UZF1, and Hydrus1D have been integrated with MODFLOW for the unsaturated-saturated water flow modeling [7,8,17].There are other quasi 3-D models adopting the fully coupled method to integrate the 1-D vertical unsaturated flow with 3-D water flow [3,18].The advantages of such quasi 3-D models based on the assumption of 1-D vertical infiltration in the soil zone are that the computational cost is reduced significantly and that the convergence speed is increased, because the solution of the non-linear 3-D Richards' equation is avoided.However, these models can only be applied to the specified conditions when the lateral transfer flow can be neglected.Otherwise, their simulations are inaccurate.
The vertical/horizontal splitting (VHS) concept is to fill the gap between computational cost and simulation accuracy.The key to the VHS concept is the decomposition of the 3-D Richards' equation into a 1-D equation for vertical flow and a 2-D equation for depth-averaged flow.The decomposition makes it possible to obtain a faster and more flexible solution.The VHS concept was first used by Lardner and Cekirge [19] to solve a 3-D equation of tidal and storm surge.They adopted a standard 2-D algorithm for the shallow-water equations to compute the depth-averaged velocity and then to compute the vertical velocity profile.This VHS-based modeling was found to be more efficient than the 3-D modeling.Then, Tsai et al. [20] developed a layered and regional land subsidence model based on the VHS concept.Similar models have been developed to simulate the sediment transport in the coastal region by coupling a 2-D model with an additional 1-D vertical profile model [21][22][23].Lin et al. [24] were the first to use the VHS concept to simplify the 3-D saturated water flow equation to a 2-D horizontal water flow equation and to integrate the horizontal flow with the vertical flux.However, the model was only suitable to a horizontal aquifer with small slopes and was invalid when the aquifer was inclined.Paulus et al. [25] proposed an innovative 3-D unsaturated flow model by decoupling the 3-D Richards' equation into a 1-D vertical equation and a 2-D depth-integrated horizontal equation.The model was demonstrated to be efficient.However, this model cannot be used for heterogeneous soils and the FDM used in the model hampers its application to irregular domains.More studies are needed to make a full use of the VHS concept for its high efficiency to solve a 2-D horizontal water flow equation and a 1-D vertical water flow equation.This is the case in comparison with the fully 3-D equation, especially the solution of the complete 3-D equation on large-scale basins, which leads to the very huge non-linear systems, whose numerical solution is very costly.Another advantage of the VHS-based modeling is that, after splitting the 3-D Richards' equation to horizontal and vertical equations, the equations can be solved by using different numerical schemes suitable to the equations.
The objective of this paper is to present a VHS-based 3-D unsaturated-saturated water flow model to split the 3-D unsaturated-saturated water flow to the 2-D horizontal flow described by a 2-D equation and the 1-D flow described by a 1-D equation.Different optional methods are used for the two equations.The 2-D horizontal equation represents the lateral flow in the horizontal plane of average head gradient, which is calculated by the water balance method.In addition, the 1-D vertical equation is discretized by the FDM to integrate the different horizontal layers.The two equations are solved simultaneously by coupling them into a unified nonlinear system with a single matrix.The efficiency and accuracy of the developed model and its computer code are evaluated by using three synthetic cases represented the 1-D infiltration flow, 2-D lateral flow, and 3-D well flow.Furthermore, the model applicability and efficiency are examined by simulating groundwater flow under complicated boundary conditions in a regional-scale case.Furthermore, the model applicability and efficiency are examined by simulating groundwater flow under complicated boundary conditions in a regional-scale case.

Discretization of the Aquifer
The model is established firstly by separating the 3-D aquifer into different layers.Each layer is divided into triangular prism elements.The horizontal plane of average head gradient is derived to represent the average horizontal flow in the triangular prism element.The aquifer domain, the divided triangular prism element and the horizontal plane of average head gradient are shown in

The Governing Equations
The general 3-D unsaturated-saturated water flow is represented as follows: where θ is the volumetric water content [-]; t is the time [T]; K is the hydraulic conductivity as a function of pressure head [L T −1 where C(h) is the soil water capacity [L −1 ] and can be calculated as

The Governing Equations
The general 3-D unsaturated-saturated water flow is represented as follows: where θ is the volumetric water content [-]; t is the time [T]; K is the hydraulic conductivity as a function of pressure head [L T ´1]; H is the hydraulic head [L]; x i , x j (i, j = x, y, z) are the coordinates in the x, y and z directions [L]; S is the sink or source term [T ´1].
Considering the first left item of Equation (1) as where C(h) is the soil water capacity [L ´1] and can be calculated as where µ is the elastic storage coefficient [L ´1]; n is the soil porosity [-]; h is the pressure head [L].
In the saturated zone, dθ/dh can be ignored, and θ/n equals 1.0.The hydraulic conductivity is the saturated hydraulic conductivity K s .Therefore, for the saturated water flow, the water flow Equation (1) can be simplified as In the unsaturated zone, θ/n ˆµ can be ignored.Then, C(h) equals to dθ/dh, and the hydraulic conductivity is a function of pressure head.These two items can be calculated according to the modified van Genuchten model [26], shown as Equations ( 6) and (7).Then, in the unsaturated zone, the water flow Equation (1) can be rewritten as Therefore, the saturated water flow Equation ( 4) and the unsaturated water flow Equation ( 5) have the same equation form but with their specified coefficients as

Depth-Averaged Water Flow Equation
Considering one horizontal layer, the depth-average 2-D equation can be obtained by integrating the 3-D flow equation over the vertical direction, which yields where z i (x,y) and z i+1 (x,y) are the bottom and upper boundary, respectively, of the aquifer.By assuming B = z i+1 (x,y)-z i (x,y) and employing Leibnitz rule, we rewrite Equation (8) as where C, H and S denote the property defined in the xy-plane (herein we call this special plane the horizontal plane of average head gradient, and the derivation can be found in Section 2.4.).The V z | z i`1 and V z | z i denote the water flux vectors at the upper and bottom aquifer, respectively.

The Horizontal Plane of Average Head Gradient in the Triangular Prism Element
The horizontal plane of average head gradient in an irregular triangular prism element means the average lateral head gradient in this irregular triangular prism element, and it is used to approximately calculate the net lateral flux.There are two assumptions made to derive the horizontal plane of average head gradient, (1) the head varies linearly between the upper and bottom nodes within the element; Water 2016, 8, 195 5 of 22 (2) the head in the horizontal plane is linearly interpolated by the head of three nodes I, J, K, which lead to H px, y, tq " 1 2∆ rpa i x `bi y `di q H I ``a j x `bj y `dj ˘HJ `pa k x `bk y `dk q H K s (11) where H P (P = I, J, K) are the head of three nodes in the horizontal plane of triangular prism element IJK [L]; H p (p = i, j, k) and H p' (p' = i', j', k') are the head of six nodes in the triangular prism element; ∆ is the area of the triangle IJK; a i " y j ´yk ; b i " x k ´xj ; d i " x j y k ´xk y j ; a j " y k ´yi ; b j " x i ´xk ; d j " x k y i ´xi y k ; a k " y i ´yj ; b k " x j ´xi ; d k " x i y j ´xj y i ; x i , y i (i = i, j, k) are the x and y coordinates of node i, j and k [L]; β P (P=I, J, K) are calculated as where z i (i = i, i', j, j', k, k') are the z-coordinates of six nodes in the triangular prism element [L]; and z is the z-coordinate of the horizontal plane of average head gradient [L].Derived from Equation ( 11), the head gradients along the xand y-directions can be calculated.Since they have a similar form, only the head gradient along the x-direction in the plane IJK is shown here in detail, where, In the triangular prism element ijki'j'k', the horizontal plane of average head gradient can be calculated as By combined consideration of Equation (11) and Equation ( 16), the vertical coordinate of the horizontal plane of average head gradient IJK in the triangular prism element ijki'j'k' is given as Water 2016, 8, 195 6 of 22

The Discretization of the Lateral flow Equation
As shown in Equation ( 9 The water balance analysis method is adopted, and the net lateral flux from the plane IJK to node I is expressed as follows [27]: where ∆ is the area of triangle IJK; and a p , b p (p = i, j, k) have the same meaning as mentioned above.
It should be noted that the x and y coordinates of the three nodes in the plane IJK equal those in the upper or bottom surface of the triangular prism element.In the triangular prism element ijki'j'k', the lateral flux to node i can be calculated by the lateral flux to node I and the thickness of the control volume of node i, which is one half the average thickness of the element.Then, by comprehensive consideration with the assumption that the head varies linearly with the depth in each triangular prism element shown in Equation (10) and Equation ( 18), the lateral flux to node i in the element ijki'j'k' can be expressed as where K is the mean hydraulic conductivity of the triangular prism element; a i 2 `bi 2 " P ii , a i a j `bi b j "

The Discretization of the Water Mass Change Item
The first item in Equation ( 9) characterizes the capacity of an aquifer to release or store water in response to a unit hydraulic head change.The water flow change in the control volume of node i in the triangular prism element can be expressed as The water mass change in the control volume of node i in the triangular prism element from the source/sink items can be calculated by

The Element Matrices
Equations ( 19), ( 20) and ( 21) provide, respectively, the lateral flow, water mass change from a unit head change and source/sink items in the control volume of node i in the element ijki'j'k'.Similarly, these items in the control volume for the other five nodes can be calculated.Then, we obtain the discretized Richards' equation in terms of matrices in the element ijki'j'k' as follows (excluding the vertical flux, which will be calculated in the next section), P ii β I P ij β J P ik β K P ii p1 ´βI q P ij p1 ´βJ q P ik p1 ´βK q P ji β I P jj β J P jk β K P ji p1 ´βI q P jj p1 ´βJ q P jk p1 ´βK q P ki β I P kj β J P kk β K P ki p1 ´βI q P kj p1 ´βJ q P kk p1 ´βK q P ii β I P ij β J P ik β K P ii p1 ´βI q P ij p1 ´βJ q P ik p1 ´βK q P ji β I P jj β J P jk β K P ji p1 ´βI q P jj p1 ´βJ q P jk p1 ´βK q P ki β I P kj β J P kk β K P ki p1 ´βI q P kj p1 ´βJ q P kk p1 ´βK q

Coupling Horizontal Flow with the Vertical Fluxes
The vertical flow flux is represented by 9).As specified in the element ijki'j'k', we calculate the net vertical flow to the control volume of node i as an example (see Figure 2), by using FDM.When the node i is located between the ground surface and the bottom of the aquifer, the vertical flux to node i can be calculated as, where q i+1/2 and q i´1/2 are the upper and bottom fluxes, respectively; A i is the control area of node i in the triangular prism element [L 2 ].Node i´1 and node i+1 are the adjacent nodes to the node i in the vertical direction; K i,i´1 is the geometric average hydraulic conductivities of the nodes i and i´1; and B i,i´1 is the distance between the nodes i and i´1.K i,i`1 and B i,i+1 are defined in the same way.
Water 2016, 8, 195 8 of 22 When the node i is located at the bottom of the aquifer, qi+1/2 becomes zero and then the vertical flux will be ( )

Illustrative Example to Elaborate the Coupling Method of Horizontal Layers
The spatial discretized Equation ( 22) is further discretized by the implicit method temporally to form the calculation matrix; then, the unified matrix of the unsaturated-saturated water flow is developed by assembling the vertical fluxes represented by Equations ( 23)-( 25) between each layer.A two layer aquifer with 3n nodes is presented to show how the coupling procedure implements to form the global matrix.For this case, the global matrix is (3n) × (3n) as shown in Equation ( 26).When the node i is located at the ground surface, q i+1/2 becomes zero and the vertical flux can be calculated as When the node i is located at the bottom of the aquifer, q i+1/2 becomes zero and then the vertical flux will be

Illustrative Example to Elaborate the Coupling Method of Horizontal Layers
The spatial discretized Equation ( 22) is further discretized by the implicit method temporally to form the calculation matrix; then, the unified matrix of the unsaturated-saturated water flow is developed by assembling the vertical fluxes represented by Equations ( 23)-( 25) between each layer.A two layer aquifer with 3n nodes is presented to show how the coupling procedure implements to form the global matrix.For this case, the global matrix is (3n) ˆ(3n) as shown in Equation (26).The small letters w j i , b i (i, j=1, 2, . . ., 3n) represent the items formed by temporally discretizing Equation ( 22) with the implicit method, and the capital letters W j i (i, j = n+1, n+2, . . ., 2n+1) represent the modified items according to Equations ( 23), ( 24) and ( 25), According to Equation ( 23), W j j pj " n `1, ¨¨¨, 2nq items for the nodes between the ground surface and the aquifer bottom are modified by According to Equation ( 24), W j j pj " 2n `1, ¨¨¨, 3nq items for the ground surface nodes are modified by Water 2016, 8, 195 According to Equation ( 25), W j j pj " 1, ¨¨¨, nq items for the bottom nodes are modified by W n`j j " w n`j j ´Aj K j,j`n B j,j`n (33)

Boundary Conditions and Source/Sink Items
The model can handle various kinds of boundary conditions and source/sink items, such as the first-and second-kind boundary conditions, root uptake, pumping well, and atmosphere boundary conditions (e.g., precipitation and evaporation or evapotranspiration).The details are as follows.
(1) First-kind boundary condition: where ϕ is the prescribed head in the boundary [L].(2) Second-kind boundary condition: where q is the prescribed flow flux in the boundary [L 2 T ´1].
(3) Pumping well: where Q p,i and Q p,i' are the pumping rate of the nodes i and i', respectively [L 3 T ´1]; L i is the length of the pumping well filtration in the i-th layer [L]; and Q w is the total pumping rate [L 3 T ´1]; β I is calculated by Equation (12).If the drawdown caused by pumping is very large to disconnect the nodes with the groundwater system, the pumping rate will be set as zero without extraction.(4) Root uptakeThe sink term of root uptake, S, is calculated by [28] as where α(h) is a pressure response function of root uptake [-], and S p is the potential water uptake by plant [T ´1].(5) Atmosphere boundary conditions The model can handle daily changing atmosphere boundary conditions such as precipitation rate and evapotranspiration rate.The precipitation rate is input according to the climate information.There are two ways to handle the evapotranspiration boundary.One is pre-processing the evapotranspiration rate as input information.The input evapotranspiration rate (ET) is calculated by multiplying the daily evaporation from water surface E 0 with a coefficient α, which includes the impacts of soil water conditions, hydraulic engineering conditions and agricultural management policies.Another way is to input the climate information needed by the Penman-Monteith equation [29] to calculate the reference evapotranspiration rate (ET 0 ) and then the model will calculate the evapotranspiration rate by subsequently multiplying ET 0 with the crop coefficient.

The Flow Chart of the Model
The model is written in FORTRAN.It runs from processing the input information, then calculating the interpolation table of the soil water parameters of θ(h), K(h) and C(h).Judgement is made on whether the nodes are in the saturated or unsaturated zone.If the nodes are in the unsaturated zone, it is necessary to calculate the soil water capacity and hydraulic conductivity according to Equations ( 6) and (7).Then, the model moves to form the element matrices, considering the boundary conditions.The implicit method of temporal discretization is adopted to form the unified global matrix.The global matrix is then modified with the vertical fluxes according to Equations ( 23)- (25).The final global matrix can be solved by solvers such as the ORTHOFEM package [30].If the simulation results converge, the model will move to the next time step.Otherwise, the iteration will be continued in the current time step.The flow chart of the model can be seen in Figure 3.

Case Study
In order to test the validity, accuracy and reliability of the model code, three hypothetical cases with 1-D infiltration flow, 2-D lateral flow and 3-D well flow are implemented.Results for the hydraulic head, soil moisture content and water table contours from the proposed model are compared with the 1-D unsaturated-saturated model Hydrus1D [31], 2-D model SWMS2D [32] and 3-D model FEFLOW.Furthermore, the model is applied to a regional scale irrigation district located in Hetao irrigation district to simulate the groundwater fluctuations to assess the model applicability and efficiency.The model calibration and prediction results of water table are assessed by measurements.

Case Study
In order to test the validity, accuracy and reliability of the model code, three hypothetical cases with 1-D infiltration flow, 2-D lateral flow and 3-D well flow are implemented.Results for the hydraulic head, soil moisture content and water table contours from the proposed model are compared with the 1-D unsaturated-saturated model Hydrus1D [31], 2-D model SWMS2D [32] and 3-D model FEFLOW.Furthermore, the model is applied to a regional scale irrigation district located in Hetao irrigation district to simulate the groundwater fluctuations to assess the model applicability and efficiency.The model calibration and prediction results of water table are assessed by measurements.

Case 1: 1-D Infiltration Flow
Case 1 considers 1-D infiltration water flow in a 3 m soil column, subjected to a second-kind upper boundary condition with a flux rate of 0.0005 m¨d ´1.The soil water flow parameters are, θ m = θ s = 0.35, θ r = θ a = 0.057, K s = 0.6 m¨d ´1, α = 4.1 m ´1, n =2.28.The column bottom is prescribed as a no-flow boundary condition.The initial water table is located at a depth of 1.3 m from the soil surface.The pressure head and soil water moisture from the proposed model are compared with those from Hydrus1D, as shown in Figures 4 and 5.The comparison results indicate that the absolute errors for the head are less than 0.004 m, and less than 0.0003 for soil moisture.The results indicate that the proposed model is able to precisely capture the flow information in the vertical profile.In this case, since 3-D meshes are used in the proposed model and consequently many more nodes are calculated, the computational cost of our model is larger than that of Hydrus1D.The simulation time of the proposed model is 42 s whereas that of Hydrus1D is less than 1 s.

Case 2: 2-D Water Flow
The model is used to simulate a typical 2-D lateral flow case with two rivers 40 m apart and with precipitation infiltration from the upper boundary.The infiltration rate is 0.002 m¨d ´1.The soil profile is homogenous and 3 m deep.The soil water parameters are θ m = θ s = 0.30, θ r = θ a = 0.02, K s = 0.5 m¨d ´1, α = 4.1 m ´1, and n = 1.964.When the water flow reaches a steady state, the analytical solution of hydraulic head for this case is given as where H is the hydraulic head at the location with a distance of x from the left river [L]; H The 2-D unsaturated-saturated model SWMS2D is also adopted to simulate this case to obtain the numerical solution.The domain is discretized into 21 layers in the simulation of the two models.The steady-state water table, where the zero location is at the bottom, i.e., z = 0 m, obtained with the proposed model, the analytical solution and SWMS2D are plotted in Figure 6.The absolute errors of the proposed model related to the analytical solution and SWMS2D are 0.0075 m and ´0.0072 m, respectively.In order to check the model accuracy in the unsaturated zone, the simulated pressure head profiles at different locations are compared with those from SWMS2D as shown in Figure 7.The results show that the model can capture the water flow information in both the saturated and unsaturated zones under the lateral flow condition.the proposed model related to the analytical solution and SWMS2D are 0.0075 m and −0.0072 m, respectively.In order to check the model accuracy in the unsaturated zone, the simulated pressure head profiles at different locations are compared with those from SWMS2D as shown in Figure 7.
The results show that the model can capture the water flow information in both the saturated and unsaturated zones under the lateral flow condition.In this example, the node number of the proposed model and FEFLOW are the same.The computational time is 90 s for the proposed model versus 137 s for FEFLOW, which indicates that the proposed model is efficient while still providing satisfactory simulation accuracy.

Model Application to a Regional Scale Irrigation District
After testing the validity of the proposed model, we further apply it to a regional-scale case, the Yonglian irrigation area in Hetao irrigation district, Inner Mongolia, China (Figure 9), with an area of approximately 29.1 km 2 .The ground surface is very flat as the elevation changes only from 1028.9 m to 1025.4 m from the southwest to the northeast.This is a typical arid area with an average annual precipitation of 169 mm, but an average annual potential evaporation of 2065 mm.Therefore, In this example, the node number of the proposed model and FEFLOW are the same.The computational time is 90 s for the proposed model versus 137 s for FEFLOW, which indicates that the proposed model is efficient while still providing satisfactory simulation accuracy.

Model Application to a Regional Scale Irrigation District
After testing the validity of the proposed model, we further apply it to a regional-scale case, the Yonglian irrigation area in Hetao irrigation district, Inner Mongolia, China (Figure 9), with an area of approximately 29.1 km 2 .The ground surface is very flat as the elevation changes only from 1028.9 m to 1025.4 m from the southwest to the northeast.This is a typical arid area with an average annual precipitation of 169 mm, but an average annual potential evaporation of 2065 mm.Therefore, irrigation is very important to ensure the water requirements of this region during the crop growing season.The irrigation and drainage canal systems clearly delineate the hydro-geological borders in this area, which are formed by the No. 6 Drainage Ditch and Yongshen Ditch in the east and the Zhaohuo Trunk Canal and Naiyong Ditch in the southwest (Figure 9).
First-kind boundary conditions are applied to these two canal segments, which are filled most of the time with irrigation or drainage water.The other segments are considered as the seepage face boundary condition.When the nodes become saturated, the pressure heads are set as zero and then the fluxes can be calculated.Otherwise, the fluxes are set as zero and the pressure heads can be calculated.The daily precipitation rate is measured in this area and the daily amount of irrigation water from Zhaohuo Trunk Canal is also measured.The irrigation rate is calculated by dividing the daily irrigation water amount diverted from Yellow river by the area of the farmland.The daily reference evapotranspiration rate is calculated by multiplying the daily evaporation from the water surface with a coefficient α, which is 0.68 in our study.From Figure it can be seen that the major land use types are farmland, bare saline soil, and villages.The assumptions here are that the area of villages is classified as bare soil, and the crop types in the farmland are not further distinguished.Therefore, the domain is divided into two sub-areas horizontally according to land use types.In each sub-area, specified upper boundary conditions (mainly with different irrigation and evaporation rates) are imposed.The details will be found in the model parameter calibration and model predictions in the following paragraphs.The aquifer geometry and hydro-geological information are obtained from the boreholes and pumping From Figure 9, it can be seen that the major land use types are farmland, bare saline soil, and villages.The assumptions here are that the area of villages is classified as bare soil, and the crop types in the farmland are not further distinguished.Therefore, the domain is divided into two sub-areas horizontally according to land use types.In each sub-area, specified upper boundary conditions (mainly with different irrigation and evaporation rates) are imposed.The details will be found in the model parameter calibration and model predictions in the following paragraphs.The aquifer geometry and hydro-geological information are obtained from the boreholes and pumping tests [33].Provided by the Geological Department of Inner Mongolia, there is a 1-m thick impervious clay layer at the depth of 53 m from the ground surface, which is, consequently, as a non-flow bottom boundary in the model.Thus, the simulated aquifer has a depth of 53 m, with the top 7 m of the aquifer consisting of loamy sand with low hydraulic conductivity, and an underlying 46 m deep sand aquifer with high permeability.The whole domain is divided into 12 layers vertically and 1336 triangular prism elements in each horizontal layer to guarantee simulation accuracy.
Ten observation wells without pumping are located in this area.The water table of these 10 wells and the soil moisture contents at the soil depths of 10 cm, 30 cm, 50 cm, 70 cm, and 140 cm are measured.The water table measurements of these 10 observation wells from 1 May to 1 November 2007 are used to calibrate the model parameters.The upper boundary conditions with the precipitation in the whole domain, the irrigation rate in the farmland, the evapotranspiration in the farmland and the evaporation in the bare soil in 2007 are shown in Figure 10a.The observed water tables at the 10 observation wells on 1 May 2007 are used to interpolate the initial hydraulic head for the transient simulations.6) and ( 7) are listed in Table 1.The groundwater level data from the 10 observation wells are used to evaluate the model performance.Observed and simulated groundwater levels are compared through linear regression, and the respective coefficient of determination for the groundwater levels is 0.8147 as plotted in Figure 11.The simulated groundwater level results of the randomly selected four observation wells, which are wells 1, 3, 5, 7, are presented in Figure 12.It can be seen that the water table dynamics have a very similar trend mainly due to the intensive irrigation and evaporation rates in this area, which result in water table fluctuations occurring mainly in the vertical direction.The same soil characteristics in the horizontal direction according to the hydro-geological data also contribute to the similar water table fluctuation trend in these four locations.The results show that the model can represent the fluctuations of the observed groundwater table well.Therefore, the parameters used in the calibrated model can be considered to carry out the model prediction the subsequent year 2008.The calibrated soil water parameters of Equations ( 6) and ( 7) are listed in Table 1.         13 shows the simulated and measured groundwater levels over time for the four randomly chosen wells.From Figure 13, it can be seen that the model maintains the same change trend as the observations and can capture the fluctuations of water tables.However, at simulation time t = 108 d, there is an abrupt rise in the water table observed at all observation wells, which is caused by an intensive precipitation rate (Figure 10b).The model fails to simulate this strong change in the water table, resulting in large errors in the water table at this time, as shown in Figure 13.The differences may be attributed to the high precipitation rate causing water ponding at the surface, which cannot be included in the proposed model.The uncertainties in the model parameters and model inputs as well as the setup of the boundary conditions would also affect the simulation results.
Then, the 2007 calibrated model is applied to simulate the unsaturated-saturated water flow in 2008, with the upper boundary conditions shown in Figure 10b.The simulation runs from 1 May 2008 to 11 November 2008, lasting 195 days.Figure 13 shows the simulated and measured groundwater levels over time for the four randomly chosen wells.From Figure 13, it can be seen that the model maintains the same change trend as the observations and can capture the fluctuations of water tables.However, at simulation time t = 108 d, there is an abrupt rise in the water table observed at all observation wells, which is caused by an intensive precipitation rate (Figure 10b).The model fails to simulate this strong change in the water table, resulting in large errors in the water table at this time, as shown in Figure 13.The differences may be attributed to the precipitation rate causing water ponding at the surface, which cannot be included in the proposed model.The uncertainties in the model parameters and model inputs as well as the setup of the boundary conditions would also affect the simulation results.
Figure 14 shows the water table distributions at two different times (t = 6 d and t = 180 d).It can be seen that the flow direction is mainly from southwest to northeast.The water table at time 180 d is higher than that of day 6 because of the large amount of irrigation water applied in the autumn season to wash out the salt in the unsaturated zone.These results demonstrate that the model works well as the variations of the hydraulic head are well simulated, even with these complicated boundary conditions and changing water supply conditions.The computational time needed to perform the calculation is of paramount importance, because it seriously affects the model's applicability to regional-scale problems.The domain is divided into  The computational time needed to perform the calculation is of paramount importance, because it seriously affects the model's applicability to regional-scale problems.The domain is divided into Figure 14 shows the water table distributions at two different times (t = 6 d and t = 180 d).It can be seen that the flow direction is mainly from southwest to northeast.The water table at time 180 d is higher than that of day 6 because of the large amount of irrigation water applied in the autumn season to wash out the salt in the unsaturated zone.These results demonstrate that the model works well as the variations of the hydraulic head are well simulated, even with these complicated boundary conditions and changing water supply conditions.
The computational time needed to perform the calculation is of paramount importance, because it seriously affects the model's applicability to regional-scale problems.The domain is divided into 8760 nodes and 16030 elements.The model runs on an 8 GB RAM, eight 3.70 GHz Intel Xeon CPU-based computer, and it takes the coupling model 642 s to perform this case with 1967 time steps.Thus, the model can be considered as highly efficient to simulate a regional scale unsaturated-saturated water flow problem.

Conclusions
In this paper, the VHS concept is used to simplify 3-D unsaturated-saturated water flow models.This paper presents not only the theoretical basis of the VHS-based modeling but also three synthetic cases for methodology evaluation and a real world application to a regional-scale problem.The model can capture the details of water flow in both the soil zone and groundwater system as shown in the three synthetic cases, which demonstrate the accuracy of the model code.The real world application to an irrigation district with complex boundary conditions shows great applicability to complicated large-scale regional groundwater flow modeling with relatively low computational cost.It means that the VHS-based modeling has a great potential for simulating water dynamics under changing management conditions of water use reduction, which is important for water sustainability in arid and semi-arid areas [34,35].
As this study mainly focuses on the model development, the uncertainty in model calibration, various scenarios of measurements, and boundary conditions are not discussed [36].Further, the current model treats rivers as the Dirichlet boundary condition by assigning hydraulic heads to the river nodes.However, a stream package for complex river conditions is needed.These two areas of study are warranted in a future investigation.

2. 1 .
Discretization of the Aquifer The model is established firstly by separating the 3-D aquifer into different layers.Each layer is divided into triangular prism elements.The horizontal plane of average head gradient is derived to represent the average horizontal flow in the triangular prism element.The aquifer domain, the divided triangular prism element and the horizontal plane of average head gradient are shown in Figure 1.The governing 3-D Richards' equation is split into a 2-D horizontal water flow equation based on the horizontal plane of average head gradient in the triangular prism element and a 1-D vertical water equation.The water balance method is used to calculate the horizontal flow flux and the FDM is adopted for calculating the vertical flow flux.The horizontal flow and the 1-D vertical flow are integrated into one unified matrix and then solved simultaneously.The details are shown in the next sections.Water 2016, 8, 195 3 of 22 vertical equation is discretized by the FDM to integrate the different horizontal layers.The two equations are solved simultaneously by coupling them into a unified nonlinear system with a single matrix.The efficiency and accuracy of the developed model and its computer code are evaluated by using three synthetic cases represented the 1-D infiltration flow, 2-D lateral flow, and 3-D well flow.

Figure 1 .
The governing 3-D Richards' equation is split into a 2-D horizontal water flow equation based on the horizontal plane of average head gradient in the triangular prism element and a 1-D vertical water equation.The water balance method is used to calculate the horizontal flow flux and the FDM is adopted for calculating the vertical flow flux.The horizontal flow and the 1-D vertical flow are integrated into one unified matrix and then solved simultaneously.The details are shown in the next sections.

Figure 1 .
Figure 1.The subsurface system and the triangular prism element as well as the horizontal plane of average head gradient in the triangular prism element.

Figure 1 .
Figure 1.The subsurface system and the triangular prism element as well as the horizontal plane of average head gradient in the triangular prism element.
˙means the lateral flow.Since B is the thickness of the triangular prism element and can be considered as constant, the key to calculate the lateral flow is to discretize the

Figure 2 .
Figure 2. The water balance volume of node i. Qi means the lateral flux shown in Equation (19).

Figure 2 .
Figure 2. The water balance volume of node i. Q i means the lateral flux shown in Equation (19).

5. 1 . 1 .
Case 1: 1-D Infiltration Flow Case 1 considers 1-D infiltration water flow in a 3 m soil column, subjected to a second-kind upper boundary condition with a flux rate of 0.0005 m•d −1 .The soil water flow parameters are, θm= θs = 0.35, θr = θa = 0.057, Ks = 0.6 m•d −1 , α = 4.1 m −1 , n =2.28.The column bottom is prescribed as a no-flow boundary condition.The initial water table is located at a depth of 1.3 m from the soil

Figure 3 .
Figure 3.The flow chart of the model.

Figure 4 .
Figure 4. Pressure head profiles from the proposed model and Hydrus1D.

Figure 5 .
Figure 5. Soil moisture profiles from the proposed model and Hydrus1D.

Figure 4 .
Figure 4. Pressure head profiles from the proposed model and Hydrus1D.

Figure 4 .
Figure 4. Pressure head profiles from the proposed model and Hydrus1D.

Figure 5 .Figure 5 .
Figure 5. Soil moisture profiles from the proposed model and Hydrus1D.

1 and H 2
are the constant hydraulic heads in the left and right rivers [L], respectively; l is the distance between the two rivers [L]; W is the precipitation rate [L T ´1]; and K s is the saturated hydraulic conductivity [L T ´1].

Figure 6 .
Figure 6.Steady-state water tables obtained with the proposed model, the analytical solution and SWMS2D.

Figure 6 .
Figure 6.Steady-state water tables obtained with the proposed model, the analytical solution and SWMS2D.

Figure 6 .
Figure 6.Steady-state water tables obtained with the proposed model, the analytical solution and SWMS2D.

Figure 7 .
Figure 7. Steady-state pressure head profiles obtained with the proposed model and SWMS2D.
Figure 8 shows the head contours of the two models at the end of the simulation.The two head contour categories of the figure indicate that the heads of the proposed model match the solution of FEFLOW well.

Figure 7 .
Figure 7. Steady-state pressure head profiles obtained with the proposed model and SWMS2D.5.1.3.Case 3: 3-D Well Flow An aquifer with the size of 200 m ˆ200 m ˆ20 m is considered.There are two pumping wells located at (100, 48) m and (100, 152) m with a pumping rate of 500 m 3 ¨d´1 .A constant head of 18 m is imposed around the domain, and the bottom is specified as a no-flow boundary.The water flow parameters are θ m = θ s = 0.30, θ r = θ a = 0.02, K s = 6.4 m¨d ´1, α = 4.1 m ´1, and n = 1.964.The simulation time is 200 days.The 3-D unsaturated-saturated water flow model FEFLOW is adopted to provide the reference solution.Figure 8 shows the head contours of the two models at the end of the simulation.The two head contour categories of the figure indicate that the heads of the proposed model match the solution of FEFLOW well.

22 Figure 8 .
Figure 8. Head contours from the proposed model and FEFLOW.

Figure 8 .
Figure 8. Head contours from the proposed model and FEFLOW.

Figure 9 .
Figure 9. Location of the Yonglian irrigation area in Hetao irrigation district.The black border line of the Yonglian irrigation area is specified with the first-kind boundary condition, and the red border line denotes the seepage face boundary condition.

Figure 9 .
Figure 9. Location of the Yonglian irrigation area in Hetao irrigation district.The black border line of the Yonglian irrigation area is specified with the first-kind boundary condition, and the red border line denotes the seepage face boundary condition.

Figure 10 .
Figure 10.Upper boundary conditions applied to the various sub-areas of the domain in (a) 2007 and (b) 2008.

Figure 10 .
Figure 10.Upper boundary conditions applied to the various sub-areas of the domain in (a) 2007 and (b) 2008.

Figure 10 .
Figure 10.Upper boundary conditions applied to the various sub-areas of the domain in (a) 2007 and (b) 2008.

Figure 11 .
Figure 11.Comparison of observed and simulated groundwater levels for the 10 observation wells in 2007.

Figure 11 .
Figure 11.Comparison of observed and simulated groundwater levels for the 10 observation wells in 2007.
water table fluctuations occurring mainly in the vertical direction.The same soil characteristics in the horizontal direction according to the hydro-geological data also contribute to the similar water table fluctuation trend in these four locations.The results show that the model can represent the fluctuations of the observed groundwater table well.Therefore, the parameters used in the calibrated model can be considered to carry out the model prediction for the subsequent year 2008.The calibrated soil water parameters of Equations (

Table 1 .
Calibrated soil water parameters of the Yonglian irrigation area.Then, the 2007 calibrated model is applied to simulate the unsaturated-saturated water flow in 2008, with the upper boundary conditions shown in Figure 10b.The simulation runs from 1 May 2008 to 11 November 2008, lasting 195 days.Figure