Environmental Modelling & Software 26 (2011) 1298e1308

Contents lists available at ScienceDirect

Environmental Modelling & Software
journal homepage: www.elsevier.com/locate/envsoft

Implementation of quadratic upstream interpolation schemes for solutetransport into HYDRUS-1D 
L.E. Neumann a, *, J. Sim
unek b, F.J. Cook c, d, e, fa

CSIRO Land and Water, 37 Graham Rd, Highett, VIC 3190, Australia
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USAc
CSIRO Land and Water, Indooroopilly, QLD 4068, Australiad
The University of Queensland, St Lucia, QLD 4068, AustraliaeeWater Cooperative Research Centre, AustraliafIrrigation Futures Cooperative Research Centre, Australiaba r t i c l e i n f o

a b s t r a c t

Article history:
Received 24 November 2010
Received in revised form
12 May 2011
Accepted 14 May 2011
Available online 24 June 2011

Numerical solution of the advectionedispersion equation, used to evaluate transport of solutes in porous
media, requires discretization schemes for space and time stepping. We examine use of quadratic
upstream interpolation schemes QUICK, QUICKEST, and the total variation diminution scheme ULTIMATE,
and compare these with UPSTREAM and CENTRAL schemes in the HYDRUS-1D model. Results for purely
convective transport show that quadratic schemes can reduce the oscillations compared to the CENTRAL
scheme and numerical dispersion compared to the UPSTREAM scheme. When dispersion is introduced
all schemes give similar results for Peclet number Pe < 2. All schemes show similar behavior for nonuniform grids that become ﬁner in the direction of ﬂow. When grids become coarser in the directionof ﬂow, some schemes produce considerable oscillations, with all schemes showing signiﬁcant clippingof the peak, but quadratic schemes extending the range of stability tenfold to Pe < 20. Similar results
were also obtained for transport of a non-linear retarded solute transport (except the QUICK scheme) andfor reactive transport (except the UPSTREAM scheme). Analysis of transient solute transport show that all
schemes produce similar results for the position of the inﬁltration front for Pe ¼ 2. When Pe ¼ 10, the
CENTRAL scheme produced signiﬁcant oscillations near the inﬁltration front, compared to only minor
oscillations for QUICKEST and no oscillations for the ULTIMATE scheme. These comparisons show that
quadratic schemes have promise for extending the range of stability in numerical solutions of solutetransport in porous media and allowing coarser grids.
Ó 2011 Elsevier Ltd. All rights reserved.

Keywords:
Quadratic interpolation
TVD scheme
Solute transport
HYDRUS

1. Introduction
The transport of solutes in soils is usually described using theadvectionedispersion equation, with extra terms to account for 
adsorption and reactions (Simunek and van Genuchten, 2006;
Zheng and Bennett, 2002). The form of these extra terms varies
according to the constituents and soil types and can take forms
such as linear or non-linear adsorption, or ﬁrst- and second-order
reactions. The one-dimensional convectionedispersion equationof a solute is described by:



vC
vðuCÞ
v
vC
D
þS
¼ 
þ
vt
vx
vx
vx

* Corresponding author. Tel.: þ61 3 9252 6224; fax: þ61 3 9252 6288.
E-mail address: luis.neumann@csiro.au (L.E. Neumann).
1364-8152/$ e see front matter Ó 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/j.envsoft.2011.05.010(1)where t is time [T], C is a solute concentration [M L3], u is thevelocity of the ﬂuid [L T1], D is the dispersion coefﬁcient [L2 T1], x
is the space dimension [L], and S(C, x, t) is a source/sink term,
accounting for reactions [M L3 T1]. Although in the examples
described below more complex forms of the convectionedispersion equation are used (e.g. with variable water contents
and/or linear/non-linear retardation factors), the simpler formshown in Eq. (1) is used here to simplify the description of thenumerical methods.
For most problems of interest, when water ﬂow is transient, andsoil properties and initial conditions are non-uniform, analyticalsolutions are generally not available and/or cannot be derived, so
numerical methods must be employed for the solution of thegoverning equations (Vanderborght et al., 2005). Numericalmethods are in general more suitable for solving practical problems
involving complicated geometries that reﬂect complex naturalgeologic and hydrologic conditions, spatially and temporarily

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

variable ﬂow and transport parameters, more realistic initial andboundary conditions, and/or non-linear constitutive relationships
 (Simunek and van Genuchten, 2006). However, numericalmodeling of advection dominated problems in the presence of
discontinuities and steep concentration fronts is not trivial(Leonard, 1979, 1991) and robust numerical methods are still beingsought.
A large number of methods are available to solve Eq. (1)numerically, but they can be broadly grouped as Eulerian,
Lagrangian, and mixed LagrangianeEulerian methods. In the
Eulerian methods, the transport equation is discretized by meansof a ﬁnite difference or ﬁnite element method on a ﬁxed grid
system, while for the Lagrangian approach the mesh moves alongwith the ﬂow or remains ﬁxed in a deforming coordinate system.
In the Lagrangian-Eulerian approach, a two-step method is used.
The ﬁrst step is to estimate convective transport using particle
trajectories in a Lagrangian approach, with all other processesmodeled with an Eulerian approach in the following step (Simunek, 2005). While Lagrangian methods are widely used
when simulating transport in groundwater, the most common
methods used to solve Eq. (1) for transport in unsaturated zone
are Eulerian, including methods of ﬁnite differences, ﬁnite
elements, and/or ﬁnite volumes, and these are the methods solely
discussed here. Many groundwater ﬂow problems are either
steady-state problems or problems with gradually changing ﬂow
ﬁelds in time, and use relatively coarse spatial discretization. The
Lagrangian methods are either well suited for such conditions(former) or even required (latter), because of the stability
constraints on Eulerian methods. On the other hand, the ﬂow
problems in the vadose zone are inherently transient with quickly
changing water contents and velocities in time and space. Eulerian methods are often considered suitable for such conditions,since the numerical solution of the highly non-linear Richards
equation often requires stricter constrains on temporal andspatial discretization than groundwater ﬂow problems.
The ﬁnite difference method can be used to solve the advectionedispersion equation either using backward, forward, or
central temporal differencing. However, for certain problems suchas convection-dominated transport or the transport of steep fronts
this method can lead to artiﬁcial oscillations (under or overshooting) or numerical dispersion due to truncation errors of thediscretization. The use of forward differencing in the temporal
discretization (explicit method) can often lead to non-convergent
schemes. If a backward scheme (implicit formulation) is chosen,
the method is always convergent but can result in numericaloscillations in the form of “wiggles” with under and over shootingof the solution (Leonard, 1979). The use of central differencing (suchas CrankeNicholson schemes) improves accuracy but is morecomputational expensive than explicit methods and is still prone tonumerical oscillation. Numerical oscillation can be minimized by
the use of upstream weighting, but this can lead to considerable
numerical dispersion owing to truncation errors (e.g. Zheng andBennett, 2002).
One alternative is the introduction of an apparent numerical
diffusion coefﬁcient, to artiﬁcially damp the numerical dispersion
introduced by the discretization (Huang et al., 1997; Moldrup et al.,
1994; van Genuchten and Gray, 1978), but this can only be used forsimple problems where the numerical dispersion coefﬁcient can be
estimated (Zheng and Bennett, 2002). The introduction of artiﬁcial
damping also affects the accuracy of the method (Leonard, 1979).
Another solution to problems with artiﬁcial oscillations is the use ofﬁner grids, with a choice based on the dimensionless Peclet
number:Pe ¼ u$Dx=D

