Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamicsin shallow water table areas

summary

Reasonable estimates of groundwater recharge and discharge through evapotranspiration is critical for sustainable water resources management in shallow water table areas. The hydrologic processes are highly interactive between the vadose zone and groundwater under shallow water table conditions. In traditional groundwater flow models, the recharge and evapotranspiration fluxes are often oversimplified as a simple sink/source term. However, the recharge and evapotranspiration are observed to vary with topography, soil type, land use, and water management practices. Additionally, they are known to vary temporally and spatially and are difficult to estimate, especially in arid and semi-arid regions. Thus, it is important to devise an appropriate method to estimate the recharge and evapotranspiration fluxes in groundwater modeling. In this study, a Soil–Water–Atmosphere–Plant (SWAP) package was integrated into a groundwater flow model (MODFLOW) in such a way that the SWAP package calculates vertical flux for MODFLOW, while MODFLOW provides averaged water table depth to determine the bottom boundary condition for SWAP zones. The SWAP zones in MODFLOW are derived from a combination of topology, soil type, land use, water management practices using geographic information systems (GIS). Then the MODFLOW with SWAP package was tested using a two-dimensional saturated–unsaturated water table recharge experiment. Results showed that the simulated water table elevations matched well with the observed ones except at the early period during which they were slightly higher than the observed ones, probably due to neglecting lateral diffusion in the unsaturated zone. Finally, we applied MODFLOW with SWAP package to simulate a regional groundwater flow problem in Hetao Irrigation District, upper Yel-low River basin of North China. The simulation results validated the applicability of the developed MOD-FLOW with SWAP package for practical regional groundwater modeling.

        合理估计通过蒸发蒸腾进行的地下水补给和排放对于浅表地区的可持续水资源管理至关重要。在浅地下水位条件下，渗流带和地下水之间的水文过程具有高度的相互作用。在传统的地下水流量模型中，补给和蒸散流量通常过于简单地定义为一个简单的汇/源项。然而，据观察，补给和蒸散量随地形、土壤类型、土地利用和水管理实践而变化。此外，它们在时间和空间上都有变化，很难估计，尤其是在干旱和半干旱地区。因此，在地下水建模中，设计一种适当的方法来估计补给和蒸散流量是很重要的。在本研究中，土壤-水-大气-植物（SWAP）软件包被集成到地下水流动模型（MODFLOW）中，使得SWAP软件包计算MODFLOW的垂直流动，而MODFLOW提供平均地下水位深度，以确定SWAP区域的底部边界条件。MODFLOW中的SWAP区域是利用地理信息系统（GIS）结合拓扑结构、土壤类型、土地利用和水管理实践得出的。然后，使用二维饱和-非饱和地下水位补给实验对带SWAP程序包的MODFLOW进行了测试。结果表明，模拟的地下水位高程与观测的水位高程匹配良好，但早期的地下水位略高于观测的水位，这可能是由于忽略了非饱和带中的横向扩散。最后，我们将MODFLOW和SWAP软件包应用于华北叶洛河上游河套灌区的区域地下水流动问题模拟。模拟结果验证了所开发的MOD-FLOW with SWAP软件包在实际区域地下水建模中的适用性。

1. 简介

Groundwater recharge and evapotranspiration (discharged by direct evaporation and crop root uptake through capillary rise)through the vadose zone are of great importance for sustainable groundwater use and control of salinity and water-logging in arid
and semi-arid regions with shallow water tables (Lerner et al., 1990; Arnold et al., 1993). Salinity has affected at least 20% of
the world’s arable land and more than 40% of all irrigated land to various degrees (Rhoades and Loveday, 1990), especially in arid and semi-arid regions such as North China (Xu et al., 2010) and Southwest Australia (Petheram et al., 2003). Thus, reasonable esti-mates of groundwater recharge and evapotranspiration is signifi-cant for sustainable groundwater management.

However, the groundwater recharge and evapotranspiration are influenced by a range of factors such as topography, soil type, land use, and water management practices (Petheram et al., 2003). The rates of re-charge and evapotranspiration are known to be the most difficult and uncertain components to estimate in groundwater budget, and they often vary spatially and temporally, especially in arid and semi-arid regions (Hendrickx and Walker, 1997; Sophocleous, 2004). The traditional physical methods for recharge and evapo-transpiration estimation such as direct measurement, water bal-ances, Darcian approaches and empirical models often require large in situ measurement data, and are still limited to small scales (Sophocleous, 2004). The tracer techniques are successfully used to estimate the recharge in arid and semi-arid regions, but they also only provide point or local scale information and do not directly quantify field scale recharge over a large area (Lerner et al., 1990). Recent developments in remote sensing techniques allow to map the groundwater recharge in a spatially distributed man-ner, as demonstrated by Brunner et al. (2004) or Tweed et al.(2007).

The groundwater recharge and evapotranspiration can also be estimated using physically-based numerical vadose zone model of which the principal advantages are that the models can allow reasonable forecast of the future recharge regimes under changing hydrological conditions such as changes in surface water manage-ment (Sophocleous, 2004; Allison et al., 1994). Most widely used groundwater flow models always oversimplify the estimates of groundwater recharge and evapotranspiration fluxes through va-dose zone for regional groundwater flow simulation. They often treat the groundwater recharge and evapotranspiration loss as simple source terms with no interface, which underestimates the effects of vadose zone on groundwater flow system. For example, in some traditional groundwater flow models such as FEFLOW (Diersch, 1996) and PLASM (Prickett and Lonnquist, 1971), these fluxes are often imposed by the modeler and estimated indepen-dently from the models when used for regional groundwater flow simulation.

In MODFLOW-2000 (Harbaugh et al., 2000), the re-charge rate is often directly specified by users as flux, and evapo-transpiration is treated as a head-dependant flux boundary, using a linear or a piecewise linear relationship between water table depth and evapotranspiration rate. Sometimes, one also uses a con-tinuous recharge–discharge function to represent groundwater discharge and recharge fluxes as a head-dependent continuous process (Doble et al., 2009). However, such an effort still may not appropriately estimate groundwater recharge and evapotranspira-tion for models due to the complex hydrological and environmen-tal conditions in the vadose zone.

Although some variably saturated models such as MODFLOW-SURFACT (HydroGeoLogic Inc., 1996), VSF package (Thoms et al., 2006), HydroGeoSphere (Therrien et al., 2007) and HYDRUS (2D/ 3D) (Šimu˚ nek et al., 2006) are able to describe the vadose zone water flow processes and estimate the recharge and dis-charge, they are all based on the use of three-dimensional Richards equations. Thus, it may impound heavily available computational resources to use these models due to the requirement of much finer temporal and spatial discretization (van Walsum and Groenendijk, 2008). In addition, some models do not reasonably consider the effects of vegetation growth and meteorological factors on groundwater recharge and evapotranspiration (such as VSF package), which have significant impacts on water flow in va-dose zone. While others (e.g., HYDRUS-2D/3D) are not appropriate for regional groundwater flow modeling due to the limitation of dealing with the boundary conditions and source/sink terms. Con-sidering the fact that water flow in vadose zone is usually assumed to be vertically one-dimensional (1-D) in field (Kroes and van Dam, 2003; Singh et al., 2006), therefore, it may be an alternative meth-od to couple a 1-D vadose zone flow model into groundwater flow models through the exchange of information between these two models. This method has also been successfully applied in some previous studies. For example, Hollanders et al. (2005) applied a loose coupling, while van Walsum and Groenendijk (2008) used an online computational scheme for coupling SWAP and MOD-FLOW. The loose coupling is relatively easy to be realized, but it cannot characterize the timely interaction between groundwater flow and vadose zone flow. Some researchers (Stoppelenburg et al., 2005; van Dam et al., 2008) have employed an offline cou-pling which can iteratively solve the groundwater flow equation and Richards equation until their solutions converge. It may be very expensive when the iteration is involved on the equation level (Shen and Phanikumar, 2010). In addition, their coupling is based on a specific assumption that a lower permeable layer underlies the unconfined aquifer, which is not always possible. Twarakavi et al. (2008) has fully coupled the simplified HYDRUS-1D and MODFLOW without considering the effects of vegetation and meteorological factors. Shen and Phanikumar (2010) have devel-oped a coupling method which lowers the dimensions of 3-D Rich-ards equation by separating it into 1-D equation for vadose zone flow and 2-D equation for groundwater flow.

