Next Article in Journal
Coliphages as Model Organisms in the Characterization and Management of Water Resources
Previous Article in Journal
The Status of Domestic Water Demand: Supply Deficit in the Kathmandu Valley, Nepal
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
2
Department of Scientific Computing, Florida State University, Tallahassee, FL 32306, USA
*
Author to whom correspondence should be addressed.
Water 2016, 8(5), 195; https://doi.org/10.3390/w8050195
Submission received: 16 October 2015 / Revised: 24 March 2016 / Accepted: 12 April 2016 / Published: 12 May 2016

Abstract

:
Due to the high nonlinearity of the three-dimensional (3-D) unsaturated-saturated water flow equation, using a fully 3-D numerical model is computationally expensive for large scale applications. A new unsaturated-saturated water flow model is developed in this paper based on the vertical/horizontal splitting (VHS) concept to split the 3-D unsaturated-saturated Richards’ equation into a two-dimensional (2-D) horizontal equation and a one-dimensional (1-D) vertical equation. The horizontal plane of average head gradient in the triangular prism element is derived to split the 3-D equation into the 2-D equation. The lateral flow in the horizontal plane of average head gradient represented by the 2-D equation is then calculated by the water balance method. The 1-D vertical equation is discretized by the finite difference method. The two equations are solved simultaneously by coupling them into a unified nonlinear system with a single matrix. Three synthetic cases are used to evaluate the developed model code by comparing the modeling results with those of Hydrus1D, SWMS2D and FEFLOW. We further apply the model to regional-scale modeling to simulate groundwater table fluctuations for assessing the model applicability in complex conditions. The proposed modeling method is found to be accurate with respect to measurements.

1. 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.

2. Mathematical and Numerical Method

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.

2.2. The Governing Equations

