Implementation of Solute Transport in the
Vadose Zone into the ‘‘HYDRUS Package for
MODFLOW’’by Sahila Beegum1,2 , Jiří Šimůnek1 , Adam Szymkiewicz3 , K. P. Sudheer4,5 , and Indumathi M. Nambi4

Abstract
The ‘‘HYDRUS package for MODFLOW’’ is an existing MODFLOW package that allows MODFLOW to simultaneously evaluate
transient water flow in both unsaturated and saturated zones. The package is based on incorporating parts of the HYDRUS-1D model(to simulate unsaturated water flow in the vadose zone) into MODFLOW (to simulate saturated groundwater flow). The coupled modelis effective in addressing spatially variable saturated-unsaturated hydrological processes at the regional scale. However, one of the
major limitations of this coupled model is that it does not have the capability to simulate solute transport along with water flow andtherefore, the model cannot be employed for evaluating groundwater contamination. In this work, a modified unsaturated flow andtransport package (modified HYDRUS package for MODFLOW and MT3DMS) has been developed and linked to the three-dimensional(3D) groundwater flow model MODFLOW and the 3D groundwater solute transport model MT3DMS. The new package can simulate,
in addition to water flow in the vadose zone, also solute transport involving many biogeochemical processes and reactions, including
first-order degradation, volatilization, linear or nonlinear sorption, one-site kinetic sorption, two-site sorption, and two-kinetic sites
sorption. Due to complex interactions at the groundwater table, certain modifications of the pressure head (compared to the original
coupling) and solute concentration profiles were incorporated into the modified HYDRUS package. The performance of the newly
developed model is evaluated using HYDRUS (2D/3D), and the results indicate that the new model is effective in simulating the
movement of water and contaminants in the saturated-unsaturated flow domains.

Introduction
The Importance and History of Integrated Modeling
The modeling of water flow and solute transport
through saturated and unsaturated soil zones requires
proper consideration of the interactions between thesetwo zones. Accurate assessment of the movement ofwater and solute through these zones is important forplanning and management of water resources (e.g.,
1 Departmentof Environmental Sciences, University ofCalifornia, Riverside, CA.
2 Corresponding author: Department of Civil Engineering,
Indian Institute of Technology Madras, Chennai 600036, India;
sahilabeegum666@gmail.com
3
Gdańsk University of Technology, Faculty of Civil andEnvironmental Engineering, Gdańsk, Poland.
4
Department of Civil Engineering, Indian Institute ofTechnology Madras, Chennai 600036, India.
5
Department of Agricultural and Biological Engineering,Purdue University, West Lafayette, IN 47906.
Article impact statement: Solute transport is integrated intothe ‘‘HYDRUS Package for MODFLOW’’ and the modified model
performance is evaluated using HYDRUS (2D/3D).
Received December 2017, accepted July 2018.
© 2018, National Ground Water Association.
doi: 10.1111/gwat.12815

392

groundwater and surface water management, irrigation
planning, groundwater pollution remediation, and artificial
recharge to groundwater) (Zhu et al. 2011, 2016). There
are various models, from simple analytical to complex
numerical, that can independently simulate a multitudeof processes in the unsaturated and saturated zones. To
consider various processes and interactions in and between
these zones, integrated models of varying degrees ofcomplexity and dimensionality have been developed in
the past by coupling models developed independently
for unsaturated and saturated soil zones (Twarakavi et al.
2008; Ragab and Bromley 2010; Kalbacher et al. 2012).
The major difficulty in integrating the existing models
is due to the differences in various factors that affectflow and transport in the two zones. This becomes morechallenging when one has to consider the dynamic natureof the water table that is dependent on external factors
such as rainfall, evapotranspiration, surface runoff, andon the soil hydraulic properties. The temporal and spatial
scales of various processes in these two zones also vary.
Flow in the saturated zone occurs at a slower rate andrequires more time to produce significant changes in water
fluxes and total heads as compared to the unsaturated
zone. The spatial scale of groundwater flow applications

Vol. 57, No. 3–Groundwater–May-June 2019 (pages 392–408)

NGWA.org

is also significantly different (larger) than the spatial
scale of the majority of vadose zone flow applications.
A three-dimensional (3D) nature of flow and transport in
integrated models results in increased computational timerequirements along with an increased demand for model
input data and physical parameters for both unsaturated
and saturated zones.
Considering that flow and transport in the unsaturated zone are predominantly vertical, many integrated
models have the unsaturated soil zone modeled as one
dimensional along with a 3D modeling of processes inthe saturated zone (e.g., Niswonger et al. 2006; Twarakavi
et al. 2008; Zhu et al. 2011). This simplification of thedimensionality of the problem has resulted in a significant reduction in the computational complexity. Someof the existing integrated models are also limited to waterflow and cannot consider solute transport (e.g., Niswonger
et al. 2006; Twarakavi et al. 2008; Ragab and Bromley2010).
Various integrated models have been developed at
different levels of temporal, spatial, and dimensional
complexities and with different considered processes (e.g.,
Refsgaard and Storm 1995; Niswonger et al. 2006; Thoms
et al. 2006; Panday and Huyakorn 2008; Twarakavi
et al. 2008; Ragab and Bromley 2010; Brunner andSimmons 2012; Kalbacher et al. 2012; Morway et al.
2013). For example, MIKE SHE (Refsgaard and Storm
1995) is a fully distributed, physically-based hydrologicalmodel that considers evapotranspiration, unsaturated flow,
overland flow, groundwater flow, channel flow, and theirinteractions, as well as solute transport. MIKE SHE
requires extensive input data and physical parameters,
which makes it rather difficult to set up the model.
There have been various approaches to integrating an
unsaturated zone model into MODFLOW (Harbaugh et al.
2000), a widely used 3D groundwater flow model, whichin its basic implementation does not consider flow inthe vadose zone. For example, MODFLOW-VSF (Thoms
et al. 2006) models 3D variably saturated flow processes
by solving the 3D Richards equation. MODFLOW-UZF(Niswonger et al. 2006) uses a one-dimensional (1D)
kinematic wave equation to describe flow processes in theunsaturated zone. Morway et al. (2013) used MOFLOWUZF to first simulate variably saturated subsurface water
flow, and then, in combination with MT3DMS, also solutetransport. MODFLOW-SURFACT (Panday and Huyakorn
2008) can simulate unsaturated flow and recharge, fracture
flow, and contaminant transport. The HYDRUS-based
flow package for MODFLOW (HPM) was developed by
Seo et al. (2007) and Twarakavi et al. (2008) to simulate
flow in both unsaturated and saturated zones. The package,
which is based on the HYDRUS-1D model (Šimůnek
et al. 2016) simulating unsaturated water flow in thevadose zone, was linked to MODFLOW (Harbaugh et al.
2000) simulating saturated groundwater flow.
Additional examples of integrated models are Integrated Hydrological Modelling System (IHMS) (Ragab
and Bromley 2010), HydroGeosphere (Brunner and Simmons 2012), and IWAS (Kalbacher et al. 2012). The

IHMS developed by Ragab and Bromley (2010) estimates
recharge from the unsaturated soil zone as a component of
the water balance and provides it as input to MODFLOW.
HydroGeoSpere (Brunner and Simmons 2012) solves 3D
saturated-unsaturated flow and transport equations for thevadose and saturated zones using the control volume finite
element method. The IWAS toolbox (Kalbacher et al.
2012) integrates OpenGeoSys (Kolditz et al. 2012), whichsolves the diffusive-wave overland flow equation and thepressure head form of the 3D Richards equation, withthe mesoscale Hydrologic Model (Samaniego et al. 2011),
which evaluates land surface hydrological processes, andthe Root (Schneider et al. 2010) model, which calculates
root water uptake.

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

The HYDRUS Package for MODFLOW: The Current Status
The HPM was developed by Seo et al. (2007) andTwarakavi et al. (2008) to simultaneously evaluate transient water flow in both unsaturated and saturated zones.
The package, which is based on the HYDRUS-1D model(Šimůnek et al. 2016) that simulates unsaturated waterflow in the vadose zone, was linked to MODFLOW(Harbaugh et al. 2000) that simulates saturated groundwater flow. Even though HYDRUS-1D has the capabilityof simulating solute transport in the vadose zone, only
the water flow part of HYDRUS-1D has been coupledwith MODFLOW in the HPM. The HPM is effective insimulating spatially variable saturated-unsaturated hydrological processes at the regional scale (Twarakavi et al.
2008). The HYDRUS-1D’s capabilities of considering
complex layering of the vadose zone and varying fluxes(both spatial and temporal) help to provide realistic input
fluxes to the MODFLOW to simulate the saturated zone
water flow. Applications of the HPM for regional scale
groundwater flow and comparisons with the UZF1 package were described in Leterme et al. (2013, 2015). The
most recent version of the HPM, which is available to theusers, is compatible with MODFLOW-2005 (Harbaugh2005). Currently, the package is being further developedin cooperation between research groups from the University of California Riverside, the Belgian Nuclear ResearchCentre (SCK•CEN), and the Gdańsk University ofTechnology.
Limitations of the HPM
The current version of the HPM can simulate waterflow in unsaturated and saturated soil zones but is unable
to simulate solute transport. Solute transport is highly
dependent on water table fluctuations due to temporaland spatial variations in groundwater recharge. This is
an important concern when the coupled model is used
for analyzing groundwater contamination due to transport
through the unsaturated zone.
Additionally, the current numerical coupling of
HYDRUS-1D and MODFLOW produces unrealistic sudden inflow and outflow fluxes at the interface when the
groundwater table depth is updated in HYDRUS-1D as aresult of its change in MODFLOW (Twarakavi et al. 2009;
Kuznetsov et al. 2012). These numerical instabilities need393