Keeping the above in view, a 1-D vadose zone flow process, which comprehensively considers the effects of soil–plant–atmo-sphere continuum on groundwater dynamics, should be fully embedded into 3-D groundwater flow models. Therefore, the objectives of this study were to adapt the simplified Soil–Water–Atmosphere–Plant (SWAP) package from the original SWAP model (Kroes and van Dam, 2003) to simulate 1-D vertical flow in the vadose zone, and to couple the SWAP package with MODFLOW-2000 for appropriate estimation of groundwater recharge and evapotranspiration. They were coupled through the exchange of water table depth and net recharge, i.e., the difference between groundwater recharge and evapotranspiration. The coupled MODFLOW-2000 with the SWAP package was first tested using a well documented two-dimensional saturated–unsaturated water table recharge experiment (Vauclin et al., 1979). Then it was further validated by applying it to simulate the groundwater dynamics in an arid irrigation district on a regional scale. Consequently, this study was expected to extend the capabilities for MODFLOW to simulate groundwater dynamics.

        通过包气带的地下水补给和蒸发蒸腾（通过直接蒸发和通过毛细上升吸收作物根系而排放）对于地下水的可持续利用以及控制地下水位较浅的干旱和半干旱地区的盐度和水涝具有重要意义（Lerner et al.，1990；Arnold et al.，1993）。盐度在不同程度上影响了世界上至少20%的耕地和40%以上的灌溉土地（Rhoades和Loveday，1990），尤其是在中国北部（Xu et al.，2010）和澳大利亚西南部（Petheram et al.，2003）等干旱和半干旱地区。因此，合理估计地下水补给和蒸发蒸腾量对可持续地下水管理至关重要。

        然而，地下水补给和蒸发蒸腾受到一系列因素的影响，如地形、土壤类型、土地利用和水管理实践（Petheram et al.，2003）。众所周知，补给率和蒸发蒸腾率是地下水预算中最难估计和最不确定的组成部分，它们往往在空间和时间上有所不同，尤其是在干旱和半干旱地区（Hendrickx和Walker，1997；Sophocleous，2004年）。补给和蒸发蒸腾估算的传统物理方法，如直接测量、水量平衡、Darcian方法和经验模型，通常需要大量的现场测量数据，并且仍然局限于小尺度（Sophocleous，2004）。示踪技术被成功地用于估计干旱和半干旱地区的补给，但它们也只提供点或局部尺度的信息，而不能直接量化大面积的野外尺度补给（Lerner等人，1990）。正如Brunner等人（2004）或Tweed等人（2007）所证明的那样，遥感技术的最新发展使人们能够绘制空间分布的地下水补给图。

        地下水补给和蒸散也可以使用基于物理的数值渗流带模型进行估计，该模型的主要优点是，该模型可以在不断变化的水文条件下（如地表水管理的变化）合理预测未来的补给状况（Sophocleous，2004；Allison等人，1994年）。大多数广泛使用的地下水流量模型总是过于简化了区域地下水流量模拟中通过va剂量区的地下水补给和蒸发蒸腾流量的估计。他们通常将地下水补给和蒸发蒸腾损失视为没有界面的简单源项，这低估了渗流带对地下水流动系统的影响。例如，在一些传统的地下水流动模型中，如FEFLOW（Diersch，1996）和PLASM（Prickett和Lonnquist，1971），当用于区域地下水流动模拟时，这些流动通常是由建模者施加的，并从模型中独立估计。

        在MODFLOW-2000（Harbaugh等人，2000）中，用户通常直接将再充电率指定为流量，蒸发蒸腾作用被视为水头相关流量边界，使用地下水位深度和蒸发蒸腾率之间的线性或分段线性关系。有时，人们也会使用连续的补给-排泄函数来表示地下水的排泄和补给流量，将其视为水头相关的连续过程（Doble等人，2009）。然而，由于渗流带复杂的水文和环境条件，这种努力可能仍然不能适当地估计模型的地下水补给和排泄。

        尽管一些可变饱和模型，如MODFLOW-SURFACT（HydroGeoLogic股份有限公司，1996年）、VSF软件包（Thoms等人，2006年）、HydroGeoSphere（Therrien等人，2007年）和HYDRUS（2D/3D）（Šimu˚nek等人，2006）能够描述渗流带水流过程并估计补给和排泄，但它们都是基于三维Richards方程的使用。因此，由于需要更精细的时间和空间离散化，它可能会占用大量可用的计算资源来使用这些模型（van Walsum和Groenendijk，2008）。此外，一些模型没有合理考虑植被生长和气象因素对地下水补给和蒸发蒸腾的影响（如VSF包），这对va剂量区的水流有显著影响。而其他（例如，HYDRUS-2D/3D）由于处理边界条件和源/汇项的限制，不适用于区域地下水流动建模。考虑到渗流带中的水流通常被假设为垂直一维（1-D）（Kroes和van Dam，2003；Singh等人，2006），因此，通过这两个模型之间的信息交换，将一维渗流带水流模型耦合到地下水水流模型中可能是一种替代方法。该方法也已成功应用于先前的一些研究中。例如，Hollanders等人（2005）应用了松散耦合，而van Walsum和Groenendijk（2008）使用了耦合SWAP和MOD-FLOW的在线计算方案。松散耦合相对容易实现，但它不能表征地下水流动和渗流带流动之间的及时相互作用。一些研究人员（Stoplenburg等人，2005；van Dam等人，2008）采用了一种迭代法，可以迭代求解地下水流量方程和Richards方程，直到它们的解收敛。当在方程层面进行迭代时，这可能会非常昂贵（Shen和Phanikumar，2010）。此外，它们的耦合是基于一个特定的假设，即非承压含水层下方有一个较低的渗透层，这并不总是可能的。Twarakavi等人（2008）在不考虑植被和气象因素影响的情况下，将简化的HYDRUS-1D和MODFLOW完全耦合。Shen和Phanikumar（2010）提出了一种耦合方法，通过将三维Rich ards方程分为渗流带渗流的一维方程和地下水渗流的二维方程，降低了三维Rich ards方程的维数。

        考虑到上述情况，应将综合考虑土壤-植物-大气连续体对地下水动力学影响的一维渗流带渗流过程完全嵌入三维地下水渗流模型中。因此，本研究的目的是采用原始土壤-水-大气-植物（SWAP）模型（Kroes和van Dam，2003）中的简化土壤-水–大气-植物包来模拟渗流带中的一维垂直流动，并将SWAP包与MODFLOW-2000相结合，以适当估计地下水补给和蒸散。它们通过地下水位深度和净补给的交换而耦合，即地下水补给和蒸散之间的差异。首先使用有充分记录的二维饱和-非饱和地下水位补给实验对MODFLOW-2000与SWAP包的耦合进行了测试（Vauclin等人，1979）。然后将其应用于区域尺度的干旱灌溉区地下水动力学模拟，进一步验证了该方法的有效性。因此，本研究有望扩展MODFLOW模拟地下水动力学的能力。

2. Methodology

2.1. Description of MODFLOW-2000

MODFLOW is a computer program that numerically solves the three-dimensional groundwater flow equation for a porous med-ium using a finite-difference method (McDonald and Harbaugh, 1988) as follows:

MODFLOW是一个计算机程序，使用有限差分法（McDonald和Harbaugh，1988）数值求解多孔介质的三维地下水流动方程，如下所示：

where Kxx, Kyy and Kzz are values of principal hydraulic conductivity along the x, y, and z coordinate axes (m d􀀁1), respectively; h is the potentiometric head (m); W is a volumetric flux per unit volume representing sources and/or sinks of water (m3 d􀀁1); Ss is the spe-cific storage of the porous material (d􀀁1); and t is time (d).

The MODFLOW-2000 is an updated version of MODFLOW-88 and MODFLOW-96. It consists of the global, groundwater flow, observation, sensitivity and parameter-estimation processes (Harbaugh et al., 2000). All these processes are further divided into independent subroutines or modules. These modules are grouped into packages which deal with specific aspects of the simulation. The modular structure of MODFLOW facilitates the development of other packages to freely handle the impacts of boundary condi-tions and sinks/sources. Visual MODFLOW which includes MOD-FLOW-2000, is a commonly used MODFLOW pre/post processor (Waterloo Hydrogeologic Inc., 2006), whose version 4.2 was em-ployed in this study.
In MODFLOW, the Recharge (RCH) and Evapotranspiration (EVT) or Evapotranspiration Segments (ETS) packages (RCH-EVT or RCH-ETS) are widely used to characterize the vadose zone flow processes (Xu et al., 2011). The recharge rate is often directly spec-ified by users as flux in the RCH package. While evapotranspiration is treated as a head-dependant flux boundary using a linear or several linear segments relationships between water table depth and evapotranspiration rate in MODFLOW-2000. The user-specified extinction depth and the maximum evapotranspiration rate (at ground surface) are used in simulation. However, ground-water recharge and evapotranspiration through the vadose zone have both been observed to vary with water table depth, soil type, land use and water management practices. Therefore, the current MODFLOW-2000 still may not appropriately describe the effect of the vadose zone on groundwater flow system, due to the complex field conditions. The applicability of above mentioned methods is apparently questionable for arid and semiarid regions where soil capillary pressures play a dominant role in vadose zone flow (Lerner et al., 1990; Hendrickx and Walker, 1997).
2.2. SWAP package used for MODFLOW-2000
The groundwater recharge and evapotranspiration depend on topography, climate, soil, land use and water management practices which vary spatially and temporally on the regional scale. The SWAP model is a 1-D eco-hydrological model for simulating water flow, salt and heat transport in close interaction with crop growth in the vadose zone, which has been widely used for irrigation water management, salinity control and crop production prediction, etc.(Kroes and van Dam, 2003; Singh et al., 2006; Droogers et al., 2000). The SWAP package used here was simplified with removal of solute and heat transport modules, and then this package was adapted to simulate 1-D vertical transient saturated–unsaturated flow with considering the impacts of irrigation, rainfall, soil evaporation, runoff, drainage/sub-infiltration and crops root water uptake. The hysteresis, scaling of soil hydraulic properties, preferen-tial flow and mobile/immobile flow were not considered in current version of SWAP package. When using the SWAP package, the sim-ulation domain was divided into sub-regions (i.e., SWAP zones) analogous to the zones used in RCH, EVT or ETS packages of MODFLOW-2000 (Xu et al., 2011). The approach used for defining the land use systems (LUS) (FAO, 1976) was adopted to determine the SWAP zones. Some results also indicated that the topography has significant influences on the estimation of groundwater evapotranspiration (Li et al., 2008). In this study, each SWAP zone was denoted as a spatially homogeneous region, which was ob-tained through combination of topography, soil types, land use, irrigation, climatic conditions and water table depth conditions (Fig. 1). The SWAP zones in polygon format can be constructed by overlapping the vector maps using ArcGIS software, and subse-quently transforming the maps into an ASCII format matrix for MODFLOW-2000 use as described by Xu et al. (2009).

The calculation of vadose zone water flow was applied to each SWAP zone when using the SWAP package. In each SWAP zone, the water movement is described by a vertically 1-D Richards equation (see Eq. (2)) which is subsequently solved using an impli-cit finite-difference scheme.

where C is the differential soil water capacity (cm􀀁1), h is the soil water pressure head (cm), t is time (d), z is the vertical coordinate (cm, positive upward), K is the hydraulic conductivity (cm d􀀁1) and Sa is the soil water extraction rate by plant roots (cm3 cm􀀁3 d􀀁1). The sink term Sa(z) refers to water stress described by a reduction function proposed by Feddes et al. (1978).
The soil hydraulic properties can be described using van Genuchten (1980) and Mualem (1976) functions, respectively:

where hr is the residual water content (cm3 cm􀀁3), hs is the satu-rated water content (cm3 cm􀀁3), h is the actual soil water content (cm3 cm􀀁3), a (cm􀀁1) and n (dimensionless) are empirical shape fac-tors, Ks is the saturated hydraulic conductivity (cm d􀀁1), Se =(h 􀀁 hr)/(hs 􀀁 hr) is the relative saturation, and k is a shape parameter (dimensionless). A common method, which is the impli-cit, backward, finite-difference scheme with explicit linearization as described by Haverkamp et al. (1977) and Belmans et al. (1983),is adopted to solve Richards equation. In addition, adaptations relative to treatments of differential water capacity and hydraulic conduc-tivity functions are made to obtain the numerical scheme by Kroes and van Dam (2003). The soil profile is vertically discretized into compartments with various thickness i.e., a 1-D finite-difference grid. A finer discretization should be used at locations where sharp pressure head gradients are expected, and this can ensure the con-vergence of the numerical solution (Kroes and van Dam, 2003; Twarakavi et al., 2008). For example, smaller thickness of compart-ment are usually needed when close to the soil surface where rapid changes in water content and pressure head gradients often occur due to the effects of meteorological factors (Šimu˚ nek et al., 2005).
The top boundary condition can be determined by the actual evaporation and transpiration rates and the irrigation and precipi-tation fluxes. For appropriate estimation of the actual evaporation and transpiration rates, the potential evapotranspiration (ETp, mm d􀀁1) is estimated first using the Penman–Monteith equation (Allen et al., 1998):

where Dm is the slope of the vapor pressure curve (kPa C􀀁1), kw is the latent heat of vaporization (J kg􀀁1), Rn is the net radiation flux at the canopy surface (J m􀀁2 d􀀁1), G is the soil heat flux (J m􀀁2 d􀀁1), p1 accounts for time unit conversion (=86,400 s d􀀁1), qair is the air density (kg m􀀁3), Cair is the heat capacity of moist air (J kg􀀁1 C􀀁1), esat is the saturation vapor pressure (kPa), ea is the actual vapor pressure (kPa), cair is the psychrometric constant (kPa C􀀁1), rcrop is the crop resistance (s m􀀁1) and rair is the aerodynamic resistance (s m􀀁1). The ETp can also be directly calculated by the product of the user-specified reference evapotranspiration (ETref) and the crop fac-tors in absence of climatic or crop data. Potential evapotranspira-tion is then partitioned into potential soil evaporation and crop transpiration by using the leaf area index or soil cover fraction. Root water extraction at various depths in the root zone is calculated from the potential transpiration, root-length density and possible reduction from wet or dry conditions. In dry soil conditions, the maximum evaporation rate, Emax (cm d􀀁1), is calculated according to Darcy’s law (van Dam et al., 1997). As the actual soil evaporation may be overestimated using Darcy’s law, two empirical exponential functions of Black et al. (1969) and Boesten and Stroosnijder (1986) are alternatively applied to calculate the actual evaporation, Eaem (cm d􀀁1). These two empirical functions are both used to estimate the cumulative actual evaporation during a dry cycle and then to obtain the daily actual evaporation. Finally, SWAP determines the actual evaporation rate by taking the minimum value of Emax and Eaem.
Surface runoff is simulated when the infiltration capacity of the soil is not sufficient to infiltrate all the water. The excess water on the soil surface builds up as a ponded reservoir until the pond water level exceeds a certain threshold ponding level (hpond), which leads to the following surface runoff as:

where qrunoff is the surface runoff (cm d􀀁1), hpond is the ponding depth of surface water (cm), Zsill is the height of the sill which is equal to the maximum ponding height specified by user (cm), csill is the runoff/inundation resistance (d) and bsill is an exponent (–).
The SWAP package provided options to prescribe the lower boundary conditions of soil profile as either the Dirichlet type (gi-ven pressure head), or Neumann type (given flux). When a shallow water depth condition existed, the bottom boundary pressure head value of soil profile can be determined according to averaged water table depth of the SWAP zone as follows:

