the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A comprehensive quasi3D model for regionalscale unsaturated–saturated water flow
Wei Mao
Heng Dai
Ming Ye
Jinzhong Yang
Jingwei Wu
For computationally efficient modeling of unsaturated–saturated flow in regional scales, the quasithreedimensional (3D) scheme that considers onedimensional (1D) soil water flow and 3D groundwater flow is an alternative method. However, it is still practically challenging for regionalscale 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 quasi3D model coupled to the UBMOD 1D soil water balance model with the MODFLOW 3D 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 postprocessing 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 realworld study. The measured groundwater table and soil water content were used to calibrate the model parameters, and the groundwater recharge data from a 2year 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.
While groundwater resources are important for domestic, agricultural, and industrial uses, groundwater is vulnerable due to overexploitation, 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 threedimensional (3D) 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 3D models have a solid theoretical foundation and have been used for regionalscale 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 largescale 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 3D numerical models.
To address the computational challenges discussed above, a variety of simplifications have been introduced for the soil water flow for regionalscale problems. One simplification is to treat the hydrological processes (e.g., infiltration, evapotranspiration, and deep percolation) occurring in the unsaturated zone as onedimensional (1D) 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 1D simplification leads to the quasi3D scheme, which ignores the lateral flow in the unsaturated zone but considers groundwater flow as a 3D problem. The quasi3D scheme avoids solving the 3D Richards equation for the unsaturated zone and thus improves computational efficiency and model stability. The quasi3D scheme is an efficient solution for largescale 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 quasi3D models for a practical regionalscale problem. Two concerns arise as follows.
The first concern is the unsaturated modeling method. Although the quasi3D scheme is computationally efficient, the numerical solutions of the 1D 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 regionalscale problems under complicated geological and climate conditions as well as anthropogenic activities. As an alternative to the numerical solutions of the 1D 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 largescale 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), SVATMODFLOW (Facchi et al., 2004), UZF1MODFLOW (Niswonger et al., 2006), HYDRUSMODFLOW (Seo et al., 2007), and SWAPMODFLOW (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 tradeoff 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 quasi3D model is developed. The UBMOD 1D 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 largescale 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 realworld application to study the regional net groundwater recharge is presented in Sect. 4.
In the new coupled model, the unsaturated–saturated domain is partitioned into a number of subareas in the horizontal direction mainly according to the spatially distributed inputs (soil types, atmosphere boundary conditions, land usage types, and crop types). A 1D soil column is used to characterize the average soil water flow in each subarea, and UBMOD is used to simulate the 1D soil water flow. MODFLOW is used to simulate the 3D 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 subarea 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 selfcontained, 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
where i indicates the vertical soil layer, i=1, …, j; q_{i} is the amount of allocated water per unit area of layer i (L); M_{i} is the thickness of layer i (L); θ_{i} is the initial soil water content of layer i (L^{3} L^{−3}); θ_{s,i} is the saturated soil water content of layer i (L^{3} 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
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.
Thirdly, the source/sink terms are used to account for soil evaporation and crop transpiration. The governing equation is as follows:
where W is the source/sink term (T^{−1}). The Penman–Monteith formula and Beer's law (also known as a Ritchietype equation) are adopted to estimate the potential soil evaporation E_{p} and potential crop transpiration T_{p}. Then E_{p} and T_{p} 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 E_{p} and T_{p} with the soil water stress coefficient.
Lastly, we calculate the diffusive movement driven by the matric potential. The governing equation is
where D(θ) is the hydraulic diffusivity (L^{2} T^{−1}). The finitedifference method is used to solve the equation. An empirical formula with four parameters (saturated hydraulic conductivity K_{s}, 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 righthand 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 3D groundwater flow equation for a porous medium using a blockcentered finitedifference method (Harbaugh, 2005). The governing equation solved by MODFLOW is
where $i,j=\mathrm{1}$–3 indicate the x, y, and z directions, respectively; K_{ij} 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 (L^{3} T^{−1}); S_{s} is the specific storage of the porous material (L^{−1}); and t is the time (T).
FloPy and ArcPy are the two peripheral tools used in the model development. FloPy is a Python package for creating, running, and postprocessing MODFLOWbased models. Unlike the common graphical user interfaces (GUIs) method, FloPy facilitates users to write a Python script to construct and postprocess 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 preprocessing, the coupled model, and the postprocessing. The preparation of geographic input information of the model shown in Fig. 1a is the major component of preprocessing. The geographic information includes the domain area, boundary conditions, subareas, 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 subareas is joined with the domain file, represented in the subarea array with different numbers specified as different subareas. 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 postprocessing 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 subareas in the horizontal direction mainly according to the spatially distributed inputs (each subarea is considered to be homogeneous in the horizontal). There are l subareas and j layers for a specific soil column shown in Fig. 1b. Soil water flow of each subarea is simulated by using one 1D 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 subarea receive the same recharge from the soil zone calculated by the representative 1D 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 1D flow, the surrounding boundary conditions for the unsaturated zone are noflux 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 subarea, and scalar d_{u} 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:
where q_{I}, q_{A}, q_{S}, and q_{D} 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 q_{I}. Then the advective flow q_{A} across the water table driven by gravitational potential is calculated by Eq. (2). The direction of these two terms is downward. The q_{S} 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 q_{S} 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, M_{j+1} (L), of the layer is set as
where z_{j+1} is the bottom depth of layer j+1 (L); d_{u} is the thickness of the unsaturated zone (L). The amount of the upward flux between the virtual layer and layer j is denoted as q_{D}. 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 Δt_{s}, and the calculation time step for UBMOD Δt_{u}. 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. Δt_{u} is a priori value and cannot be changed during the calculation. UBMOD can give acceptable results when Δt_{u} is shorter than 10 d for assumed cases and 1 d for a realworld case (Mao et al., 2018). Δt_{s} 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 t+ΔT. At the time t, the saturated hydraulic head is known, marked as H^{t} (m×n dimension). When the model runs from t to t+ΔT, firstly, the initial saturated hydraulic head H^{t+ΔT} at t+ΔT is set to be equal to H^{t}, and then the average unsaturated depth from t to t+ΔT is calculated according to H^{t+ΔT}, marked as $\mathit{D}{\mathit{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:
where $\stackrel{\mathrm{\u203e}}{D}$ is the average depth from the soil surface to the impermeable layer of the controlling domain of the soil column (L); $\stackrel{\mathrm{\u203e}}{{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 Δt_{u} to obtain the vertical recharge at each time step (marked as r_{t}) until the time comes to be t+ΔT. The total recharge during the stress period ΔT (from t to t+ΔT) R_{ΔT} can be obtained by summarizing the recharge at each unsaturated time step, as follows:
The average recharge R from t to t+ΔT can be obtained by
Then the average recharge from all 1D soil columns can be obtained, represented as ${\mathit{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 Δt_{s} to obtain the saturated hydraulic head until the time comes to t+ΔT. The hydraulic head at the time t+ΔT 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 H^{t+ΔT}. The convergence criterion is
where ε_{H} is a userspecified tolerance (L). If the criterion is met, the iteration stops, and ${\mathbf{H}}^{t+\mathrm{\Delta}T,p}$ is the convergent results at time t+ΔT, 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.
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 HYDRUS1D (Š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 (R^{2}) are used to quantitatively evaluate the misfit between the simulated results of the coupled model and reference values, which are calculated as
where the subscript i represents the serial number of the results; x represents the total number of the results; y_{i} is the simulated result of the coupled model and Y_{i} is the reference result; $\stackrel{\mathrm{\u203e}}{y}$ is the average simulated result and $\stackrel{\mathrm{\u203e}}{Y}$ is the average reference result.
3.1 Two test cases
3.1.1 Case 1: 1D 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 HYDRUS1D 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 HYDRUS1D were used as the reference of this test case. The stress period ΔT was set as 5 d, and the MODFLOW time step Δt_{s} and the UBMOD time step Δt_{u} were both set to be 1 d. The spatial discretization was 0.1 m.
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: twodimensional (2D) water table recharge experiment
This test case was used for model validation in a 2D 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 2D 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 2D domain is a rectangular sandy soil slab of $\mathrm{6.0}\times \mathrm{2.0}\times \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 noflux boundary. Because of the symmetry of the flow system, only onehalf 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. Noflow 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 subareas defined to represent the unsaturated zone, which were numbered from left to right. The first and last subareas covered 0.2 and 0.4 m in the x direction, respectively, and each of the remaining subareas covered 0.3 m in the x direction. The first and second subareas were used to define the recharge boundary, while the other subareas were used to define the norecharge boundary. The stress period ΔT was set as 1 h, and the initial MODFLOW time step Δt_{s} and UBMOD time steps Δt_{u} 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.
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 HYDRUS1D 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 HYDRUS1D. The ARE, RMSE, IA, and R^{2} 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 R^{2} were 2.2 %, 0.008 cm^{3} 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 HYDRUS1D. In general, these results indicate that the coupled model can capture the flow information under the upward flux and the heterogeneous condition.
The deviations of groundwater table depth and soil water content from the coupled model and HYDRUS1D can also be observed in Fig. 5. The deviations are caused by the different model structures of the coupled model and HYDRUS1D. HYDRUS1D 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 parameterspecific 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 R^{2} 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.
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 HYDRUS1D 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 R^{2} 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 HYDRUS1D 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 R^{2} 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 HYDRUS1D 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 quasi3D 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 subareas should move in the lateral direction, instead of moving downward to the saturated zone as in the quasi3D 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 subareas 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 quasi3D 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.
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 HYDRUS1D 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 CPUbased 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 HYDRUS1D. 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 2D 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.
4.1 Study site and input data
The coupled model was used to calculate the regionalscale groundwater recharge in a realworld 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 (107^{∘}37^{′}19^{′′}–108^{∘}51^{′}04^{′′} E, 40^{∘}45^{′}57^{′′}–41^{∘}17^{′}58^{′′} 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 km^{2}. The ground surface elevation decreases from 1028.9 to 1025.4 m from the southwest to the northeast. A 2year 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 welldefined 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 firstkind boundary condition was applied to the two segments. The nonflow 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 subareas 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 subareas. 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 ET_{0} was calculated by the measured evaporation data from the 20 cm pan (ET_{20}), multiplying by an empirical conversion coefficient. The empirical coefficient is 0.55, which was recommended by Hao (2016) by comparing monthly ET_{0} and ET_{20} with 8 years' data in this area. The ET_{0} during 2014 to 2016 was calculated by using the Penman–Monteith equation. The precipitation, irrigation, and ET_{0} are shown in Fig. 9. The crop growing season is from May to October, and the rest of the months are nocrop 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 Δt_{s} and UBMOD time step Δt_{u} were set as 1 d.
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 R^{2} 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.
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 shortterm results indicate the validity of the simulation results.
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 cm^{3} cm^{−3}, the IA values 0.621–0.775, and the R^{2} values 0.54–0.689. The corresponding values for the bared soil are 10.8 %–19.8 %, 0.038–0.052 cm^{3} 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.
Note: P is the annual precipitation; I is the irrigation water; R is the annual recharge; and R_{c} is the recharge coefficient, ${R}_{\mathrm{c}}=R/(P+I)$.
The computational cost of the realworld 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 nocrop 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 nocrop 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 nocrop 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.
This study developed a new quasi3D coupled model for the purpose of practical modeling of unsaturated–saturated flow at the regional scale. The UBMOD 1D 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 postprocessing tools. The model was evaluated by using both synthetic numerical examples and realworld experimental data. The major conclusions drawn from this research are as follows.

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.

The satisfactory results in the two testing examples demonstrate the effectiveness of the new quasi3D model with an acceptable calculative efficiency and wellmaintained saturated zone and unsaturated zone mass balance.

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.

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.

The coupled model should not be used for problems with substantial lateral flow in the unsaturated zone because of the quasi3D assumptions used in the model.

The coupled model could lead to a higher groundwater table depth since it ignores the capillary fringe.
All the data and codes used in this study can be requested by email to the corresponding author Yan Zhu at zyan0701@163.com.
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
where i indicates the vertical soil layer, i=1, …, j; q_{i} is the amount of allocated water per unit area of layer i (L); M_{i} is the thickness of layer i (L); θ_{i} is the initial soil water content of layer i (L^{3} L^{−3}); θ_{s,i} is the saturated soil water content of layer i (L^{3} 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
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.
Thirdly, the source/sink terms are used to account for soil evaporation and crop transpiration. The governing equation is as follows:
where W is the source/sink term (T^{−1}). The Penman–Monteith formula and Beer's law (also known as a Ritchietype equation) are adopted in UBMOD to estimate the potential soil evaporation E_{p} and potential crop transpiration T_{p}. Then E_{p} and T_{p} 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 E_{p} and T_{p} with the soil water stress coefficient.
Lastly, we calculate the diffusive movement driven by the matric potential. The governing equation is
where D(θ) is the hydraulic diffusivity (L^{2} T^{−1}), and $D\left(\mathit{\theta}\right)=K\left(\mathit{\theta}\right)\times \partial \phantom{\rule{0.125em}{0ex}}h\left(\right)open="/">\partial \mathit{\theta}$, where h is the matric potential (L). The finitedifference 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
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 K_{s} as
By following the steps above, the hydraulic diffusivity D(θ) of a specific soil type can be estimated with three physical meaning parameters (K_{s}, θ_{s}, and θ_{r}).
Note: the parameters K_{s}, θ_{s}, and θ_{f} are the saturated hydraulic conductivity (L T^{−1}), the saturated water content (L^{3} L^{−3}), and the field capacity (L^{3} L^{−3}); θ_{w} is the soil water content at wilting point (L^{3} L^{−3}); α in the exponential equation is a sitespecific parameter determined mainly from soil characteristics, and it has an inverse relationship with K_{s}, 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
where α and n are two parameters in the van Genuchten model. The two terms $\partial \mathit{\theta}\left./\right)\partial {\mathit{\theta}}_{\mathrm{s}}$ and $\partial \mathit{\theta}\left(\right)open="/">\partial {\mathit{\theta}}_{\mathrm{r}}$ can be easily calculated. The value of $\partial \mathit{\theta}\left(\right)open="/">\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 K_{s}. The specific equations are shown as follows:
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.
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 realworld application. All the authors contributed to model validation. WM prepared the manuscript with contributions from all the coauthors.
The authors declare that they have no conflict of interest.
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.
This research has been supported by the Natura Science Foundation of China (grant nos. 51790532, 51779178, and 51629901).
This paper was edited by Philippe Ackerer and reviewed by Gerrit H. de Rooij and one anonymous referee.
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.: CERESmaize 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/S03783774(00)000925, 2000.
Brunner, P. and Simmons, C.: HydroGeoSphere: a fully integrated, physically based hydrological model, Groundwater, 50, 170–176, https://doi.org/10.1111/j.17456584.2011.00882.x, 2012.
Clement, T., Wise, W., and Molz, F.: A physically based, twodimensional, finitedifference algorithm for modeling variably saturated flow, J. Hydrol., 161, 71–90, https://doi.org/10.1016/00221694(94)90121X, 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 twodimensional watershedscale 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 distributedparameter watershed model to estimate net infiltration below the zone, U. S. Geological Survey, Virginia, US, 2008.
Furman, A.: Modeling coupled surfacesubsurface 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 watersalt 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.: MODFLOW2005, the U.S. Geological Survey modular groundwater model – the GroundWater 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.: APSIMEvolution 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 DarziNaftchali, A.: Application of virtual water trade to evaluate cropping pattern in arid regions, Water Resour. Manage., 29, 4061–4074, https://doi.org/10.1007/s1126901510454, 2015.
Kendy, E., GérardMarchant, P., Walter, M. T., Zhang, Y. Q., Liu, C. M., and Steenhuis, T. S.: A soilwaterbalance 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 3D groundwater flow model in a semiarid region using MODFLOW and GIS tools: The Zéramdine–Béni Hassen Miocene aquifer system (eastcentral 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.: GSFLOWCoupled Groundwater and Surfacewater FLOW model based on the integration of the PrecipitationRunoff Modeling System (PRMS) and the Modular GroundWater Flow Model (MODFLOW2005), 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/S13648152(03)000781, 2003.
Niswonger, R., Prudic, D., and Regan, R.: Documentation of the UnsaturatedZone Flow (UZF1) Package for modeling unsaturated flow between the land surface and the water table with MODFLOW2005, 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 waterheatsolute 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 MODFLOW2000, the us geological survey modular groundwater model, IGWMCInternational Ground Water Modeling Center, US, 2007.
Shen, C. and Phanikumar, M.: A processbased, distributed hydrologic model based on a largescale method for surfacesubsurface 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 twodimensional variably saturated media, Research Report, California, US, 1994.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., and van Genuchten, M. T.: The HYDRUS1D 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/s1004000403856, 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łaKawęcka, A., Šimůnek, J., Leterme, B., Beegum, S., JaworskaSzulc, B., PruszkowskaCaceres, M., GorczewskaLangner, W., AnguloJaramillo, 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/johh20180005, 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 6A18, Virginia, US, 2006.
Tian, Y., Zheng, Y., Wu, B., Wu, X., Liu, J., and Zheng, C.: Modeling surface watergroundwater interaction in arid and semiarid 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 ArcGISGeospatial 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 HYDRUSbased 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 groundwater 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 steadystate 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.: HydrologicResponse simulations for the R5 catchment with a comprehensive physicsbased 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, twodimensional unsaturatedsaturated 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 APSIMSoilWat 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 watersaving 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/s1126901197932, 2011.
Xu, X., Huang, G., Zhan, H., Qu, Z., and Huang, Q.: Integration of SWAP and MODFLOW2000 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 threedimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi twodimensional case, Adv. Water Resour., 21, 679–689, https://doi.org/10.1016/S03091708(97)000316, 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 shortterm 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 unsaturatedsaturated water flow, J. Hydrol., 475, 188–203, https://doi.org/10.1016/j.jhydrol.2012.09.048, 2012.