(2)

1299

where Dx is the grid spacing. Spatial discretization resulting ina Pe number smaller than 2 can eliminate numerical oscillations,while a Pe number smaller than 10 can greatly reduce such
oscillations. However, the associated computational cost due toexcessively ﬁne grids may become impractical in some applications (Huyakorn and Pinder, 1986) usually associated withgroundwater ﬂow, multi-dimensional calculations, or multicomponent solute transport.
Two methods that potentially can be used to overcome theseoscillation and truncation problems were proposed by Leonard(1979), using quadratic upstream interpolation to solve theadvection diffusion equation. These explicit methods are known as
QUICK (Quadratic Upstream Interpolation for Convective Kinematics) and QUICKEST method (QUICK with Estimated Upstream
Terms). Both schemes have little numerical dispersion, and thelater scheme has a large (in comparison to other methods) stability
region and shows minimal oscillations.
In many applications the presence of even minimal oscillations(such as negative concentrations) can corrupt the solution. As
noted by Leonard (1991), there exists a large family of schemes that
aim to suppress such oscillations, commonly referred to as Total
Variation Diminution (TVD) schemes (e.g. MUSCL (Van Leer, 1979);
Superbee (Roe, 1985)). Leonard (1991) derived a scheme to improve
the solution near steep gradients and remove under and over shootproblems by preserving local monotonicity. This scheme is known
as ULTIMATE (the Universal Limiter for Transient Interpolative
Modeling of Advective Transport Equation). While these threemethods by Leonard (1979, 1991) have been used widely in a broad
range of convectionedispersion problems in hydrodynamics andcomputational ﬂuid mechanics (e.g. Cole and Wells, 2006; Lin andFalconer, 1997; Romero et al., 2004), the use of TVD schemes in thesubsurface hydrology has been limited to only few studies. In
porous media the non-monotonic QUICK and QUICKEST schemes
have been used in groundwater (Carter et al., 1984; Lin and Medina,
2003) and QUICK in soil (Chu and Marino, 2007), with the ULTIMATE scheme used for groundwater (Zheng and Wang, 1999) andsoils (Neumann et al., 2009).
The objective of this study is to demonstrate the implementation of the explicit QUICK, QUICKEST, and ULTIMATE 
schemes into HYDRUS-1D (Simunek et al., 2008a), one of the mostwidely used ﬂow and transport programs in vadose zonehydrology. The implementation of the quadratic schemes isveriﬁed for a series of problems of increasing complexity, against
analytical solutions for purely convective transport problems andagainst standard implicit schemes in HYDRUS-1D for non-linear
transport or transient ﬂow problems. Performance of thenumerical schemes for uniform and non-uniform grids, uniform
and non-uniform velocities, reactive transport and other conditions will be discussed. The results will demonstrate theimprovements provided by some of the quadratic schemes forsituations of larger Peclet numbers or steep fronts, where nonoscillatory behavior is attained.
2. Methods
In this section the QUICK, QUICKEST, and ULTIMATE schemes, and the main ideas
behind them, are brieﬂy introduced. Complete descriptions are presented by
Leonard (1979, 1991), while the other methods implemented in HYDRUS-1D are also 
discussed in detail elsewhere (Simunek et al., 2008a). The QUICK and QUICKEST
schemes, both explicit methods, for the solution of Eq. (1) are derived using
a control-volume approach for spatial discretization and ﬁnite differences approachfor temporal discretization of Eq. (1). This results in the requirement to calculate thediscrete ﬂuxes entering and leaving the control-volume, or the so-called face values,
associated with a node within the grid. The control-volume for node i, and upstream
and downstream nodes are deﬁned as shown in Fig. 1.
Using deﬁnitions from Fig. 1, the concentration at a nodal point i for a new timelevel (j þ 1) can be calculated as:

1300

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

Fig. 1. Schematic showing a control-volume for node i (indicated by the dashed lines),
concentrations C at upstream and downstream nodes, grid spacing and face velocities(ul, ur). Subscripts FL, L, C, R, and FR refer to the far left, left, central, right, and far rightnodes. Subscripts l and r refer to the left and right faces.

Cijþ1 ¼ Cij þ



Dt
vC
u C *  ur Cr*  Dl
DxC l l
vx

*
l


þ Dr

vC
vx

* !

þ DtS

(3)

r

where superscripts j and j þ 1 refer to previous and new time levels, Dt is the time
step, Cl* and Cr* are concentrations at the left (l) and right (r) faces of the ﬁnite
volume i, S is the average sink term for the control-volume during the time step, Dl
and Dr are corresponding dispersion coefﬁcients, and ðvC=vxÞ*l and ðvC=vxÞ*r areassociated concentration gradients. The estimation of the face and gradient values isdiscussed below.
2.1. The QUICK method
The QUICK method uses a quadratic upstream interpolation to obtain face valuesof concentrations as shown in Fig. 2. The method can be interpreted as a linear
interpolation, given by the ﬁrst term on the right hand side of Eq. (4), corrected bya term proportional to the upstream-weighted curvature, given by the second termon the right hand side.
For ﬂow to the right and for uniform grids the face values are estimated as:
11
ðC þ CC Þ  ðCFL þ CC  2CL Þ
2 L
8
1
1
Cr* ¼ ðCC þ CR Þ  ðCL þ CR  2CC Þ
2
8
Cl* ¼

2.2. The QUICKEST method(4)and for non-uniform grids we obtain:!
Dx2 1 CC  CL CL  CFL1
ðCL þ CC Þ  l

Dxl
Dxfl
2
8 DxL

2 1 C  C
D
1
x
C
 CL
R
C
Cr* ¼ ðCC þ CR Þ  r
 C
Dxr
Dxl
2
8 DxC

Cl* ¼

(5)

The general formula for the QUICK scheme is shown in Eq. (6). Note that the
QUICK scheme uses the concentrations at the beginning of the time step to calculate
the face values.

Dx2
Dx2 1 11
ðCL þ CC Þ  l Curvl ¼ ðCL þ CC Þ  lGradl  Gradfl
2
2
8
8 DxL


Dx2r
Dx2r 1
1
1*
Curvr ¼ ðCC þ CR Þ ðGradr  Gradl ÞCr ¼ ðCC þ CR Þ 
2
2
8
8 DxC

Fig. 2. Quadratic upstream interpolation for face values of concentrations Cl (a) and Cr(b), and deﬁnition of face concentration gradients (vC/vx)l and (vC/vx)r used by the
QUICK method (Leonard, 1979).

Cl* ¼

(6)

where Curv represents (upstream-weighted) concentration curvatures and Grad
concentration gradients (compare Eq. (5) and (6) for their deﬁnitions). Similar
equations can be derived for ﬂow in the opposite direction (to the left). For the
QUICK scheme, the concentration gradients in Eq. (3) are identical to the ones usedin Eq. (6). As discussed by Leonard (1979), the stability criteria for QUICK are more
restrictive than in the case of central differencing, but in theory QUICK can be used
for cases when Pe > 2. However, as will be shown in the next section, the QUICK
method is still prone to oscillations.

Rather than using the quadratic upstream interpolation only at the old-timelevel to obtain face concentration values as in the QUICK scheme, the QUICKEST
method attempts to obtain these values as an average over a time increment Dtbetween the old and new time levels. The average face concentrations are estimatedassuming that the concentration proﬁle approximated using the quadratic upstream
interpolation at the old-time level is “swept” downstream unchanged, as shown in
Fig. 3. The QUICKEST formulation uses the same approach of averaging over a time
step for the physical diffusion terms and the effect of diffusion on wall values.
This leads to the following deﬁnition of the right face concentration andconcentration gradient values, respectively:


Dxr
Dx2
1
1
ðC þ CR Þ 
Cor Gradr þ r cr  1  Co2r Curvr
2 C
3
2
2

Cr* ¼

 *
Dxr
vC
Cor Curvr¼ Gradr 
vx r
2

(7)

(8)

where Grad and Curv are deﬁned as above. Similar equations can be derived for theleft face values and for ﬂow in the opposite direction (to the left).
Eq. (7) and (8) contain two dimensionless terms, the Courant number and the
diffusion parameter c, deﬁned respectively as:
ur Dt
Dx

(9)

Dr Dt
Dx2r

(10)

Cor ¼

cr ¼

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

Fig. 3. Control-volume face values Cl and Cr are estimated in the QUICKEST method as
an average of C(j) and C(jþ1). C(tþDt) is obtained by simply translating C(t) to the rightby a distance uDt, assuming pure (retarded, when adsorption is considered)convection.
Leonard (1979) demonstrated that the QUICKEST method is stable in the (Co,c)plane if 0  Co  1 and 0 < c  0.5. The stability region is extended to parts of the(Co,c) plane if Co > 1 and c < 0.5, and to c > 0.5 if Co < 1. The method is also stable
for the region where pure convective ﬂow dominates (c ¼ 0) if the Courant number
does not exceed unity.
The adoption of the estimated upstream terms using the new time level in the
QUICKEST scheme greatly reduces the amount of over and under shooting andnumerical dispersion. However, as demonstrated in the next sections and also by
Leonard (1991) the QUICKEST method is still prone to small oscillations and/or
overshooting, especially near sharp gradients. The presence of such oscillations canlead to negative concentration values. The presence of these oscillations can beeliminated by using a scheme that preserves local monotonic resolution.
2.3. The ULTIMATE method
The ULTIMATE method adopted here aims to suppress oscillation by using a TVD
scheme to preserve local monotonicity of the solution. The ULTIMATE scheme isderived using the Normalized Variation Diagram (NVD) concept, based on
normalized variables deﬁned as:
CCj  CUj

CCnorm ¼

j

(11)

j

CD  CUCr*  CUj

Cfnorm ¼(12)

CDj  CUjjjwhere CU and CD are the concentration upstream and downstream of the central
node “C” at time step “j”. The ULTIMATE scheme uses the normalized variables fromthe NVD concept to set the boundaries, between which each nodal value mustremain to suppress oscillations. Fig. 4 shows the necessary conditions imposed on
CCnorm used in the ULTIMATE scheme to preserve monotonicity, based on the fact that
CDnorm ¼ 1 and CUnorm ¼ 0.
In order to preserve monotonicity, the following limits for the maximum
normalized face values are:CCnorm  Cfnorm  1
Cfnorm 

CCnormCo

If Ccnorm < 0

for

for

0  CCnorm  1

0 < CCnorm  1
Ccnorm > 1, it is simple to just set Cfnorm

(13)

(14)
¼ CCnorm , which is theorequivalent of setting the wall face value to the node value. The adoption of the
ULTIMATE strategy involves minimal extra computational load and as will be shownin the following section, it avoids oscillations that could lead to appearance ofnegative concentrations near steep fronts.
2.4. Implicit schemes used in HYDRUS-1D
The governing ﬂow and transport equations are solved in HYDRUS-1D numerically using Galerkin-type linear ﬁnite element schemes. Mass lumping is invoked by

1301

Fig. 4. Normalized node values in the case of locally monotonic behavior. Hatchingshows necessary conditions on the face value of interest Cfnorm , and on the corresponding upstream control-volume face value Cunorm (Leonard, 1991).

redeﬁning the nodal values of the time derivative as weighted averages over theentire ﬂow region. The Galerkin method is used only for approximating the spatial
derivatives while the time derivatives are discretized by means of ﬁnite differences.
Different ﬁnite difference schemes can be selected by users depending upon thevalue of the time-weighting coefﬁcient e (¼0: explicit scheme, ¼ 0.5:
CrankeNicholson scheme, ¼ 1: fully implicit scheme). Higher-order approximations
for the time derivative in the transport equation as derived by van Genuchten andGray (1978) are implemented.
Upstream weighing is provided as an option in HYDRUS-1D to minimize some ofthe problems with numerical oscillations when relatively steep concentration fronts
are being simulated. For this purpose the ﬂux term of Eq. (1) is not weighted using
regular linear basis functions, but instead with the special non-linear functions that
ensure that relatively more weight is placed on the ﬂow velocities of nodes locatedat the upstream side of an element. The extent of increased weight on upstream
nodes is evaluated using the equation of Christie et al. (2005); a larger weight isplaced on upstream nodes for larger Peclet numbers and a standard weight is placedon them for small Peclet numbers.
The acronym CENTRAL is used below in the Results Section for the
CrankeNicholson implicit scheme with central spatial differencing. The acronym
UPSTREAM is used below for the CrankeNicholson implicit scheme with upstream
spatial weighting. The results obtained using the quadratic explicit schemes (i.e.,
QUICK, QUICKEST, and ULTIMATE) are compared only with the results obtainedusing the CrankeNicholson implicit schemes, since these schemes are consideredeither less prone to numerical instabilities or superior in view of solution preci sion than explicit and fully explicit schemes, respectively (Simunek et al., 2008a).
Note that the explicit scheme, although in general available in HYDRUS-1D, isdisabled in its graphical user interface, in order to discourage HYDRUS-1D users
from its use.

3. Results
Implementation into HYDRUS-1D of various quadratic upstream
interpolation schemes, i.e., QUICK, QUICKEST, and ULTIMATE, andtheir performance is veriﬁed in this section on multiple problems ofincreasing complexity. The ﬁrst set of problems involves challenging purely convective transport with uniform grids andvelocities, with the results compared against analytical solutions.
The second set of problems introduces dispersion and nonuniformity of grids, with results compared against the results
obtained using the standard implicit ﬁnite element numerical
schemes of HYDRUS-1D that have been veriﬁed many times before (e.g. Scanlon et al., 2002; Sim
unek et al., 2008a; Vanderborght et al.,2005). In particular, the central spatial weighting and the
CrankeNicholson temporal differencing are used for this purpose.

1302

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

The ﬁnal set of more complex problems involves cases withtransient water ﬂow, non-uniform velocities and reactions. All
results discussed below, other than those in which negative
concentrations were set to zero, were mass conservative andproduced mass balance errors signiﬁcantly less than 1% for theoverall simulations over the entire domain.
Boundary conditions for the quadratic schemes are prescribedby either assigning a node concentration value (Dirichlet) or a cell
face concentration value as described in Tkalich (2006, 2007). Fromthe point of view of numerical stability, both Dirichlet and Cauchy
boundary conditions behaved similarly and did not present any
problems. All problems presented here used a Dirichlet boundary
condition, with the exception of the last problem where we haveused a Cauchy boundary condition by prescribing the ﬂux at theboundary and used a zero concentration gradient at the bottom.
3.1. Purely convective transport
As outlined in previous sections, standard numerical schemes,such as CrankeNicholson Finite Element (or Finite Difference)schemes with central spatial differencing are not suited to deal
with purely convective ﬂow, producing spurious oscillations. The
quadratic upstream interpolation schemes presented here were
designed to overcome partly or fully this oscillatory behavior. We
compare below these new numerical schemes using two test
examples involving (a) a unit step change at the boundary(Leonard, 1991), and (b) the transport of an isolated sine-squared
wave (Sweby, 1984). Results obtained using standard implicit
schemes are shown for comparison.
3.1.1. A unit step change at the boundary
A unit step change in concentration at the boundary at time zero
represents a basic test for most numerical schemes. The following
physical parameters are used in this example. The transportdomain is 1 m long and is discretized using 101 nodes (Dx ¼ 0:01 ).
Since dispersion, D, is assumed to be zero, the Peclet number(Pe ¼ vDx=D) is equal to inﬁnity. Mean pore water velocity, v, isassumed to be equal to 1 m/d. A constant time step (Dt) of 0.001 d isused and the ﬁnal time (tmax) is 0.5 d. The corresponding Courant
number (Co ¼ vDt=Dx) is equal to 0.1.
Fig. 5 compares results obtained using the implicit schemes (i.e.,
CENTRAL and UPSTREAM) with the explicit schemes (QUICK,

Fig. 5. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) against the analytical solution for a unit step change
example.

QUICKEST and ULTIMATE) and the analytical solution. The CENTRAL
scheme shows, as expected, clear oscillations behind the concentration step, with a maximum oscillation wave on the order of 0.2.
Oscillations, propagating all the way toward the inﬂow boundary,would keep increasing for longer simulation times, eventuallycorrupting the entire solution. On the other hand, the UPSTREAM
scheme shows, again as expected, a very gradual change ofconcentrations around the front, rather than the abrupt change as
predicted by the analytical solution. While this scheme is monotonic, i.e., it does not display any oscillations, it clearly introduces
large artiﬁcial numerical dispersion, which is an inherent propertyof this numerical scheme.
The QUICK scheme shows the largest oscillations both in front ofand behind the traveling concentration step (the largest oscillation
is about 0.4). Compared to the CENTRAL scheme, oscillations
quickly dissipate toward the inﬂow boundary. The sharp change inconcentrations is better approximated than by the CENTRAL
scheme. Leonard (1979) suggests that oscillations in the QUICK
method are less likely to corrupt the solution.
The QUICKEST scheme signiﬁcantly reduces oscillations toward
the inﬂow (there is actually only a single overshoot wave of about0.05) compared to the QUICK scheme, while preserving a good
approximation of the step increase due to small numericaldispersion. Finally, the ULTIMATE scheme entirely removes theoscillations both in front of and behind the concentration step,
while still preserving a good description of the concentration step.
Since all these results correspond closely with results of Leonard(1991), we can conclude that all three new explicit numerical
schemes were likely implemented into HYDRUS-1D correctly.
3.1.2. An isolated sine-squared wave
While in the previous test example, the concentration change
was imposed using the boundary condition, in this example
concentrations are speciﬁed using the initial condition. This test
follows Sweby (1984) by specifying an isolated sine-squared wave
of width 20Dx as an initial condition:

Cðx; t ¼ 0Þ ¼ sin2
¼ 0

 px 
20Dx

for0  x  20Dx

(15)otherwiseThis function represents a relatively smooth proﬁle witha continuously changing gradient on both sides of a single local
maximum (Leonard, 1991). All other physical conditions, i.e., spatial
and temporal discretizations and velocity, are the same as for theunit step change example.
Fig. 6 compares results obtained using the implicit schemeswith the explicit schemes and the analytical solution. Two solutions
are shown with the CENTRAL scheme, one where negative
concentrations are not allowed (by resetting them to zero; denoted
CENTRAL1), and a second denoted CENTRAL2 without any resetting. Apart from minor oscillations behind the traveling wave, thestandard CENTRAL scheme shows a very good approximation of thetraveling sine-squared wave. When negative concentrations are notallowed (by resetting them to zero), the CENTRAL1 scheme showsan almost perfect description of the traveling wave. The UPSTREAM
scheme shows, again as expected, an unacceptably large clipping(underestimation of maximum concentrations) of the concentration peak due to numerical dispersion, even though the scheme ismonotonic and free of any oscillations.
The worst results were obtained using the QUICK scheme thatquite dramatically overestimated the concentration peak, as well as
displayed signiﬁcant oscillations behind the traveling wave. Onceagain, the QUICKEST scheme improves on the QUICK results, withan almost perfect description of the traveling wave. It shows a small
amount of clipping as it slightly underestimated the concentration

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

1303

Fig. 6. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) against the analytical solution for an isolated sine-squared
wave and purely convective transport.

Fig. 7. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) for an isolated sine-squared wave and a convectivedispersive transport with Pe ¼ 2.

peak (Cmax ¼ 0.982). The scheme also shows small undershoot bothin front of and behind the traveling wave (<0.018). Minor oscillations produced by the QUICKEST scheme were completely removedby the ULTIMATE scheme, although at the expense of a small
increase in clipping (Cmax ¼ 0.96).
All schemes, other than the UPSTREAM scheme, described wella continuously changing gradient on both sides of a single local
maximum of the sine-squared traveling wave.

made gradually ﬁner from left to right using a multiplication factorof 1/1.01. In the second case, the smallest element was on the left
side of the grid, which was made gradually coarser from left to rightusing a multiplication factor of 1.01. In both cases the smallest
element was 0.00181818 m and the largest 0.0181818 m. In bothcases an isolated sine-squared wave was speciﬁed as an initialcondition over the same number (i.e., 20) of grid nodes. Consequently, different solute mass was initially in the transport domain
in these two cases since the initial sine-squared wave was spreadover different parts of the domain (grids have different sizes).
As shown in Fig. 8, all numerical schemes, except for UPSTREAM
which again produced signiﬁcant numerical dispersion, performed
reasonably well when the grid becomes ﬁner in the direction of
ﬂow. However, none of the numerical schemes performed wellwhen the grid becomes increasingly coarser in the direction of ﬂow,
as shown in Fig. 9. In this case both oscillatory schemes, i.e.,
CENTRAL and QUICK, suffered from large (larger than for theuniform grid) oscillations behind the traveling wave, and as in theuniform grid case, oscillations of the CENTRAL scheme propagated

3.2. Convective-dispersive transport
3.2.1. Uniform grids
As described earlier, numerical oscillations in implicit numerical
schemes (e.g. CENTRAL) can be virtually eliminated when local
Peclet numbers do not exceed about 2. Although acceptably small
oscillations may be obtained with local Peclet numbers as high as
10 (Huyakorn and Pinder, 1986) unacceptably large oscillations are
obtained when local Peclet numbers are larger than 10 (as shown
above in Fig. 5). As such, for cases where Pe equals 2 all numerical
schemes should converge to the true solution.
The previous example with a traveling isolated sine-squared
wave was modiﬁed by introducing a diffusion coefﬁcient D equal
to 0.005 m2/d and keeping all other parameters the same, so the
local Peclet number is equal to 2 (P ¼ vDx/D ¼ 1  0.01/0.005 ¼ 2).
The result is shown in Fig. 7 and, in this particular example, resultsof all numerical schemes, other than QUICK and UPSTREAM,
converged to the same solution. While the QUICK scheme slightly
overestimated the concentration peak (by about 0.02), the
UPSTREAM scheme underestimated the concentration peak moresigniﬁcantly (by about 0.05). Even for the local Peclet number of 2,
the UPSTREAM scheme introduced some, although limited,
numerical dispersion.
3.2.2. Non-uniform grids
All numerical schemes were further tested using spatially nonuniform grids, with the grid spacing either gradually increasing or
decreasing. All other physical parameters other than the initialcondition and spatial discretization are identical to the example fora purely convective transport shown in Fig. 6. In both cases thetransport domain was discretized using 101 nodes; in the ﬁrst case,the largest element was on the left side of the grid, and the grid wasFig. 8. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) against the analytical solution for an isolated sine-squarewave, a purely convective transport, and a gradually ﬁner grid.

1304

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

Fig. 9. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) against the analytical solution for an isolated sine-squared
wave, a purely convective transport, and a gradually coarser grid.

much further upstream than for the QUICK scheme. All nonoscillatory schemes produce signiﬁcant clipping of the peakconcentration and the QUICKEST scheme again produced onlya single undershoot at both sides of the traveling wave, with noundershoot in the ULTIMATE case. Once more, the UPSTREAM
scheme produced much higher levels of numerical dispersion,
although in this case, all non-oscillatory regimes are affected bynumerical dispersion.
Although this is clearly a disturbing ﬁnding, it is important toemphasize that these results are for purely convective transport, forwhich the standard numerical schemes, such as CENTRAL and
UPSTREAM, are not expected to be either stable or accurate. In mostsoil applications (if not all), such conditions are not encountered as
the transport is likely to also always contain dispersion.
It is thus necessary to investigate, for what conditions with
gradually coarser grids these quadratic upstream weighting
schemes (mainly QUICKEST and ULTIMATE) provide reasonable
solutions. Fig. 10 presents the same problem with additional
dispersion so that the three scenarios shown have Peclet numbers
2, 20 and 200. To provide a benchmark, these three additional runs
are compared to the CENTRAL scheme (denoted below as
CENTRAL#) applied on a much ﬁner uniform grid, for which thestability conditions were satisﬁed, i.e., Pe < 2.
As expected all numerical schemes studied, i.e., CENTRAL,QUICKEST and ULTIMATE converged to a true solution provided by
CENTRAL# when the Peclet number was equal to 2. Relatively goodresults were still obtained for a Peclet number equal to 20 by the
QUICKEST and ULTIMATE schemes, which only slightly underestimated the concentration peak. The CENTRAL scheme alreadyshows relatively minor oscillations beyond the sine-squared wave.
Finally, for a Peclet number of 200, all three schemes produced
signiﬁcant clipping, similarly as for purely convective transport,
while the CENTRAL scheme additionally suffered from signiﬁcantoscillations and the QUICKEST scheme shows undershooting before
and after the sine wave. Signiﬁcant oscillations encountered with
the CENTRAL scheme caused slight lagging of the sine wave front
behind the analytical solutions.
Although for increasingly coarser non-uniform grids the
QUICKEST and ULTIMATE schemes cannot be used for inﬁnite Peclet
numbers (i.e., purely convective ﬂow), they nevertheless quite
dramatically expand the interval of stability (Pe ¼ 20, i.e., at least
ten times) compared to standard implicit schemes.

Fig. 10. Comparison of various numerical schemes (CENTRAL, QUICKEST, and ULTIMATE) for a convective-dispersive transport characterized by a Peclet number of 2 (a),20 (b), and 200 (c) on a gradually coarser grid. CENTRAL# represents the CENTRAL
scheme on a very ﬁne grid, representing the true solution.

3.3. Non-uniform velocities and adsorption
Non-uniform velocities cannot be generated for a onedimensional problem without compensating at the same time forchanges in velocities by changes in water contents, which would
lead to a transient water ﬂow problem (addressed below).
However, non-uniform “retarded” velocities can be generated byconsidering transport of a solute undergoing non-linear sorption
while keeping water contents constant (unity so far). In this casethe convectionedispersion equation

R

vC
vC
v2 C
¼ u
þD 2
vt
vx
vx

(16)

can be reformulated using the so-called retarded velocities (vR) andretarded dispersion coefﬁcients (DR) (Jury and Horton, 2004) as
follow:vC
u vC D v2 C
¼ þvt
R vx R vx2
vC
vC
v2 C
¼ uR þ DR 2
vt
vx
vx

(17)

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e13081305

where for non-linear sorption the retardation factor R is deﬁned as:

R ¼ 1þ

rds
qdC

(18)

and for the Freundlich adsorption isotherm:s ¼ Kd C n
R ¼ 1þ

rnKd C n1q(19)

where s is the sorbed concentration [MM1], r is the soil bulk
density [ML3], q is the water content [-], and Kd [M1  nL3n M1]and n [-] are empirical Freundlich coefﬁcients. Note that Eq. (17)with retarded velocity and dispersion has the same form as Eq.
(1), for which the quadratic numerical schemes discussed in thispaper were developed, and they can thus be used for problems
involving retardation, without any major modiﬁcations.
However, for the QUICKEST (and consequently ULTIMATE)
scheme, the method estimates the wall values as an average over
the time interval, so we solve a modiﬁed form of Eq. (16) with R
inside of the temporal derivative, and a different deﬁnition of R as
given below (Huang et al., 1997):

vRC
vC
v2 C
¼ u
þD 2
vt
vx
vx
R ¼ 1þ

rKd C n1
q

(20)

(21)

Note that the Freundlich exponent n appears as multiplicator in
the deﬁnition of the retardation factor in Eq. (19) and not in Eq. (21).
This is because the retardation factor deﬁned by Eq. (19) is outsideof the temporal partial derivative in Eq. (17) and thus the sorbed
concentration was subject to the ordinary derivative given in Eq.
(18), while the retardation factor deﬁned in Eq. (21) is inside of thetemporal partial derivative in Eq. (20). Since the explicit approach
was used for evaluating the retardation factor, i.e., the retardation
factor was evaluated using concentrations at the old-time level, no
iterative procedure was required for non-linear problems.
Thus, when there is a traveling concentration wave in the transport domain, although actual water velocities may be uniform,
retarded velocities will be non-uniform, because of spatially changing
concentrations and consequently changing retardation factors.
The problem presented here, involving steady-state water ﬂow
and non-linear solute transport in a homogeneous soil proﬁle,follows a benchmark traveling wave problem suggested byVanderborght et al. (2005) and an analytical solution by van der Zee(1990). The physical parameters for this problem were as follows.
The soil column was 200 cm deep and was discretized using 101
nodes. The water content and velocities were uniform and equal to0.004 and 0.05 m/d, respectively. The soil proﬁle was initially solute
free and the inﬂow concentration was 10 mg mL1. The soil bulk
density was 1000 mg mL1, the Freundlich exponent n was 2/3, andthe Freundlich sorption coefﬁcient Kd was 0.001 mLn mg1 mg(1  n).
The dispersivity l was equal to 0.01 m, leading to a local Peclet
number (Pe ¼ vDx/D ¼ vDx/lv ¼ Dx/l) of 2. A constant time step (Dt)of 0.04 d was used and the ﬁnal time (tmax) was 60 d. The corresponding maximum local Courant number (variable in spacebecause of non-uniform vR, Co ¼ vRDt/Dx) was about 0.05.
Notice that due to the relatively low Pe number all numerical
schemes, except for the QUICK scheme, provided as expected almost
identical solutions similar to the one obtained by the analyticalsolution (Fig. 11). The UPSTREAM schemes shows again limited
numerical dispersion because of a relatively low local Peclet number.
Similar veriﬁcation for transport with the Freundlich exponent
larger than 1 (n ¼ 4/3) shows similar agreement between theFig. 11. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) for a traveling wave problem.

models. However, since concentration fronts for n < 1 are steeper,and thus more difﬁcult to solve numerically, only results for thiscase (i.e., n ¼ 2/3) are presented here.
3.4. Convective-dispersive-reactive transport
The next problem involves, in addition to physical processes ofdispersion and retardation, also reactive chemical processes ofdegradation and production. The problem involves ﬁrst-orderdecay and zero-order production. HYDRUS-1D has been shown in 
multiple studies (e.g. Sim
unek et al., 2008a; Vanderborght et al.,2005) to correctly solve such problems, and thus in the following
example it is considered to provide a true solution that otherschemes will be compared against.
The following physical parameters are used in this example. The
transport domain is again 1 m long and is discretized using 101
nodes (Dx ¼ 0.01 m). Mean pore water velocity, v, is assumed to beequal to 1 m/d. Since dispersion, D, is assumed to be 0.001 m2 d1,
the Peclet number (Pe ¼ vDx/D) is equal to 10. The linear distribution coefﬁcient Kd is equal to 1 g1 cm3. A constant time step (Dt)of 0.005 d is used and the ﬁnal time (tmax) is 1.0 d. The corresponding Courant number (Co ¼ vRDt/Dx) is equal to 0.25. The
solute is additionally subject to the ﬁrst-order decay with a rate
constant m equal to 1 d1 and the zero-order production with a rate
constant g equal to 0.0001 g m3 d1. The reaction term of Eq. (1) isthus deﬁned as:

S ¼ mc þ g(22)Comparison of various numerical schemes for this reactive
transport problem is shown in Fig. 12. Notice that the CENTRAL,QUICKEST and ULTIMATE schemes produce almost identical results,
similarly as in the two previous examples. The QUICK and
UPSTREAM schemes, similarly as above in Figs. 6 and 7, overestimate and underestimate the concentration peak, respectively.
When the dispersion D was increased to 0.005 m2 d1, i.e., thePeclet number was reduced to 2, the correspondence between
results obtained by individual numerical schemes (not shown) was
similar as in Fig. 7.
3.5. Transient water ﬂow
An ultimate test for each numerical scheme for solute transport
is how it performs under transient water ﬂow conditions. So in the1306

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

Fig. 12. Comparison of various numerical schemes (CENTRAL, UPSTREAM, QUICK,QUICKEST, and ULTIMATE) for a transport of an isolated sine-squared wave subject toretardation, production, convection, and dispersion.

last example transient water ﬂow conditions are generated by
imposing ponding (a pressure head of zero) at the surface of a 1-m
deep loamy soil proﬁle.
The following physical parameters are used in this example. The
transport domain is discretized using 101 nodes (Dx ¼ 0.01 m).
Diffusion is not considered, but the dispersivity l is consideredusing two values, either 0.005 or 0.0005 m. Even though the local
velocities are not constant, the local Peclet numbers are constantthroughout the proﬁle when a uniform grid is used (since Pe ¼ vDx/
vl ¼ Dx/l) and equal to 2 (l ¼ 0.005 m) and 20 (l ¼ 0.0005 m). The
former case is used to demonstrate that various numerical schemes
converge to the same solutions when Pe is about 2, while the lattercase shows the improvement of the ULTIMATE scheme for problems with high Peclet numbers.
To obtain transient values of water contents and pore watervelocities required by the convectionedispersion equation fortransient water ﬂow, the Richards equation for variably-saturated
water ﬂow in porous media was solved using standard numerical
techniques of HYDRUS-1D. Soil hydraulic parameters for loam fromthe Carsel and Parrish (1988) database were used. For this soil, thesaturated hydraulic conductivity Ks is equal to 0.249 m d1, theresidual and saturated water contents are 0.078 and 0.43, respectively, and the van Genuchten shape parameters a and n (vanGenuchten, 1980) are 3.6 m1 and 1.56, respectively. Initial pressure head was equal to 100 cm and the free drainage boundary
condition was used at the bottom of the soil column.
The time step (Dt) was determined by requirements of thenumerical solution of the Richards equation and varied from theinitial value of 0.0001 d to a maximum value of 0.008 d. The ﬁnal
time (tmax) was 1.0 d. The corresponding Courant numbers(Co ¼ vDt/Dx) varied between 0.1 and 0.5. The solute was assumed
not to be subject to any chemical reactions. The governing convectionedispersion equation for transient ﬂow equation solved inthis case is:



vqc
vuqc
v
vc
qD
¼ 
þ
vt
vx
vx
vx

Fig. 13 shows generated water content proﬁles after 0.2, 0.4, 0.6,
and 1.0 d. The water content front reaches the depth of about
0.40 m after 0.2 d and the bottom of the soil proﬁle after about
0.65 d. Fig. 13 also shows concentration proﬁles after 0.2, 0.4, 0.6,
0.8, and 1.0 d obtained using the CENTRAL, QUICKEST, and ULTIMATE schemes. All three numerical schemes converged to the same
solution since the local Peclet numbers were equal to 2, i.e., nooscillations were expected to develop.
Fig. 14 shows concentration proﬁles for the same transient ﬂow
conditions, however for transport conditions with the Peclet
number equal to 20. As discussed above, this is signiﬁcantly above
the range of Peclet numbers for which the solution is expected to be
oscillation free for the CENTRAL scheme. And indeed, signiﬁcantoscillations developed behind the concentration front. In this casethe concentrations were constrained to positive values, similarly asin the example shown in Fig. 6 above. If this constraint was not
used, small negative oscillations in front of the concentration frontand signiﬁcantly larger oscillations behind the front would havedeveloped (data not shown). The QUICKEST scheme showed only
a minor undershoot and overshoot in front of and immediately
behind the concentration front, respectively. The ULTIMATE
scheme completely eliminated negative values ahead of theconcentration front and dramatically reduced the overshoot behindthe front.
4. Discussion
Test examples shown in the previous section demonstrate thatthe quadratic interpolation schemes QUICKEST and ULTIMATE,(23)

Similarly as above for the non-linear problem, the partially implicit
approach in the QUICKEST and ULTIMATE methods was applied
only to concentrations, and water contents, ﬂuxes, and dispersion
coefﬁcients were evaluated at the old-time level.

Fig. 13. Water contents (top) and concentrations (bottom) for the transient ﬂow
example and dispersivity of 0.5 cm (Pe ¼ 2). Results obtained for t ¼ 1 day using
CENTRAL (solid lines), QUICKEST (B), and ULTIMATE (D) schemes at are also shown.

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

Fig. 14. Concentrations for the transient ﬂow example and dispersivity of 0.05 cm(Pe ¼ 20). Results obtained using CENTRAL (solid lines), QUICKEST (B), and ULTIMATE(D) schemes.

implemented in HYDRUS-1D, were able to adequately solve various
presented transport problems. It is also important to emphasizethat the examples, in which the standard implicit schemes, i.e.,
CENTRAL and UPSTREAM, failed or did not produce reasonable
results, were for purely convective transport or very large local
Peclet numbers. Such transport is only encountered in few applications and is very unlikely to occur in vadose zone hydrology. As
shown above, numerical oscillations in implicit numerical schemes(e.g. CENTRAL) can be virtually eliminated when local Peclet
numbers do not exceed about 2. Furthermore with local Peclet
numbers as high as 10 (Huyakorn and Pinder, 1986) acceptably
small oscillations can still be obtained. Even for problems withlarger Peclet numbers, the oscillations can be minimized or virtually eliminated when the time step is adjusted appropriately.
Perrochet and Bérod (1993) proposed a stability criterion:

Pe$Cr  us ð ¼ 2Þ(24)

where us is the performance index [-]. This criterion indicates thatconvection-dominated transport problems having large Peclet
numbers can be safely simulated provided Co (by, in fact, reducing
the time step) is reduced according to the equation above(Perrochet and Bérod, 1993). Imposing these constraints, i.e., either
Pe<2 or us < 2, numerical oscillations can be virtually eliminatedand the CENTRAL scheme is expected to provide precise solutions (Simunek et al., 2008b). When small oscillations in the solution canbe tolerated, us can be increased to about 5 or 10.
For real-life one-dimensional problems encountered in vadosezone hydrology, it is usually not a problem to adjust the spatial
discretization grid so that these criteria are fulﬁlled, and usually thenumerical solution of the Richards equation requires ﬁner temporal
discretization than that required by solute transport. If needed,having hundreds or even several hundreds of nodes, is usually not
a problem either in terms of memory or solution time for most onedimensional applications. It is worthwhile noting that, as shown in
Figs. 6, 10 and 14, the QUICKEST and ULTIMATE schemes quite
dramatically expand the stability range, compared to standard
implicit schemes, to Pe w 20 for non-uniform grids and up to
inﬁnity for uniform grids.
There are, however, several other areas where these nonoscillatory techniques can be used with a great advantage. One
such area is in multi-dimensional applications. Larger twodimensional, and certainly three-dimensional, applications in the1307

vadose hydrology can easily have several hundreds of thousands, oreven millions, of elements (volumes), as required by various
stability criteria. Reducing the number of elements for theseapplications by relaxing the constraints on maximum grid spacing
can still produce signiﬁcant computational savings. Both QUICKEST
and ULTIMATE schemes have been used in many multidimensional codes of hydrodynamics (e.g. Cole and Wells, 2006;
Lin and Falconer, 1997; Romero et al., 2004) and it is our goal,
after this initial testing presented in this manuscript, to also
implement them in the multi-dimensional versions of HYDRUS (Simunek et al., 2008b).
The second area where these numerical schemes can be usedwith a great advantage is in codes simulating multicomponent
reactive chemistry (e.g. Jacques et al., 2008; Wissmeier and Barry,
2010, 2011). Such codes are usually developed by coupling the
ﬂow and transport codes with the reactive chemistry codes, such as
PHREEQC (Parkhurst and Appelo, 1999), and simulations with thesecodes often fail when they are provided with negative element
concentrations by the transport codes. As shown in Figs. 13 and 14,
the ULTIMATE scheme does not produce oscillations andsuppresses negative concentrations. The ULTIMATE scheme
preserves sharp concentration fronts while maintaining monotonicity of the solution, which is especially relevant in ﬂow situations with high grid Peclet numbers.
As shown above, the accuracy of the ULTIMATE scheme greatly
exceeds the accuracy of commonly applied ﬁrst- and second-orderschemes currently used in HYDRUS for problems, in whichconvective transport dominates over the dispersive transport, andother problems with sharp concentration fronts. Note that sharp
concentration fronts may not necessarily be a result of a dominantconvective transport, but can also be result of involved chemical
reactions, such a sorption (Borkovec et al., 1996; Frey, 1987), cationexchange (Voegelin et al., 2000), or precipitation/dissolution (Lakeet al., 2002). All these reactions can under certain conditions (e.g.
for n < 1 in the Freundlich adsorption isotherm) lead to selfsharpening concentration fronts, which are difﬁcult to handle bytraditional implicit schemes.
5. Conclusions
Two quadratic interpolations schemes (QUICK and QUICKEST)
and one TVD scheme (ULTIMATE) have been successfully implemented in HYDRUS-1D and tested against a series of analytical andnumerical solutions. Both quadratic schemes showed small
numerical dispersion. The QUICKEST scheme showed reduced
oscillatory behavior when compared to the QUICK scheme or
unbound central differencing. The incorporation of the TVD scheme
removed the oscillatory behavior near sharp fronts for large Peclet
numbers.
The QUICKEST and ULTIMATE schemes presented here havea larger stability region in terms of the Peclet number, compared tostandard implicit schemes. Even if this advantage is of limited
application in 1D cases, it shows a large potential to be employed in
multi-dimensional applications to relax Peclet number constraints
on grid spacing.
Finally, the ULTIMATE scheme shows improved stability nearsteep fronts due to non-oscillatory behavior and reduced numericaldispersion. This suggests an improvement in stability in reactive
transport applications, as the absence of negative values is likely to
improve the convergence of such models.
Acknowledgments
The CSIRO Computational and Simulation Sciences Transformational platform is thanked for providing the funding for Prof.

1308

L.E. Neumann et al. / Environmental Modelling & Software 26 (2011) 1298e1308

 Jirka Simunek to travel to Australia and work on this project. Drs
Luis Neumann and Freeman Cook would like to thank the eWater
CRC for funding their time on this project.
References
Borkovec, M., Bürgisser, C.S., Cerník,
M., Glättli, U., Sticher, H., 1996. Quantitativedescription of multi-component reactive transport in porous media: anempirical approach. Transport Porous Med. 25 (2), 193e204.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil
water retention characteristics. Water Resour. Res. 24 (5), 755e769.
Carter, K., Saenz, M., Wilson, D., Rosten, P., 1984. Migration of pollutants ingroundwater. II. adsorbable pollutants and numerical dispersion reduction.
Environ. Monit. Assess. 4 (2), 171e202.
Christie, I., Grifﬁths, D., Mitchell, A., Zienkiewicz, O., 2005. Finite element methods
for second order differential equations with signiﬁcant ﬁrst derivatives. Int. J.
Numer. Methods Eng. 10 (6), 1389e1396.
Chu, X., Marino, M.A., 2007. Improved compartmental modeling and application to
three-phase contaminant transport in unsaturated porous media. J. Environ.
Eng. 132 (2), 211e219.
Cole, T.M., Wells, S.A., 2006. CE-QUAL-W2: A two-dimensional, laterally averaged,hydrodynamic and water quality model, Version 3.5, Instruction Report EL06e01. US Army Engineering and Research Development Center: Vicksburg, MS.
Frey, D.D., 1987. Free energy consumption at self-sharpening concentration fronts in
isothermal ﬁxed-bed adsorption: theory for isotherms of arbitrary shape andapplications in displacement chromatography. J. Chromatogr. A 409, 1e13.
Huang, K., Simunek, J., Genuchten, M.T.V., 1997. A third-order numerical schemewith upwind weighting for solving the solute transport equation. Int. J. Numer.
Methods Eng. 40 (9), 1623e1637.
Huyakorn, P.S., Pinder, G.F., 1986. Computational Methods in Subsurface Flow.
Academic Press, Inc, Orlando, FL.
 