where hn,i is the pressure head at the bottom of the soil profile (cm), GWLi is the groundwater level (negative = below surface level, cm), zn,i is the position of bottom nodal point (negative, cm), hresis,i is the head difference between the groundwater level and hydraulic head of the bottom nodal point in the previous time step, the subscript i is the number of the SWAP zone and the subscript n represents the number of the bottom nodal point. If water table was deep and low-er than the bottom of soil profile, the free drainage boundary condi-tion at bottom of soil profile was provided for the SWAP zone.
The basic drainage module was used to describe the drainage /infiltration processes at the field scale in this study. The interac-tion between groundwater and up to five surface water systems, including drainage ditches and canals, can be simulated as follows:

where qdrain,i represents the drainage/infiltration (cm d􀀁1) to and from the surface water system i, the drainage base /drain,i is equal to the surface water level of the system i (cm below the soil surface), /gwl is the groundwater level (cm below the soil surface), and cdrain,i is the drainage or infiltration resistance of the system i (d).

2.3. Model coupling approach

Due to the rapid fluctuations of vadose zone water dynamics and slow movement of groundwater flow, the groundwater model-ing always has larger spatial and temporal scales than vadose zone modeling. Moreover, solving the non-linear Richards equation (Eq.(2)) often requires much shorter time steps compared to solving the groundwater flow equation (Eq. (1)). Thus, to couple the vadose water flow processes with groundwater flow models, one chal-lenge is to deal with these different temporal and spatial scale problems. Fig. 2 depicts the conceptualization of the coupled mod-el development, using three simulation units and three corre-sponding SWAP zones as an example. In MODFLOW, the modeling domain was discretized into a 3-D finite-difference grid. The flux components, including the vertical net recharge from the soil profile, groundwater abstraction, drainage, river recharge, leakage, etc., were imposed on the grid cells to solve the finite-dif-ference approximation of mass conservation equation. The SWAP zones were defined to represent a homogeneous sub-region of the vadose zone processes, and were specified to the correspond-ing MODFLOW grid. Each SWAP zone consists of one or more cells of MODFLOW grid (Fig. 2). It is feasible to define the SWAP zones on a cellwise basis, however, it is not always necessary for most practical situations because of the very complicated data process-ing and large computation. The soil profile for each SWAP zone was vertically discretized into a 1-D finite-difference grid, whose thick-ness was much smaller than the vertical layer thickness of the MODFLOW grid.

由于渗流带水动力学的快速波动和地下水流动的缓慢运动，地下水模型总是比渗流带模型具有更大的空间和时间尺度。此外，与求解地下水流量方程（方程（1））相比，求解非线性Richards方程（方程。（2））通常需要更短的时间步长。因此，为了将渗流过程与地下水渗流模型相结合，一个挑战是处理这些不同的时间和空间尺度问题。图2以三个模拟单元和三个相应的SWAP区域为例，描述了耦合模式开发的概念。在MODFLOW中，建模域被离散为三维有限差分网格。将流量分量，包括土壤层的垂直净补给、地下水抽取、排水、河流补给、渗漏等，施加在网格单元上，以求解质量守恒方程的有限差分近似值。SWAP区被定义为代表渗流区过程的均匀亚区，并被指定为相应的MODFLOW网格。每个SWAP区由一个或多个MODFLOW栅格单元组成（图2）。在单元基础上定义SWAP区域是可行的，但是，由于数据处理非常复杂，计算量大，在大多数实际情况下并不总是必要的。将每个SWAP区域的土壤文件垂直离散为一维差分网格，其厚度远小于MODFLOW网格的垂直层厚度。

A disparity exists with respect to temporal scales of vadose zone flow and groundwater flow. To deal with this issue, different time steps were used to solve the vadose zone flow and groundwater flow in the coupled model. Each MODFLOW time step consists of many SWAP time steps. As shown in Fig. 3, the vadose zone and groundwater flow models are interactive through the exchange of averaged water table depth and vertical net recharge flux and the calculated specific yield Sy if needed in each MODFLOW time step. The initial hydraulic heads of MODFLOW were used to calcu-late the averaged water table depth for each SWAP zone at the first MODFLOW time step of the first stress period. Initially, the bottom flux at the soil profile of each SWAP zone was calculated by the SWAP package, and was subsequently specified to the correspond-ing SWAP zone as a net recharge flux at start of each MODFLOW time step. The Sy value for phreatic aquifers is affected by soil water retention properties, water table depth and its change rate, etc. (Sophocleous, 1985), especially under shallow water table con-ditions. In the coupled model, three options were provided to determine the values of Sy when using the SWAP package. The val-ues of Sy can be specified directly by users as a constant according to the results of pumping or slug test. The Sy value can also be cal-culated for each SWAP zone at start of each new MODFLOW time step by using Eq. (9) (Stoppelenburg et al., 2005) or Eq. (10) (Duke, 1972; Said et al., 2005):

渗流带渗流和地下水渗流的时间尺度存在差异。为了解决这个问题，使用不同的时间步长来求解耦合模型中的渗流带流量和地下水流量。每个MODFLOW时间步长由许多SWAP时间步长组成。如图3所示，如果每个MODFLOW时间步长需要，渗流带和地下水流量模型通过平均地下水位深度和垂直净补给流量以及计算的特定流量Sy的交换而相互作用。在第一个应力期的第一个MODFLOW时间步长，使用MODFLOW的初始水头计算每个SWAP区的平均地下水位深度。最初，每个SWAP区土壤剖面的底部流量由SWAP程序包计算，随后指定为相应的SWAP区，作为每个MODFLOW时间步长开始时的净补给流量。潜水含水层的Sy值受土壤持水特性、地下水位深度及其变化率等的影响（Sophocleous，1985），尤其是在浅地下水位条件下。在耦合模型中，提供了三个选项来确定使用SWAP包时Sy的值。Sy的值可以由用户根据泵送或段塞试验的结果直接指定为常数。在每个新的MODFLOW时间步长开始时，也可以使用等式（9）（Stoplenburg等人，2005）或等式（10）（Duke，1972；Said等人，2005年）计算每个SWAP区域的Sy值：

where DZi is the thickness of the ith soil compartment (cm) and k is the total number of soil compartments above the water table; / is the porosity (cm3 cm􀀁3), Sr is the soil specific retention (cm3 cm􀀁3), s and ha are the pore size distribution index and the soil air entry pressure (cm) of Brooks and Corey water retention model (Brooks and Corey, 1966), respectively, and GWTD is the depth to water ta-ble (cm). The updated averaged water table depth for each SWAP zone was calculated by MODFLOW, and assigned to SWAP zones in order to determine the hydraulic head at the bottom of the soil profile for the next MODFLOW time step using Eq. (7). This coupled model ran iteratively through the exchange of the averaged water table depth and the net recharge flux in SWAP zones for each MOD-FLOW time step, as shown in Fig. 3. However, if the water table be-comes lower than the bottom elevation of the soil profile, the bottom boundary will be switched into a free drainage condition for the corresponding SWAP zone. The net recharge to water table can be optionally specified to the cells of the first layer or the upper active cells of MODFLOW by users.
Any length of time step can be selected by users as needed in MODFLOW. However, a 1-day time step was used for MODFLOW if Penman–Monteith equation was chosen for calculating ETp be-cause this equation is based on a daily calculation in the SWAP package.

The results of the vadose zone flow and groundwater flow were saved separately. MODFLOW-2000 reported the groundwater head and groundwater balance at the end of a user-specified time step.

Meanwhile, SWAP provided the moisture information of the soil profile for specified time. In this study, we have the option of using the SWAP package to replace the RCH-EVT or RCH-ETS packages of MODFLOW. However, for those sub-regions without enough soil information to use the SWAP package, the RCH, EVT or ETS pack-ages of MODFLOW can still be simultaneously used with the SWAP package. RCH package is also needed to simulate recharge in SWAP zones from other sinks or sources besides of the groundwater re-charge and evapotranspiration through the vadose zone. Such kinds of flexibility make this coupled model friendly to users for various applications.

3. Test study

3.1. Case study 1: 2-D water table recharge experiment

3.1.1. 2-D water table recharge experiment

Due to the lack of analytical solutions to such a coupled system, a 2-D water table recharge experiment, presented in detail by Vauclin et al. (1979), was chosen to test the applicability of the SWAP package. The dataset of this experiment has been widely used to test the variably saturated flow models (Clement et al., 1994), or the coupled saturated–unsaturated models with simpli-fied treatment of the unsaturated flow process (Thoms et al., 2006; Shen and Phanikumar, 2010; Twarakavi et al., 2008). The flow domain consists of a rectangular sandy soil slab of 6.0 m long, 2 m high and 0.05 m thick. The bottom of domain can be defined as the reference datum and the initial pressure head is 0.65 m above the domain bottom. At the soil surface, a constant flux of q = 3.55 m d􀀁1 was applied to over the central 1.0 m by 0.05 m area, while the rest soil surface was covered to prevent soil evaporation. Because of the symmetry of the flow system, only flow in one half of do-main (right side) with size of 3.0 m 􀀃 2.0 m needs to be modeled (Fig. 4). The setup of the coordinate system is shown in Fig. 4. No-flow boundaries were defined on the bottom and the left side of the modeled domain (water divide due to the symmetry of flow). The water table elevation on the right boundary of the flow domain is 0.65 m throughout the entire experiment. The soil hydraulic properties were obtained from Vauclin et al. (1979), and the parameter values for van Genuchten model are presented in Table 1 (Thoms et al., 2006).

由于缺乏这种耦合系统的解析解，Vauclin等人（1979）详细介绍了二维地下水位补给实验，以测试SWAP包的适用性。该实验的数据集已被广泛用于测试可变饱和流量模型（Clement等人，1994），或对不饱和流量过程进行简化处理的饱和-不饱和耦合模型（Thoms等人，2006；Shen和Phanikumar，2010；Twarakavi等人，2008）。流域由一块长6.0m、高2m、厚0.05m的矩形砂质土板组成。域底部可定义为参考基准，初始压头位于域底部上方0.65 m处。在土壤表面，q=3.55 m d的恒定流量􀀁1施用于中心1.0米乘0.05米的区域，同时覆盖其余土壤表面以防止土壤蒸发。由于流动系统的对称性，只有一半的do main（右侧）流动，尺寸为3.0 m􀀃 2.0 m需要建模（图4）。坐标系的设置如图4所示。在建模区域的底部和左侧定义了无流动边界（由于流动的对称性而造成的分水）。在整个实验过程中，流域右边界的地下水位高程为0.65 m。土壤水力特性由Vauclin等人获得。（1979），van Genuchten模型的参数值如表1所示（Thoms等人，2006）。

3.1.2. Model setup

In the simulation with MODFLOW, the flow domain was divided into 1 row, 30 columns and 1 layer. The grid was horizontally discretized into uniform rectangular cells of 5 cm height and 10 cm width, and the thickness of the layer was assigned to 200 cm. It was assumed that the domain is homogeneous and isotropic. The hydrogeological parameters were obtained from Vauclin et al. (1979). The hydraulic conductivity of 840 cm d􀀁1 was used, and the value of Sy was calculated using Eq. (9) in each MODFLOW time step.
The simulation period lasted 8 h with one stress period and the time step of 2 min. The bottom of domain was used as the refer-ence datum in MODFLOW. The initial hydraulic heads are 65 cm for all cells. The SWAP package was used to define the upper boundary conditions. Eleven SWAP zones were defined considering the recharge zones (zones 2 and 3) and non-recharge zones (zones 4–11) (Fig. 4), while zone 1 was used to define the zone with zero net recharge. The right side of flow domain was assigned as a 65 cm constant-head boundary. No-flow boundaries were defined on the bottom and the left side of the modeled domain in MOD-FLOW (Fig. 4).
With the SWAP package, identical homogeneous soil profile was considered for all SWAP zones with the same values of van Genuchten parameters as shown in Table 1. The recharge flux was simulated as the precipitation since the reference evaportran-spiration rate was assigned to zero. The bottom boundary (i.e., the hydraulic head at the bottom of the soil profile) was determined using the averaged water table depth of corresponding SWAP zone. The initial conditions of soil moisture were assigned to each com-partment of soil profile according to the measurements of Vauclin et al. (1979).

3.1.3. Simulation results

The simulated water table using MODFLOW-2000 with the SWAP package was compared with the measured water table as shown in Fig. 5. It was found that the simulated water tables at 3, 4 and 8 h matched well with the observed water tables; how-ever, the model obviously overestimates the water table at 2 h. This overestimation may result from that the model does not ac-counts for the horizontal water flow in the unsaturated zone, which can be significant during the early period due to the rela-tively low initial soil moisture close to land surface. Lateral mois-ture diffusion caused the inflow from recharge zone to be redistributed horizontally in the vadose zone before it reached the water table, and the unsaturated zone can store a portion of the inflow. However, our model assumes that lateral redistribution of water occurs only after water has entered into the saturated zone. Therefore, our model overestimated the net groundwater recharge and subsequently resulted in a higher water table at the early stage for this particularly test problem, as also evident in Fig. 5. However, as the water table rose, more groundwater will be drained out through the constant head on the right side. The flow system then gradually approached to the steady-state condi-tion. Subsequently, the simulated water tables agreed well with the observed ones at 3, 4 and 8 h. Especially at 8 h, the simulated and observed water tables were in a good agreement. These phenomena are quite similar to the result obtained by Shen and Phanikumar (2010).
The observed and simulated distributions of soil water content in soil profile for initial and 8 h are presented in Fig. 6. The simu-lated and observed soil water content agreed well at x =20cm where the vertical flow was dominant in the unsaturated zones. But neglect of lateral moisture diffusion may result in a slightly higher simulated water moisture than the observed one at 8 h (Fig. 6a). Much closer agreements were observed at x = 140 and 200 cm where the lateral water flows were so small and nearly negligible, which may be due to that this region was away from the recharge area (x = 0–50 cm) (Fig. 6d and e). However, the sim-ulated water moisture was apparently lower than the observed one in the unsaturated zone at locations of x = 60 and 80 cm, which was closer to the recharge area of (x = 0–50 cm) (Fig. 6b and c). It was known that the effect of lateral diffusion effect in vadose zone may be weaker in actual regional flow problems compared to this experiment, and the vadose zone water flow can be assumed to mainly occur in the vertical direction for many irrigated areas such as described by Singh et al. (2006) and Droogers et al. (2000). Therefore, the coupled model in this study may result in less error due to the ignorance of later diffusion effect for most of the prac-tical problems.

3.2. Case study 2: Regional groundwater flow simulation

3.2.1. Study area

In order to verify the proposed MODFLOW with SWAP package for real applications on a regional scale, the Yonglian Irrigation Sys-tem (YLIS) of Hetao Irrigation District (Hetao) in North China was selected as a study area in this case (Fig. 7). It has a typically arid continental climate with shallow water tables caused by improper agricultural irrigation practices. The mean annual precipitation is only 169 mm, with most of the rainfall occurring from July to Sep-tember in the YLIS. The mean annual temperature is 6 C, with the lowest and highest monthly averages being 􀀁13 C and 22 Cin January and July, respectively. The mean annual evaporation is about 2000 mm. Soils in the southern part include alluvial silt sed-iments with textures of sandy loam types, while in the northern part they are finely textured, such as silt loam or clay loam.
The YLIS is located east of the Naiyong sub-main drainage ditch and west of the Yongshen sub-main drainage ditch (see Fig. 7). It is bounded by the sixth sub-main drainage ditch in the north, and by the south part of Yongshen sub-main drainage ditch and the upper reaches of the Zaohuo sub-main canal on the south border. The topography is flat, with a gradient of 1/4000 from south to north.
Due to the special climatic condition in the region, irrigation is essential during the entire crop growing season. All irrigation water is diverted from the Yellow River, while surface basin irriga-tion remains the major irrigation method. Canal seepage and field percolation cause a high water table in Hetao. Subsequently, the amount of groundwater evapotranspiration is very large, which results in severe problems of soil secondary saline-alkalization (IWC-IM, 1999).
The aquifer geometry and hydrogeological parameter data were obtained from boreholes and pumping tests (Wu, 2007; IWC-IM, 1999). This region is underlain by Quaternary sediments, mainly lake sediments and alluvial deposits from the Yellow River. The unconfined aquifer is composed of two water-bearing strata (Q4 and Q3) with respective thicknesses of 5–7 m and 39–45 m. The Q4 stratum is mainly sandy loam, and becomes finer from south to north, while the Q3 stratum consists of fine-medium sand, fine sand and silt sand with clay inter-layers. The Q4 stratum has a low-er permeability and the Q3 stratum has a relatively high perme-
ability. The upper layer of Q2 (Q22), located below the Q3 and Q4, is stable muddy clay, acting as a low permeable aquitard (Wu, 2007; IWC-IM, 1999).
The flow of groundwater moves from south to north according to hydraulic head measurements at 10 observation wells (see Fig. 7). Hydrogeologists from the Hetao Administration Office mon-itor the water table manually once every 5 days. The hydraulic gra-dient of the water table is about 1/3000–1/5000, which is approximately the same as the topographic slope (IWC-IM, 1999).

3.2.2. Model setup

The data sets used for model construction in this section were primarily obtained from several previous studies (Wu, 2007; IWC-IM, 1999; Bameng Survey, 1994) and from field survey. The aquifer system was divided into two sub-aquifers, with the first sub-aquifer on top of the second sub-aquifer. The first sub-aquifer has one layer. The second sub-aquifer was assumed to have uni-form hydrogeological parameter values and was divided into three layers with uniform thickness, due to the relatively larger thick-ness compared to the first layer. Thus, the aquifer system was ver-tically discretized into four layers. A uniform rectangular grid with 100 m square cells was used. The first layer, representing the Q4 stratum, is a low-permeable layer with a thickness of 5–7 m, while layers 2–4 consist of Q3 stratum and have higher permeability, with thicknesses varying from 13 to 15 m from south to north.
It is assumed that the aquifer is horizontally isotropic with its horizontal hydraulic conductivity being five times larger than the vertical one. The hydrogeological parameters were available from the previous studies (Wu, 2007; IWC-IM, 1999). The horizontal hydraulic conductivity for the first layer varies between 0.6 and 1.0 m d􀀁1 in the north–south (N–S) direction. The value of Sy was estimated using Eq. (10) for the phreatic aquifer – layer 1. The val-ues of (/-Sr), s, ha are respectively equal to 0.09, 0.43, 23 cm for south area with sandy loam and 0.084, 0.39, 25 cm for north area with loam soil. Layers 2–4 have the same hydraulic conductivity, which range from 8 to 10 m d􀀁1 along the N–S direction; their cor-responding Sy equals to 0.20 and their specific storage is 10􀀁6 m􀀁1. The simulation period lasted from May 1, 2004 to October 31, 2004 with 1-day time step. Initial groundwater heads were obtained through interpolation of observed groundwater levels from 10 wells on May 1, 2004.
The drainage ditches were considered as drain boundaries using the Drainage (DRN) package in MODFLOW for the first layer on the north, west and east borders. In the south, the upper part of the Yongshen sub-main drainage ditch and the upper reaches of the Zaohuo canal were simulated using a drain boundary and a flux boundary for the first layer, respectively. No-flow boundary condi-tions were defined for the lower layers. The recharge due to canal seepage was a calculated average, assigned to the canal control re-gion, and computed with the RCH Package in MODFLOW.

Nineteen SWAP zones were defined for the study area through the combination of soil type, land use, irrigation system and depth to water table for vadose zone flow modeling (Fig. 8). Two soil types including sandy loam and loam were considered, corre-sponding to zones 2–8 in the south and zones 9–19 in the north, respectively. These two soil types were represented by two soil profiles with a thickness of 3 m for each. The soil profiles were di-vided into three horizontal layers and forty compartments in the SWAP modeling. The values of Mualem–van Genuchten model parameters (see in Table 2) for the horizontal layers of two soil profiles were determined according to previous studies (Wu, 2007).

The crops cultivated in the study area were simplified into two main crops: summer maize and sunflower. Maize grew in the south area (zones 2–8), while sunflower growth was in the north area (zones 9–19) (Fig. 8). The crop growth data, including devel-opment stage, leaf area, crop height and root depth were collected from Wu (2007). Other parameters required by the simple crop module were assigned according to Kroes and van Dam (2003).
The weather data obtained from the nearby Hangjinhouqi County weather station was used for simulation. The rainfall was obtained from the rain gauge measurements in the YLIS. The basic drainage module was used and its parameters were obtained from Wu (2007). The bottom boundary (i.e., the hydraulic head at the bottom of the soil profile) was determined using the water table depth.

3.2.3. Simulation results

Groundwater level data from 10 observation wells (Fig. 7) were used to evaluate the model performance. The root mean square er-ror (RMSE) and the model efficiency (EFF) were used as indicators of model fitting:

where GWobs and GWsim are respectively the ith value of the observed and calculated water table (i =1, 2, .. ., n), n is the total num-ber of observations, and GWmean is the average observed water table over the total number of observations. Good agreements between measured and simulated groundwater levels were achieved for the 10 selected observation wells (Fig. 9). The RMSE and EFF of groundwater level for a 5-day observation interval were 0.23 m and 0.93, respectively. Fig. 9 also shows that the simulated ground-water level agrees well with the observation. The observed and sim-ulated groundwater levels were compared through a linear regression forced to the origin, and the respective coefficients of regression and determination were b = 1.0 and R2 = 0.93. The simu-lated daily groundwater levels of four observation wells (B3, B4, B5 and B8) are selected and presented in Fig. 10. The observation wells B3&B4 and B5&B8 were respectively located in the south area with sandy loam and in the north area with loam soil. It showed that fluctuations of simulated groundwater levels agreed reasonably with the observation data. However, some discrepancy still exist between the simulated groundwater levels and the observations, e.g., the simulated groundwater levels at observation well B3 show slightly lower than the observed ones, while the simulated ground-water levels at B5 are higher than observed ones in June. This im-plies that uncertainties may be imposed by model parameters, model structure and model input.
The MODFLOW with SWAP package gave a reasonable descrip-tion of water table fluctuations responding to rainfall, irrigation and evapotranspiration, compared to artificially-designated recharge or discharge rates, as shown in Fig. 10. The groundwater level rose quickly after irrigation or rainfall, and declined due to direct evaporation and crop root uptake during May to late July. A continuous declining trend was observed due to direct evaporation and crop root uptake and less irrigation or rainfall from early August to late September (Fig. 11). The water table depth decreased from 0.6–1.0 m to more than 2.0 m during the same period. In October,the water table increased markedly, a result of the prescribed autumn irrigation process for leaching salt and storing water in the soil profile for crop use for the next agricultural year.
The value of actual evapotranspiration (ETa) in vadose zone was in a range of 453–536 mm from May 1st to September 20th for all SWAP zones, which is close to the results of Wu (2007) with a water balance method in YLIS. The temporal variation of ETa and net bottom flux (Qbot) in SWAP zone 13 is presented in Fig. 11.A positive Qbot represents upward groundwater flow to recharge the vadose zone, and a negative Qbot means downward flow to recharge groundwater from the vadose zone. Result showed that groundwater recharge was closely related to the irrigation and rainfall, while large groundwater capillary rise was contributed to crop water use during the crop growing period from June to middle August. Those imply the coupled model can reasonably reproduce the flow pattern in vadose zone.

4. Discussion

The two case studies showed that the coupled model can rea-sonably simulate groundwater and vadose zone water dynamics under the conditions that the vertical water flow is dominant in va-dose zone. The coupled model simulation can help us gain greater insights into the spatial and temporal pattern of groundwater re-charge and evapotranspiration, and fluctuation of water tables accounting for the spatial and temporal variations of topography, soil, land use and water management practices on a regional scale. However, in spite of these favorable advantages, it is necessary to know its limitations and improvements needed when using the SWAP package. Firstly, the current version of SWAP package does not take into account the effects of terrain changes (i.e., slope and aspect) on surface water movement and its interaction among SWAP zones. Therefore, it is not applicable for the regions with strong surface runoff such as hilly and mountain regions in humid and semi-humid regions. This may be solved by employing the concept of rainfall–runoff mode as used in the Soil and Water Assessment Tool (SWAT) model (Arnold et al., 1993) in further investigation. Secondly, the solute transport has not been em-ployed at this study. Thus, the current version of SWAP package is not suitable for regions which have severe salinization problems. In those regions, vegetation growth will be seriously affected by solute transport processes, which will significantly affect ground-water recharge and evapotranspiration pattern. Waterlogged areas in arid and semi-arid regions often have the salinization problems or the potential threat to become salinization over time (Petheram et al., 2003). Thus, integrating a solute transport module both for vadose zone and saturated zone such as MT3D (Zheng, 1990) into our developed model will be a follow-up subject to investigate. Thirdly, for regions with deep water tables (thick vadose zones), it is increasingly difficult to characterize the hydrological proper-ties of the unsaturated zones because of the heterogeneous nature of the media, which will make our model less reliable. Finally, the SWAP package of coupled model is more suitable for the areas with relatively small changes in depth to water table in space. Other-wise, a finer discretization is needed for local regions with greatly spatial variations of depth to water table, which may impound heavily on available computational resources or bring some uncer-tainties. An alternative method is to use RCH-ETS packages instead of SWAP package for estimating recharge and evportranspiration in local regions with highly spatial variations of depth to water table.

5. Summary and conclusion