to be resolved before the solute transport component canbe incorporated into the HPM.
Objectives
The primary objective of the study is to incorporate
the solute transport into the HPM. However, to overcome
the limitation of the HPM, which generated instabilities
in water fluxes at the interface between the vadose andgroundwater zones, was considered a priority as well,as these instabilities affect solute transport computations.
The incorporation of the solute transport is achieved in
two steps: (1) integrating the solute transport componentof HYDRUS-1D into the HPM (hereafter referred to as“the modified HPM”) and (2) coupling the modified HPM
with MODFLOW and linking it with MT3DMS. The performance of the modified HPM and its coupling with
MODFLOW and MT3DMS is evaluated by comparing
its predictions of pressure heads, solute concentrations,and water and solute fluxes at different locations in thetransport domain with the results obtained by HYDRUS(2D/3D). HYDRUS (2D/3D) (Šimůnek et al. 2016) simulates water flow and solute transport in both unsaturated
and saturated domains by solving 3D versions of thestandard Richards and convection–dispersion equations.
Finally, the dependence of time requirements of the coupled model and its particular parts (e.g., MODFLOW, themodified HPM) to spatial and temporal discretization isanalyzed and discussed.
Since this manuscript frequently refers to various
versions of the HYDRUS models and its modifications,as well as to the MODFLOW and MT3DMS models, thenomenclature presented in Table 1 is used to refer eachmodel in this paper.

The HYDRUS Package for MODFLOW:
Modeling Approach
Water Flow
Governing Equations
Groundwater flow is modeled in MODFLOW bysolving the mass conservation equation using the finite
difference approximation. While MODFLOW is capableof simulating 3D flow in confined and unconfined
aquifers, here we limit our analysis (without the lossof generality) to a single unconfined aquifer, represented
by one layer of grid blocks. The two-dimensional (2D)
movement of groundwater of constant density in an
unconfined aquifer is described by the partial differential
equation (derived by applying the Dupuit assumption) as
follows:




∂H
∂H
∂
∂
∂H
H
+ Ky
H
= Sy
(1)
Kx
∂x
∂x
∂y
∂y
∂t
where K x and K y are the hydraulic conductivities[LT−1 ] in the x and y directions, respectively, H is thepiezometric head [L], S y is the specific yield [−] of theporous material, and t is time [T].
394

S. Beegum et al. Groundwater 57, no. 3: 392–408

HYDRUS-1D solves water flow in the unsaturated
zone using the modified 1D Richards equation:


∂h
∂
∂θ
K (h)
(2)
=
− K (h) − S∂t∂z∂z
where θ is the volumetric water content (dimensionless),
h is the soil water pressure head [L], t is time [T], z is thevertical coordinate [L], S is the sink term [T−1 ], and K (h)is the unsaturated hydraulic conductivity [LT−1 ]. The
unsaturated hydraulic conductivity, K (h), and the watercontent, θ (h), depend on the soil water pressure head. This
makes the Richards equation a highly nonlinear equationthat needs to be solved numerically. HYDRUS permits
the use of five different analytical models to describe thesoil hydraulic properties (i.e., Brooks and Corey 1964;
van Genuchten 1980; Vogel and Cislerova 1988; Durner
1994; Kosugi 1996).
Numerical Implementation
In MODFLOW, the groundwater modeling domain
is discretized into regular grids. These grids can be
combined into multiple zones based on similarities in soilhydrology, topographical characteristics, and the depth to
the groundwater. In the HPM, each of these zones is
assigned one soil 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(hereafter referred to as “HPM profile”). Each of the HPM
profiles is divided into finite elements. The soil profile canhave an arbitrary number of soil layers with different soil
hydraulic properties. The user can specify the number andposition of nodes used to discretize each HPM profile. In
order to ensure the numerical stability and efficiency, thesize of finite elements should be small at locations where
sharp pressure head gradients are expected, especiallyclose to the soil surface. The details of the numerical
scheme used to solve the Richards equation in HYDRUS1D (on which the HPM is based) can be found in Šimůnek
et al. (2013).
In MODFLOW, the total simulation period is divided
into stress periods during which the external stresses areconstant. Each stress period is divided into time steps. For
each MODFLOW time step, a simulation of flow in theunsaturated zone is performed by the HPM for each HPM
profile. The HPM uses its own time stepping algorithm,
based on user-defined criteria. In general, HPM time steps
are smaller than MODFLOW time steps because flow
conditions in the vadose zone vary more rapidly than in
the saturated zone and because of the requirements of
the numerical solution of the nonlinear Richards equation.
After every MODFLOW time step, the flux at the bottomof the HPM profile is given 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 conditionin the HPM for the next MODFLOW time step. If a
specific HPM profile is associated with more than one
MODFLOW grid cell, the average value of the water
NGWA.org

Table 1
Nomenclature of the Models (HYDRUS and Its Modifications) Considered in the Study
No.

Name

Description

1
2

MODFLOW
MT3DMS3HYDRUS-1D

4

HYDRUS (2D/3D)

5

HPM

6

Modified HPM

A modular 3D groundwater flow model (Harbaugh et al. 2000).
A modular 3D multispecies solute transport model for groundwater (Zheng andWang 1999)
The software package for simulating the 1D movement of water, heat, and multiple
solutes in variably-saturated media (Šimůnek and Bradford 2008; Šimůnek et al.
2016).
The software package for simulating the 2D and 3D movement of water, heat, andmultiple solutes in variably-saturated media (Šimůnek and Bradford 2008;
Šimůnek et al. 2016).
The 1D unsaturated water flow package for MODFLOW based on the HYDRUS-1Dmodel (Seo et al. 2007; Twarakavi et al. 2008).
A newly-developed 1D unsaturated water flow and solute transport package for
MODFLOW and MT3DMS

table position is taken. The bottom boundary condition
thus changes stepwise in the HPM. Refer to Seo et al.
(2007) for more details.
Updates of Pressure Heads to Eliminate Spurious Fluxes
As discussed above, the pressure head bottom boundary condition changes in the HPM profile stepwise at theend of each MODFLOW time step. A sudden change inthe position of the water table at the HPM-MODFLOWinterface (i.e., the groundwater table) after every time step,and thus a change in the value of the bottom boundary
condition (a time-variable pressure head boundary condition) of the HPM profile, led to either a sudden inflow
into (a higher water table) or outflow from (a lower watertable) the soil profile and to an incorrect cumulative bottom flux over the duration of the simulation (Twarakavi
et al. 2009; Kuznetsov et al. 2012). This sudden change
in the water flux when the position of the water table
shifted was eliminated in the updated code by adjusting
the pressure head distribution in each HPM profile after
every MODFLOW time step as described below.
The updated water table level in each region at theend of the time step is obtained by the MODFLOW
simulation and is then considered as the bottom boundary
condition during the new time step in the HPM. The
bottom pressure head in each HPM profile is set according
to the new water table depth, and the pressure heads above
the soil profile bottom are adjusted so that they correspondto the bottom flux at the previous time step. While thisadjustment is straightforward in the saturated zone (dueto a linear nature of Darcy’s law), a root-finding routine
was used to adjust pressure heads in the unsaturated zone
(due to a non-linear nature of the Darcy-Buckinghamlaw). According to the discretization scheme used by
HYDRUS-1D and the HPM, the water flux between two
adjacent nodes “i ” and “i + 1” can be expressed as:q=−
NGWA.org

K (hi ) + K (hi+1 )2