Jacques, D., Sim
unek, J., Mallants, D., van Genuchten, M.T., 2008. Modeling coupledhydrologic and chemical processes: long-term uranium transport following
phosphorus fertilization. Vadose Zone J. 7 (2), 698e711.
Jury, W., Horton, R., 2004. Soil Physics, sixth ed. John Wiley and Sons, Inc., Hoboken,
New Jersey.
Lake, L.W., Bryant, S.L., Araque-Martinez, A.N., 2002. Geochemistry and Fluid Flow.
Elsevier, Amsterdam; London.
Leonard, B.P., 1979. A stable and accurate convective modelling procedure based on
quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng. 19,59e98.
Leonard, B.P., 1991. The ULTIMATE conservative difference scheme applied tounsteady one-dimensional advection. Comput. Methods Appl. Mech. Eng. 88(1), 17e74.
Lin, B., Falconer, R.A., 1997. Tidal ﬂow and transport modeling using ULTIMATEQUICKEST scheme. J. Hydraul Eng-ASCE 123 (4), 303e314.
Lin, Y., Medina, M., 2003. Incorporating transient storage in conjunctive streameaquifer modeling. Adv. Water Resour. 26 (9), 1001e1019.
Moldrup, P., Yamaguchi, T., Rolston, D.E., Vestergaard, K., Hansen, J.A., 1994.
Removing numerically induced dispersion from ﬁnite difference models forsolute and water transport in unsaturated soils. Soil Sci. 157 (3), 153e161.
Neumann, L.N., Cook, F., Western, A.W., Verburg, K., 2009. A one dimensional solutetransport model for hydrological response units. In: Anderssen, R.S.,
Braddock, R.D., Newham, L.T.H. (Eds.), 18th World IMACS Congress and MODSIM09 International Congress on Modelling and Simulation. Cairns, Australia,
Modelling and Simulation Society of Australia and New Zealand and International Association for Mathematics and Computers in Simulation, pp.
3528e3534.
Parkhurst, D.L., Appelo, C.A.J., 1999. User’s Guide to PHREEQC (Version 2):
a Computer Program for Speciation, Batch-reaction, One-dimensional Trans-

port, and Inverse Geochemical Calculations. U.S. Geological Survey, Denver,
Colorado.
Perrochet, P., Bérod, D., 1993. Stability of the standard CrankeNicolsoneGalerkinscheme applied to the diffusion-convection equation: some new insights.
Water Resour. Res. 29 (9), 3291e3297.
Roe, P., 1985. Large Scale Computations in Fluid Mechanics, Part 2. Springer, Berlin.
Romero, J., Antenucci, J., Imberger, J., 2004. One-and three-dimensional biogeochemical simulations of two differing reservoirs. Ecol. Modell 174 (1e2),
143e160.
 
Scanlon, B.R., Christman, M., Reedy, R.C., Porro, I., Sim
unek, J., Flerchinger, G.N.,
2002. Intercode comparisons for simulating water balance of surﬁcial sediments in semiarid regions. Water Resour. Res. 38 (12), 1323.
 Simunek, J., 2005. Models of water ﬂow and solute transport in the unsaturated
zone. In: Anderson, M.G., McDonnell, J.J. (Eds.), Encyclopedia of Hydrological
Sciences. John Wiley & Sons, Ltd., Chichester, England, pp. 1171e1180.
 
Sim
unek, J., van Genuchten, M.T., 2006. Contaminant transport in the unsaturated
zone theory and modeling. In: Delleur, J. (Ed.), The Handbook of Groundwater
Engineering, second ed. CRC Press, pp. 22.21e22.46.
 

Sim
unek, J., Sejna,
M.H., Saito, H.M., Sakai, M., van Genuchten, M.T., 2008a. The
HYDRUS-1D Software Package for Simulating the Movement of Water, Heat,and Multiple Solutes in Variably Saturated Media, Version 4.0. HYDRUS Software Series 3. Department of Environmental Sciences, University of California
Riverside, Riverside, California, USA, p. 315.
 
Simunek, J., van Genuchten, M.T., Sejna, M., 2008b. Development and applicationsof the HYDRUS and STANMOD software packages and related codes. Vadose
Zone J. 7 (2), 587e600.
Sweby, P., 1984. High resolution schemes using ﬂux limiters for hyperbolic
conservation laws. SIAM J. Numer. Anal. 21 (5), 995e1011.
Tkalich, P., 2006. Derivation of high-order advection-diffusion schemes.
J. Hydroinformatics 8 (3), 149e164.
Tkalich, P., 2007. Discussion of derivation of high-order advection-diffusion
schemes by Pavel Tkalich, 2006. J. Hydroinf. 8 (3), 149e164. Author’s reply.
Journal of Hydroinformatics 9(2) 159e162.
van Genuchten, M., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892e898.
van Genuchten, M.T., Gray, W., 1978. Analysis of some dispersion correctednumerical schemes for solution of the transport equation. Int. J. Numer.
Methods Eng. 12 (3), 387e404.
Van Leer, B., 1979. Towards the ultimate conservative difference scheme. V.
a second-order sequel to Godunov’s method. J. Comput. Phys. 32 (1), 101e136.
van der Zee, S.E.A.T.M., 1990. Analytical traveling wave solutions for transport withnonlinear and nonequilibrium adsorption. Water Resour. Res. 26 (10),2563e2578.
Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiery, D., Vanclooster, M.,
Mouvet, C., Vereecken, H., 2005. A set of analytical benchmarks to test
numerical models of ﬂow and transport in soils. Vadose Zone J. 4 (1), 206e221.
Voegelin, A., Vulava, V.M., Kuhnen, F., Kretzschmar, R., 2000. Multicomponent
transport of major cations predicted from binary adsorption experiments.
J. Contam. Hydrol 46 (3e4), 319e338.
Wissmeier, L., Barry, D., 2010. Implementation of variably saturated ﬂow into
PHREEQC for the simulation of biogeochemical reactions in the vadose zone.
Environ. Modell Softw. 25 (4), 526e538.
Wissmeier, L., Barry, D.A., 2011. Simulation tool for variably saturated ﬂow withcomprehensive geochemical reactions in two- and three-dimensional domains.
Environ. Modell Softw. 26 (2), 210e218.
Zheng, C., Bennett, G.D., 2002. Applied Contaminant Transport Modeling, second ed.
John Wiley & Sons, New York.
Zheng, C., Wang, P.P., 1999. MT3DMS: A Modular Three-Dimensional Multispecies
Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User’s
Guide. Contract Report SERDP-99e101. US Army Engineer Research andDevelopment Center, Vicksburg, MS.