The general 3-D unsaturated-saturated water flow is represented as follows:
θ t = x i [ K H x j ) ] S
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]; xi, xj (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
θ t = C ( h ) H t
where C(h) is the soil water capacity [L−1] and can be calculated as
C ( h ) = θ n μ + d θ d h
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 Ks. Therefore, for the saturated water flow, the water flow Equation (1) can be simplified as
μ H t = x i [ K s H x j ] S
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
C ( h ) H t = x i [ K ( h ) H x i ] S
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
C = { m n α n ( θ m θ a ) [ 1 + ( α h ) n ] m 1 h n 1 ( in the unsaturated zone ) μ ( in the saturated zone )
K = { K s ( θ θ r θ s θ r ) 1 / 2 [ 1 ( ( 1 θ θ r θ s θ r ) 1 / m ) m ] 2 ( in the unsaturated zone ) K s ( in the saturated zone )
where θm is the fictitious saturated water content in the soil water retention function [–]; θr is the residual water content [–]; θa is the parameter in the soil water retention function [–]; θs is the saturated water content [–]; Ks is the saturated hydraulic conductivity [L T−1]; α is the coefficient in the soil water retention function [L−1]; m, n are the exponents in the soil water retention function [–].

2.3. 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
z i ( x , y ) z i + 1 ( x , y ) C H t d z = z i ( x , y ) z i + 1 ( x , y ) [ x ( K H x ) + y ( K H y ) + z ( K H z ) ] d z z i ( x , y ) z i + 1 ( x , y ) S d z
where zi(x,y) and zi+1(x,y) are the bottom and upper boundary, respectively, of the aquifer. By assuming B = zi+1(x,y)- zi(x,y) and employing Leibnitz rule, we rewrite Equation (8) as
B C ¯ H ¯ t = B x ( K H ¯ x ) + B y ( K H ¯ y ) + V z | z i V z | z i + 1 B S ¯
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.

2.4. 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; (2) the head in the horizontal plane is linearly interpolated by the head of three nodes I, J, K, which lead to
H P = β P H p + ( 1 β P ) H p
H ( x , y , t ) = 1 2 Δ [ ( a i x + b i y + d i ) H I + ( a j x + b j y + d j ) H J + ( a k x + b k y + d k ) H K ]
where HP (P = I, J, K) are the head of three nodes in the horizontal plane of triangular prism element IJK [L]; Hp (p = i, j, k) and Hp’ (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 y k ; b i = x k x j ; d i = x j y k x k y j ; a j = y k y i ; b j = x i x k ; d j = x k y i x i y k ; a k = y i y j ; b k = x j x i ; d k = x i y j x j y i ; xi, yi (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
β I = ( z ¯ z i ) / ( z i z i ) β J = ( z ¯ z j ) / ( z j z j ) β K = ( z ¯ z k ) / ( z k z k )
where zi (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 x- and 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,
H x = 1 2 Δ [ ( a i A i i + a j A j j + a k A k k ) z + ( a i H i i 0 + a j H j j 0 + a k H k k 0 ) ]
where,
Δ = 1 2 ( a j b k a k b j )
A i i = ( H i H i z i z i ) , H i i 0 = H i H i H i z i z i z i z i < z < z i   A j j = H j H j z j z j , H j j 0 = H j H j H j z j z j z j z j < z < z j A k k = H k H k z k z k , H k k 0 = H k H k H k z k z k z k z k < z < z k
In the triangular prism element ijki’j’k’, the horizontal plane of average head gradient can be calculated as
z k z i ( H x S ) d z z k z i S d z = 1 2 Δ ( a i A i i + a j A j j + a k A k k ) × { ( a j b k a k b j ) 24 [ ( z i z j ) ( z i + 3 z j ) ( z i z k ) ( z i + 3 z k ) ] + Δ ( z j z i ) ( z j + z i ) 2 + Δ ( z i z k ) ( z i + z k ) 2 } ( a j b k a k b j ) 2 [ ( z i z j ) ( z i z k ) ] 3 + Δ ( z j z k ) + 1 2 Δ ( a i H i i 0 + a j H j j 0 + a k H k k 0 )
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
z ¯ = { 1 12 [ ( z i z j ) ( z i + 3 z j ) ( z i z k ) ( z i + 3 z k ) + ( z j z i ) ( z j + z i ) 2 + ( z i z k ) ( z i + z k ) 2 ] } { [ ( z i z j ) ( z i z k ) ] 3 + ( z j z k ) }

2.5. The Discretization of the Lateral flow Equation

As shown in Equation (9), the item B x ( K H ¯ x ) + B y ( K H ¯ y ) 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 x ( K H ¯ x ) + y ( K H ¯ y ) term. The water balance analysis method is adopted, and the net lateral flux from the plane IJK to node I is expressed as follows [27]:
Q I = K 4 Δ [ ( b i 2 + a i 2 ) H I + ( b j b i + a j a i ) H J + ( b k b i + a k a i ) H K ]
where Δ is the area of triangle IJK; and ap, bp (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
Q i = K ¯ × B 8 Δ [ P i i β I H i + P i i ( 1 β I ) H i + P i j β J H j + P i j ( 1 β J ) H j + P i k β K H k + P i k ( 1 β K ) H K ]
where K ¯ is the mean hydraulic conductivity of the triangular prism element; a i 2 + b i 2 = P i i , a i a j + b i b j = P i j , a i a k + b i b k = P i k .

2.6. 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
Q c = Δ 3 C ¯ B 2 Δ H i Δ t
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
Q s = B 2 × Δ 3 × S ¯ = 1 6 Δ × B × S ¯

2.7. 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),
K ¯ B 8 Δ [ P i i β I P i j β J P i k β K P i i ( 1 β I ) P i j ( 1 β J ) P i k ( 1 β K ) P j i β I P j j β J P j k β K P j i ( 1 β I ) P j j ( 1 β J ) P j k ( 1 β K ) P k i β I P k j β J P k k β K P k i ( 1 β I ) P k j ( 1 β J ) P k k ( 1 β K ) P i i β I P i j β J P i k β K P i i ( 1 β I ) P i j ( 1 β J ) P i k ( 1 β K ) P j i β I P j j β J P j k β K P j i ( 1 β I ) P j j ( 1 β J ) P j k ( 1 β K ) P k i β I P k j β J P k k β K P k i ( 1 β I ) P k j ( 1 β J ) P k k ( 1 β K ) ] [ H i H j H k H i H j H k ] + C ¯ Δ B 6 [ 1 1 1 1 1 1 ] [ d H i / d t d H j / d t d H k / d t d H i / d t d H j / d t d H k / d t ] Δ B S ¯ 6 [ 1 1 1 1 1 1 ] = 0

2.8. Coupling Horizontal Flow with the Vertical Fluxes

The vertical flow flux is represented by V z | z i V z | z i + 1 in Equation (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,
Q v = A i ( q i 1 / 2 q i + 1 / 2 ) = A i K i , i 1 ¯ B i , i 1 H i 1 ( A i K i , i + 1 ¯ B i , i + 1 + A i K i , i 1 ¯ B i , i 1 ) H i + A i K i , i + 1 ¯ B i , i + 1 H i + 1
where qi+1/2 and qi−1/2 are the upper and bottom fluxes, respectively; Ai is the control area of node i in the triangular prism element [L2]. 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 Bi,i−1 is the distance between the nodes i and i−1. K i , i + 1 ¯ and Bi,i+1 are defined in the same way.
When the node i is located at the ground surface, qi+1/2 becomes zero and the vertical flux can be calculated as
Q v _ u p p e r = A i ( q i 1 / 2 ) = A i K i , i 1 ¯ B i , i 1 H i 1 A i K i , i 1 ¯ B i , i 1 H i
When the node i is located at the bottom of the aquifer, qi+1/2 becomes zero and then the vertical flux will be
Q v _ B o t t o m = A i ( q i + 1 / 2 ) = A i K i , i + 1 ¯ B i , i + 1 H i + A i K i , i + 1 ¯ B i , i + 1 H i + 1

2.9. 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 i j , bi (i, j=1, 2, …, 3n) represent the items formed by temporally discretizing Equation (22) with the implicit method, and the capital letters W i j (i, j = n+1, n+2, …, 2n+1) represent the modified items according to Equations (23), (24) and (25),
[ W 1 1 w 1 n W 1 n + 1 w n 1 W n n W n 2 n W n + 1 1 W n + 1 n + 1 w n + 1 2 n W n + 1 2 n + 1 W 2 n n w 2 n n + 1 W 2 n 2 n W 2 n 3 n W 2 n + 1 n + 1 W 2 n + 1 2 n + 1 w 2 n + 1 3 n W 3 n 2 n w 3 n 2 n + 1 W 3 n 3 n ] [ H 1 H n H n + 1 H 2 n H 2 n + 1 H 3 n ] = [ b 1 b n b n + 1 b 2 n b 2 n + 1 b 3 n ]
According to Equation (23), W j j ( j = n + 1 , , 2 n ) items for the nodes between the ground surface and the aquifer bottom are modified by
W j j = w j j + ( A j K j , j + n ¯ B j , j + n + A j K j , j n ¯ B j , j n )
W j j n = w j j n A j K j , j n ¯ B j , j n
W j n + j = w j n + j A j K j , j + n ¯ B j , j + n
According to Equation (24), W j j ( j = 2 n + 1 , , 3 n ) items for the ground surface nodes are modified by
W j j = w j j + A j K j , j n ¯ B j , j n
W j j n = w j j n A j K j , j n ¯ B j , j n
According to Equation (25), W j j ( j = 1 , , n ) items for the bottom nodes are modified by
W j j = w j j + A j K j , j + n ¯ B j , j + n
W j n + j = w j n + j A j K j , j + n ¯ B j , j + n

3. 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:
H 1 ( x , y , z , t ) = φ ( x , y , z , t ) ( x , y , z ) Γ D
where φ is the prescribed head in the boundary [L].
(2)
Second-kind boundary condition:
K H n = q ( x , y , z , t ) ( x , y , z ) Γ N
where q is the prescribed flow flux in the boundary [L2 T−1].
(3)
Pumping well:
Q p , i = L i K i x ( L i K i x ) Q w β I , Q p , i = L i K i x ( L i K i x ) Q w ( 1 β I )
where Qp,i and Qp,i’ are the pumping rate of the nodes i and i’, respectively [L3 T−1]; Li is the length of the pumping well filtration in the i-th layer [L]; and Qw is the total pumping rate [L3 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 uptake
The sink term of root uptake, S, is calculated by [28] as
S = α ( h ) S p
where α(h) is a pressure response function of root uptake [–], and Sp 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 E0 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 (ET0) and then the model will calculate the evapotranspiration rate by subsequently multiplying ET0 with the crop coefficient.

4. 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.

5. 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.

5.1. Model Code Verification Cases

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 surface. The pressure head and soil water moisture from the proposed model are compared with those from Hydrus1D, as shown in Figure 4 and Figure 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.

5.1.2. 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, Ks = 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
H 2 = H 1 2 + H 2 2 H 1 2 l x + W K s ( l x x 2 )
where H is the hydraulic head at the location with a distance of x from the left river [L]; H1 and H2 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 Ks is the saturated hydraulic conductivity [L T−1].
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.

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 m3·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, Ks = 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.
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.

5.2. 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 km2. 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 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, considered 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.
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 for the subsequent year 2008. The calibrated soil water parameters of Equations (6) and (7) are listed in Table 1.
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 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.
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.

6. 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.

Acknowledgments

The study is in part supported by the Natural Science Foundation of China through grants 51409192, 41272270, 51479143, 51479144 and 51279141. The first author is also supported by grant (2014M560627) of the China Postdoctoral Science Foundation. We also would like to thank three anonymous reviewers and the editor for their valuable comments that significantly improved the manuscript.

Author Contributions

Yan Zhu derived the numerical method, performed the model code and implemented the validation cases and application case study; Jinzhong Yang conceived the model theory and designed the experiment; Jingwei Wu and Liangsheng Shi contributed to background, field experiment and data analysis input; Lihong Cui was responsible for checking the model code; Ming Ye contributed to the final revision of the manuscript; Yan Zhu and Jinzhong Yang interpreted the results and wrote the paper; and all authors participated in the final review and editing of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kinzelbach, W.; Barer, P.; Siegfried, T.; Brunner, P. Sustainable groundwater management-Problems and scientific tools. Episodes 2003, 26, 279–284. [Google Scholar]
  2. Petheram, C.; Dawes, W.; Grayson, R.; Bradford, A.; Walker, G. A sub-grid representation of groundwater discharge using a one-dimensional groundwater model. Hydrol. Process. 2003, 17, 2279–2295. [Google Scholar] [CrossRef]
  3. Zhu, Y.; Shi, L.S.; Yang, J.Z.; Wu, J.W.; Mao, D.Q. Coupling methodology and application of a fully integrated model for contaminant transport in the subsurface system. J. Hydrol. 2013, 501, 56–72. [Google Scholar] [CrossRef]
  4. Brunner, P. Water and Salt Management in the Yanqi Basin, China; Eidgenössische Technische Hochschule ETH: Züric, Switzerland, 2005. [Google Scholar]
  5. Pikul, M.F.; Street, R.L.; Remson, I. A numerical model based on coupled one-dimensional Richards and Boussinesq equations. Water Resour. Res. 1974, 10, 295–302. [Google Scholar] [CrossRef]
  6. Van Dam, J.C. Field-Scale Water Flow and Solute Transport. SWAP Model Concepts, Parameter Estimation, and Case Studies. Ph.D. Dissertation, Wageningen University, Wageningen, The Netherlands, 2000; p. 167. [Google Scholar]
  7. Niswonger, R.G.; Prudic, D.E.; Regan, R.S. 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; U.S. Geological Survey: Reston, VA, USA, 2006; p. 62.
  8. Twarakavi, N.K.C.; Šimůnek, J.; Sophia, S. Evaluating interactions between groundwater and vadose zone using the HYDRUS-based flow package for MODFLOW. Vadose Zone J. 2008, 7, 757–768. [Google Scholar] [CrossRef]
  9. McDonald, M.G.; Harbaugh, A.W. A Modular Three-Dimensional Finite-Difference Groundwater Flow Model; Techniques of Water-Resources Investigations of the United States Geological Survey: Reston, VA, USA, 1988.
  10. Diersch, H.J.G.; Kolditz, O. Coupled groundwater flow and transport: 2. Thermohaline and 3D convection systems. Adv. Water Resour. 1998, 21, 401–425. [Google Scholar] [CrossRef]
  11. Spitz, F.J.; Nicholson, R.S.; Daryll, A.P. A nested rediscretization method to improve pathline resolution by eliminating weak sinks representing wells. Ground Water 2001, 39, 778–785. [Google Scholar] [CrossRef] [PubMed]
  12. Mehl, S.; Hill, M.C. Three-dimensional local grid refinement method for block-centered finite-difference groundwater models using iteratively coupled shared nodes: A new method of interpolation and analysis of errors. Adv. Water Resour. 2004, 27, 899–912. [Google Scholar] [CrossRef]
  13. Di Giammarco, P.; Todini, E.; Lamberti, P. A Conservative finite elements approach to overland flow: The control volume finite element formulation. J. Hydrol. 1996, 175, 267–291. [Google Scholar] [CrossRef]
  14. Zhu, Y.; Shi, L.S.; Lin, L.; Yang, J.Z.; Ye, M. A Fully Coupled Numerical Modeling for Regional Unsaturated-Saturated Water Flow. J. Hydrol. 2012, 475, 188–203. [Google Scholar] [CrossRef]
  15. Sophocleus, M.; Perkins, S.P. Methodology and application of combined watershed and groundwater models in Kansas. J. Hydrol. 2000, 236, 185–201. [Google Scholar] [CrossRef]
  16. SWIM—Soil and Water Integrated Model. Available online: https://www.pik-potsdam.de/research/climate-impacts-and-vulnerabilities/models/swim (accessed on 20 April 2016).
  17. Facchi, A.; Ortuani, B.; Maggi, D.; Gandolfi, C. Coupled SVAT–groundwater model for water resources simulation in irrigated alluvial plains. Environ. Model. Softw. 2004, 19, 1053–1063. [Google Scholar] [CrossRef]
  18. Yakirevich, A.; Weisbrod, N.; Kuznetsov, M.; Villarreyes, C.A.R.; Benavent, I.; Chavez, A.M.; Ferrando, D. Modeling the impact of solute recycling on groundwater salinization under irrigated lands: A study of the Alto Piura aquifer, Peru. J. Hydrol. 2013, 482, 25–39. [Google Scholar] [CrossRef]
  19. Lardner, R.; Cekirge, H. A new algorithm for three-dimensional tidal and storm surge computations. Appl. Math. Model. 1988, 12, 471–481. [Google Scholar] [CrossRef]
  20. Tsai, T.L.; Yang, J.C.; Wu, C.H. Layered and regional land subsidence model. J. Chin. Inst. Civ. Hydraul. Eng. 1998, 10, 617–626. [Google Scholar]
  21. Drønen, N.; Deiggard, R. Quasi-three-dimensional modelling of the morphology of longshore bars. Coast. Eng. 2007, 54, 197–215. [Google Scholar] [CrossRef]
  22. Fernando, P.T.; Pan, S. Modelling wave of hydrodynamics around a scheme of detached leaky breakwaters. In Proceeding of the 29th International Conference on Coastal Engineering, World Scientific, Lisbon, Portugal, 19–24 September 2004; pp. 830–841.
  23. Li, M.; Fernando, P.T.; Pan, S.Q.; O’Connor, B.A.; Chen, D.Y. Development of a quasi-3d numerical model for sediment transport prediction in the coastal region. J. Hydro-Environ. Res. 2007, 1, 143–156. [Google Scholar] [CrossRef]
  24. Lin, L.; Yang, J.Z.; Zhang, B.; Zhu, Y. A simplified numerical model of 3-D groundwater and solute transport at large scale area. J. Hydrodyn. 2010, 22, 319–328. [Google Scholar] [CrossRef]
  25. Paulus, R.; Dewals, B.J.; Erpicum, S.; Pirotton, M.; Archambeau, P. Innovative modelling of 3D unsaturated flow in porous media by coupling independent models for vertical and lateral flows. J. Comput. Appl. Math. 2013, 246, 38–51. [Google Scholar] [CrossRef]
  26. Van Genuchten, M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Am. J. 1980, 44, 892–898. [Google Scholar] [CrossRef]
  27. Zhang, W. Unsteady Groundwater Flow Calculation and Groundwater Resources Evaluation; Science Press: Beijing, China, 1983. [Google Scholar]
  28. Feddes, R.A.; Kowalik, P.J.; Zaradny, H. Simulation of Field Water Use and Crop Yield; John Wiley & Sons: New York, NY, USA, 1978. [Google Scholar]
  29. Doorenbos, J.; Pruitt, W.O. Guidelines for Predicting Crop Water Requirements, 24, 2nd Edition; FAO Irrigation and Drainage Paper: Rome, Italy, 1977. [Google Scholar]
  30. Mendoza, C.A.; Therrien, R.; Sudicky, E.A. ORTHOFEM User’s Guide Version 1.02. Waterloo Center for Groundwater Research; University of Waterloo: Waterloo, ON, Canada, 1991. [Google Scholar]
  31. Vogel, T.; Huang, K.; Zhang, R. The Hydrus Code for Simulating One-Dimensional Water Flow, Solute Transport, and Heat Movement; U.S. Salinity Laboratory Agriculture Research Service, U.S. Department of Agriculture Riverside: California, CA, USA, 1996.
  32. Šimůnek, J.; Vogel, T.; van Genuchten, M.T. The SWMS_2D Code for Simulating Water Flow and Solute Transport in Two-Dimensional Variably Saturated Media, Version 1.1, Research Report No. 126; U.S. Salinity Laboratory, USDA, ARS: Riverside, CA, USA, 1992.
  33. Xu, X.; Huang, G.H.; Zhan, H.B.; Qu, Z.Y.; Huang, Q.Z. Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas. J. Hydrol. 2012, 412–413, 170–181. [Google Scholar] [CrossRef]
  34. Li, H.T.; Brunner, P.; Kinzelbach, W.; Li, W.P.; Dong, X.G. Calibration of a groundwater model using pattern information from remote sensing data. J. Hydrol. 2009, 377, 120–130. [Google Scholar] [CrossRef]
  35. Li, H.T.; Kinzelbach, W.; Brunner, P.; Li, W.P.; Dong, X.G. Topography representation methods for improving evaporation simulation in groundwater modeling. J. Hydrol. 2008, 356, 199–208. [Google Scholar] [CrossRef]
  36. Brunner, P.; Doherty, J.; Simmons, C.T. Uncertainty assessment and implications for data acquisition in support of integrated hydrologic models. Water Resour. Res. 2012, 48. [Google Scholar] [CrossRef]
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. The subsurface system and the triangular prism element as well as the horizontal plane of average head gradient in the triangular prism element.
Water 08 00195 g001
Figure 2. The water balance volume of node i. Qi means the lateral flux shown in Equation (19).
Figure 2. The water balance volume of node i. Qi means the lateral flux shown in Equation (19).
Water 08 00195 g002
Figure 3. The flow chart of the model.
Figure 3. The flow chart of the model.
Water 08 00195 g003
Figure 4. Pressure head profiles from the proposed model and Hydrus1D.
Figure 4. Pressure head profiles from the proposed model and Hydrus1D.
Water 08 00195 g004
Figure 5. Soil moisture profiles from the proposed model and Hydrus1D.
Figure 5. Soil moisture profiles from the proposed model and Hydrus1D.
Water 08 00195 g005
Figure 6. Steady-state water tables obtained with the proposed model, the analytical solution and SWMS2D.
Figure 6. Steady-state water tables obtained with the proposed model, the analytical solution and SWMS2D.
Water 08 00195 g006
Figure 7. Steady-state pressure head profiles obtained with the proposed model and SWMS2D.
Figure 7. Steady-state pressure head profiles obtained with the proposed model and SWMS2D.
Water 08 00195 g007
Figure 8. Head contours from the proposed model and FEFLOW.
Figure 8. Head contours from the proposed model and FEFLOW.
Water 08 00195 g008
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. 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.
Water 08 00195 g009
Figure 10. Upper boundary conditions applied to the various sub-areas of the domain in (a) 2007 and (b) 2008.
Figure 10. Upper boundary conditions applied to the various sub-areas of the domain in (a) 2007 and (b) 2008.
Water 08 00195 g010
Figure 11. Comparison of observed and simulated groundwater levels for the 10 observation wells in 2007.
Figure 11. Comparison of observed and simulated groundwater levels for the 10 observation wells in 2007.
Water 08 00195 g011
Figure 12. Simulated and measured water table changes over time for the four selected wells, (a) Well 1; (b) Well 3; (c) Well 5; (d) Well 7, in 2007.
Figure 12. Simulated and measured water table changes over time for the four selected wells, (a) Well 1; (b) Well 3; (c) Well 5; (d) Well 7, in 2007.
Water 08 00195 g012
Figure 13. Predicted groundwater table levels in 2008 compared with the measurements at (a) well 1; (b) well 3; (c) well 5; and (d) well 7.
Figure 13. Predicted groundwater table levels in 2008 compared with the measurements at (a) well 1; (b) well 3; (c) well 5; and (d) well 7.
Water 08 00195 g013aWater 08 00195 g013b
Figure 14. Simulated water table contours at times (a) t = 6 days; (b) t = 180 days.
Figure 14. Simulated water table contours at times (a) t = 6 days; (b) t = 180 days.
Water 08 00195 g014
Table 1. Calibrated soil water parameters of the Yonglian irrigation area.
Table 1. Calibrated soil water parameters of the Yonglian irrigation area.
Depth (m)θr (−)θs (−)α (m)n (−)Ks (m d−1)θa (−)θm (−)θk (−)
0–70.020.432.11.611.20.020.430.43
7–530.010.422.11.615.20.010.420.42

Share and Cite

MDPI and ACS Style

Zhu, Y.; Shi, L.; Wu, J.; Ye, M.; Cui, L.; Yang, J. Regional Quasi-Three-Dimensional Unsaturated-Saturated Water Flow Model Based on a Vertical-Horizontal Splitting Concept. Water 2016, 8, 195. https://doi.org/10.3390/w8050195

AMA Style

Zhu Y, Shi L, Wu J, Ye M, Cui L, Yang J. Regional Quasi-Three-Dimensional Unsaturated-Saturated Water Flow Model Based on a Vertical-Horizontal Splitting Concept. Water. 2016; 8(5):195. https://doi.org/10.3390/w8050195

Chicago/Turabian Style

Zhu, Yan, Liangsheng Shi, Jingwei Wu, Ming Ye, Lihong Cui, and Jinzhong Yang. 2016. "Regional Quasi-Three-Dimensional Unsaturated-Saturated Water Flow Model Based on a Vertical-Horizontal Splitting Concept" Water 8, no. 5: 195. https://doi.org/10.3390/w8050195

APA Style

Zhu, Y., Shi, L., Wu, J., Ye, M., Cui, L., & Yang, J. (2016). Regional Quasi-Three-Dimensional Unsaturated-Saturated Water Flow Model Based on a Vertical-Horizontal Splitting Concept. Water, 8(5), 195. https://doi.org/10.3390/w8050195

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop