Journal cover Journal topic
Hydrology and Earth System Sciences An interactive open-access journal of the European Geosciences Union
Journal topic
Hydrol. Earth Syst. Sci., 23, 3481–3502, 2019
https://doi.org/10.5194/hess-23-3481-2019
Hydrol. Earth Syst. Sci., 23, 3481–3502, 2019
https://doi.org/10.5194/hess-23-3481-2019

Research article 28 Aug 2019

Research article | 28 Aug 2019

# A comprehensive quasi-3-D model for regional-scale unsaturated–saturated water flow

A comprehensive quasi-3-D model for regional-scale unsaturated–saturated water flow
Wei Mao1, Yan Zhu1, Heng Dai2, Ming Ye3, Jinzhong Yang1, and Jingwei Wu1 Wei Mao et al.
• 1State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan, Hubei 430072, China
• 2Institute of Groundwater and Earth Sciences, Jinan University, Guangzhou, Guangdong 32306, China
• 3Department of Earth, Ocean, and Atmospheric Science, Florida State University, Tallahassee, FL 32306, USA

Correspondence: Yan Zhu (zyan0701@163.com)

Abstract

For computationally efficient modeling of unsaturated–saturated flow in regional scales, the quasi-three-dimensional (3-D) scheme that considers one-dimensional (1-D) soil water flow and 3-D groundwater flow is an alternative method. However, it is still practically challenging for regional-scale problems due to the highly nonlinear and intensive input data needed for soil water modeling and the reliability of the coupling scheme. This study developed a new quasi-3-D model coupled to the UBMOD 1-D soil water balance model with the MODFLOW 3-D hydrodynamic model. A new implementation method of the iterative scheme was developed in which the vertical net recharge and unsaturated zone depth were used as the exchange information. A modeling framework was developed to organize the coupling scheme of the soil water model and the groundwater model and to handle the pre- and post-processing information. The strength and weakness of the coupled model were evaluated by using two published studies. The comparison results show that the coupled model is satisfactory in terms of computational accuracy and mass balance error. The influences of spatial and temporal discretization as well as the stress period on the model accuracy were discussed. Additionally, the coupled model was used to evaluate groundwater recharge in a real-world study. The measured groundwater table and soil water content were used to calibrate the model parameters, and the groundwater recharge data from a 2-year tracer experiment were used to evaluate the recharge estimation. The field application further shows the practicability of the model. The developed model and the modeling framework provide a convenient and flexible tool for evaluating unsaturated–saturated flow systems at the regional scale.

1 Introduction

While groundwater resources are important for domestic, agricultural, and industrial uses, groundwater is vulnerable due to over-exploitation, climate change, and biochemical pollution (Bouwer, 2000; Sophocleous, 2005; Evans and Sadler, 2008; Karandish et al., 2015; Zhang et al., 2018). To protect or exploit groundwater resources, understanding soil water flow systems is necessary as soil water is the major source of groundwater recharge and destination of phreatic consumption (Yang et al., 2016; Wang et al., 2017). The Richards equation is usually used to describe the soil water flow and groundwater flow. Many numerical schemes have been developed to solve the three-dimensional (3-D) Richards equation (Weill et al., 2009) in computer codes, such as HYDRUS (Šimůnek et al., 2012), FEFLOW (Diersch, 2013), HydroGeoSphere (Brunner and Simmons, 2012), InHM (VanderKwaak and Loague, 2001), and MODHMS (Tian et al., 2015). These fully 3-D models have a solid theoretical foundation and have been used for regional-scale unsaturated–saturated water flow simulation. However, since the soil water flow is highly nonlinear in nature and sensitive to atmospheric changes, soil utilizations, and human activities, the numerical schemes require use of fine discretization in vertical space and time for accurate numerical solutions (Downer and Ogden, 2004; Varado et al., 2006). This makes the numerical solutions computationally expensive, especially for large-scale modeling (Van Walsum and Groenendijk, 2008; Shen and Phanikumar, 2010; Yang et al., 2016; Szymkiewicz et al., 2018). There are also many conceptual unsaturated–saturated water flow models, e.g., SWAT (Arnold et al., 2012), INFIL 3.0 (Fill, 2008), HSPF (Duda et al., 2012), and SALTMOD (Oosterbaan, 1998), which show advantages in mass balance and computational cost. However, these models usually adopt many empirical equations which result in poor performance compared with the fully 3-D numerical models.

To address the computational challenges discussed above, a variety of simplifications have been introduced for the soil water flow for regional-scale problems. One simplification is to treat the hydrological processes (e.g., infiltration, evapotranspiration, and deep percolation) occurring in the unsaturated zone as one-dimensional (1-D) processes in the vertical direction. Field experiments at the regional scale also show that, in the unsaturated zone, the lateral hydraulic gradient is usually significantly smaller than the vertical gradient (Sherlock et al., 2002). This 1-D simplification leads to the quasi-3-D scheme, which ignores the lateral flow in the unsaturated zone but considers groundwater flow as a 3-D problem. The quasi-3-D scheme avoids solving the 3-D Richards equation for the unsaturated zone and thus improves computational efficiency and model stability. The quasi-3-D scheme is an efficient solution for large-scale unsaturated–saturated flow modeling (Twarakavi, et al., 2008; Yang et al., 2016) and is popular among groundwater modelers (Havard et al., 1995; Harter and Hopmans, 2004; Graham and Butts, 2005; Stoppelenburg et al., 2005; Seo et al., 2007; Markstrom et al., 2008; Ranatunga et al., 2008; Kuznetsov et al., 2012; Xu et al., 2012; Zhu et al., 2012; Leterme et al., 2015). However, it is still challenging when using the quasi-3-D models for a practical regional-scale problem. Two concerns arise as follows.

The first concern is the unsaturated modeling method. Although the quasi-3-D scheme is computationally efficient, the numerical solutions of the 1-D Richards equation still require intensive input data and face numerical instability and mass balance errors under some specific situations (Zha et al., 2017). These problems limit their practical application for simulating regional-scale problems under complicated geological and climate conditions as well as anthropogenic activities. As an alternative to the numerical solutions of the 1-D Richards equation, water balance models have been used to describe soil water movements, which not only reduce the amount of input data, but also improve computational efficiency and stability. The water balance models can be coupled with groundwater models. Facchi et al. (2004) coupled a SVAT conceptual soil water movement model with MODFLOW to simulate the hydrologically relevant processes in the alluvial irrigated plains. Kim et al. (2008) integrated SWAT with MODFLOW to describe the exchange between hydrologic response units in the SWAT model and MODFLOW cells. The traditional water balance models however may oversimplify soil water movement, and thus cannot accurately represent certain important features of soil water flow, e.g., the upward flux and soil heterogeneity. To extend the application of the water balance model for more complicated conditions, Mao et al. (2018) developed a soil water balance model (called the UBMOD model), which can simulate both upward and downward soil water movement in a heterogeneous situation. And the model can be used with a coarse discretization in space and time, all of which make it suitable for the large-scale modeling.

Another concern is the scheme when coupling saturated models with unsaturated models. There are three different numerical coupling schemes categorized by Furman (2008): uncoupled, iterative coupled, and fully coupled. The uncoupled scheme is widely used when using soil water flow packages with MODFLOW, such as LINKFLOW (Havard et al., 1995), SVAT-MODFLOW (Facchi et al., 2004), UZF1-MODFLOW (Niswonger et al., 2006), HYDRUS-MODFLOW (Seo et al., 2007), and SWAP-MODFLOW (Xu et al., 2012). While this scheme is easy to implement, its results may not be reliable when recharge from the unsaturated zone causes substantial changes to the water table. Additionally, this scheme may result in the mass balance error (Shen and Phanikumar, 2010; Kuznetsov et al., 2012). The fully coupled scheme is mathematically and computationally rigorous, because it solves unsaturated and saturated flows simultaneously with internal boundary conditions of the two flows (Zhu et al., 2012). However, the fully coupled scheme is computationally expensive (Furman, 2008). The iterative coupled scheme offers a trade-off between model accuracy and computational cost (Yakirevich et al., 1998; Liang et al., 2003). And it has been widely used to couple two hydrodynamic models, both of which calculate the hydraulic head and use the hydraulic head as the exchange information (Stoppelenburg et al., 2005; Kuznetsov et al., 2012). However, the soil water content is the variable calculated by soil water balance models other than the hydraulic head. Therefore, the traditional implementation method of the iterative scheme is inapplicable, and a specific implementation method of the iterative scheme should be developed to couple the soil water balance model and the hydrodynamic groundwater model.

In this study, a new quasi-3-D model is developed. The UBMOD 1-D water balance model developed by Mao et al. (2018) is integrated with MODFLOW (Harbaugh, 2005). A new implementation method of the iterative scheme is established for numerical solutions, and the net groundwater recharge and the depth of the unsaturated zone (which is equal to the groundwater table depth) are chosen as the exchange information. The coupled model can achieve mass balance and keep numerical stability well, and it is suitable for large-scale modeling based on the characteristics of MODFLOW and UBMOD. Moreover, instead of developing a new package for MODFLOW, a framework of organizing the modeling procedures is developed. This paper elaborates the methodology of coupling the unsaturated and saturated water flow and the modeling framework in Sect. 2. Two published studies are used to test the performance of the coupled model when handling different water flow conditions in Sect. 3. A real-world application to study the regional net groundwater recharge is presented in Sect. 4.

2 Methodology and model development

In the new coupled model, the unsaturated–saturated domain is partitioned into a number of sub-areas in the horizontal direction mainly according to the spatially distributed inputs (soil types, atmosphere boundary conditions, land usage types, and crop types). A 1-D soil column is used to characterize the average soil water flow in each sub-area, and UBMOD is used to simulate the 1-D soil water flow. MODFLOW is used to simulate the 3-D groundwater flow of the whole domain. It is assumed that the flow in the unsaturated zone is in the vertical direction and that there is only vertical exchange flux between the unsaturated and saturated zones. It is further assumed that using the vertical column can reasonably simulate the unsaturated flow in each sub-area while ignoring the horizontal heterogeneity. In this section, UBMOD is first presented, followed by a brief introduction to MODFLOW and two peripheral tools (FloPy, Bakker et al., 2016, and ArcPy, Toms, 2015) used in the model. The procedures of the new model and the modeling framework are described in Sect. 2.3, and the specific implementation method of the unsaturated and saturated coupling scheme is described in Sect. 2.4.

## 2.1 The UBMOD soil water balance model

This section describes the UBMOD soil water balance model to make this paper self-contained, and more details of UBMOD are referred to in Mao et al. (2018) or the Appendix. UBMOD is a water balance model based on a hybrid of numerical and statistical methods. The model can effectively and efficiently simulate both downward and upward soil water movement with only four physically meaningful parameters, which makes it suitable for practical application.

There are four major components to describe the soil water movement in UBMOD. Firstly, the vertical soil column is divided into a cascade of “buckets” and each “bucket” corresponds to a soil layer. The “buckets” will be filled to saturation from the top layer to the bottom layer if there is infiltration, which is referred to as the allocation of infiltration water. Specifically speaking, the infiltration water first fills the top “bucket”, and then the excessive infiltration water moves downward to the next “bucket”, until all the infiltration water is allocated in the “buckets”. The governing equation of layer i is

$\begin{array}{}\text{(1)}& {q}_{i}=min\left({M}_{i}×\left({\mathit{\theta }}_{\mathrm{s},i}-{\mathit{\theta }}_{i}\right),I-{I}_{d,i-\mathrm{1}}\right),\end{array}$

where i indicates the vertical soil layer, i=1, …, j; qi is the amount of allocated water per unit area of layer i (L); Mi is the thickness of layer i (L); θi is the initial soil water content of layer i (L3 L−3); θs,i is the saturated soil water content of layer i (L3 L−3); I is the quantity of infiltration rate (L); ${I}_{d,i-\mathrm{1}}$ is the consumed infiltration water per unit area by all upper layers above layer i (L). The infiltration rate I is an input data in the model, and the partitioning of rainfall between infiltration and runoff has not been considered by now.

Secondly, when the soil water content exceeds the field capacity, the soil water will move downward driven by the gravitational potential. The governing equation is

$\begin{array}{}\text{(2)}& \frac{\partial \mathit{\theta }}{\partial t}=-\frac{\partial K\left(\mathit{\theta }\right)}{\partial z},\end{array}$

where t is the time (T); z is the elevation in the vertical direction (L). The vertical coordinate is positive downward. K(θ) is the unsaturated hydraulic conductivity (L T−1) as a function of soil water content, which is characterized by empirical formulas referred to as drainage functions. The commonly used equations can be found in Mao et al. (2018) and the Appendix.

Figure 1(a) The procedures of geographic input information preparation. (b) The spatial scheme of the coupled model.

Thirdly, the source/sink terms are used to account for soil evaporation and crop transpiration. The governing equation is as follows:

$\begin{array}{}\text{(3)}& \frac{\partial \mathit{\theta }}{\partial t}=-W,\end{array}$

where W is the source/sink term (T−1). The Penman–Monteith formula and Beer's law (also known as a Ritchie-type equation) are adopted to estimate the potential soil evaporation Ep and potential crop transpiration Tp. Then Ep and Tp are distributed to each layer based on the evaporation cumulative distribution function and the root density function. The actual soil evaporation and crop transpiration are obtained by discounting Ep and Tp with the soil water stress coefficient.

Lastly, we calculate the diffusive movement driven by the matric potential. The governing equation is

$\begin{array}{}\text{(4)}& \frac{\partial \mathit{\theta }}{\partial t}=\frac{\partial }{\partial z}\left(D\left(\mathit{\theta }\right)\frac{\partial \mathit{\theta }}{\partial z}\right),\end{array}$

where D(θ) is the hydraulic diffusivity (L2 T−1). The finite-difference method is used to solve the equation. An empirical formula with four parameters (saturated hydraulic conductivity Ks, saturated water content θs, field capacity θf, and residual water content θr) is used to describe the hydraulic diffusivity D(θ). The heterogeneity of soils is also taken into account by adding a correction item on the right-hand side, which makes the model applicable to heterogeneous situations. With the help of the diffusive term, UBMOD can consider upward soil water movement, which is ignored by most water balance models. The details of D(θ) are shown in the Appendix.

The original UBMOD is a soil water balance model, which cannot consider groundwater table. For the purpose of saturated–unsaturated coupling, the model has been improved to calculate the groundwater recharge, which is expatiated in Sect. 2.4.

## 2.2 Brief introduction to MODFLOW and two peripheral tools

MODFLOW is a computer program that numerically solves the 3-D groundwater flow equation for a porous medium using a block-centered finite-difference method (Harbaugh, 2005). The governing equation solved by MODFLOW is

$\begin{array}{}\text{(5)}& \frac{\partial }{\partial {x}_{i}}\left({K}_{ij}\frac{\partial H}{\partial {x}_{j}}\right)+W={S}_{\mathrm{s}}\frac{\partial H}{\partial t},\end{array}$

where $i,j=\mathrm{1}$–3 indicate the x, y, and z directions, respectively; Kij is the saturated hydraulic conductivity (L T−1); H is the hydraulic head (L); W is the volumetric flux per unit volume representing sources and/or sinks of water (L3 T−1); Ss is the specific storage of the porous material (L−1); and t is the time (T).

Figure 2(a) The spatial coupling scheme for one saturated cell and one unsaturated soil column. (b) The temporal coupling scheme and the relationship between the stress period (ΔT) and the time steps for UBMOD (Δtu) and MODFLOW (Δts). (c) The specific implementation of the iterative coupling scheme from t to tT. Note: zi (i=1, …, j) is the vertical elevation of layer i; du is the thickness of the unsaturated zone; Mi is the thickness of layer i; R is the groundwater recharge for one cell; qI, qA, qS, and qD are the fluxes across the water table caused by allocation of the infiltration water, the advective movement, source/sink terms, and the water diffusion per unit area, respectively; Du is the thickness of the unsaturated zone (l dimension); R is the vertical net recharge for the regional scale (m×n dimension); t is the time; p is the iteration level; H is the saturated hydraulic head (m×n dimension); εH is a user-specified tolerance.

FloPy and ArcPy are the two peripheral tools used in the model development. FloPy is a Python package for creating, running, and post-processing MODFLOW-based models. Unlike the common graphical user interfaces (GUIs) method, FloPy facilitates users to write a Python script to construct and post-process MODFLOW models, and it has been shown as a convenient and powerful tool by Bakker et al. (2016). The geographic information system (GIS) is a helpful tool for groundwater modeling by providing geospatial database and results presentation (Xu et al., 2011; Lachaal et al., 2012). ArcPy is an application program interface (API) of ArcGIS for Python (Toms, 2015), which provides a useful and productive way to perform geographic data analysis, data conversion, data management, and map automation with Python.

## 2.3 The process of geographic input information

The procedures of the modeling framework are composed of three major parts, including the pre-processing, the coupled model, and the post-processing. The preparation of geographic input information of the model shown in Fig. 1a is the major component of pre-processing. The geographic information includes the domain area, boundary conditions, sub-areas, digital elevation model (DEM), hydraulic conductivity, and porosity. The shapefile of the domain area (usually irregular in shape) is first discretized by a regular boundary with both active and inactive cells. The discretized domain can be joined with the shapefile of the boundary condition to generate the “ibound” array of MODFLOW as shown in Fig. 1a, which is used to specify which cells are active, inactive, or fixed head in MODFLOW. The shapefile of sub-areas is joined with the domain file, represented in the sub-area array with different numbers specified as different sub-areas. The raster files of DEM, hydraulic conductivity, and porosity are further joined, and the values of these variables are listed in the arrays shown in Fig. 1a. The unsaturated–saturated flow model coupling scheme will be described in the next section. The result presentations are accomplished by the post-processing tool, which contains a series of utilities developed based on Python packages.

## 2.4 Coupling scheme of UBMOD and MODFLOW

Figure 1b demonstrates the sketch map of the specific implementation method of the unsaturated and saturated coupling schemes. The unsaturated–saturated domain is partitioned into a number of sub-areas in the horizontal direction mainly according to the spatially distributed inputs (each sub-area is considered to be homogeneous in the horizontal). There are l sub-areas and j layers for a specific soil column shown in Fig. 1b. Soil water flow of each sub-area is simulated by using one 1-D soil column. The recharge at the bottom boundary calculated by UBMOD is treated as the upper boundary condition of MODFLOW. The whole saturated zone is discretized into a grid with cells, and there are m row and n column cells of the saturated zone as shown in Fig. 1b. All cells in the same sub-area receive the same recharge from the soil zone calculated by the representative 1-D soil column. In the vertical direction, both the saturated domain and the soil columns are discretized into different layers based on available data and information, and the layer discretization remains unchanged during the simulation. The lower boundary condition of the whole region is set in MODFLOW. As the soil water movement is reduced to 1-D flow, the surrounding boundary conditions for the unsaturated zone are no-flux boundary, while the surrounding boundary conditions for the saturated zone are set in MODFLOW as practical. Note that the saturated zone and the unsaturated zone are independent, but some layers may transform between the saturated zone and the unsaturated zone, which are referred to as the overlap region. Fine vertical discretization of UBMOD in the overlap region is needed to improve the simulation accuracy.

Since the independent variable of UBMOD is the soil water content and the independent variable of MODFLOW is the hydraulic head, this study uses the vertical net recharge and the unsaturated zone depth to couple the unsaturated zone and saturated zone. The domain shown in Fig. 1b is used as an example to illustrate the spatial and temporal coupling methods in the study. The vertical net recharge is represented by matrix R with m×n elements, and the unsaturated zone depth by vector Du with l elements, as illustrated in Fig. 1b. Scalar R is used to denote the specific net recharge of a soil column to the corresponding saturated sub-area, and scalar du denotes the depth of the soil column. Figure 2a shows the spatial coupling method of a soil column connected with the groundwater system. The water table is located in the jth layer. The net recharge R from the soil zone is calculated by UBMOD as follows:

$\begin{array}{}\text{(6)}& R={q}_{\mathrm{I}}+{q}_{\mathrm{A}}+{q}_{\mathrm{S}}+{q}_{\mathrm{D}},\end{array}$

where qI, qA, qS, and qD are the fluxes across the water table caused by allocation of the infiltration water, the advective movement driven by the gravitational potential, source/sink terms, and the water diffusion driven by the matric potential per unit area, respectively (L).

These four terms correspond to the four major components in UBMOD, as described in Sect. 2.1. Specifically, the infiltration water is allocated first according to Eq. (1) if there is precipitation or irrigation. When there is residual infiltration water across the water table in the jth layer, the amount of residual infiltration is denoted as qI. Then the advective flow qA across the water table driven by gravitational potential is calculated by Eq. (2). The direction of these two terms is downward. The qS term is caused by evapotranspiration. When the critical depth of evapotranspiration is shallower than the groundwater table depth, the groundwater can be consumed by evapotranspiration and it causes an upward qS term. A virtual layer is needed when calculating the diffusive movement driven by matric potential across the water table based on Eq. (4). As shown in Fig. 2a, the virtual layer will be added under the water table, numbered as layer j+1. The thickness, Mj+1 (L), of the layer is set as

$\begin{array}{}\text{(7)}& {M}_{j+\mathrm{1}}={z}_{j+\mathrm{1}}-{d}_{\mathrm{u}},\end{array}$

where zj+1 is the bottom depth of layer j+1 (L); du is the thickness of the unsaturated zone (L). The amount of the upward flux between the virtual layer and layer j is denoted as qD. Then, the net recharge matrix R for the whole area is obtained and used for the Recharge (RCH) package of MODFLOW.

The time coupling method is shown in Fig. 2b. There are three levels of time discretization in the coupled model as follows: the stress period ΔT used in MODFLOW, the calculation time step for MODFLOW Δts, and the calculation time step for UBMOD Δtu. The stress time step (ΔT) is also used in the iterative process, and the UBMOD unsaturated model and MODFLOW saturated model exchange information at the end of each stress period. Δtu is a priori value and cannot be changed during the calculation. UBMOD can give acceptable results when Δtu is shorter than 10 d for assumed cases and 1 d for a real-world case (Mao et al., 2018). Δts is set as the technical report described by Harbaugh (2005) and can be changed during the calculation.

The implementation of an iterative coupling scheme is shown in Fig. 2c, which shows the calculation period from t to tT. At the time t, the saturated hydraulic head is known, marked as Ht (m×n dimension). When the model runs from t to tT, firstly, the initial saturated hydraulic head HtT at tT is set to be equal to Ht, and then the average unsaturated depth from t to tT is calculated according to HtT, marked as $\mathbit{D}{\mathbit{u}}^{t+\mathrm{\Delta }T,p}$ (l elements). p is the iteration level. The ${d}_{\mathrm{u}}^{t+\mathrm{\Delta }T,p}$ for one soil column is calculated as follows:

$\begin{array}{}\text{(8)}& {d}_{\mathrm{u}}^{t+\mathrm{\Delta }T,p}=\stackrel{\mathrm{‾}}{D}-\stackrel{\mathrm{‾}}{{H}^{t+\mathrm{\Delta }T}},\end{array}$

where $\stackrel{\mathrm{‾}}{D}$ is the average depth from the soil surface to the impermeable layer of the controlling domain of the soil column (L); $\stackrel{\mathrm{‾}}{{H}^{t+\mathrm{\Delta }T}}$ is the average thickness of the controlling saturated domain of the soil column (L).

Secondly, the model runs UBMOD with the unsaturated time step Δtu to obtain the vertical recharge at each time step (marked as rt) until the time comes to be tT. The total recharge during the stress period ΔT (from t to tT) RΔT can be obtained by summarizing the recharge at each unsaturated time step, as follows:

$\begin{array}{}\text{(9)}& {R}_{\mathrm{\Delta }T}=\sum _{t}^{t+\mathrm{\Delta }T}{r}_{t}.\end{array}$

The average recharge R from t to tT can be obtained by

$\begin{array}{}\text{(10)}& R={R}_{\mathrm{\Delta }T}/\mathrm{\Delta }T.\end{array}$

Then the average recharge from all 1-D soil columns can be obtained, represented as ${\mathbit{R}}^{t+\mathrm{\Delta }T,p}$, which is then used by the MODFLOW RCH package. Subsequently, the model runs the MODFLOW model with the saturated time step Δts to obtain the saturated hydraulic head until the time comes to tT. The hydraulic head at the time tT is marked as ${\mathbf{H}}^{t+\mathrm{\Delta }T,p}$ (m×n dimension). The convergence of the iteration is determined by using the difference of the hydraulic head between the present ${\mathbf{H}}^{t+\mathrm{\Delta }T,p}$ and the initial HtT. The convergence criterion is

$\begin{array}{}\text{(11)}& \mathrm{if}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}max\left(\left|{\mathbf{H}}^{t+\mathrm{\Delta }T,p}-{\mathbf{H}}^{t+\mathrm{\Delta }T}\right|\right)<{\mathit{\epsilon }}_{H},\end{array}$

where εH is a user-specified tolerance (L). If the criterion is met, the iteration stops, and ${\mathbf{H}}^{t+\mathrm{\Delta }T,p}$ is the convergent results at time tT, and the model proceeds to the next stress period. Otherwise, the iteration continues to p+1 and ${\mathbf{H}}^{t+\mathrm{\Delta }T,p}$ will be used to calculate the average unsaturated depth shown in Eq. (8). The above procedures will be repeated until the convergence criterion of Eq. (11) is met.

Table 1The hydraulic parameters of case 1 and case 2.

Note: θr is the residual water content; θs is the saturated water content; Ks is the saturated hydraulic conductivity; α and n are parameters depending on the pore size distribution; θf is the field capacity and μ is the specific yield.

Figure 3The values of actual precipitation and potential transpiration rates of case 1.

3 Model evaluation

In this section, two test cases were designed to evaluate the model accuracy and the performance of the numerical coupling scheme under complicated soil and boundary conditions. The simulation results were compared with numerical results obtained using HYDRUS-1D (Šimůnek et al., 2008) and SWMS2D (Šimůnek et al., 1994), and with published experimental data. For these cases, the mean absolute relative error (ARE), the root mean squared error (RMSE), the index of agreement (IA), and the determination coefficient (R2) are used to quantitatively evaluate the misfit between the simulated results of the coupled model and reference values, which are calculated as

$\begin{array}{}\text{(12)}& \mathrm{ARE}=\frac{\mathrm{1}}{x}\sum _{i=\mathrm{1}}^{x}\frac{\left|{y}_{i}-{Y}_{i}\right|}{{Y}_{i}}×\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathit{%},\text{(13)}& \mathrm{RMSE}=\sqrt{\frac{\mathrm{1}}{x}\sum _{i=\mathrm{1}}^{x}{\left({Y}_{i}-{y}_{i}\right)}^{\mathrm{2}}},\text{(14)}& \mathrm{IA}=\mathrm{1}-\frac{\sum _{i=\mathrm{1}}^{x}{\left({y}_{i}-{Y}_{i}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^{x}{\left[\left|{y}_{i}-\stackrel{\mathrm{‾}}{y}\right|+\left|{Y}_{i}-\stackrel{\mathrm{‾}}{Y}\right|\right]}^{\mathrm{2}}},\text{(15)}& {R}^{\mathrm{2}}=\mathrm{1}-\frac{\sum _{i=\mathrm{1}}^{x}{\left({y}_{i}-{Y}_{i}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^{x}{\left({Y}_{i}-\stackrel{\mathrm{‾}}{Y}\right)}^{\mathrm{2}}},\end{array}$

where the subscript i represents the serial number of the results; x represents the total number of the results; yi is the simulated result of the coupled model and Yi is the reference result; $\stackrel{\mathrm{‾}}{y}$ is the average simulated result and $\stackrel{\mathrm{‾}}{Y}$ is the average reference result.

## 3.1 Two test cases

### 3.1.1 Case 1: 1-D upward flux with atmospheric condition

This case was to test the performance of the coupling scheme explained in Sect. 2.4. The case simulated a single field soil profile of the Hupselse Beek watershed in the Netherlands, which was used as a demo in the HYDRUS-1D technical manual (Šimůnek et al., 2008). The soil profile consists of a 0.4 m thick upper layer and a 1.9 m thick bottom layer. The depth of the root zone is 0.3 m. The hydraulic parameters of the two soil layers are presented in Table 1. The surface boundary condition involves actual precipitation and potential transpiration rate as shown in Fig. 3. The groundwater level was initially set at 0.55 m below the soil surface. Only one vertical soil column and one MODFLOW cell were used in the coupled model. The parameters used in the coupled model are also listed in Table 1. The results from HYDRUS-1D were used as the reference of this test case. The stress period ΔT was set as 5 d, and the MODFLOW time step Δts and the UBMOD time step Δtu were both set to be 1 d. The spatial discretization was 0.1 m.

Figure 4(a) The sketch of the 2-D recharge experiment of case 2. (b) The comparison of the water table between simulated results by the coupled model, SWMS2D, and observation data of case 2.

To figure out the influence of the temporal and spatial discretization as well as the stress period on the simulation results, scenarios with different temporal and spatial resolutions and the stress period of the coupled model were performed. Scenario 1 was set as the same as the above case. The UBMOD time steps of scenario 2 and scenario 3 were 0.5 and 2 d, while other inputs were the same as scenario 1. The spatial discretizations of scenario 4 and scenario 5 were set as 0.05 and 0.2 m, while other inputs were the same as scenario 1. The stress periods of scenario 6, scenario 7, and scenario 8 were set as 8, 10, and 15 d, while other inputs were the same as scenario 1. The eight scenarios were marked as S1–S8.

### 3.1.2 Case 2: two-dimensional (2-D) water table recharge experiment

This test case was used for model validation in a 2-D unsaturated–saturated flow system. The purpose of the case is to discuss the performance of the model under the condition with large lateral flux in the unsaturated zone. The numerical simulation of our model was compared with the data of a 2-D water table recharge experiment conducted by Vauclin et al. (1979). The experimental data have been used to test the variably saturated flow models (Clement et al., 1994) and coupled unsaturated–saturated flow models (Thoms et al., 2006; Twarakavi et al., 2008; Shen and Phanikumar, 2010; Xu et al., 2012). The 2-D domain is a rectangular sandy soil slab of $\mathrm{6.0}×\mathrm{2.0}×\mathrm{0.05}$ m. The initial pressure head is 0.65 m at the domain bottom. At the soil surface, a constant flux of q=3.55 m d−1 is applied at the central 1.0 m, and the remaining soil surface is the no-flux boundary. Because of the symmetry of the flow system, only one-half of the domain (right side) with a size of 3.0 m × 2.0 m × 0.05 was simulated. The setup of the simulation is shown in Fig. 4a. No-flow boundaries were defined on the bottom and the left side, and the specified hydraulic head boundary of 0.65 m was set on the right side. The values of soil hydraulic parameters are listed in Table 1. The simulation period is 8 h. In our coupled model, there were 30 uniform rectangular cells used by MODFLOW, and there were 10 sub-areas defined to represent the unsaturated zone, which were numbered from left to right. The first and last sub-areas covered 0.2 and 0.4 m in the x direction, respectively, and each of the remaining sub-areas covered 0.3 m in the x direction. The first and second sub-areas were used to define the recharge boundary, while the other sub-areas were used to define the no-recharge boundary. The stress period ΔT was set as 1 h, and the initial MODFLOW time step Δts and UBMOD time steps Δtu were set as 0.167 h. The spatial discretization of UBMOD was uniformly 0.1 m. The experiment was also simulated by using SWMS2D, which considered the lateral flow. The mean time step of SWMS2D was set to be 0.0225 h, and 20 200 finite elements were used.

Figure 5The comparison of the results calculated by HYDRUS-1D and the coupled model of case 1.

Table 2The statistical index values of the coupled model of case 1.

## 3.2 Results and discussions of model performance

### 3.2.1 Computational accuracy of the coupling scheme

Figure 5 shows the comparison of the results simulated by HYDRUS-1D and the coupled model of case 1. The statistical indexes are listed in Table 2. Figure 5a demonstrates that the groundwater table depth calculated by the coupled model has a similar pattern to that of HYDRUS-1D. The ARE, RMSE, IA, and R2 values were 17.0 %, 0.171 m, 0.976, and 0.977. The soil water contents at the depth of 1.15 m over time from the two models are compared in Fig. 5b. The ARE, RMSE, IA, and R2 were 2.2 %, 0.008 cm3 cm−3, 0.991, and 0.976. The simulated soil water content profiles at different times are shown in Fig. 5c–e and the evaluation indexes demonstrate the satisfactory performance of the model. Moreover, the net groundwater consumption at the end of the simulation period was compared, which is 0.132 m calculated by the coupled model, and it is the same as that from HYDRUS-1D. In general, these results indicate that the coupled model can capture the flow information under the upward flux and the heterogeneous condition.

Table 3The statistical index values of SWMS2D and the coupled model of case 2.

The deviations of groundwater table depth and soil water content from the coupled model and HYDRUS-1D can also be observed in Fig. 5. The deviations are caused by the different model structures of the coupled model and HYDRUS-1D. HYDRUS-1D solves the saturated–unsaturated flow together, and the groundwater table is determined at the depth with the matric potential equaling zero. The soil water content of the capillary fringe above the groundwater table is almost saturated. However, the UBMOD model cannot simulate the capillary fringe. And there is a parameter, the field capacity used to calculate the downward movement of soil water, which is defined under a free drainage condition. So, the coupled model could lead to the lower soil water content in the capillary fringe and higher groundwater table as shown in Fig. 5. And there is another parameter-specific yield used in the coupled model to determine the groundwater table, which is also attributed to the deviation of the groundwater table.

Figure 4b shows the comparison of simulated water tables at four different times using the coupled model and SWMS2D and the observation data in case 2. The index values are listed in Table 3. The coupled model matched the observation data well at the simulation times of 3, 4, and 8 h, with the ARE values smaller than 3 %, the RMSE values smaller than 0.03 m, and the IA and R2 values close to 1. The observed and simulated soil water content profiles for the initial and ending times are presented in Fig. 6. The statistical index values are also listed in Table 3. The simulations by the coupled model agree well with the observations at the locations of x=0.2 m, x=1.4 m, and x=2 m (Fig. 6a, d, and e) where the lateral water flow is negligible. The calculated recharge is 3.55 m d−1 per unit area when the flow becomes steady, which equals the input flux. These results demonstrate the accuracy of the coupled model and the reliability of the coupling scheme shown in Sect. 2.4.

Figure 6Comparison of soil water content profiles between the simulations from the coupled model, SWMS2D, and the observations at different locations: (a) x=0.2 m; (b) x=0.6 m; (c) x=0.8 m; (d) x=1.4 m; (e) x=2 m.

Figure 7The influence of (a) temporal and (b) spatial discretization and (c) stress period on simulation results.

### 3.2.2 Influence of the temporal and spatial discretization as well as the stress period on simulation results

The groundwater table depths calculated by scenarios with different temporal discretizations (S1–S3) are compared with those from HYDRUS-1D in Fig. 7a. The statistical index values are shown in Table 2. It can be found that the water table depths calculated by different scenarios have the same variation trend. The ARE values of the three scenarios are smaller than 20 %, and the maximum RMSE value is 0.171 m. The IA and R2 values are larger than 0.95. The groundwater table depths calculated by scenarios with different spatial discretizations (S1, S4, and S5) are compared with those from HYDRUS-1D in Fig. 7b. The ARE values of the three scenarios are smaller than 25 %. The maximum RMSE value is 0.215 m. The IA and R2 values are larger than 0.95. The water table depths calculated by scenarios with different stress periods (S1, S6, S7, and S8) are compared with those from HYDRUS-1D in Fig. 7c. It should be noted that the model collapsed at the time of 227 d when the stress period is 15 d (S8). The statistical index values for S1, S6, and S7 are shown in Table 2. The ARE and RMSE values of the three scenarios are very similar. Considering the water balance method and empirical formulas adopted in the coupled model, the results calculated by all the scenarios except S8 are acceptable. These results indicate that the temporal and spatial discretizations have slight influence on the modeling results. It should be noted that the impact of the stress period on a certain scale (< 10 d in this case) has no significant impact on the simulation results. However, too large a stress period will cause improper results.

### 3.2.3 Limitations of the coupled model

Although the coupled model had a sufficient computational accuracy as shown above, there were limitations because of the quasi-3-D assumptions. The coupled model overestimates the water table at the time of 2 h in case 2 as shown in Fig. 4b. This is caused by a significant lateral flow in the unsaturated zone during the early period due to the relatively low initial soil water content condition. Therefore, a portion of the infiltration water in the first and second sub-areas should move in the lateral direction, instead of moving downward to the saturated zone as in the quasi-3-D model. The coupled model thus overestimates the recharge flux and results in a higher water table at the early period. Additionally, the simulated soil water content by the coupled model has poor performance at the locations of x=0.6 m and x=0.8 m (Fig. 6b and c). These two sub-areas are close to the recharge zone and affected by the lateral flow, which is ignored in the coupled model. These phenomena are similar to the results calculated by other quasi-3-D models (Xu et al., 2012; Shen and Phanikumar, 2010). Therefore, the coupled model overestimates the recharge and underestimates the soil water content when the lateral flow cannot be ignored. Its application should be limited to cases in which the soil flow mainly occurs in the vertical direction.

Figure 8(a) The geographic location of the Yonglian irrigation area. (b) The land use map. (c) The surface DEM.

### 3.2.4 Water mass balance and computational cost

The mass balance error of the coupled model is small, with the maximum values 0.012 % in case 1 and 0.004 % in case 2, while they are 1.6 % for the HYDRUS-1D model and 0.133 % for the SWMS2D model. The cases were run on a 6 GB RAM, double 2.93 GHz intel Core (TM) 2 Duo CPU-based personal computer. The computational cost of different scenarios in case 1 of the coupled model ranges from 49 to 63 s as listed in Table 2. It is 1.4 s by HYDRUS-1D. The temporal and spatial discretization has slight influence on computational cost, while the stress period has significant influence on the computation cost. The iteration and information exchange are responsible for the high computational cost. For case 2, the computational costs of the coupled model and the SWMS2D model are 46 and 95 s, respectively. The coupled model has a better efficiency compared with the complete 2-D model due to its simpler numerical solutions and coarse discretization in space and time. The advantage of decreasing computational cost will be more obvious when the application scale becomes larger. Generally speaking, the coupled model provides satisfactory mass balance and good computational efficiency.

Figure 9Daily climate data in the Yonglian irrigation area.

Figure 10Comparison between simulated and observed water table depth of the real-world application.

4 Real-world application

## 4.1 Study site and input data

The coupled model was used to calculate the regional-scale groundwater recharge in a real-world case, where the shallow groundwater has significant impact on the soil water movement. Figure 8a shows the location of the study site, the Yonglian irrigation area (1073719′′–1085104′′ E, 404557′′–411758′′ N) in Inner Mongolia, China. The irrigation area is 12 km long from north to south and 3 km wide from east to west. The whole domain size is 29.75 km2. The ground surface elevation decreases from 1028.9 to 1025.4 m from the southwest to the northeast. A 2-year tracer experiment from 2014 to 2016 was conducted to obtain the groundwater recharge (Yang, 2018), and the experimental locations are shown in Fig. 8a. This irrigation area has well-defined hydrogeological borders by the channel network. Since the Zaohuo Trunk Canal and No. 6 Drainage Ditch are filled with water over the simulation time, the first-kind boundary condition was applied to the two segments. The non-flow boundary condition was used for the other segments. The irrigation water of this area is diverted from the Renmin Canal. This irrigation area was divided into three sub-areas according to the land usage since they own significantly different upper boundary conditions, which are farm land, villages, and bared soil, as shown in Fig. 8b. The crop types in the farm land were not considered for determining the sub-areas. The surface digital elevation model (DEM) is shown in Fig. 8c.

The measured soil water content and groundwater table in the crop growing season from May to October of 2004 were used to calibrate the hydraulic parameters, and the tracer experiment from 2014 to 2016 was used for the groundwater recharge evaluation. A uniform daily rainfall rate was applied to the whole domain. The irrigation water was only applied to the farm land. Due to lack of the weather data in 2004, the potential evapotranspiration ET0 was calculated by the measured evaporation data from the 20 cm pan (ET20), multiplying by an empirical conversion coefficient. The empirical coefficient is 0.55, which was recommended by Hao (2016) by comparing monthly ET0 and ET20 with 8 years' data in this area. The ET0 during 2014 to 2016 was calculated by using the Penman–Monteith equation. The precipitation, irrigation, and ET0 are shown in Fig. 9. The crop growing season is from May to October, and the rest of the months are no-crop growing season. Based on the hydrogeological characteristics of the study area provided by the Geological Department of Inner Mongolia, the top aquifer within the depth of 7 m is loamy sand and loam with small hydraulic conductivity; an underlying sand aquifer with a thickness of 46 m has high permeability, and the sand aquifer lies on an impervious 1 m thick clay layer. The clay layer was used as the bottom of the simulation domain, and seven different geological layers were used in the MODFLOW model. The first layer was set to be the top aquifer, and the second aquifer was divided into six layers for numerical simulation. Ten groundwater monitoring wells were set in this district, and the groundwater tables were observed every 6 d. Well 1, well 2, well 3, well 5, and well 6 are located in farm land areas, well 4 and well 8 in villages, and well 7, well 9, and well 10 in bared soil areas. Additionally, there are five soil water content monitoring points in the farm land and two points in the bared soil area, as shown in Fig. 8a. Soil water contents within 1 m depth were observed one to three times every month from May to October in 2004.

Five GIS files are prepared as the shapefile files of the study domain, the land usage types, the boundary conditions, and raster files of the surface DEM and initial hydraulic head. There were 150 rows and 50 columns used in the MODFLOW model. The spatial discretization of UBMOD was set to be 0.1 m. The stress period ΔT was set as 5 d, and the MODFLOW time step Δts and UBMOD time step Δtu were set as 1 d.

Table 4The unsaturated hydraulic parameters of the real-world application.

Figure 11Spatial simulated water table depth at different output times of the real-world application.

## 4.2 Model calibration results

There are two soil types in the first layer as loamy sand and loam. The unsaturated hydraulic parameters of the two soils are listed in Table 4. The hydraulic conductivity of the top aquifer in MODFLOW was set as the same as the unsaturated layer, the hydraulic conductivity of the bottom sand aquifer was set as 3.5 m d−1 during the calibration, and the specific yields of the top and bottom were set as 0.08 and 0.1, respectively. Figure 10 shows the comparison of the simulated and observed water table depths for the whole area and locations of different monitoring wells. The statistical index values are listed in Table 5. It can be found that the ARE, RMSE, IA, and R2 values are 9.9 %, 0.203 m, 0.869, and 0.71 for the regional average water table depth. Larger deviations of simulated water table depth can be found for the locations of monitoring wells, with RMSE values ranging from 0.25 to 0.39 m. Figure 11 further shows the spatial distribution of the simulated water table depth at different output times. The increasing trend is obviously found in Fig. 11a to c in the crop growing season, during which the groundwater was consumed by crop transpiration and strong soil evaporation. When the intensive autumn irrigation happened after the 160th day, the water table depth in the farm land decreased rapidly, as shown in Fig. 11d. These results indicate that our model can reasonably simulate the water table depth trend in space and time.

Table 5The statistical index values of the real-world application.

The recharge during the short term was calculated for further checking the results by comparing the results with those from reference papers. The calculated recharge in farm land during the autumn irrigation (from 16 to 31 October) is 93.3 mm, and the coefficient of recharge from the autumn irrigation is 0.37. Zhang (2011) suggested that the coefficient of recharge from the autumn irrigation is approximately 0.3. Yang (2016) proposed that the coefficient of the recharge from the autumn irrigation is between 0.36 and 0.4. Yu (2017) used the coefficient of recharge from autumn irrigation as 0.33 for the district. The calculated result is consistent with the previous studies. The phreatic evaporation coefficient was estimated during the period from 15 to 30 September with no precipitation or irrigation. The quantity of the recharge from saturated zone to unsaturated zone is 10.1 mm during the period in the farm land. The phreatic evaporation coefficient is 0.179, and the averaged water table depth is 1.51 m during the period. The phreatic evaporation coefficient measured by Wang (2002) is 0.172 at the depth of 1.5 m. The short-term results indicate the validity of the simulation results.

Figure 12Comparison between simulated and observed regional average soil water content profiles of the real-world application.

Figure 12 shows the comparison between the simulated and average observed soil water content profiles of the farm land and bared soil at different times. The statistical index values are listed in Table 5. The ARE values of the farm land at the times of 40, 85, 125, and 166 d are 15.3 %–24.9 %, the RMSE values 0.044–0.066 cm3 cm−3, the IA values 0.621–0.775, and the R2 values 0.54–0.689. The corresponding values for the bared soil are 10.8 %–19.8 %, 0.038–0.052 cm3 cm−3, 0.823–0.905, and 0.620–0.813, respectively. The larger measured soil water content in the root zone for the farm land can be observed than the simulations, while the simulated soil water content profiles in the bared soil agree well with the observations, as shown in Fig. 12. The reason may be that the sampling locations are at the border of fields, which leads to an overestimation of the soil water content in the root zone due to smaller crop root uptake.

Table 6The recharge sources and results of the tracer experiment.

Note: P is the annual precipitation; I is the irrigation water; R is the annual recharge; and Rc is the recharge coefficient, ${R}_{\mathrm{c}}=R/\left(P+I\right)$.

Figure 13The recharge rate in the farm land calculated by the coupled model.

The computational cost of the real-world application is 120 s, which is efficient considering the scale of the problem.

## 4.3 Regional groundwater recharge

In the tracer experiment, bromide (Br) was used as the tracer for calculating groundwater recharge. The tracer was injected at 1 m depth at two locations shown in Fig. 8a in October 2014. Based on two sampling locations in October 2015 and 2016, the downward recharge is estimated according to the movement distance of the tracer peak and the average water content from the initial position of the tracer to the final position (Tan et al., 2014). The soil water content at the depth of 1 m is relatively stable according to the measurements and the results of Peng (2015), which ensures the reliability of the experiment. As shown in Table 6, the annual average recharge R is 33.8 mm yr−1, and the recharge coefficient is 0.055 during the period of 2014–2016.

The calibrated coupled model was used to estimate the groundwater recharge from 1 October 2014 to 30 September 2016. Figure 13 shows the time series of simulated recharge rate in the farm land, and Table 6 lists the simulation results. The simulation results indicate that groundwater is recharged in the no-crop growing season and consumed in the crop growing season. The two peak values of groundwater recharge in Fig. 13 are due to the autumn irrigation after harvest for washing salt out. The no-crop growing season provides 92.30 mm yr−1 groundwater recharge over a year and the average recharge coefficient is 0.346, which indicates that the autumn irrigation in the no-crop growing season provides the primary groundwater recharge in the year. In the crop growth season, the recharge is negative, which means that groundwater is consumed by crop transpiration and soil evaporation. As calculated by the coupled model, the annual groundwater recharge is 36.21 mm yr−1 during the period from 1 October 2014 to 30 September 2016 in the farm land, which is similar to the result of the tracer experiment. The results confirm the coupled model for groundwater recharge evaluation, which is helpful for scheduling the irrigation amount in the crop growing season under the water saving policies.

5 Conclusions

This study developed a new quasi-3-D coupled model for the purpose of practical modeling of unsaturated–saturated flow at the regional scale. The UBMOD 1-D water balance model describing the unsaturated soil water flow was integrated with MODFLOW iteratively. A developed framework implemented the modeling procedures and provided the pre- and post-processing tools. The model was evaluated by using both synthetic numerical examples and real-world experimental data. The major conclusions drawn from this research are as follows.

1. The new iteration coupling scheme iteratively integrating a hydrodynamic model with a water balance model is reliable. The vertical net recharge and the depth of the unsaturated zone are effective to be used as the exchange information to couple the unsaturated zone and saturated zone.

2. The satisfactory results in the two testing examples demonstrate the effectiveness of the new quasi-3-D model with an acceptable calculative efficiency and well-maintained saturated zone and unsaturated zone mass balance.

3. The spatial and temporal discretization has slight impact on the simulation results. The stress period should be not too large and it also has slight impact on the simulation results in a certain range.

4. The model gives a satisfactory performance for calculating the groundwater recharge measured from the tracer experiment. The calculated annual groundwater recharge is 36.21 mm yr−1 and the recharge coefficient is 0.059 in the study area.

5. The coupled model should not be used for problems with substantial lateral flow in the unsaturated zone because of the quasi-3-D assumptions used in the model.

6. The coupled model could lead to a higher groundwater table depth since it ignores the capillary fringe.

Data availability
Data availability.

All the data and codes used in this study can be requested by email to the corresponding author Yan Zhu at zyan0701@163.com.

Appendix A

UBMOD is a water balance model based on a hybrid of numerical and statistical methods. Compared with the traditional ones, the model can simulate upward soil water movement in heterogeneous situations.

There are four major components to describe the soil water movement in the UBMOD model, as shown in Fig. A1. Firstly, the vertical soil column is divided into a cascade of “buckets” and each “bucket” corresponds to a soil layer. The “buckets” will be filled to saturation from the top layer to the bottom layer if there is infiltration, which is referred to as the allocation of infiltration water. Specifically speaking, the infiltration water first fills the top “bucket”, then the excessive infiltration water moves downward to the next “bucket”, until all the infiltration water is allocated in the “buckets”, as shown in Fig. A1a. The governing equation of layer i is

$\begin{array}{}\text{(A1)}& {q}_{i}=min\left({M}_{i}×\left({\mathit{\theta }}_{\mathrm{s},i}-{\mathit{\theta }}_{i}\right),I-{I}_{\mathrm{d},i-\mathrm{1}}\right),\end{array}$

where i indicates the vertical soil layer, i=1, …, j; qi is the amount of allocated water per unit area of layer i (L); Mi is the thickness of layer i (L); θi is the initial soil water content of layer i (L3 L−3); θs,i is the saturated soil water content of layer i (L3 L−3); I is the quantity of infiltration water per unit area (L); ${I}_{\mathrm{d},i-\mathrm{1}}$ is the consumed infiltration water per unit area by upper layers for the specific layer i (L); and $I-{I}_{\mathrm{d},i-\mathrm{1}}$ is the infiltration water applied to the upper boundary of soil layer i (L). The infiltration rate I is input data in the model, and the partitioning of rainfall between infiltration and runoff has not been considered by now. As shown in Fig. A1a, the first three layers are filled to saturation, and the fourth layer is filled with the residual infiltration water.

Secondly, when the soil water content exceeds the field capacity, the soil water will move downward driven by the gravitational potential. The governing equation is

$\begin{array}{}\text{(A2)}& \frac{\partial \mathit{\theta }}{\partial t}=-\frac{\partial K\left(\mathit{\theta }\right)}{\partial z},\end{array}$

where t is the time (T); K(θ) is the unsaturated hydraulic conductivity (L T−1) as a function of soil water content; z is the elevation in the vertical direction (L). The vertical coordinate is positive downward. The unsaturated hydraulic conductivity K(θ) in Eq. (A2) is a function of soil water content θ. The relationship between K(θ) and θ is characterized by empirical formulas for the purpose of simplifying calculation and eliminating the soil hydraulic parameters. These empirical formulas are referred to as drainage functions. The commonly used equations are listed in Table A1.

Figure A1The schematic procedures of UBMOD with (a) the allocation of the infiltration water, (b) the soil water advective movement driven by the gravitational potential, (c) the source/sink terms, and (d) the soil water diffusive movement driven by the matric potential.

Thirdly, the source/sink terms are used to account for soil evaporation and crop transpiration. The governing equation is as follows:

$\begin{array}{}\text{(A3)}& \frac{\partial \mathit{\theta }}{\partial t}=-W,\end{array}$

where W is the source/sink term (T−1). The Penman–Monteith formula and Beer's law (also known as a Ritchie-type equation) are adopted in UBMOD to estimate the potential soil evaporation Ep and potential crop transpiration Tp. Then Ep and Tp are distributed to each layer based on the evaporation cumulative distribution function and the root density function. The actual soil evaporation and crop transpiration are obtained by discounting Ep and Tp with the soil water stress coefficient.

Lastly, we calculate the diffusive movement driven by the matric potential. The governing equation is

$\begin{array}{}\text{(A4)}& \frac{\partial \mathit{\theta }}{\partial t}=\frac{\partial }{\partial z}\left(D\left(\mathit{\theta }\right)\frac{\partial \mathit{\theta }}{\partial z}\right),\end{array}$

where D(θ) is the hydraulic diffusivity (L2 T−1), and $D\left(\mathit{\theta }\right)=K\left(\mathit{\theta }\right)×\partial \phantom{\rule{0.125em}{0ex}}h/\partial \mathit{\theta }$, where h is the matric potential (L). The finite-difference method is used to solve the governing equation. A new empirical formula is presented to describe the hydraulic diffusivity D(θ). The expression formula of D(θ) has an exponential form as

$\begin{array}{}\text{(A5)}& D\left(\mathit{\theta }\right)={\mathrm{10}}^{a×S\left(\mathit{\theta }\right)+b},\end{array}$

where S(θ) is the effective saturation (–); a and b are two intermediate parameters. In order to eliminate the parameters, we calculate the hydraulic diffusivity D(θ) of different soils by the van Genuchten model firstly, and then fit the hydraulic diffusivity D(θ) by Eq. (A5). Furthermore, we establish the relationship between the two intermediate parameters (a and b) and the saturated hydraulic conductivity Ks as

$\begin{array}{}\text{(A6)}& \left\{\begin{array}{l}b=-\mathrm{3.55}+\mathrm{0.55}×{\mathrm{log}}_{\mathrm{10}}\left({K}_{\mathrm{s}}\right)-\mathrm{1.36}×{\mathrm{log}}_{\mathrm{10}}{\left({K}_{\mathrm{s}}\right)}^{\mathrm{2}}\\ a=\mathrm{3.72}+\mathrm{0.61}×{\mathrm{log}}_{\mathrm{10}}\left({K}_{\mathrm{s}}\right)+\mathrm{1.52}×{\mathrm{log}}_{\mathrm{10}}{\left({K}_{\mathrm{s}}\right)}^{\mathrm{2}}\end{array}\right\\end{array}$

By following the steps above, the hydraulic diffusivity D(θ) of a specific soil type can be estimated with three physical meaning parameters (Ks, θs, and θr).

Table A1Empirical drainage functions representing the relationship between the unsaturated hydraulic conductivity K(θ) and soil water content θ.

Note: the parameters Ks, θs, and θf are the saturated hydraulic conductivity (L T−1), the saturated water content (L3 L−3), and the field capacity (L3 L−3); θw is the soil water content at wilting point (L3 L−3); α in the exponential equation is a site-specific parameter determined mainly from soil characteristics, and it has an inverse relationship with Ks, and the value ranges between 10 and 30 (Jiang et al., 2008). The parameter β in the power equation ensures K(θ) approaches zero when θ approaches θf, and Arnold et al. (1990) proposed an empirical formula as $\mathit{\beta }=-\mathrm{2.655}/\mathrm{log}\left({\mathit{\theta }}_{\mathrm{f}}/{\mathit{\theta }}_{\mathrm{s}}\right)$.

Soil water content is discontinuous at the material interface when a soil profile is heterogeneous. When adopting the van Genuchten model to represent the soil water characteristics, Eq. (A4) can be expressed as

$\begin{array}{}\text{(A7)}& \begin{array}{rl}& \frac{\partial \mathit{\theta }}{\partial t}=\frac{\partial }{\partial z}\\ & \left(D\left(\mathit{\theta }\right)\left(\frac{\partial \mathit{\theta }}{\partial z}-\left(\frac{\partial \mathit{\theta }}{\partial {\mathit{\theta }}_{\mathrm{s}}}\frac{\partial {\mathit{\theta }}_{\mathrm{s}}}{\partial z}+\frac{\partial \mathit{\theta }}{\partial {\mathit{\theta }}_{\mathrm{r}}}\frac{\partial {\mathit{\theta }}_{\mathrm{r}}}{\partial z}+\frac{\partial \mathit{\theta }}{\partial \mathit{\alpha }}\frac{\partial \mathit{\alpha }}{\partial z}+\frac{\partial \mathit{\theta }}{\partial n}\frac{\partial n}{\partial z}\right)\right)\right),\end{array}\end{array}$

where α and n are two parameters in the van Genuchten model. The two terms $\partial \mathit{\theta }/\partial {\mathit{\theta }}_{\mathrm{s}}$ and $\partial \mathit{\theta }/\partial {\mathit{\theta }}_{\mathrm{r}}$ can be easily calculated. The value of $\partial \mathit{\theta }/\partial \mathit{\alpha }$ is close to zero, which can be ignored. A regression formula is developed to characterize the relationship between n and the saturated hydraulic conductivity Ks. The specific equations are shown as follows:

$\begin{array}{}\text{(A8)}& \left\{\begin{array}{l}\frac{\partial \mathit{\theta }}{\partial {\mathit{\theta }}_{\mathrm{s}}}=S\left(\mathit{\theta }\right),\\ \frac{\partial \mathit{\theta }}{\partial {\mathit{\theta }}_{\mathrm{r}}}=\mathrm{1}-S\left(\mathit{\theta }\right),\\ \frac{\partial \mathit{\theta }}{\partial n}=\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)S\left(\mathit{\theta }\right)\\ \left[\frac{\mathrm{ln}\left(S\left(\mathit{\theta }\right)\right)}{n\left(n-\mathrm{1}\right)}-\frac{n-\mathrm{1}}{{n}^{\mathrm{2}}}\left(S{\left(\mathit{\theta }\right)}^{n/\mathrm{1}-n}-\mathrm{1}\right)\right\\ \mathrm{ln}\left(S{\left(\mathit{\theta }\right)}^{n/\mathrm{1}-n}-\mathrm{1}\right)S{\left(\mathit{\theta }\right)}^{-n/\mathrm{1}-n}],\\ n=\mathrm{0.9505}+\mathrm{0.8883}{e}^{\mathrm{0.7751}×{\mathrm{log}}_{\mathrm{10}}\left({K}_{\mathrm{s}}\right)}.\end{array}\right\\end{array}$

Therefore, the diffusive term of the heterogeneous soil can be calculated. With the help of the diffusive term, UBMOD can consider upward soil water movement, which is ignored by most water balance models.

Author contributions
Author contributions.

WM, YZ, and JY proposed the structure and coupling method of the model, and developed the model code. YZ, JY, and JW provided data support for the real-world application. All the authors contributed to model validation. WM prepared the manuscript with contributions from all the co-authors.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The study was supported by the Natural Science Foundation of China. The author would thank the General Administration of Hetao Irrigation District and the members of the 311 laboratory for their help in field investigation.

Financial support
Financial support.

This research has been supported by the Natura Science Foundation of China (grant nos. 51790532, 51779178, and 51629901).

Review statement
Review statement.

This paper was edited by Philippe Ackerer and reviewed by Gerrit H. de Rooij and one anonymous referee.

References

Abedinpour, M., Sarangi, A., Rajput, T. B. S., and Singh, M.: Prediction of maize yield under future water availability scenarios using the AquaCrop model, J. Agric. Sci., 152, 558–574, https://doi.org/10.1017/S0021859614000094, 2014.

Adnan, A. A., Jibrin, J. M., Kamara, A. Y., Abdulrahman, B. L., Shaibu, A. S., and Garba, I. I.: CERES-maize model for determining the optimum planting dates of early maturing maize varieties in Northern Nigeria, Front. Plant Sci., 8, 1118, https://doi.org/10.3389/fpls.2017.01118, 2017.

Arnold, J. G., Williams, J. R., Griggs, R. H., and Sammons, N. B.: SWRRB – A basin scale simulation model for soil and water resources management, M Press, Texas A, 1990.

Arnold, J., Moriasi, D., Gassman, P., Abbaspour, K., White, M., Srinivasan, R., Santhi, C., Harmel, R., van Griensven, A., Van Liew, M., Kannan, N., and Jha, M.: SWAT: model use, calibration, and validation, T. ASABE, 55, 1491–1508, https://doi.org/10.13031/2013.42256, 2012.

Bakker, M., Post, V., Langevin, C., Hughes, J., White, J., Starn, J., and Fienen, M.: Scripting MODFLOW model development using Python and FloPy, Groundwater, 54, 733–739, https://doi.org/10.1111/gwat.12413, 2016.

Bouwer, H.: Integrated water management: emerging issues and challenges, Agric. Water Manage., 45, 217–228, https://doi.org/10.1016/S0378-3774(00)00092-5, 2000.

Brunner, P. and Simmons, C.: HydroGeoSphere: a fully integrated, physically based hydrological model, Groundwater, 50, 170–176, https://doi.org/10.1111/j.1745-6584.2011.00882.x, 2012.

Clement, T., Wise, W., and Molz, F.: A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow, J. Hydrol., 161, 71–90, https://doi.org/10.1016/0022-1694(94)90121-X, 1994.

Diersch, H.: FEFLOW: finite element modeling of flow, mass and heat transport in porous and fractured media, Springer Science & Business Media, Berlin, Germany, 2013.

Downer, C. and Ogden, F.: Appropriate vertical discretization of Richards' equation for two-dimensional watershed-scale modelling, Hydrol. Process., 18, 1–22, https://doi.org/10.1002/hyp.1306, 2004.

Duda, P., Hummel, P., Donigian Jr., A., and Imhoff, J.: BASINS/HSPF: model use, calibration, and validation, T. ASABE, 55, 1523–1547, https://doi.org/10.13031/2013.42261, 2012.

Evans, R. and Sadler, E.: Methods and technologies to improve efficiency of water use, Water Resour. Res., 44, W00E04, https://doi.org/10.1029/2007WR006200, 2008.

Facchi, A., Ortuani, B., Maggi, D., and Gandolfi, C.: Coupled SVAT–groundwater model for water resources simulation in irrigated alluvial plains, Environ. Modell. Softw., 19, 1053–1063, https://doi.org/10.1016/j.envsoft.2003.11.008, 2004.

Fill, V.: Documentation of computer program INFIL3.0 – A distributed-parameter watershed model to estimate net infiltration below the zone, U. S. Geological Survey, Virginia, US, 2008.

Furman, A.: Modeling coupled surface-subsurface flow processes: a review, Vadose Zone J., 7, 741–756, https://doi.org/10.2136/vzj2007.0065, 2008.

Graham, D. and Butts, M.: Flexible, integrated watershed modelling with MIKE SHE, in: Watershed Models, edited by: Singh, V. and Frevert, D., CRC Press, Cleveland, Ohio, US, 2005.

Hao, P.: Regional soil water-salt balance in Hetao Irrigation District with drip irrigation and combined use of surface water and groundwater, Master thesis, School of Water Resources and Hydropower Engineering, Wuhan University, China, 2016.

Harbaugh, A.: MODFLOW-2005, the U.S. Geological Survey modular ground-water model – the Ground-Water Flow Process, US Geological Survey, Virginia, US, 2005.

Harter, T. and Hopmans, J.: Role of vadose zone flow processes in regional scale hydrology: Review, opportunities and challenges, in: Unsaturated Zone Modeling: Progress, Challenges and Applications, edited by: Feddes, R., de Rooij, G., and van Dam, J., Kluwer Academic Publishers, Dordrecht, Netherlands, 179–210, 2004.

Havard, P., Prasher, S., Bonnell, R., and Madani, A.: Linkflow, a water flow computer model for water table management: Part I. Model development, T. ASABE., 38, 481–488, https://doi.org/10.13031/2013.27856, 1995.

Holzworth, D. P., Huth, N. I., deVoil, P. G., Zurcher, E. J., Herrmann, N. I., McLean, G., Chenu, K., van Oosterom, E. J., Snow, V., Murphy, C., Moore, A. D., Brown, H., Whish, J. P. M., Verrall, S., Fainges, J., Bell, L. W., Peake, A. S., Poulton, P. L., Hochman, Z., Thorburn, P. J., Gaydon, D. S., Dalgliesh, N. P., Rodriguez, D., Cox, H., Chapman, S., Doherty, A., Teixeira, E., Sharp, J., Cichota, R., Vogeler, I., Li, F. Y., Wang, E., Hammer, G. L., Robertson, M. J., Dimes, J. P., Whitbread, A. M., Hunt, J., van Rees, H., McClelland, T., Carberry, P. S., Hargreaves, J. N. G., MacLeod, N., McDonald, C., Harsdorf, J., Wedgwood, S., and Keating, B. A.: APSIM-Evolution towards a new generation of agricultural systems simulation, Environ. Modell. Softw., 62, 327–350, https://doi.org/10.1016/j.envsoft.2014.07.009, 2014.

Jiang, J., Zhang, Y., Wegehenkel, M., Yu, Q., and Xia, J.: Estimation of soil water content and evapotranspiration from irrigated cropland on the North China Plain, J. Plant Nutr. Soil Sci., 171, 751–761, https://doi.org/10.1002/jpln.200625179, 2008.

Karandish, F., Salari, S., and Darzi-Naftchali, A.: Application of virtual water trade to evaluate cropping pattern in arid regions, Water Resour. Manage., 29, 4061–4074, https://doi.org/10.1007/s11269-015-1045-4, 2015.

Kendy, E., Gérard-Marchant, P., Walter, M. T., Zhang, Y. Q., Liu, C. M., and Steenhuis, T. S.: A soil-water-balance approach to quantify groundwater recharge from irrigated cropland in the North China Plain, Hydrol. Process., 17, 2011–2031, https://doi.org/10.1002/hyp.1240, 2003.

Kim, N., Chung, I., Won, Y., and Arnold, J.: Development and application of the integrated SWAT–MODFLOW model, J. Hydrol., 356, 1–16, https://doi.org/10.1016/j.jhydrol.2008.02.024, 2008.

Kuznetsov, M., Yakirevich, A., Pachepsky, Y., Sorek, S., and Weisbrod, N.: Quasi 3D modeling of water flow in vadose zone and groundwater, J. Hydrol., 450, 140–149, https://doi.org/10.1016/j.jhydrol.2012.05.025, 2012.

Lachaal, F., Mlayah, A., Bédir, M., Tarhouni, J., and Leduc, C.: Implementation of a 3-D groundwater flow model in a semi-arid region using MODFLOW and GIS tools: The Zéramdine–Béni Hassen Miocene aquifer system (east-central Tunisia), Comput. Geosci.-UK., 48, 187–198, https://doi.org/10.1016/j.cageo.2012.05.007, 2012.

Leterme, B., Gedeon, M., Laloy, E., and Rogiers, B.: Unsaturated flow modeling with HYDRUS and UZF: calibration and intercomparison, In Proc. MODFLOW and More 2015, Integrated GroundWater Modeling Center, 31 May– 3 June, Golden, CO, 2015.

Liang, X., Xie, Z., and Huang, M.: A new parameterization for surface and groundwater interactions and its impact on water budgets with the variable infiltration capacity (VIC) land surface model, J. Geophys. Res.-Atmos., 108, 8613, https://doi.org/10.1029/2002JD003090, 2003.

Mao, W., Yang, J., Zhu, Y., Ye, M., Liu, Z., and Wu, J.: An efficient soil water balance model based on hybrid numerical and statistical methods, J. Hydrol., 559, 721–735, https://doi.org/10.1016/j.jhydrol.2018.02.074, 2018.

Markstrom, S., Niswonger, R., Regan, R., Prudic, D., and Barlow, P.: GSFLOW-Coupled Ground-water and Surface-water FLOW model based on the integration of the Precipitation-Runoff Modeling System (PRMS) and the Modular Ground-Water Flow Model (MODFLOW-2005), US Geological Survey, Virginia, US, 2008.

Merritt, W. S., Letcher, R. A., and Jakeman, A. J.: A review of erosion and sediment transport models, Environ. Modell. Softw., 18, 761–799, https://doi.org/10.1016/S1364-8152(03)00078-1, 2003.

Niswonger, R., Prudic, D., and Regan, R.: Documentation of the Unsaturated-Zone Flow (UZF1) Package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005, US Geological Survey, Virginia, US, 2006.

Oosterbaan, R.: SALTMOD: Description of Principles and Application, ILRI, Wageningen, 1998.

Peng, Z. Y.: Mechanism and modeling of coupled water-heat-solute movement in unidirectional freezing soils, Doctor thesis, School of Water Resources and Hydropower Engineering, Wuhan University, China, 2015.

Raes, D.: BUDGET: a soil water and salt balance model, Reference manual, version 5, available at: http://iupware.be/?page_id=_820 (last access: 13 August 2019), 2002.

Ranatunga, K., Nation, E., and Barratt, D.: Review of soil water models and their applications in Australia, Environ. Modell. Softw., 23, 1182–1206, https://doi.org/10.1016/j.envsoft.2008.02.003, 2008.

Seo, H., Šimůnek, J., and Poeter, E.: Documentation of the hydrus package for MODFLOW-2000, the us geological survey modular ground-water model, IGWMC-International Ground Water Modeling Center, US, 2007.

Shen, C. and Phanikumar, M.: A process-based, distributed hydrologic model based on a large-scale method for surface-subsurface coupling, Adv. Water Resour., 33, 1524–1541, https://doi.org/10.1016/j.advwatres.2010.09.002, 2010.

Sheikh, V., Visser, S., and Stroosnijder, L.: A simple model to predict soil moisture: bridging event and continuous hydrological (BEACH) modelling, Environ. Modell. Softw., 24, 542–556, https://doi.org/10.1016/j.envsoft.2008.10.005, 2009.

Sherlock, M., McDonnell, J., Curry, D., and Zumbuhl, A.: Physical controls on septic leachate movement in the vadose zone at the hillslope scale, Putnam County, New York, USA, Hydrol. Preocess., 16, 2559–2575, https://doi.org/10.1002/hyp.1048, 2002.

Sluiter, R.: Desertification and Grazing On south Crete, A model approach. Department of Physical Geography, Faculty of Geographical Sciences, University of Utrecht, Utrecht, 1998.

Šimůnek, J., Vogel, T., and van Genuchten, M. T.: The SWMS_2D code for simulating water flow and solute transport in two-dimensional variably saturated media, Research Report, California, US, 1994.

Šimůnek, J., Šejna, M., Saito, H., Sakai, M., and van Genuchten, M. T.: The HYDRUS-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 4.0, HYDRUS Software Series 3, Department of Environmental Sciences, University of California Riverside, Riverside, California, US, 2008.

Šimůnek, J., van Genuchten, M. T., and Šejna, M.: HYDRUS: Model use, calibration and validation, T. ASABE., 55, 1261–1274, https://doi.org/10.13031/2013.42239, 2012.

Sophocleous, M.: Groundwater recharge and sustainability in the High Plains aquifer in Kansas, USA, Hydrogeol. J., 13, 351–365, https://doi.org/10.1007/s10040-004-0385-6, 2005.

Stoppelenburg, F., Kovar, K., Pastoors, M., and Tiktak, A.: Modelling the interactions between transient saturated and unsaturated groundwater flow, RIVM report 500026001, 2005.

Szymkiewicz, A., Gumuła-Kawęcka, A., Šimůnek, J., Leterme, B., Beegum, S., Jaworska-Szulc, B., Pruszkowska-Caceres, M., Gorczewska-Langner, W., Angulo-Jaramillo, R., and Jacques, D.: Simulations of freshwater lens recharge and salt/freshwater interfaces using the HYDRUS and SWI2 packages for MODFLOW, J. Hydrol. Hydromech., 66, 246–256, https://doi.org/10.2478/johh-2018-0005, 2018.

Tan, X., Wu, J., Cai, S., and Yang, J.: Characteristics of groundwater recharge on the North China Plain. Groundwater, 52, 798–807, https://doi.org/10.1111/gwat.12114, 2014.

Thoms, R., Johnson, R., and Healy, R.: User's guide to the variably saturated flow (VSF) process for MODFLOW, US Geological Survey Techniques and Methods 6-A18, Virginia, US, 2006.

Tian, Y., Zheng, Y., Wu, B., Wu, X., Liu, J., and Zheng, C.: Modeling surface water-groundwater interaction in arid and semi-arid regions with intensive agriculture, Environ. Modell. Softw., 63, 170–184, https://doi.org/10.1016/j.envsoft.2014.10.011, 2015.

Toms, S.: ArcPy and ArcGIS-Geospatial Analysis with Python, Packt Publishing Ltd, Birmingham, UK, 2015.

Twarakavi, N., Šimůnek, J., and Seo, H.: Evaluating interactions between groundwater and vadose zone using the HYDRUS-based flow package for MODFLOW, Vadose Zone J., 7, 757–768, https://doi.org/10.2136/vzj2007.0082, 2008.

Vaccaro, J.: A deep percolation model for estimating ground-water recharge, Documentation of modules for the modular modeling system of the US Geological Survey, available at: https://pubs.er.usgs.gov/publication/sir20065318 (last access: 13 August 2019), 2007.

Van Walsum, P. and Groenendijk, P.: Quasi steady-state simulation of the unsaturated zone in groundwater modeling of lowland regions, Vadose Zone J., 7, 769–781, https://doi.org/10.2136/vzj2007.0146, 2008.

VanderKwaak, J. and Loague, K.: Hydrologic-Response simulations for the R-5 catchment with a comprehensive physics-based model, Water Resour. Res., 37, 999–1013, https://doi.org/10.1029/2000WR900272, 2001.

Varado, N., Ross, P., and Haverkamp, R.: Assessment of an efficient numerical solution of the 1D Richards equation on bare soil, J. Hydrol., 323, 244–257, https://doi.org/10.1016/j.jhydrol.2005.07.052, 2006.

Vauclin, M., Khanji, J., and Vachaud, G.: Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem, Water Resour. Res., 15, 1089–1101, https://doi.org/10.1029/WR015i005p01089, 1979.

Verburg, K.: Methodology in soil water and solute balance modelling: an evaluation of APSIM-SoilWat and SWIMv2 models, Report of APSRU, CSIRO Division of Soils Workshop, Brisbane, Qld, 1995.

Wang, W., Zhang, Z., Yeh, T., Qiao, G., Wang, W., Duan, L., Huang, S., and Wen, J.: Flow dynamics in vadose zones with and without vegeration in an arid region, Adv. Water Resour., 106, 68–79, https://doi.org/10.1016/j.advwatres.2017.03.011, 2017.

Wang, Y.: The analysis of the regional scale groundwater table variation before and after the water-saving transformation in Hetao Irrigation District, Water Saving Irrigation, 1, 15–17, 2002.

Weill, S., Mouche, E., and Patin, J.: A generalized Richards equation for surface/subsurface flow modeling, J. Hydrol., 336, 9–20, https://doi.org/10.1016/j.jhydrol.2008.12.007, 2009.

Xu, X., Huang, G., Qu, Z., and Pereira, L.: Using MODFLOW and GIS to access changes in groundwater dynamics in response to water saving measures in irrigation districts of the upper Yellow River basin, Water Resour. Manage., 25, 2035–2059, https://doi.org/10.1007/s11269-011-9793-2, 2011.

Xu, X., Huang, G., Zhan, H., Qu, Z., and Huang, Q.: Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas, J. Hydrol., 412, 170–181, https://doi.org/10.1016/j.jhydrol.2011.07.002, 2012.

Yakirevich, A., Borisov, V., and Sorek, S.: A quasi three-dimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi two-dimensional case, Adv. Water Resour., 21, 679–689, https://doi.org/10.1016/S0309-1708(97)00031-6, 1998.

Yang, J., Zhu, Y., Zha, Y., and Cai, S.: Mathematical model and numerical method of groundwater and soil water movement, Science press, Beijing, China, 2016.

Yang, W.: Numerical simulation of conjunctive use of groundwater and surface water in Yongji Irrigation Field of Hetao Irrigation District, Master thesis, School of Water Resources and Hydropower Engineering, Wuhan University, China, 2016.

Yang, X.: Soil salt balance in Hetao Irrigation District based on the SaltMod and tracer experiment, Master thesis, School of Water Resources and Hydropower Engineering, Wuhan University, China, 2018.

Yu, L.: Numerical simulation of conjunctive use of groundwater and surface water in Hetao Irrigation District and water resources forecase, Master thesis, School of Water Resouces and Hydropower Engineering, Wuhan Univerisity, China, 2017.

Wang, X., Williams, J. R., Gassman, P. W., Baffaut, C., Izaurralde, R. C., Jeong, J., and Kiniry, R.: EPIC and APEX: model use, calibration, and validation, T. ASABE, 55, 1447–1462, https://doi.org/10.13031/2013.42253, 2012.

Zha, Y., Yang, J., Yin, L., Zhang, Y., Zeng, W., and Shi, L.: A modified Picard iteration scheme for overcoming numerical difficulties of simulating infiltration into dry soil, J. Hydrol. 551, 56–69, https://doi.org/10.1016/j.jhydrol.2017.05.053, 2017.

Zhang, J., Zhu, Y., Zhang, X., Ye, M., and Yang, J.: Developing a long short-term memory (LSTM) based model for predicting water table depth in agricultural areas, J. Hydrol., 561, 918–929, https://doi.org/10.1016/j.jhydrol.2018.04.065, 2018.

Zhang, Z.: Irrigation infiltration and recharge coefficient in Hetao Irrigation District and the primary study on threshold value of water in different diversion, Master thesis, Inner Mongolia Agricultural University, China, 2011.

Zhu, Y., Shi, L., Lin, L., Yang, J., and Ye, M.: A fully coupled numerical modeling for regional unsaturated-saturated water flow, J. Hydrol., 475, 188–203, https://doi.org/10.1016/j.jhydrol.2012.09.048, 2012.