hi+1 − hi+1
zi+1 − zi(3)The above equation has to be solved for h i + 1 , whilethe value h i is known and q is equal to the bottom flux.
Since K depends on h i + 1 in a strongly nonlinear manner,the equation is solved iteratively using the “false position”
method.
The steady-state nodal fluxes in the adjusted pressure
head profile are compared with the nodal fluxes obtainedat the end of the previous time step. The node, in which
the relative difference between these two fluxes is larger
than 0.1% of the flux (plus a small round off value of
10−12 ), is set as the uppermost node, in which the pressure
head profile is updated. The pressure head values below
this node are set equal to the pressure heads obtained by
the steady-state profile calculations for a given flux and agiven position of the water table. The pressure head values
above this node remain the same as the pressure heads
at the end of the previous time step. This procedure isfollowed after every MODFLOW time step in each HPM
profile to eliminate sudden fluxes due to abrupt changes
in the water table elevation.
Example Illustrating the Effects of the Pressure Head Updates
The impact of updating (or not) the pressure headvalues at the bottom of the HPM profiles when theposition of the groundwater table changes is illustrated
in the following example. A domain with a length anda width both equal to 0.75 m and a depth of 10 m wasconsidered. The MODFLOW domain was divided intonine equal grid cells (three rows and three columns) withno flow boundaries so that the position of the water table
could change only due to recharge. A constant inflow
of 0.002 m day−1 was applied at the soil surface for aduration of 600 days. The soil profile was considered tohave a hydrostatic initial pressure head distribution with awater table at a 5-m depth. The duration of the simulationof 600 days was divided into six equal MODFLOW
time steps. The HPM column was divided into 100
finite elements of equal dimensions. The maximum
considered allowable time step in HPM simulations was
0.1 day.
S. Beegum et al. Groundwater 57, no. 3: 392–408

395

(b)

Water table elevation (m)1098With pressure head modification
Without pressure head modification

7
6

5
4
3

Cumulative bottom flux (m)

(a)

0.55
0.45
0.35
0.25
0.15
0.05
-0.05

0

100 200 300 400 500 600
Time (days)

Cumulative bottom flux (m)(c)

0

100 200 300 400 500 600
Time (days)

0.00100.00050.0000-0.0005-0.0010-0.0015-0.0020-0.0025-0.0030
0 10 20 30 40 50 60 70 80 90
Time (days)

Figure 1. Water table elevations with and without pressure head modifications at each MODFLOW time step (a), the flux atthe bottom of the HPM profile without pressure head modifications (b), and the close-up view of the cumulative bottom flux
up to 90 days (c).

The van Genuchten analytical model (van Genuchten1980) was used to describe the soil hydraulic properties
with the following parameters: the residual water content,θ r = 0.01, the saturated water content, θ s = 0.3, the saturated hydraulic conductivity, K s = 8.4 m day−1 , the poreconnectivity parameter, l = 0.5, and the shape parameters,
α = 3.3 m−1 and n = 4.1. The specific yield was set equal
to 0.29 and the specific storage to 0.0015 m−1 . Note thatthe value of the specific yield in MODFLOW is definedindependently from the soil hydraulic functions used in the
HPM. While in principle, one could use the vadose zone
flow simulations to improve specific yield estimations, nosuch attempt was undertaken in this study.
In the HPM, the water table level obtained from
MODFLOW after every time step was considered as thebottom boundary condition for the vadose zone profile
simulations. The change in the water table level after
every MODFLOW time step is shown in Figure 1a.
Figure 1b shows the water flux at the bottom of the
HPM profile when the pressure heads were not modified.
The bottom flux was zero during approximately the first
45 days of the simulation, reflecting the initial hydrostatic
boundary condition. A close up of the Figure 1b up
to 90 days is shown in Figure 1c. Once the moisture
front reached the groundwater table, the bottom flux
equilibrated with the surface flux of 0.002 m day−1 .
396

S. Beegum et al. Groundwater 57, no. 3: 392–408

However, the bottom flux displayed a sudden increase(Figure 1b) after every MODFLOW time step because of
the sudden change in the bottom pressure head boundary
condition in the HPM profiles, reflecting the water tablelevels calculated by MODFLOW. This sudden change inthe bottom flux can be eliminated by adjusting the pressure
head profile after every time step using the method
described above. Figure 2a shows the flux at the bottomof the HPM profile when the pressure head profile wasmodified after each MODFLOW time step. The bottom
flux was still initially zero for about 45 days, but thenit increased to 0.002 m day−1 and remained constant atthat level during the rest of the simulation. The spurious
fluxes that occurred when the bottom boundary condition
was changed were fully eliminated by modifying thepressure head profile using the method described above.
Note that upward fluxes are positive and downward fluxes
are negative in HYDRUS-1D and HPM.
Figure 2b shows the cumulative bottom flux when thepressure head distribution either was or was not modified
after each MODFLOW time step. While the cumulativeflux for the case with the pressure head modifications
was a smooth line corresponding to the applied surface
flux of 0.002 m day−1 , the cumulative flux for the casewithout the pressure head modifications exhibited sudden
jumps corresponding to sudden influxes when the bottom
NGWA.org

(a)(b)
Cumulative bottom flux (m)0.0010Bottom flux (m)0.00050.0000-0.0005-0.0010-0.0015-0.0020-0.0025-0.0030

0
0.00
-0.20
-0.40
-0.60

-0.80
-1.00
-1.20

0

Time (days)
100 200 300 400 500 600

With pressure head
modification
Without pressure
head modification

100 200 300 400 500 600
Time (days)

Depth (m)

(c)
0.00
0
1
2
3
4
5
6
7
8
9
10

Water content (-)
0.10
0.20
0.30

Initial water content
profileFinal water content
profileFigure 2. The flux at the bottom of the HPM profile for the simulation with the pressure head modifications (a), cumulative
fluxes at the bottom of the HPM profile for simulations with and without pressure head modifications (b), and the initial andfinal water content profiles for the simulation with the pressure head modification (c).

pressure head was changed. Overall, this line has a smaller
average slope than the line corresponding to the applied
surface flux and is thus incorrect. The cumulative flux
for the simulation with the pressure head modifications
was larger than when there were no pressure head
modifications (Figure 2b). This is because of the removalof the sudden upward fluxes after every time step withan increase in the water table, which is generated when
the pressure head distribution is not modified. Figure 1ashows the water table level with and without pressure
head modifications. The water table level was found tobe higher when the pressure head profile was modified
compared to when it was not, which was also because of
the removal of upward inflow after every time step.
The initial and final water content profiles for thesimulation with the pressure head modification are shownin Figure 2c. The total inflow is 1.20 m, the total
outflow is 1.11 m, and the change in storage is equal to
0.09 m. The mass balance was considered separately for
MODFLOW and HYDRUS-1D simulations in the HPM.
The mass balance calculation in the original HYDRUS1D source code was modified to account for the amountof water added or removed from the soil profile due toupdating of the pressure head profile. The mass balancein the groundwater domain is calculated by including thecoupling flux as the recharge to the groundwater domain(Szymkiewicz et al. 2018).
In order to test the effectiveness of the updated
coupling algorithm in the case of a falling water table,a case study similar to the one above, but considering
a pumping rate of 0.002 m3 /day in the middle of thedomain, was performed. This case study considered thesame discretization and initial conditions as before.
The cumulative bottom flux for the case with thepressure head modifications was a smooth line whilesudden jumps were observed (corresponding to sudden
outflow fluxes when the bottom pressure head waschanged) when the pressure head profile was not updated.
The cumulative flux for the simulation with the pressure
head modifications was smaller than when there were nopressure head modifications (Figure 3b). This is becauseof the removal of sudden downward fluxes after every
time step with a decrease in the water table, whichis generated when the pressure head distribution is not
modified. A corresponding change in the water tablelevel can also be observed when the pressure headprofile is not modified. Water table elevations (Figure 3a)and cumulative bottom fluxes (Figure 3b) show theeffectiveness of the implementation of the new couplingalgorithm.

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

397

Time (days)

(b)
6.5
6.0
5.5
5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.0

With pressure head modification
Without pressure head modification

0

Cumulative bottom flux (m)Water table elevation (m)

(a)

0.0
-0.2
-0.4
-0.6
-0.8
-1.0
-1.2
-1.4
-1.6
-1.8

0

100 200 300 400 500 600

With pressure head
modification
Without pressure
head modification

100 200 300 400 500 600
Time (days)

Figure 3. Water table elevations (a) and cumulative bottom fluxes (b) in the HPM profile with and without pressure head
modifications after each MODFLOW time step.

The Effects of Time Steps on the Cumulative Bottom Flux
The method discussed above for eliminating sudden
fluxes at the bottom of the HPM profile when the positionof the water table changes at the HYDRUS-MODFLOWinterface was further tested for its functionality for
different MODFLOW time steps. This was done using
two case studies. Case studies 1 and 2 consider a constantand variable recharge at the soil surface, respectively.
Case Study 1: Constant Recharge: The same domain, thesame discretization, and the same initial and boundary
conditions as used above were considered again. The
effect of the time step on the cumulative bottom flux andthe water table elevation was studied by considering threedifferent numbers of MODFLOW time steps: 6, 60, and600.
Figure 4a shows the cumulative bottom flux when
the pressure head profile in the HPM either was or was
not modified for different MODFLOW time steps. Whilethe cumulative flux for all three cases with the pressure
head modifications was a smooth line corresponding to
the applied surface flux of 0.002 m day−1 , the cumulativeflux for the cases without the pressure head modifications
exhibited sudden jumps corresponding to sudden influxes
when the bottom pressure head was changed. Overall,these lines have a smaller average slope than the line
corresponding to the applied surface flux and are thus
incorrect. Cumulative fluxes for simulations with thepressure head modifications were larger (in absolutevalues) compared to those without the pressure head
modifications for all time steps. This is because of theremoval of upward inflow fluxes after every time stepwith an increase in the water table which was generated
when the pressure heads were not modified.
Figure 4b shows the water table elevation when
the pressure head distributions either were or were notmodified for different MODFLOW time steps. Water table
elevations were higher when the pressure head distribution
was modified for all three time steps, compared to when
it was not. This is because of the removal of the sudden398

S. Beegum et al. Groundwater 57, no. 3: 392–408

upward flux after every time step, which resulted in alower water table rise. The water table started rising
after about 45 days. Simulations with 60 and 600 time
steps predicted the initiation of the water table rise after
approximately the same time. In the case of 6 time steps,
the water table elevation (the bottom boundary condition)in MODFLOW remained constant for 100 days, which
created the delayed increase in the elevation of the watertable level.
This indicates the importance of choosing optimal
time steps within each MODFLOW stress period inorder to describe the water table dynamics accurately.
Guidelines for choosing the time step in case of a transient
simulation is explained in Reilly and Harbaugh (2004) andHarbaugh et al. (2000).
Case Study 2: Variable Recharge: The same domain as inprevious examples was used in this example as well. The
length and width were both equal to 0.75 m and a depthof 10 m was considered. The MODFLOW domain was
defined as such so that the position of the water table couldchange only due to recharge. Daily values of precipitationand potential evapotranspiration were considered for aduration of 365 days (Figure 5). The data were obtained
from the weather station at the Gdańsk University ofTechnology, Poland, and correspond to the year 2015. The
potential evapotranspiration values were divided into thepotential evaporation and transpiration values using theleaf area index, which was considered to follow a typicalwinter wheat growth function during the entire growth
period. Root water uptake for the winter wheat crop
was considered using the Feddes water uptake reduction
model (Feddes et al. 1978). The following values wereused for the stress response function: the pressure head
below which roots start extracting water from the soil,
P 0 = 0 m, the pressure head below which roots extract
water at the maximum possible rate, POpt = −0.005 m,
the limiting pressure head below which roots cannot
extract water at the maximum rate, P 2 = −3.0 m, andthe pressure head below which root water uptake ceases,
NGWA.org

(a)

(b)
10.0

Time (days)
100 200 300 400 500 600

Cumulative bottom flux (m)

-0.2
-0.4
-0.6
-0.8
-1.0
-1.2

6 TS, modified
6 TS, unmodified
60 TS, modified
60 TS, unmodified
600 TS, modified
600 TS, unmodified

Water table elevation (m)00.0

6 TS, modified
6 TS, unmodified
60 TS, modified
60 TS, unmodified
600 TS, modified
600 TS, unmodified

9.5
9.0
8.5
8.0
7.5
7.0

6.5
6.0
5.5
5.0
0

100 200 300 400 500 600
Time (days)

Figure 4. Cumulative fluxes at the bottom of the HPM profile (a) and water table elevations (b) obtained using the HPM fordifferent MODFLOW time steps (TS) with and without pressure head modifications.

P 3 = −5.0 m. Root water uptake was considered to belinearly decreasing from the soil surface to a depth of
1 m. Hysteresis in the soil hydraulic functions was notconsidered in this study.
The HPM profile had a depth of 5 m and a constantinitial pressure head of −0.283 m down to a depth of
3.5 m and a hydrostatic equilibrium below this depth.
The HPM profile was divided into 201 finite elements.
Relatively small finite elements (the top element was
0.05 cm) were used at the top of the column. This isdone to ensure convergence of the numerical solution
and proper estimation of evaporation fluxes since themeteorological factors are expected to cause rapid changes
in water content and pressure head gradients near thesurface.
Two soil layers were considered: a sandy soil down
to a depth of 2.5 m from the surface and loamy sandbelow this depth. The van Genuchten analytical model(van Genuchten 1980) was used to describe the soil
hydraulic properties with the following parameters: theresidual water content, θ r = 0.045 (0.057), the saturated
water content, θ s = 0.43 (0.41), the saturated hydraulicconductivity, K s = 7.128 m day−1 (3.502 m day−1 ), thepore-connectivity parameter, l = 0.5 (0.5), and the shapeparameters, α = 14.5 m−1 (12.4 m−1 ) and n = 2.68 (2.28)
for the sandy soil (loamy sand). The specific yield was
set to 0.353 and the specific storage to 0.0015 m−1 .
Cumulative bottom fluxes and water table elevations
are analyzed for different MODFLOW time steps of 10,
20, 30, and 365 days. A maximum HPM time step of
0.1 days was considered in all simulations. Figure 6aand 6b show cumulative bottom fluxes and corresponding water table elevations, respectively, obtained using
the HPM for different MODFLOW time steps when
the bottom pressure head distributions either were orwere not modified for each considered MODFLOW
time step.
Cumulative fluxes for the simulation with the pressure head modifications were higher compared to when
there were no pressure head modifications (Figure 6a).
Cumulative bottom fluxes were higher for smaller MODFLOW time steps when compared to larger MODFLOW
time steps.
For all considered MODFLOW time steps, watertable elevations were higher when the pressure headprofiles were modified compared to when they were not(Figure 6b). This is because of the removal of the upward
flux after every time step when the groundwater table rose,
which resulted in a lower water table rise. As the numberof MODFLOW time steps increased the water table also
increased.
The relationship between modified and unmodified
solutions is similar to those shown in Figure 4, when asteady infiltration rate was considered.

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

Precipitation / Potential
evapotranspiration (m/day)

0.035
Precipitation
Potential evapotranspiration

0.030
0.025
0.020
0.015
0.010
0.005

0.000
1

61

121 181 241 301 361
Time (day)Figure 5. Daily precipitation and potential evapotranspiration for the variable recharge case study.

399

(a)

Time (days)
0

(b)
1.8

50 100 150 200 250 300 350

0.00

Cumulative bottom flux (m)

-0.05

-0.10

-0.15

-0.20

10 TS, modified
10 TS, unmodified
20 TS, modified
20 TS, unmodified
30 TS, modified
30 TS, unmodified
365 TS, modified
365 TS, unmodified-0.25

Water table elevation (m)

1.6
1.4
1.2
1.0
10 TS, modified
10 TS, unmodified
20 TS, modified
20 TS, unmodified
30 TS, modified
30 TS, unmodified
365 TS, modified
365 TS, unmodified

0.8
0.6
0.4

0.2
0.0
0

50 100 150 200 250 300 350
Time (days)

Figure 6. Cumulative fluxes at the bottom of the HPM profile (a) and water table elevations (b) obtained using the HPM fordifferent MODFLOW time steps (TS) with and without pressure head modifications.

The method developed above to eliminate sudden
inflow into (for an increasing water table) or outflow
from (for a decreasing water table) the soil profile when
the position of the water table changes at the HYDRUSMODFLOW interface has been found to be functioning
well for different MODFLOW time steps and for different
surface boundary conditions.
Solute Transport
Governing Equations: Solute Transport in the Unsaturated
Soil Zone
The HYDRUS-1D model (Šimůnek et al. 2016)
simulates solute transport in variably saturated porous
media using the standard advection-dispersion transport
equation of the form:


∂s
∂
∂c
∂qc
∂θ c
+ρ
=
θD
−
−ϕ
(4)
∂t
∂t
∂z
∂z
∂z
where c is the solution concentration [ML−3 ], s is
the sorbed concentration [MM−1 ], D is the dispersion
coefficient [L2 T−1 ], ρ is the bulk density of the porous
medium [ML−3 ], q is the volumetric flux density [LT−1 ],which is obtained using the Darcy-Buckingham law, andϕ is a sink-source term accounting for various zero- andfirst-order or other reactions [ML−3 T−1 ] (Šimůnek andvan Genuchten 2008).
The solute transport subroutines from the standard
computational module of HYDRUS-1D were integrated
into the HPM to make it capable of additionally simulating
the solute transport in the unsaturated zone. Once waterflow is simulated in the unsaturated zone at each
HYDRUS-1D time step, the nodal values of the fluid flux
are determined from the nodal values of the pressure head
and the hydraulic conductivity by applying Darcy’s law.
The solute transport part of the code then uses these values
400

S. Beegum et al. Groundwater 57, no. 3: 392–408

of water contents and fluid fluxes to solve the transport
equations.
Similar to HYDRUS-1D, the solute transport subroutines incorporated into the modified HPM can consider
linear and nonlinear sorption, zero- and first-order production and degradation in liquid and solid phases, the onesite kinetic sorption model, the two-site sorption modelwith equilibrium and kinetic sorption, and the two-kinetic
site sorption model (Šimůnek and van Genuchten 2008).
Governing Equations: Solute Transport in the Saturated Soil
Zone
The Modular Three-Dimensional Multispecies Transport Model (MT3DMS) (Zheng and Wang 1999) was used
to simulate solute transport in the saturated soil zone.
MT3DMS can simulate advection, dispersion/diffusion,and chemical reactions of contaminants in the saturated
zone under general hydro-geologic conditions using thefollowing partial differential equation:



∂
∂
∂C
∂ (θ C )
θ Dij
−
=
Rn
(θ vi C) + qs Cs +
∂t
∂xi
∂xj
∂xi
(5)
where C is the dissolved solute concentration [ML−3 ], θis the porosity of the subsurface medium (dimensionless),t is time [T], x i is the distance along the respective
Cartesian coordinate axis [L], D ij is the hydrodynamic
dispersion coefficient tensor [L2 T−1 ], v i is the pore watervelocity [LT−1 ], q s is the volumetric flow rate per unit
volume of the aquifer representing fluid sources (positive)
and sinks (negative) [T−1 ], C s is theconcentration of the
R n is the chemicalsource or sink flux [ML−3 ], andreaction term [ML−3 T−1 ].
NGWA.org

Figure 7. Flowchart describing the integration of the modified HPM (including solute transport component) into MODFLOW.

Coupling Procedure
Figure 7 shows the flowchart describing the integration of the HPM (flow and solute transport) with MODFLOW and the linking of MT3DMS with this model. The
top section of the flow chart depicts the implementationof the HPM, and the bottom section depicts MODFLOW
and the linking of MT3DMS, as well as their interactions
at the water table. In MODFLOW, the entire simulation time is divided into stress periods, during which theexternal stresses are constant. These stress periods arefurther divided into smaller “MODFLOW” time steps.
For example, in Figure 7, the total simulation period is
divided into three stress periods and these stress periods
are divided into two, one, and three time steps, respectively. The 3D solute transport in the saturated soil zone
can be further simulated using the MT3DMS package.
The solution of the transport equation using MT3DMS
requires much smaller time steps compared to time steps
used in the MODFLOW simulation because of the stability criteria and the accuracy requirements (Zheng 2009).
For example, the MODFLOW time Step 2 in the stressperiod 1 is divided into four smaller transport time steps(Figure 7). The HPM performs multiple time steps withineach MODFLOW time step since the solution of the
Richards equation requires smaller time steps comparedto the time steps used by MODFLOW.
Twenty-four distinct stages are required to carry out
the simulation illustrated in the flow chart (Figure 7).
Initially, average water table levels from the MODFLOW
grids are used as bottom boundary conditions in the HPM
profiles (Stage 1). The HPM performs the water flow andsolute transport simulation in the unsaturated zone for aduration of the MODFLOW time Step 1 and calculates
the total bottom flux and the solute concentration flux
reaching the groundwater (Stage 2). The total bottom fluxis then given as recharge for the groundwater simulation
carried out by MODFLOW (Stage 3). MODFLOW solves
the groundwater flow equation for time Step 1 andcalculates the updated water table level (Stage 4). Thiswater table level is assigned as the bottom boundary
condition in the HPM (Stage 5). The solute flux can beprovided as the concentration flux input to MT3DMS tosimulate the solute transport in the saturated soil zone.
Before moving to Stage 6, the pressure head distribution at the bottom of the HPM profiles at the end of
the first time step is modified to eliminate sudden fluxes
using the method discussed above. This is followed by
the modifications of solute concentrations corresponding
to the modifications of pressure heads using the method
discussed below. The HPM again solves the solute transport and water flow for a duration of MODFLOW time
Step 2 (Stage 6). Stages 7, 11, 15, 19, and 23 perform thesame task as Stage 3. Stages 6, 10, 14, 18, and 22 performtasks similar to Stage 2, and Stages 8, 12, 16, 20, and 24
perform steps similar to Stage 4.

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

Updates of Solute Concentrations
Since the pressure head (water content) distributions
are updated after each MODFLOW time step (to eliminate
the effects of sudden fluxes at the HYDRUS-MODFLOWinterface and to fix the overall mass balance errors),
solute concentrations need to be correspondingly adjusted
to preserve the mass balance. This was done by either
lowering (diluting) or increasing (concentrating) soluteconcentrations, depending on the change in water contents
due to the modification of the pressure heads. The newconcentration was obtained by equating the total solute
mass (the mass of solute per volume of soil) before andafter the adjustment of the water content profile at each401

node where water contents were updated:ρKd cβ+ θc
1 + ηcβ




=

old

ρKd cβ+ θc
1 + ηcβ(6)
new

In Equation 6, K d , β, and η are empirical coefficientsof the adsorption isotherm that is used in HYDRUS-1D
and the modified HPM, θ is the volumetric water content,
c new is the updated solution concentration, c old is theoriginal solution concentration, and ρ is the bulk densityof the porous medium.
For linear sorption (β = 1 and η = 0), the solutionof Equation 6 is trivial, and c new is obtained by solving
Equation 7:
cnew = cold

ρKd + θold
(ρKd + θ )old
= coldρKd + θnew(ρKd + θ )new(7)

For nonlinear sorption (i.e., Langmuir or Freundlich,η > 0 or β = 1, respectively), c new is obtained by the Brent
method for finding roots of nonlinear equations (Press
et al. 1992).
Example Illustrating the Effects of the Concentration
Updates
The same example discussed above in Section“Example Illustrating the Effects of the Pressure Head
Updates” is considered here. A solute concentration of
5 mg l−1 was considered along with a constant inflow of
0.002 m day−1 . The bulk density of the soil was assumed
to be 1.5 g cm−3 with a longitudinal dispersivity 0.1 m.
Since the entire simulation for 600 days was divided
into six equal MODFLOW time steps, the concentrations after every MODFLOW time step were updated
correspondingly to the modifications of the pressure
heads and water contents. Figure 8 shows the change inthe water content profile (corresponding to the pressure
head modifications) after the fifth time step and corresponding changes in the concentration profile. The modified pressure head profile resulted in higher water contents
in the soil immediately above the water table. Hence, theupdating method lowers (dilutes) the solute concentration
correspondingly to an increase in the water content.
Performance Evaluation of the Modified HPM Using
HYDRUS (2D/3D)
The performance of the newly developed modified HPM was evaluated by comparing its resultswith the results of the 2D simulation with HYDRUS(2D/3D) (referred to below as HYDRUS-2D). Similar to
HYDRUS-1D, HYDRUS-2D simulates water flow andsolute transport in both unsaturated and saturated domains
by solving the 2D versions of the standard Richards
and convection–dispersion equations (Equations 2 and 4,
respectively).
Case Study
A hypothetical 2D domain was considered for the
performance evaluation of the modified HPM with
402

S. Beegum et al. Groundwater 57, no. 3: 392–408

Depth (m)

Relative solute concentraion (-)
Relative water content (-)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

10
9
8
7
6
5
4
3
2
1
0
Relative water content before*
Relative water content after*
Relative solute conentration before*
Relative solute conentration after*

Figure 8. The original and updated water contents andsolute concentrations in response to the modification of thepressure heads (before* and after* refers to before and after
the pressure head profile updating).

HYDRUS-2D. The domain has a width of 100 m anda depth of 10 m. This domain was considered to have
20 MODFLOW grids and 20 HPM (HP-1 through HP20) soil profiles. The initial pressure heads in the HPM
profiles were assumed to be in hydrostatic equilibrium.
The initial water table in the entire domain was set at
a depth varying from 4.0 to 6.0 m between the left andright sides of the domain. While four HPM profiles (HP-9,
HP-10, HP-11, HP-12) were assumed to have a constantrecharge of 0.04 m day−1 with a solute concentrationof 5 mg/L at the soil surface, there was no recharge
considered in the remaining HPM profiles throughout thesimulation period. Such a setting represents infiltrationfrom a pond or effluent lagoon, or artificial rechargeusing treated wastewater. The total simulation time was
100 days. The number of stress periods (MODFLOW
time steps) was also 100, each having a time step
duration of 1 day. A constant head of 6.0 and 4.0 m wasassumed at left and right ends of the domain, respectively.
Figure 9 shows the transport domain considered in ourstudy.
The van Genuchten analytical model was used todescribe the soil hydraulic properties with the following parameters: θ r = 0.01, θ s = 0.3, K s = 8.4 m day−1 ,l = 0.5, α = 3.3 m−1 , and n = 4.1. The specific yield wasset equal to 0.29 and the specific storage to 0.0015 m−1 . A
bulk density of 1.5 g/cm3 and a longitudinal dispersivityof 1 m was assumed. Solute transport without adsorption was considered here. Each 10 m deep HPM profile
was divided into 100 finite elements of equal widths. A
maximum allowable time step of 1 day was considered in
the simulation.
NGWA.org

Figure 9. A schematic showing the flow and transport domain. Cases 1, 2, 3, and 4 indicate evaluated scenarios with different
numbers of considered HPM profiles.

The same domain was modeled using HYDRUS2D. In HYDRUS-2D, the entire domain was discretizedspatially into a finite element mesh, with relatively small
elements at locations where large hydraulic gradients were
expected. An average value of 0.05 m was consideredas the size of the finite element toward the soil surface
whereas a size of 0.3 m was considered toward the bottomof the domain. The same water flow and solute transportparameters used in the modified HPM were consideredalong with a transverse dispersivity of 0.1 m. A maximum
time step of 1 day was considered here as well.
Figure 10 shows a comparison of water table variations, solute concentrations reaching the water table, andwater and solute fluxes reaching the water table calculatedwith the modified HPM (HPM-20 HP) and HYDRUS2D. Figure 10a shows a comparison of water table elevations at the center of the transport domain calculatedby HYDRUS-2D and the modified HPM during the simulation. An increase in the water table elevation is almost
the same for both models. Similarly, solute concentrations(Figure 10c) and water and solute fluxes (Figure 10b and10d) reaching the water table at the center of the domain
are also comparable in both models. This comparison
shows that the solute transport component incorporated
into the HPM is functioning well and can simulate thesolute transport in the unsaturated zone accurately.
The Effects of Spatial Discretization on the Computational
Time and Accuracy of the Results
A significant saving of computational time can beachieved by reducing the number of HPM profiles considered in the unsaturated soil zone if this can be achieved
without compromising the accuracy of the results. The
same example problem discussed in the previous section
was further analyzed for the computational time requirement and the accuracy of the results by varying the number of HPM profiles. Four cases were considered in theanalysis. In Case 1, the domain had 20 HPM profiles(HP-1, HP-2, HP-3, . . . , up to HP-20) (Figure 9). Case 2
considered 10 HPM profiles (HP-1, HP-2, HP-3, . . . , up
to HP-10), Case 3 considered only 5 HPM profiles (HP1, HP-2, HP-3, HP-4 and HP-5), and Case 4 considered

only 1 HPM profile (HP-1). The HP-9, HP-10, HP-11,
and HP-12 profiles in Case 1, the HP-5 and HP-6 profiles
in Case 2, the HP-3 profile in Case 3, and the HP-1 profile in Case 4 were assumed to have a constant rechargeof 0.04 m day−1 with a solute concentration of 5 mg/L
at the soil surface. No recharge was considered in other
HPM profiles during the entire simulation period in all
three cases. In Case 4, no HPM profiles were requiredin the part of the domain, which did not receive any
recharge.
Figure 10 shows a comparison of water table variations, solute concentrations reaching the water table, andwater and solute fluxes reaching the water table in the center of the domain obtained in four scenarios analyzed by
the modified HPM and HYDRUS-2D.Water table elevations, bottom fluxes, solute concentrations, and concentration fluxes simulated in all four modified HPM scenariosand by HYDRUS-2D are almost the same. In Figure 10a,the water table elevation in the center of the domain was
found to be slightly lower when a smaller number of HPM
profiles was considered. This is because of the averaging of the water table elevation in the MODFLOW gridcells assigned to the same HPM profiles. The averagingof the water table results in a lower level of the watertable, which is given as a bottom boundary condition
to the HPM profile. For example, in Case 1, in which
each MODFLOW grid was assigned one HPM profile,
the water table level obtained by MODFLOW is given to
the corresponding HPM profile without any averaging. In
Case 2, the water table averaged after every time step forseveral MODFLOW grid cells (i.e., 1, 2; 3, 4; 5, 6; 7,
8; 9, 10; 11, 12; 13, 14; 15, 16; 17, 18; and 19, 20) are
given as the bottom boundary condition to the HPM profiles (HP-1, HP-2, HP-3, HP-4, HP-5, HP-6, HP-7, HP-8,
HP-9, and HP-10, respectively). In Case 3, the averagewater table in MODFLOW grid cells (1, 2, 3, 4; 5, 6, 7,
8; 9, 11, 12; 13, 14, 15, 16; and 17, 18, 19, 20) are given
as the bottom boundary condition to the HP-1, HP-2, HP3, HP-4 and HP-5 profiles, respectively. In Case 4, theaverage water table from MODFLOW grids (9, 10, 11,
and 12) are given as the bottom boundary condition to
the HP-1 profile.

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

403

(a)

Time (days)(b)

0 10 20 30 40 50 60 70 80 90 100
0.00

5.4
5.3
5.2
5.1

5.0

HYDRUS-2D
HPM-20 HP
HPM-10 HP
HPM-5 HP
HPM-1 HP

Water flux (m/day)

Water table elevation (m)

5.5

-0.01
-0.02

HYDRUS-2D
HPM-20 HP
HPM-10 HP
HPM-5 HP
HPM-1 HP

-0.03
-0.04
-0.05

0 10 20 30 40 50 60 70 80 90 100
Time (days)

(d)

5
4
3
2
1
0

HYDRUS-2D
HPM-20 HP
HPM-10 HP
HPM-5 HP
HPM-1 HP
0 10 20 30 40 50 60 70 80 90 100
Time (days)

Solute flux (mg/m2/day)

Solute concentration (mg/l)

(c)

0.20
0.15
0.10
0.05
0.00

HYDRUS-2D
HPM-20 HP
HPM-10 HP
HPM-5 HP
HPM-1 HP
0 10 20 30 40 50 60 70 80 90 100
Time (days)

Figure 10. Water table elevations (a), water fluxes (b), solute concentrations (c), and solute fluxes (d) at the water table inthe center of the domain obtained using HYDRUS-2D and the modified HPM with different numbers of profiles.

The total computational time (on a PC with theIntel i7 processor and the Windows 7 operating system)
required in Case 1 was 17.36 s, whereas in Cases 2, 3,
and 4 it was 9.31, 4.42, and 2.01 s, respectively. The
total number of nodes in the HPM profiles in Cases 1,
2, 3, and 4 were (20 profiles*101 nodes) 2020, (10*101)
1010, (5*101) 505, and (1*101)101, respectively, alongwith 20 MODFLOW grids. This has resulted in asignificant decrease in the computational time requirementof the simulations. This shows that we can reduce thecomputational time of the modified HPM using a proper
spatial discretization and the allocation of HPM profiles
to each MODFLOW grid.
The Use of MT3DMS for the Saturated Zone Solute
Transport Modeling
Recharge and Solute Flux from the Modified HPM as Input
to MT3DMS
The solute transport in the saturated zone is simulated using MT3DMS, which is not fully coupled with
MODFLOW and the modified HPM. Since the water andsolute concentration fluxes reaching the groundwater are
obtained from the modified HPM, this information can be
fed into MT3DMS to simulate the solute transport in thesaturated soil zone. The MT3DMS code is a standalone
404

S. Beegum et al. Groundwater 57, no. 3: 392–408

transport simulation model, which can be used with any
finite difference flow model. MT3DMS is often used inconjunction with the widely accepted groundwater flow
model MODFLOW. LMT (Zheng et al. 2001) is an addon package for MT3DMS that has been developed tosave the flow solution (cell by cell heads and fluxes, andlocations and flow rates from sources and sinks solved
by MODFLOW) that is required for the simulation of
the solute transport in the saturated soil zone. The LMT
package is used for linking the results of MODFLOW
with the modified HPM and MT3DMS. The solute concentration flux reaching the water table calculated by themodified HPM is provided as the concentration flux input
to MT3DMS to simulate the solute transport in the saturated soil zone. Since MT3DMS is not fully integrated
into the MODFLOW-HYDRUS program, the solute transport simulation can be carried out using the LMT package(Zheng et al. 2001), which allows MODFLOW to store
the flow information (cell by cell heads and fluxes, andlocations and flow rates from sources and sinks solved
by MODFLOW) for each time step in the saturated zone,
which is then further used by MT3DMS.
Case Study
The case study discussed in the previous section
is further analyzed in combination with MT3DMS. The
NGWA.org

(a)

(b)

Figure 11. Concentration profiles in the transport domain at the end of the simulation obtained using MT3DMS (left) andHYDRUS-2D (right). Two vertical profiles (1 and 2) in the left figure are used in Figure 12.

4.52
5.07
5.62
6.16
6.71
7.26
7.81
8.36
8.90
9.45
10.00

Solute concentration (mg/l)
0
1
2
3
4

(b)
0

5

Section 1HYDRUS-2D
Section 1MT3DMS

Depth below the water table (m)Depth below the water table (m)

(a)

5.03
5.53
6.03
6.52
7.02
7.52
8.01
8.51
9.01
9.50
10.00

Solute concentration (mg/l)
1
2
3
4
5
Section 2HYDRUS-2D
Section 2MT3DMS

Figure 12. A comparison of solute concentrations in the saturated zone obtained by MT3DMS and HYDRUS-2D on verticalprofiles defined in Figure 11.

solute concentration flux at the water table at each time
step obtained using the modified HPM is considered as
the concentration flux input for the MT3DMS simulationsof the saturated soil zone. The same solute transportparameters as used in HYDRUS-2D were also usedin MT3DMS. An effective porosity of 0.3 and alongitudinal dispersivity of 1 m was considered. The
ratio of the horizontal transverse dispersivity to thelongitudinal dispersivity was assumed as 0.1, and the ratioof the vertical transverse dispersivity to the longitudinal
dispersivity was assumed as 0.01.
Figure 11 shows the concentration profiles in thetransport domain obtained by MT3DMS and HYDRUS2D at the end of the simulation. Two cross sections(Figure 11, left), one at the center of the domain (1) andone at a distance of 20 m from the center of the domain(2), were considered in the further analysis. Figure 12a
and 12b shows the comparison of concentrations along
these two cross sections obtained using MT3DMS andHYDRUS-2D at the end of the simulation.
The solute concentration below the water table (thewater table is 4.52 m below the soil surface in Section
Introduction and 5.033 m in Section The HYDRUS Package for MODFLOW: Modeling Approach) is comparable in both models. The concentration in the saturated

zone obtained using MT3DMS is influenced by various
factors such as horizontal and vertical transverse dispersivity, the spatial and temporal discretization, and thewetting factor used in MODFLOW for converting thedry cell to wet cell. The differences in the concentration profiles obtained using MT3DMS and HYDRUS2D are mainly caused by the differences in the spatial
discretization used. While rectangular grids were considered in the MT3DMS simulation, the finite element
mesh in HYDRUS-2D was constructed by dividing theflow region into quadrilateral and/or triangular elements.
These discretizations are not the same due to different
solution techniques and convergence criteria used in bothmodels.

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

Computational Time Requirements of the Coupled Modeland Its Components
The computational efficiency of the developed model
was analyzed by comparing the total computational timeused by (1) MODFLOW with the modified HPM, (2) the
modified HPM itself, and (3) the HYDRUS-2D modelusing the case study discussed above.
The total simulation time was 100 days. The computer used for running the model was Intel i7 processor
with Windows 7 operating system. When the number of405

Table 2
Computational Time of the Entire MODFLOW
Model with the Modified HPM and the Vadose
Zone Simulations in the Modified HPM for
Different MODFLOW Time Steps
Computational Time (Seconds)
No. of MODFLOWTime Steps
5
10
20
40
60
80
100

The
Modified
HPM

MODFLOW
with theModified HPM

4.92
5.12
5.53
6.52
9.15
12.80
14.95

5.67
5.89
6.52
7.82
11.10
14.98
17.36

considered time steps was 100 (with each MODFLOW
time step duration of 1 day), the total computational timerequired for MODFLOW with the modified HPM was
17.36 s, from which the modified HPM required 14.95 s.
This shows that most computational time required by themodel was used for simulating processes in the unsaturated zone. This is because of the increased computational time requirements for the numerical solution of the
Richards equation for simulating water flow in the unsaturated soil zone when compared to the solution of thegroundwater flow equation for the saturated soil zone.
Note that in this case study, the total number of nodes
in the HPM profiles in the unsaturated zone was 2020
(20 profiles * 101 nodes), whereas the saturated zone had
only 20 MODFLOW cells, which well explains the higher
computational time requirement of the modified HPM.
When the same case study was solved using
HYDRUS-2D (with a maximum allowable time step of
1 day), the computational time was 295.14 s. HYDRUS2D required more computational time when compared to
MODFLOW with the modified HPM.
The total computational time of the MODFLOW
model coupled with the modified HPM and the HPM
depends on the number of MODFLOW time steps(Table 2). In this study, the total simulation time isdivided into MODFLOW time steps of equal duration.
The computational time of MODFLOW with the modified
HPM is found to be increasing with an increase in
the number of time steps. It can also be observed thatthe vadose zone simulation time in the modified HPM
constitutes a major percentage of the total computationaltime for different time steps.
It can also be observed in Table 2 that the overall
computational time, and especially the modified HPM
computation time, is increasing only slowly when there
are fewer MODFLOW time steps (5, 10, 20, and 40).
However, with a further increase in the number of
MODFLOW time steps, the computational time starts
increasing exponentially. This is likely due to temporal
406

S. Beegum et al. Groundwater 57, no. 3: 392–408

discretization criteria used in the HPM. The numberof time steps in the modified HPM is almost the same
when the number of the MODFLOW time steps is small(less than 40 in this case) and larger MODFLOW time
steps have only a limited impact on HPM time steps.
However, when the number of MODFLOW time steps
is large, i.e., MODFLOW time steps are small, the HPM
is required to take smaller time steps as well and their
number increases. More information on the time step
control can be found in Šimůnek et al. (2013).

Summary and Conclusions
The major limitation of the HPM was that it generated
instabilities in water fluxes at the saturated-unsaturated
interface when the position of the water table changed.
This limitation is addressed by updating the couplingalgorithm. The new algorithm has been tested on a large
number of examples and found to be functioning well fordifferent MODFLOW time steps. The modified HPM does
not create spurious fluxes at the bottom of the soil profiles
for any temporal discretization used in MODFLOW.
Since the original HPM did not have the capabilities
of simulating solute transport, the HPM was modified
by incorporating the solute transport component of
HYDRUS-1D and linking it with MT3DMS. The linkingof the HPM with MT3DMS requires additional processingof information from MODFLOW and the HPM before
MT3DMS can be run.
The performance of the entire modeling system
involving MODFLOW, the modified HPM, and MT3
DMS, proposed in this study, was evaluated by comparing
its predictions of pressure heads, solute concentrations,and water and solute fluxes at different locations in
the transport domain with the results obtained using
HYDRUS (2D/3D). Since the simulation results were inclose agreement, the proposed modeling system can beused to comprehensively simulate water flow and solutetransport processes in the saturated-unsaturated zones. The
future research can focus on tighter coupling of MODFLOW with the modified HPM and MT3DMS so thatuser interventions can be fully eliminated, The analysis
of computational time requirements of MODFLOW withthe modified HPM shows that most CPU time is used by
the HPM to simulate processes in the unsaturated zone andonly a small fraction of time is used by MODFLOW itself.
Therefore, by choosing appropriate MODFLOW time
steps and spatial discretizations, one can achieve a considerable reduction in the computational time requirement.

Acknowledgments
The first author S.B. would like to acknowledge
the support from the United States India Education
Foundation (USIEF) through Fulbright-Nehru Doctoral
Research Fellowship for the research stay at University ofCalifornia Riverside for pursuing this work. A.S. has been
supported by National Science Centre, Poland, through the
NGWA.orgproject 2015/17/B/ST10/03233 “Groundwater recharge on
outwash plain.”