We have presented the development of a SWAP package for the MODFLOW-2000 model to simulate the vadose zone flow pro-cesses and estimate the groundwater recharge and evapotranspira-tion for groundwater modeling in relation to the shallow water problems. The SWAP package was simplified from the original SWAP model by only considering the flow process and neglecting the solute and heat transport processes. This simplified SWAP package was adapted to describe the vertically 1-D vadose zone water movement. It considered the effects of infiltration, soil evap-oration, crop root uptake, soil moisture storage, drainage in field scale, crop growth, and surface runoff in the vadose zone. There-fore, it includes most important processes related to groundwater recharge and evapotranspiration estimation. The MODFLOW-2000 model was used to simulate three-dimensional groundwater flow, interacting with the SWAP package through an exchange of net re-charge flux and averaged water table depth in each SWAP zone. The MODFLOW-2000 coupled with the SWAP package was then tested using a two-dimensional saturated–unsaturated water table recharge experiment of Vauclin et al. (1979) and a regional ground-water flow simulation in an arid irrigation district of North China. The agreement of simulated and observed water table validates the practical applicability of this coupled model. Finally, the limitation and improvement needed for the current version of the coupled model were analyzed. This calls for a follow-up investigation on integration of the solute transport, surface water movement to the current model as well as a combination with the remote sens-ing technologies.
Acknowledgements
This research was jointly supported by the National Science Foundation (NSF) for Distinguished Young Scholars of China to Prof. Huang G and other three NSF Projects of China (Grant Numbers: 52039007, 51069006, 50979106), the 12th Five-year Research Program of the Chinese Ministry of Science and Technol-ogy, Chinese Universities Scientific Fund (2011JS122) and the Program for Changjiang Scholars and Innovative Research Team in University. A scholarship to Dr. Xu from China Scholarship Coun-
cil to support his one year stay at Texas A&M University was also
acknowledged. The constructive comments from the reviewers
are greatly appreciated.
References
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop Evapotranspiration:
Guidelines for Computing Crop Water Requirements. Irrigation and Drainage Paper 56, Rome, Italy.
Allison, G.B., Gee, G.W., Tyler, S.W., 1994. Vadose-zone techniques for estimating
groundwater recharge in arid and semiarid regions. Soil Sci. Soc. Am. J. 58, 6–14. Arnold, J.G., Allen, P.M., Bernhardt, G., 1993. A comprehensive surface-groundwater
flow model. J. Hydrol. 142, 47–69.
Bameng Survey, 1994. Researches on the Water Resources and Salinization Control System Management Model in Hetao Irrigation District of Inner Mongolia.
Bameng Hydrology and Water Resources Survey Teams of Inner Mongolia,
Hohhot (in Chinese).
Belmans, C., Wesseling, J.G., Feddes, R.A., 1983. Simulation model of the water
balance of a cropped soil: SWATRE. J. Hydrol. 63, 271–286.
Black, T.A., Gardner, W.R., Thurtell, G.W., 1969. The prediction of evaporation,
drainage, and soil water storage for a bare soil. Soil Sci. Soc. Am. J. 33 (5), 655–
660.
Boesten, J.J.T.I., Stroosnijder, L., 1986. Simple model for daily evaporation from
fallow tilled soil under spring conditions in a temperate climate. Neth. J. Agric.
Sci. 34, 75–90.
Brooks, R.H., Corey, A.T., 1966. Properties of porous media affecting fluid flow. J.
Irrig. Drain. Eng. 92 (2), 61–90.
Brunner, P., Bauer, P., Eugster, M., Kinzelbach, W., 2004. Using remote sensing to
regionalize local precipitation recharge rates obtained from the Chloride
Method. J. Hydrol. 294 (4), 241–250.
Clement, T.P., Wise, W.R., Molz, F.J., 1994. A physically based, two-dimensional,
finite-difference algorithm for modeling variably saturated flow. J. Hydrol. 161,
71–90.
Diersch, H.J.G., 1996. Interactive, Graphics-based Finite-element Simulation System
FEFLOW for Modeling Groundwater Flow, Contaminant Mass and Heat
Transport Processes, User’s Manual Version 4.50. WASY Institute for Water
Resources Planning and Systems Research Ltd., Berlin.
Doble, R.C., Simmons, C.T., Walker, G.R., 2009. Using MODFLOW 2000 to model ET
and recharge for shallow ground water problems. Ground Water 47 (1), 129–
135.
Droogers, P., Bastiaanssen, W.G.M., Beyazg, M., Kayam, Y., Kite, G.W., Murray-Rust,
H., 2000. Distributed agro-hydrological modeling of an irrigation system in
western Turkey. Agric. Water Manage. 43 (2), 183–202.
Duke, H.R., 1972. Capillary properties of soils—influence upon specific yield. Trans.
ASAE 15 (4), 688–691.
FAO, 1976. A Framework for Land Evaluation. FAO Soils Bulletin 32. Rome, Italy.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. Simulation Monographs, Pudoc, Wageningen, The Netherlands.
Harbaugh, A.W., Banta, E.R., Hill, M.C., McDonald, M.G., 2000. MODFLOW-2000, the
US Geological Survey Modular Groundwater Mode—User Guide to
Modularization Concepts and the Ground-water Flow Process. Open-File
Report 00-92, Reston, Virginia.
Haverkamp, R., Vauclin, M., Touma, P.J., Vachaud, G., 1977. A comparison of
numerical simulation models for one-dimensional infiltration. Soil Sci. Soc. Am.
J. 41, 285–294.
Hendrickx, J.M.H., Walker, G.R., 1997. Recharge from precipitation. In: Simmers, I.,
Balkema, A.A. (Eds.), Recharge of Phreatic Aquifers in (Semi-) Arid Areas.
Rotterdam, The Netherlands, pp. 19–111.
Hollanders, P., Schultz, B., Wang, S.L., Cai, L.G., 2005. Drainage and salinity
assessment in the Huinong Canal Irrigation District, Ningxia, China. Irrig.
Drain. 54 (2), 155–173.
HydroGeoLogic Inc., 1996. MODFLOW-SURFACT Software (Version 2.2) Overview:
Installation, Registration, and Running Procedures. HydroGeoLogic Inc.,
Herndon, Virginia.
IWC-IM, 1999. Construction and Rehabilitation Planning Project for Water-saving in
Hetao Irrigation District of the Yellow River Basin. Inner Mongolia. Institute of
Water Conservancy and Hydropower of Inner Mongolia, Hohhot (in Chinese).
Kroes, J.G., van Dam, J.C., 2003. Reference Manual SWAP version 3.03. Alterra-report
773. Alterra, Green World Research, Wageningen. ISSN 1566-7197.
Lerner, D.N., Issar, A.S., Simmers, I., 1990. Groundwater Recharge: A Guide to
Understanding and Estimating Natural Recharge. Report 8. Int. Assoc.
Hydrogeol., Kenilworth, UK.
Li, H.T., Kinzelbach, W., Brunner, P.A., Li, W.P., Dong, X.G., 2008. Topography
representation methods for improving evaporation simulation in groundwater
modelling. J. Hydrol. 356, 199–208.
McDonald, M.G., Harbaugh, A.W., 1988. A Modular Three-dimensional Finitedifference
Ground-water Flow Model. Techniques of Water-Resources
Investigations of US Geological Survey. Book 6, CH.A1, Reston, Virginia.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of
unsaturated porous media. Water Resour. Res. 12 (3), 513–522.
Petheram, C., Dawes, W., Grayson, R., Bradford, A., Walker, G., 2003. A sub-grid
representation of groundwater discharge using a one-dimensional groundwater
model. Hydrol. Process. 17, 2279–2295.
Prickett, T.A., Lonnquist, C.G., 1971. Selected Digital Computer Techniques for
Groundwater Resource Evaluation. Illinois State Water Survey, Bulletin 55.
Rhoades, J.D., Loveday, J., 1990. Salinity in irrigated agriculture. In: Stewart, B.A.,
Nielson, D.R. (Eds.), American Society of Civil Engineers, Irrigation of
Agricultural Crops. Am. Soc. Agronomists, Madison, WI, pp. 1089–1142.
Said, A., Nachabe, M., Ross, M., Vomacka, J., 2005. Methodology for estimating
specific yield in shallow water environment using continuous soil moisture
data. J. Irrig. Drain. Eng. 131, 533–538.
Shen, C., Phanikumar, M.S., 2010. A process-based, distributed hydrologic model
based on a large-scale method for surface-subsurface coupling. Adv. Water Res..
doi:10.1016/j.advwatres.2010.09.002.
Šimu˚ nek, J., van Genuchten, M.T., Šejna, M., 2005. The HYDRUS-1D Software
Package for Simulating the One-Dimensional Movement of Water, Heat, and
Multiple Solutes in Variably-Saturated Media, Version 3.0, HYDRUS Software
Series 1, Department of Environmental Sciences, University of California
Riverside, Riverside, California, USA, p. 240.
Šimu˚ nek, J., van Genuchten, M.T., Šejna, M., 2006. The HYDRUS Software Package for
Simulating Two-and Three-dimensional Movement of Water, Heat, and
Multiple Solutes in Variably-saturated Media. Version1.0. PC-Progress, Prague,
Czech Republic.
Singh, R., Kroes, J.G., van Dam, J.C., Feddes, R.A., 2006. Distributed ecohydrological
modelling to evaluate the performance of irrigation system in Sirsa district,
India: I. Current water management and its productivity. J. Hydrol. 329 (3–4),
692–713.
Sophocleous, M., 1985. The role of specific yield in ground-water recharge
estimations: a numerical study. Ground Water 23, 52–58.
Sophocleous, M., 2004. Groundwater Recharge. In: Silveira, Luis, Wohnlich, Stefan,
Usunoff, Eduardo J. (Eds.), Groundwater. Encyclopedia of Life Support Systems
(EOLSS), Developed under the Auspices of the UNESCO. Eolss Publishers, Oxford,
UK. <http://www.eolss.net>.
Stoppelenburg, F.J., Kovar, K., Pastoors, M.J.H., Tiktak, A., 2005. Modeling the
Interactions between Transient Saturated and Unsaturated Groundwater. Offline
Coupling of LGM and SWAP. RIVM In Report 500026001/2005, Bilthoven,
Netherlands.
Therrien, R., Mclaren, R.G., Sudicky, E.A., Panday, S.M., 2007. HydroGeoSphere: A
Three-dimensional Numerical Model Describing Fully-integrated Subsurface
and Surface Flow and Solute Transport.
Thoms, R.B., Johnson, R.L., Healy, R.W., 2006. User’s Guide to the Variably Saturated
Flow (VSF) Process for MODFLOW. Techniques and Methods 6-A18, USGS,
Reston, VA. <http://pubs.usgs.gov/tm/2006/tm6a18/> (verified 15.04.08).
Twarakavi, N.K.C., Šimu˚ nek, J., Seo, S., 2008. Evaluating interactions between
groundwater and vadose zone using the HYDRUS-based flow package for
MODFLOW. Vadose Zone J. 7 (2), 757–768.
Tweed, S., Leblanc, M., Webb, J., Lubczynski, M., 2007. Remote sensing and GIS for
mapping groundwater recharge and discharge areas in salinity prone
catchments, southeastern Australia. Hydrogeol. J. 15 (1), 75–96.
van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabat, P., van Walsum, P.E.V.,
Groenendijk, P., van Diepen, C.A., 1997. Theory of SWAP Version 2.0. Simulation
of Water Flow, Solute Transport and Plant Growth in the Soil–Water–
Atmosphere–Plant Environment. Report 71, Sub Department of Water
Resources, Wageningen University, Technical document 45. Alterra Green
World Research, Wageningen, The Netherlands.
van Dam, J.C., Groenendijk, P., Hendriks, R.F.A., Kroes, J.G., 2008. Advances of
modeling water flow in variably saturated soils with SWAP. Vadose Zone J. 7,
640–653.
van Genuchten, M.T.H., 1980. A closed-form equation for predicting the hydraulic
conductivity of unsaturated Soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Walsum, P.E.V., Groenendijk, P., 2008. Quasi steady-state simulation of the
unsaturated zone in groundwater modeling of lowland regions. Vadose Zone J. 7
(2), 769–781.
Vauclin, M., Khanji, J., Vachaud, G., 1979. Experimental and numerical study of a
transient, two-dimensional unsaturated-saturated water table recharge
problem. Water Resour. Res. 15 (5), 1089–1101.
Waterloo Hydrogeologic Inc., 2006. Visual MODFLOW Pro User’s Manual. Waterloo
Hydrogeologic Inc., Waterloo, Ontario, Canada.
Wu, X.N., 2007. Influence of restrictive irrigation on water balance of Hetao
Irrigation District in Inner Mongolia. PhD dissertation. Wuhan University,
Wuhan (in Chinese).
Xu, X., Huang, G.H., Qu, Z.Y., 2009. Integrating MODFLOW and GIS technologies for
assessing impacts of irrigation management and groundwater use in the Hetao
Irrigation District, Yellow River basin. Sci. China Ser. E – Technol. Sci. 52 (11),
3257–3263.
Xu, X., Huang, G.H., Qu, Z.Y., Pereira, L.S., 2010. Assessing the groundwater dynamics
and impacts of water saving in the Hetao Irrigation District, Yellow River basin.
Agric. Water Manage. 98 (2), 301–313.
Xu, X., Huang, G.H., Qu, Z.Y., Pereira, L.S., 2011. Using MODFLOW and GIS to assess
changes in groundwater dynamics in response to water saving measures in
irrigation districts of the upper Yellow River basin. Water Resour. Manage..
doi:10.1007/s11269-011-9793-2.
Zheng, C.M., 1990. MT3D: A Modular Three-Dimensional Transport Model for
Simulation of Advection, Dispersion and Chemical Reaction of Contaminants in
Groundwater Systems. US Environmental Protection Agency, Ada, Oklahoma.