Published online July 19, 2018

Technical Notes
Core Ideas• The coupling algorithm between
HYDRUS-1D and MODFLOW is
updated.
• Pressure head profiles are updated
based on the groundwater table andbottom flux after each time step.
• This eliminates sudden variationsin inflow and outflow fluxes withgroundwater table changes.

Updating the Coupling Algorithm
between HYDRUS and MODFLOW
in the HYDRUS Package for MODFLOW
Sahila Beegum,* Jiří Šimůnek, Adam Szymkiewicz,
K.P. Sudheer, and Indumathi M. Nambi
The HYDRUS-based flow package for MODFLOW (the HPM or the HYDRUS package) is an existing unsaturated zone flow package for MODFLOW. In MODFLOW
with the HPM, the groundwater modeling domain is discretized into regular gridsthat can be combined into multiple zones based on similarities in soil hydrology,
topographical characteristics, and the depth to the groundwater. Each of thesezones is assigned one unsaturated soil profile (the HPM profile). In this model,
after every MODFLOW time step, the flux at the bottom of the HPM profile isgiven as an input recharge flux to MODFLOW. MODFLOW simulates groundwater
flow, and the water table depth at the end of the MODFLOW time step is assigned
as the bottom boundary condition in the HPM profile. The current coupling algorithm assumes that the groundwater table in the HPM profile remains constantthroughout the entire MODFLOW time step. This results in unrealistic sudden
inflow and/or outflow fluxes at the bottom of the HPM profile after every time
step. The objective of this study was to develop a methodology to eliminate theerror in the determination of the recharge flux at the bottom of the HPM profile.
This was achieved by updating or modifying the pressure head profile in the HPM
profile after every MODFLOW time step. The effectiveness and the applicability of
the new coupling algorithm were evaluated using different case studies. The newcoupling algorithm is effective in eliminating unrealistic sudden variations in thebottom flux in the HPM profiles.
Abbreviations: 1D, one-dimensional; 3D, three-dimensional; HPM, HYDRUS Package for MODFLOW.

S. Beegum and J. Šimůnek, Dep. of Environmental Sciences, Univ. of California, Riverside,
CA 92521; A. Szymkiewicz, Faculty of Civil
and Environmental Engineering, Gdańsk
Univ. of Technology, 80-233 Gdańsk, Poland;
S. Beegum, K.P. Sudheer, and I.M. Nambi,
Dep. of Civil Engineering, Indian Instituteof Technology, Madras, Chennai 600 036,
India; K.P. Sudheer, Dep. of Agricultural andBiological Engineering, Purdue Univ., West
Lafayette, IN 47906. *Corresponding author(sahilabeegum666@gmail.com).
Received 20 Feb. 2018.
Accepted 28 May 2018.
Citation: Beegum, S., J. Šimůnek, A. Szymkiewicz, K.P. Sudheer, and I.M. Nambi. 2018.
Updating the coupling algorithm between
HYDRUS and MODFLOW in the HYDRUS
Package for MODFLOW. Vadose Zone J.
17:180034. doi:10.2136/vzj2018.02.0034© Soil Science Society of America.
This is an open access article distributed under
the CC BY-NC-ND license(http://creativecommons.org/licenses/by-ncnd/4.0/).

In the past few decades, the coupling of one-dimensional (1D) unsaturated zone
models with three-dimensional (3D) groundwater flow models to simulate various processes in and interactions between unsaturated and saturated soil zones has received a
lot of attention. This is mainly because of the computational complexity and requirements of fully 3D variably saturated flow models when modeling larger domains with
unsaturated soil zones (e.g., SHE model [Abbott et al., 1986], MODFLOW-SURFACT[HydroGeoLogic, 1996], FEFLOW [Diersch and Kolditz, 1998], TOUGH2 [Pruess etal., 1999], VSF [Thoms et al., 2006], HydroGeoSphere [Therrien et al., 2010], HYDRUS(2D/3D) [Šimůnek et al., 2016], MIN3P [Mayer et al., 2012], and PARFLOW [Maxwell
et al., 2016]).
MODFLOW (Harbaugh et al., 2000) is a widely accepted 3D groundwater flow
model. There have been many attempts to incorporate unsaturated zone flow models
into MODFLOW (e.g., Havard et al., 1995; Facchi et al., 2004; Niswonger et al., 2006;
Twarakavi et al., 2008; Lin et al., 2010; Zhu et al., 2011; Xu et al., 2012). The basic principle behind linking independent models for unsaturated and saturated soil zones is theexchange of information regarding the recharge flux from the unsaturated zone to thesaturated zone and the elevation of the water table at appropriate time steps.
The HYDRUS-based flow package for MODFLOW (further referred to as the
HPM or the HYDRUS package) was developed by Seo et al. (2007) and Twarakavi etal. (2008) to simultaneously evaluate transient water flow in unsaturated and saturated
zones. In this package, the subroutines from the computational module of HYDRUS-1D
Vadose Zone Journal | Advancing Critical Zone Science(Šimůnek et al., 2016) simulating unsaturated water flow in thevadose zone were coupled to MODFLOW (Harbaugh et al., 2000)
simulating saturated groundwater flow. The HYDRUS package
can represent the effects of unsaturated zone processes, such as
infiltration, evaporation, root water uptake, capillary rise, andrecharge in homogeneous or layered soil profiles. The coupledmodel is effective in addressing spatially variable saturated–unsaturated hydrological processes at the regional scale (Twarakavi et al.,2008). The most recent version of the HYDRUS package, whichis available to the public, is compatible with MODFLOW-2005(Harbaugh, 2005).

Governing Equations in MODFLOW
and the HYDRUS Package

Current Coupling Algorithm

Groundwater flow is modeled in MODFLOW by solving themass conservation equation using the finite difference approximation. The two-dimensional movement of groundwater of constantdensity in an unconfined aquifer is described by the partial differential equation (derived by applying the Dupuit assumption) as
Kx¶ æç ¶ H ö÷¶ æ ¶ H ö÷
¶H÷÷ = S y÷÷ + K y ççç H
ççè H
÷
¶x
¶x ø
¶ y çè ¶ y ø¶t[1]

where K x and Ky are the hydraulic conductivities [L T−1] in the x
and y directions, respectively, H is the piezometric head [L], Sy isthe specific yield (dimensionless) of the porous material, and t is
time [T].
The HYDRUS package simulates water flow in the unsaturated zone using the modified one-dimensional Richards equation:ù¶q ¶ é¶h=ê K ( h ) - K ( h )ú - Súû
¶ t ¶ z êë¶z[2]where q is the volumetric water content (dimensionless), h is the soil
water pressure head [L], t is time [T], z is the vertical coordinate [L],S is a sink term [T−1], and K(h) is the unsaturated hydraulic conductivity [L T−1]. The unsaturated hydraulic conductivity, K(h),and the water content, q(h), depend on the soil water pressure head.
This makes the Richards equation a highly nonlinear equation thatneeds to be solved numerically. The HPM permits the use of fivedifferent analytical models to describe the soil hydraulic properties(i.e., Brooks and Corey, 1964; van Genuchten, 1980; Vogel andCíslerova, 1988; Kosugi, 1996; Durner, 1994).

Spatial and Temporal Discretizationin the HYDRUS Package
In MODFLOW, the groundwater modeling domain is discretized into regular grids. These grids are combined into multiple
zones based on similarities in soil hydrology, depth to the groundwater, and topographical characteristics. Each of these zones is
assigned one vertical soil profile (the HPM profile) that extends
from the soil surface down to a depth, which is below the deepest possible water table level that can occur during the simulation.

VZJ | Advancing Critical Zone ScienceThe total simulation period in MODFLOW is divided into
stress periods (during which the external stresses such as pumping,recharge, or river stage remain constant), and each stress period isfurther divided into time steps. For each MODFLOW time step,
the HYDRUS package simulates water flow in the unsaturated
zone for each HPM profile. The HYDRUS package uses its own
time stepping algorithm, based on user-defined criteria. Since theflow conditions in the vadose zone vary more rapidly than in thesaturated zone, the HPM time steps are generally smaller than the
MODFLOW time steps. The numerical solution of the nonlinear
Richards equation also requires smaller time steps for simulating
unsaturated zone flow.

The HPM simulates water flow in the unsaturated zone
during a particular MODFLOW time step by considering thewater table depth calculated by MODFLOW during the previous time step as the bottom boundary condition. When the HPM
profile is associated with more than one MODFLOW grid cell, theaverage value of the water table position in these cells is used as thebottom boundary condition. The bottom flux calculated by the
HPM during the current time step is passed to MODFLOW as arecharge flux to simulate the groundwater flow and to calculate
the groundwater level change during the current time step. The
bottom boundary condition thus changes stepwise in the HPM.
See Seo et al. (2007) for more details.

Limitations of the HYDRUS Package
The current coupling algorithm implicitly assumes that thepressure head at the bottom of the HPM profile remains constantthroughout the MODFLOW time step. This is less realistic when
there is a large change in the water table level during this time step
due to recharge from the unsaturated zone or due to some otherlateral inflow into the saturated soil zone. In reality, the groundwater table varies continuously with time in response to vadose zone
flow. A sudden change in the position of the groundwater tablein the HPM profile after every MODFLOW time step leads to achange in the bottom boundary condition (a time-variable pressure
head boundary condition). Since the bottom boundary condition
does not at this moment correspond with the flow profile at theprevious time step in the HPM profile, this results in a sudden
inflow into the soil profile when the water table level is increased
or a sudden outflow when the water table level is decreased. This
results in the inaccurate estimation of the cumulative bottom flux
throughout the duration of the simulation (Twarakavi et al., 2009;
Kuznetsov et al., 2012).

Objectives
The main objective of this study was to eliminate sudden
inflow or outflow fluxes at the bottom of the HPM profile
when the groundwater table depth changes. This was achievedby updating the coupling algorithm between MODFLOW andthe HPM. The second objective was to evaluate the applicability
p. 2 of 8

of the updating algorithm by comparing its predictions of watertable elevations obtained using MODFLOW with the HYDRUS
package with the results obtained using the HYDRUS-1D and
HYDRUS (2D/3D) models.

Updates of Bottom Pressure Head in the HYDRUS
Package Profile to Eliminate Spurious Fluxes
Sudden inflow or outflow fluxes at the bottom of the HPM
profile when the groundwater table depth changes can be eliminated if we update pressure heads at the bottom of the HPM profileat the beginning of the MODFLOW time step to an equivalentpressure head profile that would exist if there was a continuouschange in the water table elevation. A pressure head profile updating algorithm was developed to achieve this objective (Fig. 1).
Figure 1 shows how the pressure head profile in the HPM profile is updated after every MODFLOW time step. Let T1 and T2
be two consecutive MODFLOW time steps. Initially, the averagewater table level (WT1) from the MODFLOW grids is given as thebottom boundary condition (a constant bottom pressure head, h1) tothe corresponding HPM profile. The HYDRUS package simulates
unsaturated flow using the 1D Richards equation until T1. Figure
1a shows the pressure head profile in the HPM profile at time T1.
The bottom flux q1 calculated using this pressure head profile is given as the recharge flux to MODFLOW for simulating
groundwater flow and the new water table level (WT2) at time T1.
During the next time step, T2 , the bottom boundary condition in
the HPM profile will be h2 , i.e., the pressure head corresponding
to WT2 . However, imposing this pressure head boundary condition at the bottom of the HPM profile while keeping the samepressure head profile (from T1) would produce a sudden inflow
or outflow from the soil profile depending on whether the watertable increased or decreased, respectively, during the time step T1.
To eliminate this sudden inflow or outflow flux, a new algorithm
is proposed that adjusts pressure heads at the bottom of the HPM
profile while taking into account the new groundwater table position WT2 and the final recharge flux q1 from the previous stressperiod T1.
This can be done in a straightforward way in the saturated zone (due to the linear nature of Darcy’s law), whereas aroot-finding routine is used to adjust pressure heads in the unsaturated zone (due to the nonlinear nature of the Darcy–Buckinghamlaw). According to the discretization scheme used in the HYDRUS
package, the water flux between two adjacent finite element nodes
i and i + 1 in the HPM profile is expressed asq =-öK ( hi ) + K ( hi +1 )æç hi +1 - hi
+1÷÷÷
çç
÷ø
çè z i +1 - z i2[3]

where q is the bottom flux [L T−1], K is the unsaturated hydraulicconductivity [L T−1], hi and hi+1 are the pressure heads in the two
nodes at the bottom of the soil profile [L], and zi and zi+1 are thevertical coordinates of these two nodes [L]. Equation [3] has to besolved for hi+1, in which the value hi is known and q is equal tothe bottom flux. Since K depends on hi+1 in a strongly nonlinear
manner, the equation is solved iteratively using the “false position”
method (e.g., Press et al., 2007). The steady-state pressure headprofile corresponding to a flux equal to q1 and the bottom pressure
head equal to h2 is shown in Fig. 1b.
The nodal fluxes in the HPM profile at time T1 (Fig. 1a) arethen compared with the bottom recharge flux q1, and a node isfound where the relative difference is >0.1% (plus a small roundoff value of 10−12). The nodal pressure heads below this node arethen set equal to the steady-state pressure head profile for q1 (Fig.
1b), and the nodal pressure heads above this node are kept equal
to those at T1. This procedure is performed at every MODFLOW
time step in each of the HPM profiles to eliminate sudden fluxes
due to abrupt changes in the water table elevation. Figure 1c (redline) shows the pressure head profile in the HPM profile updated
at the end of T1 before moving to T2 . Note that the transition
between two parts of the pressure head profile is not smooth, but
it improves the consistency of the model as shown in the illustrating examples.
MODFLOW with the HPM does not represent a fully integrated modeling of flow processes in the unsaturated and saturated
soil zones. Because of the stepwise change in the bottom boundary condition, calculations of a combined mass balance for both
unsaturated and saturated soils are not possible. Therefore, themass balance is considered separately for MODFLOW and theFig. 1. Illustration of the pressure head profile updating algorithm incorporated in the
HYDRUS package (HPM) to eliminate sudden variations in the bottom flux in the HPM
profiles: (a) pressure head profile in the HPM
profile at the end of time step T1; (b) steadystate pressure head profile for water table level
WT2; and (c) pressure head profile in the
HPM profile T1 updated before moving to T2(red line).

VZJ | Advancing Critical Zone Sciencep. 3 of 8

HYDRUS package simulations. The mass balance in the original
HYDRUS-1D source code is modified to account for the amountof water added or removed from the soil profile due to updatingof the pressure head profile. The mass balance in the groundwater
domain is calculated by including the coupling flux as the recharge
to the groundwater domain (Szymkiewicz et al., 2018).

Illustrating Examples
The performance of the pressure head profile updating algorithm is illustrated using two case studies. The first case study
illustrates the elimination of the sudden change in the bottom fluxin the HPM profile. This case study also examines the effects ofdifferent MODFLOW time steps on the cumulative bottom fluxin the HPM profile and the water table elevation obtained using
MODFLOW with the HYDRUS package. The performance of
the algorithm is checked by comparing the results obtained using
either MODFLOW with the HPM or HYDRUS-1D.
The second case study highlights the effectiveness of
MODFLOW with the HYDRUS package (which is a 3D saturated flow model coupled with a 1D unsaturated flow model)
by comparing its results with HYDRUS (2D/3D) (a fully 3Dmodel). HYDRUS (2D/3D) (Šimůnek et al., 2016) simulates
water flow in both unsaturated and saturated domains by solving
a three-dimensional version of the standard Richards equation.
This comparison is performed to show that MODFLOW with
the HYDRUS package can be used for modeling larger domains
without compromising the accuracy of the results.

Case Study 1: Evaluating the Effectiveness of MODFLOW
with the HYDRUS Package Using HYDRUS-1D
A domain with a length and a width both equal to 1 m and a
depth of 10 m was considered in this case study. The MODFLOW
domain was divided into four equal grid cells with no-flow boundaries so that the position of the water table could change only
due to recharge. An atmospheric boundary condition with daily
varying precipitation and potential evapotranspiration was considered at the soil surface (Fig. 2). This meteorological dataset was
obtained from the weather station in Gdańsk, Poland, for 2011.
The HPM profiles were considered with a depth of 5 m. The
soil profile was considered to have a constant pressure head equal
to −0.283 m down to a depth of 3.5 m from the soil surface and ahydrostatic initial pressure head distribution below this depth. The

Fig. 2. Daily precipitation and potential evapotranspiration used inCase Study 1.

groundwater table was considered to be initially at a depth of 3.95
m. The HPM profile was divided into 200 finite elements withrelatively smaller elements near the top to ensure the convergenceof the numerical solution.
The HPM profile was considered to have two soil layers: (i) asandy soil down to a depth of 2.5 m and (ii) a loamy sand soil below
this depth. The van Genuchten–Mualem analytical model (vanGenuchten, 1980) was used to describe the soil hydraulic propertiesof the sandy and loamy sand soils. The following van Genuchten–
Mualem parameters were considered for the sandy and loamy
sand soils: the residual water content, q r = 0.045 and 0.057; thesaturated water content, q s = 0.43 and 0.41; the saturated hydraulic conductivity, Ks = 7.13 and 3.50 m d−1; the pore-connectivity
parameter, l = 0.5 and 0.5; and the shape parameters, a = 14.5 and12.4 m−1, and n = 2.68 and 2.28. The simulations were performedfor 1 yr. The maximum allowed time step in the HPM simulations
was 0.1 d.
The effect of different MODFLOW time steps was analyzed
by dividing the duration of the simulation of 365 d into a different number of MODFLOW time steps (20, 40, and 365). In
the MODFLOW domain, we assumed a constant and isotropic
hydraulic conductivity equal to 3.50 m d−1, the specific yield equal
to 0.255, and the specific storage coefficient of 0.0015 m−1.
Figure 3 shows the cumulative fluxes at the bottom of the
HPM profile obtained using MODFLOW with the HYDRUS
package for different numbers of MODFLOW time steps when
the pressure head distribution either was or was not modified. For

Fig. 3. Cumulative fluxes at the bottom of the HYDRUS
package profile obtained using MODFLOW with the
HYDRUS package for a different number of MODFLOW time steps (TS) with (modified) and without(unmodified) pressure head modifications. Negativevalues indicate downward flow.

VZJ | Advancing Critical Zone Sciencep. 4 of 8all time steps, it can be observed that the cumulative bottom flux
(in absolute values) increased with time. Cumulative fluxes for thesimulation with the pressure head modifications were higher thanwhen there were no pressure head modifications. This is because of
the removal of upward inflow fluxes after every time step with anincrease in the water table, which were generated when the pressure heads were not modified. The difference between the two setsof lines (for simulations with and without pressure head modifications) corresponds to the amount of water added to the HPM
profiles as a result of the pressure head modifications.
Figure 4 shows the water table elevation obtained using
MODFLOW with the HYDRUS package with and withoutpressure head modifications for different MODFLOW time
steps. Water table elevations were found to be at a higher level
when the bottom pressure head distributions were modified for
all three considered MODFLOW time steps than when they were
not modified. This is because of the removal of sudden upward
fluxes after every MODFLOW time step, which resulted in a lower
water table rise.
When the number of MODFLOW time steps is 365, the watertable level is calculated by MODFLOW and provided as the bottom
boundary condition to the HPM each day. When the number of
MODFLOW time steps is 40 and 20, the water table elevation (thebottom boundary condition) is updated approximately only after
every 9 and 18 d, respectively. As a result, the initiation of the watertable rise is delayed in the latter two cases by 9 and 18 d, respectively,
and the overall increase in the water table level was lower than whena smaller MODFLOW time step was used.

A constant value of the specific yield was considered in
MODFLOW calculations with the HPM. The specific yield canbe considered to be constant only when the aquifer response islinear, i.e., when the volume of water released is linearly proportional to the water table fluctuation. However, this is not always
true because this parameter is highly dependent on the unsaturated
zone properties. The specific yield determines the change in the
groundwater table position in response to the simulated recharge
flux and thus it determines the updated boundary condition at
the bottom of each HPM profile. There have been many studies
evaluating the value of the specific yield as a function of various
factors that affect the specific yield, such as the transient natureof the water release from the unsaturated soil zone, soil hydraulicproperties, depth of the groundwater table, or hysteresis (Nachabe,
2002; Said et al., 2005). The use of the proper value of the specificyield, taking into account unsaturated flow, is the subject of ongoing research (Stoppelenburg et al., 2005; Xu et al., 2012).
The water table elevation obtained in this case study was further compared with the results obtained using the HYDRUS-1Dmodel by assuming a zero flux as the bottom boundary condition. The same domain and initial conditions were used in
HYDRUS-1D. Figure 5 shows the comparison of the water table
elevation obtained using HYDRUS-1D and MODFLOW with
the HYDRUS package for different specific yield values. The
temporal evolution of the water table elevation obtained with a
specific yield of 0.255 was found to be the same as that obtainedusing HYDRUS-1D, whereas a lower value of the specific yield wasfound to give a larger increase in the water table elevation and vice

Fig. 4. Water table elevations obtained using MODFLOW with the HYDRUS package for a different
number of MODFLOW time steps (TS) with (modified) and without (unmodified) pressure head
modifications.

Fig. 5. The water table elevation obtained using
MODFLOW with the HYDRUS package, with 365
MODFLOW time steps (TS) and different values ofthe specific yield (Sy), and the HYDRUS-1D model.

VZJ | Advancing Critical Zone Science

p. 5 of 8

Fig. 7. Monthly precipitation considered in Case Study 2.

Fig. 6. Hypothetical domain and the initial water table level considered in Case Study 2.

versa. The temporal evolution of the water table elevation is highly
dependent on the value of the specific yield.

Case Study 2: Evaluating the effectiveness of MODFLOW
with the HYDRUS Package Using HYDRUS (2D/3D)
A hypothetical domain used in a similar study by Morway
et al. (2013) was considered in this case study. The domain had
dimensions of 8000 by 4000 by 15 m, with a slope of 0.001 m
along the 4000-m side. The initial condition considered a groundwater table level varying between 7 m above the bottom on one sideof the domain and 0.9 m on the opposite side of the domain (Fig.
6). A constant pressure head was considered along the 8000-m
sides throughout the simulation. A no-flow boundary condition
was assigned on the other two sides as well as at the bottom of
the domain. An atmospheric boundary condition was applied at
the soil surface with a monthly varying rainfall as shown in Fig. 7.

The same atmospheric boundary condition was repeated during
a 5-yr simulation.
In MODFLOW with the HYDRUS package, the entire
domain was divided into 80 rows and 40 columns (each grid had
a dimension of 100 by 100 m). The MODFLOW grid cells weredivided into 40 zones, and each of these zones was assigned an HPM
profile (numbered from 1 to 40 in Fig. 6). The total simulation timeof 1825 d was divided into 60 stress periods, with time steps equal to
the number of days in the month. A specific yield of 0.28 was used.
Each HPM profile was divided into 60 finite elements of equalsize. The maximum allowable time step in the HYDRUS package is 1 d. The soil was considered to be homogenous throughout
the domain. The van Genuchten–Mualem analytical model (vanGenuchten, 1980) was used to describe the soil hydraulic properties
with the following parameters: q r = 0.1, q s = 0.45, Ks = 50 m d−1,
l = 0.5, a = 1.65 m−1, and n = 2. The same domain was modeledusing HYDRUS (2D/3D) with the same soil hydraulic properties
and with the same maximum allowable time step. Figures 8a and8b show the initial and final pressure head distributions, respectively, in the domain obtained using HYDRUS (2D/3D).

Fig. 8. (a) The initial pressure head distribution and (b) the final pressure head distribution obtained using HYDRUS (2D/3D).

VZJ | Advancing Critical Zone Science

p. 6 of 8

suing this work. Adam Szymkiewicz has been supported by the National ScienceCentre, Poland, through the project 2015/17/B/ST10/03233 “Groundwater recharge on outwash plain.”

References
Abbott, M.B., J.C. Bathurst, J.A. Cunge, P.E. O’Connell, and J. Rasmussen.
1986. An introduction to the European Hydrological System—Systeme Hydrologique Europeen, “SHE”: 1. History and philosophy of aphysically-based, distributed modelling system. J. Hydrol. 87:45–59.
doi:10.1016/0022-1694(86)90114-9
Brooks, R., and A. Corey. 1964. Hydraulic properties of porous media.
Hydrol. Pap. 3. Colorado State Univ., Fort Collins.
Diersch, H.J., and O. Kolditz. 1998. Coupled groundwater flow and transport: 2. Thermohaline and 3D convection systems. Adv. Water Resour.
21:401–425. doi:10.1016/S0309-1708(97)00003-1
Durner, W. 1994. Hydraulic conductivity estimation for soils withheterogeneous pore structure. Water Resour. Res. 30:211–223.
doi:10.1029/93WR02676
Facchi, A., B. Ortuani, D. Maggi, and C. Gandolfi. 2004. Coupled
SVAT–groundwater model for water resources simulation in irrigated alluvial plains. Environ. Model. Softw. 19:1053–1063.
doi:10.1016/j.envsoft.2003.11.008
Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald. 2000. MODFLOW-2000, the U.S. Geological Survey modular ground-water model:
User guide to modularization concepts and the ground-water flow
process. Open-File Rep. 00-92. USGS, Reston, VA.

Fig. 9. Water table elevations (a) at the center of the domain and (b)at the end of the simulation along the width of the domain obtainedusing MODFLOW with the HYDRUS package and HYDRUS(2D/3D).

Harbaugh, A.W. 2005. MODFLOW-2005, The U.S. Geological Survey Modular Ground-Water Model: The ground-water flow process. Techniquesand Methods 6-A16. USGS, Reston, VA.

The temporal evolution of the water table elevation at thecenter of the domain throughout the simulation period and thefinal water table level along the width of the domain obtainedusing MODFLOW with the HYDRUS package and HYDRUS(2D/3D) are shown in Fig. 9a and 9b, respectively. The results
obtained by the two models were found to be similar, which indicates that MODFLOW with the HYDRUS package is able tosimulate the water table elevation accurately.

HydroGeoLogic. 1996. MODFLOW-SURFACT software (Version 2.2) overview: Installation, registration, and running procedures. HydroGeoLogic Inc., Herndon, VA.

Concluding Remarks

Maxwell, R.M., S.J. Kollet, S.G. Smith, C.S. Woodward, R.D. Falgout, I.M. Ferguson, et al. 2016. ParFlow user’s manual. Rep. GWMI 2016-01. Integr.
Groundw. Model. Ctr., Colorado School of Mines, Golden.

The coupling algorithm between MODFLOW and the
HYDRUS package was modified to eliminate the spurious fluxes
at the HYDRUS–MODFLOW interface. The performance of the
proposed method was evaluated by comparing the results obtainedusing MODFLOW with the HYDRUS package with the results
obtained using either the HYDRUS-1D or HYDRUS (2D/3D)
models. The current version of the HYDRUS package neglects
the effects of unsaturated zone soil properties on the specific yield.
Future research should incorporate the influence of vadose zone
properties on the specific yield in the HYDRUS package, which
would lead to improved modeling of water flow in the saturated–unsaturated soil zone and more accurate description of the watertable dynamics.

Acknowledgments

Sahila Beegum would like to acknowledge support from the United States India
Education Foundation (USIEF) through the Fulbright–Nehru Doctoral Research
Fellowship for the research stay at the University of California–Riverside for pur-

VZJ | Advancing Critical Zone ScienceHavard, P.L., S.O. Prasher, R.B. Bonnell, and A. Madani. 1995. Linkflow, awater flow computer model for water table management: I. Modeldevelopment. Trans. ASAE 38:481–488. doi:10.13031/2013.27856

Kosugi, K.I. 1996. Lognormal distribution model for unsaturated soil
hydraulic properties. Water Resour. Res. 32:2697–2703.
Kuznetsov, M., A. Yakirevich, Y.A. Pachepsky, S. Sorek, and N. Weisbrod.
2012. Quasi 3D modeling of water flow in vadose zone and groundwater. J. Hydrol. 450–451:140–149. doi:10.1016/j.jhydrol.2012.05.025
Lin, L., J.-Z. Yang, B. Zhang, and Y. Zhu. 2010. A simplified numerical modelof 3-D groundwater and solute transport at large scale area. J. Hydrodyn., Ser. B 22:319–328. doi:10.1016/S1001-6058(09)60061-5

Mayer, K., R. Amos, S. Molins, and F. Gerard. 2012. Reactive transport modeling in variably saturated media with MIN3P: Basic model formulationand model enhancements. In: F. Zhang et al., editors, Groundwater reactive transport models. Bentham Sci. Publ., Sharjah, UAE. p. 186–211.
Morway, E.D., R.G. Niswonger, C.D. Langevin, R.T. Bailey, and R.W. Healy.
2013. Modeling variably saturated subsurface solute transportwith MODFLOW-UZF and MT3DMS. Groundwater 51:237–251.
doi:10.1111/j.1745-6584.2012.00971.x
Nachabe, M.H. 2002. Analytical expressions for transient specific yieldand shallow water table drainage. Water Resour. Res. 38(10):1193.
doi:10.1029/2001WR001071
Niswonger, R.G., D.E. Prudic, and S.R. Regan. 2006. Documentation of theUnsaturated-Zone Flow (UZF1) package for modeling unsaturated
flow between the land surface and the water table with MODFLOW-2005. Techniques and Methods 6-A19. USGS, Reston, VA.
Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. 2007. Numerical recipes: The art of scientific computing. 3rd ed. Cambridge Univ.
Press, New York.
Pruess, K., C.M. Oldenburg, and G.J. Moridis. 1999. TOUGH2 user’s guideversion 2. LBNL-43134). Lawrence Berkeley Natl. Lab., Berkeley, CA.
doi:10.2172/751729

p. 7 of 8

Seo, H., J. Šimůnek, and E. Poeter. 2007. Documentation of the HYDRUS
package for MODFLOW-2000, the US Geological Survey modular
ground-water model. IGWMI 2007-01. Integr. GroundWater Model. Ctr.,
Colorado School of Mines, Golden, CO.
Said, A., M. Nachabe, M. Ross, and J. Vomacka. 2005. Methodology
for estimating specific yield in shallow water environment using
continuous soil moisture data. J. Irrig. Drain. Eng. 131:533–538.
doi:10.1061/(ASCE)0733-9437(2005)131:6(533)Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2016. Recent developments and applications of the HYDRUS computer software packages.
Vadose Zone J. 15(7). doi:10.2136/vzj2016.04.0033
Stoppelenburg, F.J., K. Kovar, M.J.H. Pastoors, and A. Tiktak. 2005. Modelling the interactions between transient saturated and unsaturated
groundwater flow: Offline coupling of LFM and SWAP. RIVM Rep.
500026001/2005. Rijksinstituut voor Volksgezondheid en Milieu,
Bilthoven, the Netherlands.
Szymkiewicz, A., A. Gumuła-Kawęcka, J. Šimůnek, B. Leterme, S. Beegum,
B. Jaworska-Szulc, et al. 2018. Simulation of the freshwater lens recharge using the HYDRUS and SWI2 packages for MODFLOW. J, Hydrol,
Hydromech, 66:246–256. doi:10.2478/johh-2018-0005
Therrien, R., R.G. McLaren, E.A. Sudicky, and S.M. Panday. 2010. HydroGeoSphere: A three-dimensional numerical model describing fully-integrated subsurface and surface flow and solute transport. Groundwater
Simulations Group, Univ. of Waterloo, Waterloo, ON, Canada.
Thoms, R.B., R.L. Johnson, and R.W. Healy. 2006. User’s guide to the Vari-VZJ | Advancing Critical Zone Science

ably Saturated Flow (VSF) process to MODFLOW. Techniques andMethods 6-A18. USGS, Reston, VA.
Twarakavi, N.K.C., J. Šimůnek, and S. Seo. 2008. Evaluating Interactions between groundwater and vadose zone using the HYDRUSbased flow package for MODFLOW. Vadose Zone J. 7:757–768.
doi:10.2136/vzj2007.0082
Twarakavi, N.K.C., J. Šimůnek, and S. Seo. 2009. Reply to “Commenton ‘Evaluating interactions between groundwater and vadosezone using the HYDRUS-based flow package for MODFLOW’ ” by
N.K.C. Twarakavi, J. Šimůnek, and S. Seo. Vadose Zone J. 8:820–821.
doi:10.2136/vzj2008.0004L
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44:892–898. doi:10.2136/sssaj1980.03615995004400050002x
Vogel, T., and M. Císlerova. 1988. On the reliability of unsaturated hydraulic conductivity calculated from the moisture retention curve. Transp.
Porous Media 3:1–15. doi:10.1007/BF00222683
Xu, X., G. Huang, H. Zhan, Z. Qu, and Q. Huang. 2012. Integration of
SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas. J. Hydrol. 412–413:170–181.
doi:10.1016/j.jhydrol.2011.07.002
Zhu, Y., Y. Zha, J. Tong, and J. Yang. 2011. Method of coupling 1-D unsaturated flow with 3-D saturated flow on large scale. Water Sci. Eng.
4(4):357–373. doi:10.3882/j.issn.1674-2370.2011.04.001

p. 8 of 8