Brooks, R., and A. Corey. 1964. Hydraulic properties of porous
media. Hydrology Papers. Colorado State University, 37.
https://doi.org/10.13031/2013.40684
Brunner, P., and C.T. Simmons. 2012. HydroGeoSphere: A
fully integrated, physically based hydrological model.
Groundwater 50, no. 2: 170–176. https://doi.org/10.1111/
j.1745-6584.2011.00882.x
Durner, W. 1994. Hydraulic conductivity estimation
for soils with heterogeneous pore structure. Water
Resources Research 30: 211–223. https://doi.org/10.1029/
93WR02676
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978Simulation of
Field Water Use and Crop Yield . New York: John Wiley &Sons.
van Genuchten, M.T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Science Society of America Journal 44, no. 5: 892–898.
Harbaugh, A.W. 2005. MODFLOW-2005, The U.S. Geological
Survey Modular Ground-Water Model—The Ground-Water
Flow Process, 253. U.S. Geological Survey Techniques andMethods. U.S. Geological Survey Techniques and Methods,
6-A16 . Reston, Virginia: US Geological Survey.
Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald.
2000. MODFLOW-2000, The US Geological Survey Modular Groundwater Model: User Guide to Modularization
Concepts and the Groundwater Flow Process, 121. Boulder,
Colorado: U.S. Geological Survey.
Kalbacher, T., J.O. Delfs, H. Shao, W. Wang, M. Walther,
L. Samaniego, C. Schneider, R. Kumar, A. Musolff, F.
Centler, F. Sun, A. Hildebrandt, R. Liedl, D. Borchardt, P.
Krebs, and O. Kolditz. 2012. The IWAS-ToolBox: Software
coupling for an integrated water resources management.
Environmental Earth Sciences 65, no. 5: 1367–1380.
https://doi.org/10.1007/s12665-011-1270-y
Kolditz, O., S. Bauer, L. Bilke, N. Böttcher, J.O. Delfs,
T. Fischer, U.J. Görke, T. Kalbacher, G. Kosakowski,
C.I. McDermott, C.H. Park, F. Radu, K. Rink, H. Shao,
H.B. Shao, F. Sun, Y.Y. Sun, A.K. Singh, J. Taron,
M. Walther, W. Wang, N. Watanabe, Y. Wu, M. Xie,
W. Xu, and B. Zehner. 2012. OpenGeoSys: An opensource initiative for numerical simulation of thermo-hydromechanical/chemical (THM/C) processes in porous media.
Environmental Earth Sciences 67, no. 2: 589–599. https://doi.org/10.1007/s12665-012-1546-x
Kosugi, K.I. 1996. Lognormal distribution model for unsaturated
soil hydraulic properties. Water Resources Research 32, no.
9: 2697–2703.
Kuznetsov, M., A. Yakirevich, Y.A. Pachepsky, S. Sorek,
and N. Weisbrod. 2012. Quasi 3D modeling of waterflow in vadose zone and groundwater. Journal of Hydrology 450–451: 140–149. https://doi.org/10.1016/j.jhydrol.
2012.05.025
Leterme, B., M. Gedeon, and D. Jacques. 2013. Groundwater
recharge modeling of the Nete catchment (Belgium) using
the HYDRUS-1D—MODFLOW package. In Proceedingsof the 4t h International Conference HYDRUS Software
Applications to Subsurface Flow and Contaminant Transport Problems, ed. J. Šimůnek, M.Th. van Genuchten, andR. Kodešová, 235–243. Prague, Czech Republic.
Leterme, B., Gedeon, M., Laloy, E., and Rogiers, B. 2015.
Unsaturated flow modeling with HYDRUS and UZF:

Calibration and intercomparison. MODFLOW and More
2015, Golden, CO, 5 .
Morway, E.D., R.G. Niswonger, C.D. Langevin, R.T. Bailey,
and R.W. Healy. 2013. Modeling variably saturated subsurface solute transport with MODFLOW-UZF and MT3DMS.
Groundwater 51, no. 2: 237–251. https://doi.org/10.1111/
j.1745-6584.2012.00971.x
Niswonger, R.G., D.E. Prudic, and S.R. Regan. 2006. Documentation of the Unsaturated-Zone Flow (UZF1) Package
for Modeling Unsaturated Flow Between the Land Surface
and the Water Table With MODFLOW-2005 . U.S. Geological Survey Techniques and Methods, 71. Reston, Virginia:
U.S. Department of the Interior.
Panday, S., and P.S. Huyakorn. 2008. MODFLOW SURFACT:
A state-of-the-art use of vadose zone flow and transport
equations and numerical techniques for environmental
evaluations. Vadose Zone Journal 7, no. 2: 610. https://doi.org/10.2136/vzj2007.0052
Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery.
1992. Numerical Recipes: The Art of Scientific Computing,
2nd ed. New York: Cambridge University Press.
Ragab, R., and J. Bromley. 2010. IHMS-integrated hydrologicalmodelling system. Part 1. Hydrological processes andgeneral structure. Hydrological Processes 24, no. 19:
2663–2680. https://doi.org/10.1002/hyp.7681
Refsgaard, J.C., and B. Storm. 1995. MIKE SHE. In Computer
Models of Watershed Hydrology, ed. V.P. Singh, 809–846.
Colorado: Water Resources Publications.
Reilly, T.E., and A.W. Harbaugh. 2004. Guidelines for Evaluating Ground-Water Flow Models, 30. Scientific Investigations Report 2004–5038. Reston, Virginia: US Departmentof the Interior, US Geological Survey.
Samaniego, L., R. Kumar, and C. Jackisch. 2011. Predictions
in a data-sparse region using a regionalized grid-basedhydrologic model driven by remotely sensed data. Hydrology Research 42, no. 5: 338–355. https://doi.org/10.2166/
nh.2011.156
Schneider, C.L., S. Attinger, J.-O. Delfs, and A. Hildebrandt.
2010. Implementing small scale processes at the soil-plantinterface—The role of root architectures for calculatingroot water uptake profiles. Hydrology and Earth System
Sciences 14: 279–289. https://doi.org/10.5194/hess-14279-2010
Seo, H., J. Šimůnek, and E. Poeter. 2007. Documentationof the HYDRUS package for MODFLOW-2000, the US
Geological Survey modular ground-water model. GWMI
2007-01. In International Groundwater Modeling Center,
1–98. Golden, Colorado: Colorado School of Mines.
Šimůnek, J., and S.A. Bradford. 2008. Vadose zone modeling:
Introduction and importance. Vadose Zone Journal 7, no.
2: 581–586. https://doi.org/10.2136/vzj2008.0012
Šimůnek, J., and M.T. van Genuchten. 2008. Modeling nonequilibrium flow and transport processes using HYDRUS.
Vadose Zone Journal 7, no. 2: 782–797.
Šimůnek, J., M. Šejna, H. Saito, M. Sakai, and M.T. van
Genuchten. 2013. The Hydrus-1D Software Package forSimulating the Movement of Water, Heat, and Multiple
Solutes in Variably Saturated Media, Version 4.17, HYDRUS
Software Series 3 , 342. Riverside, California: Departmentof Environmental Sciences, University of California Riverside.
Šimůnek, J., M.T. van Genuchten, and M. Šejna. 2016. Recent
developments and applications of the HYDRUS computer
software packages. Vadose Zone Journal 15, no. 7: 0–25.
https://doi.org/10.2136/vzj2016.04.0033
Szymkiewicz, A., A. Gumuła-Kawȩcka, J. Šimůnek, B. Leterme,
S. Beegum, B. Jaworska-Szulc, M. Pruszkowska-Caceres,
W. Gorczewska-Langner, R. Angulo-Jaramillo, and D.
Jacques. 2018. Simulation of the freshwater lens rechargeusing the Hydrus and SWI2 packages for Modflow. Journal

NGWA.org

S. Beegum et al. Groundwater 57, no. 3: 392–408

Authors’ Note:of interest.

The author(s) does not have any conflicts

References407

of Hydrology and Hydromechanics 66, no. 2: 246–256.
https://doi.org/10.2478/johh-2018-0005
Thoms, R., R. Johnson, and R. Healy. 2006. User’s guide to
the variably saturated flow (VSF) process for MODFLOW.
Methods, no. 6-A18: 1–58.
Twarakavi, N.K.C., J. Šimůnek, and S. Seo. 2008. Evaluatinginteractions between groundwater and vadose zone using
the HYDRUS-based flow package for MODFLOW. Vadose
Zone Journal 7, no. 2: 757–768. https://doi.org/10.2136/
vzj2007.0082
Twarakavi, N.K.C., J. Šimůnek, and S. Seo. 2009. Reply to“comment on ‘evaluating interactions between groundwaterand vadose zone using the HYDRUS-based flow package
for MODFLOW’” by N. K. C. Twarakavi, J. Šimůnek, and
S. Seo. Vadose Zone Journal 8, no. 3: 820–821.
Vogel, T., and M. Cislerova. 1988. On the reliability ofunsaturated hydraulic conductivity calculated from themoisture retention curve. Transport in Porous Media 3, no.
1: 1–15. https://doi.org/10.1007/BF00222683
Zheng, C. 2009. Recent developments and future directionsfor MT3DMS and related transport codes. Groundwater
47, no. 5: 620–625. https://doi.org/10.1111/j.1745-6584.
2009.00602.x

408

S. Beegum et al. Groundwater 57, no. 3: 392–408

Zheng, C., and Wang, P. P. 1999. MT3DMS, a modular threedimensional multi-species transport model for simulationof advection, dispersion and chemical reactions of contaminants in groundwater systems—Documentation and user’s
guide. Contract Report SERDP-99-1. Vicksburg, Mississippi: U.S. Army Engineer Research and Development
Center, 202.
Zheng, C., Hill, M. C., and Hsieh, P. A. 2001. MODFLOW2000, the US Geological Survey Modular Ground-Water
Model: user guide to the LMT6 package, the linkage with
MT3DMS for multi-species mass transport modelling. No.
2001-82. U.S. Geological Survey.
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 Science and Engineering 4, no. 4: 357–373. https://doi.org/10.3882/j.issn.1674-2370.2011.04.001
Zhu, Y., L. Shi, J. Wu, M. Ye, L. Cui, and J. Yang. 2016.
Regional quasi-three-dimensional unsaturated-saturated
water flow model based on a vertical-horizontal splitting
concept. Water (Switzerland) 8, no. 5: 195. https://doi.org/
10.3390/w8050195

NGWA.org