Journal of Contaminant Hydrology 88 (2006) 197 – 218
www.elsevier.com/locate/jconhyd

Operator-splitting errors in coupled reactive transport
codes for transient variably saturated flow andcontaminant transport in layered soil profiles
D. Jacques a,⁎, J. Šimůnek b , D. Mallants a , M.Th. van Genuchten caWaste and Disposal Department, SCK•CEN, Boeretang 200, B-2400 Mol, Belgium
Department of Environmental Sciences, A135 Bourns Hall, University of California Riverside,
900 University Avenue, Riverside, CA 92521, USAc
George E. Brown, Jr. Salinity Laboratory, USDA, ARS, 450 W. Big Spring Rd, Riverside, CA 92507, USAb

Received 15 June 2005; received in revised form 20 June 2006; accepted 30 June 2006
Available online 17 August 2006

Abstract
One possible way of integrating subsurface flow and transport processes with (bio)geochemical reactions is tocouple by means of an operator-splitting approach two completely separate codes, one for variably-saturated flow
and solute transport and one for equilibrium and kinetic biogeochemical reactions. This paper evaluates theaccuracy of the operator-splitting approach for multicomponent systems for typical soil environmental problems
involving transient atmospheric boundary conditions (precipitation, evapotranspiration) and layered soil profiles.
The recently developed HP1 code was used to solve the coupled transport and chemical equations. For steady-stateflow conditions, the accuracy was found to be mainly a function of the adopted spatial discretization and to a lesser
extent of the temporal discretization. For transient flow situations, the accuracy depended in a complex manner on
grid discretization, time stepping and the main flow conditions (infiltration versus evaporation). Whereas a finer
grid size reduced the numerical errors during steady-state flow or the main infiltration periods, the errors sometimes
slightly increased (generally less than 50%) when a finer grid size was used during periods with a highevapotranspiration demand (leading to high pressure head gradients near the soil surface). This indicates that
operator-splitting errors are most significant during periods with high evaporative boundary conditions. The
operator-splitting errors could be decreased by constraining the time step using the performance index (the productof the grid Peclet and Courant numbers) during infiltration, or the maximum time step during evapotranspiration.
Several test problems were used to provide guidance for optimal spatial and temporal discretization.
© 2006 Elsevier B.V. All rights reserved.
Keywords: Coupled reactive transport; Unsaturated flow modelling; Geochemical reactions; Operator-splitting errors⁎ Corresponding author.
E-mail address: djacques@sckcen.be (D. Jacques).
0169-7722/$ - see front matter © 2006 Elsevier B.V. All rights reserved.
doi:10.1016/j.jconhyd.2006.06.008

198

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

1. Introduction
The migration of many contaminants in the subsurface is affected by a multitude of complex,
interactive physical, chemical, mineralogical, geological, and biological processes. Cycles ofprecipitation and evaporation largely determine if contaminants remain near the soil surface (e.g.,Öztürk and Özkan, 2004). Changes in the chemical composition or pH of the soil solution may
impact the retention of heavy metals, radionuclides, or other contaminants on organic matter or ironoxides (e.g., Kohler et al., 1996; Meeussen et al., 1998; Kent et al., 2000; Davis et al., 2004).
Dissolution and precipitation of minerals generally buffer the transport of a solution with a different
pH or chemical composition through the soil profile (e.g., Ayora et al., 1998; Adler, 2001). The
mobility of many elements such as arsenic and iron may be further impacted significantly by changes
in the redox state resulting from fluctuations in the phreatic surface, or because of temporary
saturation of certain soil layers during periods of heavy rainfall. Simulation of these and related
processes requires a coupled reactive transport code that integrates the physical processes of waterflow and advective–dispersive solute transport with a range of biogeochemical processes.
As discussed by Šimůnek and Valocchi (2002), coupled codes can be developed in several ways
such as using specific chemistry or a more general biogeochemistry module. Many existing codeswith general geochemistry are limited to solute transport and biogeochemical reactions, while waterflow paths have to be calculated outside of the reactive transport code (e.g., CRUNCH, Steefel,
2000; PHREEQC, Parkhurst and Appelo, 1999; PHT3D, Prommer et al., 2003). Other codes, such as
HYTEC (van der Lee et al., 2003), allow the velocity field to be calculated internally from thepermeability field. Still other models have coupled transport and biochemical reactions withsaturated flow (TBC, Schäfer et al., 1998; PHAST, Parkhurst et al., 2004) or variably saturated flow
(HYDRUS-2D-CW2D, Langergraber and Šimůnek, 2005). The more advanced codes integrate
general (bio)geochemistry models with transient unsaturated flow and transport models. These
include the multi-dimensional codes 3DHYDROGEOCHEM (Yeh and Cheng, 1999), CORE2D(Samper et al., 2000), MIN3P (Mayer et al., 2002), and RETRASO (Saaltink et al., 2004).
Two approaches exist to solve the coupled equations for transient flow, transport and geochemistry(e.g., see Steefel and MacQuarrie, 1996, for an overview): (1) the operator-splitting approach in which
the transport and chemistry equations are solved separately (an overview of operator-splitting schemes
is given in Carrayrou et al., 2004), and (2) the one-step global implicit approach in which the transportand chemistry equations are solved simultaneously. The main advantages of the global implicit
approach are that (a) all processes are treated simultaneously in the temporal, spatial and chemical
domains, (b) the method is mass-conservative, and (c) the approach may lead to better convergence forlarge time steps, which may be an important consideration for multi-dimensional flow problems(Mayer, 1999; van der Lee and De Windt, 2001). However, the global implicit method is generally
quite complex mathematically and its memory requirements can be large.
The main advantages of the operator-splitting approach are that: (a) the method is easy to
implement when different codes or modules are coupled (i.e., well-documented and tested codes canbe used, while the modules later still be updated), (b) the approach allows one to use different andspecific methods for solving the transport and reaction equations (Carrayrou et al., 2004), and (c) the
resulting code can be readily modified for parallel computations (Hundsdorfer and Verwer, 1995). A
disadvantage is the occurrence of operator-splitting errors resulting from both the operator-splitting
approach itself and from the time and spatial discretizations (Carrayrou et al., 2004). Operatorsplitting errors have been studied for a range of geochemical problems involving first-orderdegradation reactions (Kaluarachchi and Morshed, 1995), Monod kinetics (Morshed andKaluarachchi, 1995; Steefel and MacQuarrie, 1996), non-equilibrium sorption and equilibrium

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

199

linear adsorption (Barry et al., 1996), equilibrium nonlinear adsorption (Barry et al., 1997),
irreversible and reversible kinetic reactions (Carrayrou et al., 2004) and multicomponent transportwith surface complexation (Steefel and MacQuarrie, 1996). All of these studies were carried out forsteady-state flow conditions; almost no information about the operator-splitting errors is currently
available for transient flow conditions.
The objective of this study was to analyze numerical errors associated with the operator-splitting
approach in coupled reactive codes for typical soil-related flow and transport problems, especially for
applications involving transient flow and transport in heterogeneous soil profiles. Two benchmark
problems were considered for this purpose: (i) a first-order decay chain of linearly or nonlinearlysorbing contaminants during precipitation and evapotranspiration, and (ii) a multicomponent
transport problem involving cation exchange reactions in a multi-layered soil profile. Simulations
were carried out with the recently developed HP1 code (Jacques and Šimůnek, 2005) for variablysaturated water flow and multicomponent solute transport problems in soils. Results are comparedwith simulations using HYDRUS-1D (Šimůnek et al., 1998) and CRUNCH (Steefel, 2000, based on
the earlier GIMRT/OS3D package of Steefel and Yabusaki, 1996) for the first and second benchmark
problem, respectively.
2. Theoretical basis
2.1. Numerical model
The HP1 code (Jacques and Šimůnek, 2005) was used to evaluate the operator-splitting approachfor different flow and transport problems. The simulator is the result of coupling the variablysaturated water flow and solute transport model HYDRUS-1D (Version 2.0, Šimůnek et al., 1998)
and the geochemical model PHREEQC (Version 2.4, Parkhurst and Appelo, 1999). HP1 combines
the unique features of the two original models in one comprehensive tool by containing modules tosimulate (1) one-dimensional transient water flow in variably-saturated media, including root water
uptake and evaporation, (2) the transport of multiple components with provisions for physicalnonequilibrium (dual-porosity type) transport, (3) mixed equilibrium/kinetic geochemical reactions(aqueous speciation, ion exchange, mineral equilibrium), and (4) heat transport. This coupling of
HYDRUS-1D and PHREEQC makes HP1 a very flexible tool for addressing a large number of flow
and multicomponent transport problems in soils, including problems exhibiting heterogeneity inphysical and chemical properties. A full description of the processes, constitutive equations andnumerical solution strategies are found in the manuals of HYDRUS-1D (Šimůnek et al., 1998),
PHREEQC (Parkhurst and Appelo, 1999) and HP1 (Jacques and Šimůnek, 2005; www.sckcen.be/hp1). The main limitations of the current version 1.0 of HP1 are the absence of diffusion ofcomponents in the gas phase, and applicability to only one-dimensional flow and transport problems.
For completeness, the mathematical equations used in the examples discussed in this paper are
briefly summarized here. One-dimensional water flow in variably-saturated porous media is
described using the Richards equation:



∂hðhÞ
∂
∂h
¼
KðhÞ
þ 1 −SðhÞ
ð1Þ
∂t
∂x
∂x
where h is the water pressure head [L], θ is the volumetric water content [L3L− 3], t is time [T], x is
the spatial coordinate [L] (positive upward), S is a sink term [L3L− 3T− 1], and K is the unsaturated
hydraulic conductivity [LT− 1]. The water content and hydraulic conductivity are both nonlinear

200

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

functions of the pressure head. One possible way of describing these soil hydraulic properties iswith the closed-form expressions of van Genuchten (1980):hðhÞ ¼ hr þ

hs −hrð1 þ jahjn ÞmKðhÞ ¼ Ks Sel ½1−ð1−Se1=m Þm 2

ð2Þ
ð3Þ

where θr is the residual water content [L3L− 3], θs is the saturated water content [L3L− 3], α [L− 1],n [−] and m (= 1 − 1 / n) [−] are shape parameters, l is the pore connectivity parameter [−], Ks is thesaturated hydraulic conductivity [LT− 1], and Se = (θ − θr) / (θs − θr) is the effective saturation.
One-dimensional multicomponent transport in HP1 is described by means of aqueous masterspecies (also called components or primary species), which reduces the number of transport
equations to be solved. It is possible to write a given set of aqueous equilibrium reaction equationsin terms of a limited set of independent master species or components (Morel and Hering, 1993):
NmX

mji Amj ¼ Aið4Þ

j¼1

where Nm is the number of aqueous master species, i = 1, …, Ns, Ns is the number of aqueous
secondary species, Ajm and Ai are the chemical equations for the master and secondary species,
respectively, and νji are the stoichiometric coefficients in the reaction. The number of aqueous masterspecies equals the total number of aqueous species minus the number of aqueous equilibrium
reactions (if the reactions are written stoichiometrically independent). For each aqueous mastercomponent, solute transport is defined as (e.g., Lichtner, 1996; Mayer, 1999):


∂hCj
∂Cj
∂qCj
∂
hDw
¼
−SCr; j þ Ro; j
−
∂x
∂t
∂x
∂x

ð5Þ

where j = 1, ….., Nm, q is the volumetric flux density [LT− 1], S is the sink term in the flow equation(Eq. (1)), Dw is the dispersion coefficient in the liquid phase [L2T− 1], and Cj is the total liquid
concentration of master species j [ML− 3]:
Cj ¼ cj þNsXmji cið6Þ

i¼1

Additionally, Cr,j in Eq. (5) is the total concentration of the sink term (defined similarly as Cj), andRo, j is a reaction term that accounts for removal of a master species from the liquid phase due toaqueous kinetic reactions, and equilibrium or kinetic heterogeneous reactions [ML− 3T− 1] (e.g.,
sorption, ion exchange, mineral dissolution/precipitation). Kinetic reactions of virtually any form
can be given by means of small BASIC programs in the input file, thus providing great flexibility foraddressing a broad range of applications.
HP1 uses the operator-splitting approach without any iterations during a single time step (the noniterative sequential approach). Among others, Walter et al. (1994) used this approach in a simple twostep procedure to couple physical and chemical processes. HP1 uses a more complex four-step
approach (Fig. 1). Starting from the variables related to water flow (Hn, e.g., the pressure heads), heatD. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

201

Fig. 1. Schematic flow chart of the operator-splitting approach used in HP1.

flow (Tn, e.g., temperatures), solute transport (Cjn, e.g., the concentrations of master species j) and thesolid phase (Gn, e.g., the mineralogical composition, exchange sites) at a given time step n, thesequence of solving the governing equation in the next time step Δt is:(1) Calculate the new water flow properties pw(T n) and solve the water flow Eq. (1) using H n
to obtain H n + 1;

202

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218(2) Calculate the new heat transport properties ph(H n + 1) and solve the heat transport equation(Eq. 4.1 in Šimůnek et al., 1998; not further given here since only isothermal flow andtransport problems are considered here) using H n + 1 and T n to obtain T n + 1;
(3) Calculate the new solute transport properties ps(T n + 1,H n + 1) and solve the advection–
dispersion equation for each master species j without any chemical reactions (Eq. (5) without
the reaction term Ro,j; this source-sink term is solved in step 4 for all master species
simultaneously) using H n + 1 and Cjn to obtain Cjn,phys (i.e., the transported concentrations
without any chemical source or sink reactions). The advection–dispersion equation during
each time step is solved Nm times (i.e., separately for each master species); and(4) Calculate new geochemical properties pgc(T n + 1) and solve the system of equations for geochemical speciation using Cjn,phys ( j = 1, …, Nm) and G n to obtain Cjn + 1 (i.e., the concentrations
at the end of the time step resulting from equilibrium and kinetic reactions) and G n + 1. Thissystem is solved using the Newton–Raphson method for equilibrium reactions, together with theRunga–Kutta method for mixed equilibrium-kinetic reactions (see Parkhurst and Appelo, 1999).
The first three steps are solved with the HYDRUS-1D part of HP1, which thus acts as thetransport solver, whereas the PHREEQC part of HP1 is used for the fourth (geochemical) step.
A critical aspect of the solution scheme is optimal control of the time stepping process. Time
stepping in HP1 is controlled by means of three criteria within the HYDRUS-1D transport solver.
HYDRUS-1D has a self-adjusting time stepping scheme which changes the size of the time step
depending upon the number of iterations needed to solve the Richards equation (Šimůnek et al., 1998).
The time step, Δt, calculated in this manner is limited by two user-defined constraints: (i) the maximum
allowed time step, and (2) the performance index ωs (Perrochet and Berod, 1993), which is defined as
the product of the grid Peclet number, Pe (=qΔx/θDw) and the Courant number, Cr (=qΔt/θΔx). At
the beginning of a new time step, the product of Pe and Cr is calculated for each grid cell using the
updated Δt, but with q and θ from the previous time step. If the maximum value of this product islarger than the defined ωs, Δt is decreased to fulfil the performance index requirement. Finally, Δtcannot be larger than a predefined maximum allowed time step, Δt. This means that for relatively highwater fluxes, Δt is constrained mainly by ωs, whereas the maximum allowed time step constrains Δtmostly during periods with low water fluxes. Šimůnek et al. (1998) suggested values less than 5 for Pe
and less than 2 for ωs (Perrochet and Berod, 1993) when using HYDRUS-1D separately from HP1.
2.2. Test criteria for numerical accuracy
The accuracy of the operator-splitting approach was assessed by comparing concentration
profiles obtained with HP1 against concentration profiles obtained with codes that use a global
implicit approach (which is considered here to give the correct reference solution). Errors will bequantified using the L2 norm defined as (Steefel and MacQuarrie, 1996):
vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
u Ny
uX
y2i
L2 ¼ tð7Þi¼1

where y is the vector of concentration differences [ML− 3] between the HP1 results and the reference
solution and Ny is the number of data points. Note that L2 has the dimension of concentration, but thiswill not be further indicated in this paper. Depending upon the example problem discussed here, thereference solution is obtained by using (1) a global implicit approach that is implemented inD. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

203

HYDRUS-1D for the transport of a first-order decay chain of adsorbing contaminants duringtransient variably saturated flow (example 1), and (2) a global implicit approach that is used in
CRUNCH (Steefel, 2000) for the transport of major cations and heavy metals during steady-stateflow in a soil profile having a pH-dependent cation exchange complex (example 2).
3. Results and discussion
3.1. Transport of a first-order decay chain of adsorbing contaminants during transient flow
3.1.1. Problem definition
This example considers the transport of three hypothetical contaminants (A1, A2, A3) in a
homogeneous soil profile covered with grass subject to atmospheric boundary conditions defined bya time series of precipitation P and potential evapotranspiration ETP rates, and with free drainage as
the lower boundary condition. The time series for P, potential evaporation Ep, and potential
transpiration Tp were derived from climatological data for 3 years (1975–1977) from the Campaign
region (Belgium) using a crop factor of 0.8 for grass (Vanclooster et al., 1994) and a time-variableleaf area index (see e.g. Hupet et al., 2002 for estimating Ep and Tp from ETP). Fig. 2 shows the
cumulative P, Ep, and Tp rates during the 3 year simulation period. Roots were assumed to beuniformly distributed within a 20-cm deep root zone, thus implying a constant potential root water
uptake distribution in the root zone. The actual root water uptake rate may be lower due to soil
moisture stress, which was described with the Feddes water uptake reduction model (Feddes et al.,
1978; Table 1). The actual evaporation rate can similarly be lower than Ep when the hydraulicconductivity of the soil is too low to move enough water to the soil surface, thus causing a shortage ofmoisture near the soil surface. The actual evaporation rate was calculated from the potential
evaporation rate by limiting the surface pressure head to the value determined from equilibrium
conditions between soil water and atmospheric water vapour, and evaluating the actual surface flux
for this limiting pressure head. Geochemical reactions involve both equilibrium (adsorption) andkinetic (first-order degradation) reactions. The adsorption reactions in this example were assumed tobe linear for A1 and nonlinear (Freundlich sorption) for A2 and A3. Degradation reactions wereassumed to occur only in the liquid phase for A1, but in both the liquid and solid phases for A2 andA3. Note that both the equilibrium adsorption and degradation reactions were solved with theFig. 2. Potential and simulated actual fluxes during three growing seasons (1975–1977) for the transport of a first-order decay
chain of adsorbing contaminants. P — precipitation; Ep, Ea— potential and actual evaporation; respectively; Tp, Ta—potential and actual transpiration, respectively.

204

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

chemical solver of HP1 (i.e. PHREEQC), whereas inert transport was solved with the transport
solver of HP1 (i.e., HYDRUS-1D). Table 1 lists the various hydraulic, solute transport, root water
uptake and geochemical parameters of the assumed 1-m deep uniform medium-textured soil profile.
The numerical accuracy of the operator-splitting technique used in HP1 was assessed bycomparisons against a numerical reference solution obtained with HYDRUS-1D. The reference
solution was generated with a very small uniform spatial grid of 0.5 cm, a Peclet number of 0.5, and amaximum time step of 0.1 day. The maximum Courant number during the simulation was 0.197. The
0.5-cm discretization produced accurate water fluxes near the soil surface during both infiltration andevaporation (Van Dam and Feddes, 2000). The selected Pe and Cr values should also lead to a very
accurate solution of the solute transport equation for this problem (Šimůnek et al., 1998).
Fig. 2 shows the cumulative actual transpiration rate (root water uptake) Ta, the cumulative actual
evaporation rate Ea, and the cumulative actual surface flux P − Ea. Since the second growing season(year 1976) was very dry during summer, not enough soil water was available for potential
transpiration and, to a lesser extent, for potential evaporation. Because of the transpiration reduction,
the actual root water uptake rate hence was smaller than the potential root water uptake during thisperiod. The smaller Ea produced a slightly larger cumulative actual surface flux compared to theTable 1
Soil hydraulic, solute transport, root water uptake and geochemical parameters for simulating the transport of a first-orderdecay chain of adsorbing contaminantsParameter

SymbolValuesUnit

1

Soil hydraulic parameters
Saturated water content
Residual water content
Shape parameter
Shape parameter
Saturated hydraulic conductivity

Ks

0.43
0.078
0.036
1.56
24.96

cm3 cm− 3
cm3 cm− 3
cm− 1–
cm d− 1

Solute transport parameters
Bulk density
Dispersivity

ρb
λ

1.5
1

g cm− 3cm

Root water uptake parameters
Pressure head below which water uptake starts
Pressure head below which water uptake is maximal
Pressure head below which water uptake is not more at its maximum
Pressure head below which water uptake ceases

P0
Popt
P2
P3

Biogeochemical parameters
Adsorption isotherm coefficient for A1
Adsorption isotherm coefficient for A2
Freundlich exponent for A2
Adsorption isotherm coefficient for A3
Freundlich exponent A4First-order degradation rate for A1 in aqueous phase
First-order degradation rate for A2 in aqueous phase
First-order degradation rate for A3 in aqueous phase
First-order degradation rate for A2 in solid phase
First-order degradation rate for A3 in solid phase

Ka1
Ka2
n2
Ka3
n3
ka1
ka2
ka3
ks2
ks3

1

Parameters representative of a loamy soil (from Šimůnek et al., 1998).

θs
θr
α

−10
−25
− 1000− 8000

1
2.5
0.9
5
0.8
0.005
0.06
0.02
0.006
0.002

cm
cm
cm
cm
cm3 g− 1
cm3 g− 1–
cm3 g− 1–day− 1day− 1day− 1day− 1day− 1

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

205

potential surface flux. Note that differences between the potential and actual fluxes were almost
negligible during the first and third growing seasons (years 1975 and 1977, respectively).
Prior to analyzing the accuracy of HP1, the effect of time stepping and discretization on inert solutetransport was evaluated using HYDRUS-1D. We subsequently evaluated the accuracy of the HP1
results using three cases: (1) only equilibrium adsorption reactions; no degradation, (2) only degradation reactions; no sorption, and (3) mixed equilibrium sorption and kinetic degradation reactions.
3.1.2. Inert transportThis problem was simulated with HYDRUS-1D using three grid sizes (Δx = 0.5, 1, and 2 cm)
assuming a maximum time step of 1 day, as well as with the reference solution scheme as discussed
earlier (i.e., with Δx = 0.5 cm and a maximum Δt = 0.1 day). The L2 norm for the simulated
concentration profiles was calculated every 15 day. Fig. 3 shows L2 values as a function of time.
Larger grid sizes were found to produce larger values of L2. Notice that L2 varied considerably intime, with the largest values occurring during periods with high Ep, even when we used the same
spatial discretization but a different maximum time step. During periods with high P, L2 values werequite similar for the simulations with Δx = 0.5 cm and Δx = 1 cm (i.e., between 0.9 and 1.25 year andbetween 1.75 and 2.4 year). L2 values decreased slightly when a smaller maximum time step wasused (results not shown). The different spatial and time discretizations also produced differences in
the water content and flux profiles (results not shown). Van Dam and Feddes (2000) showed that
calculated evaporation rates sometimes change with increasing grid size, leading to errors in thewater content and the fluxes. These errors in the flow field will induce numerical errors in the soluteconcentration, and hence may interfere with an accurate assessment of the operator-splitting errors
discussed in the next sections. However, the concentration profiles obtained with the different grid
sizes were consistent with each other and showed only limited enhanced numerical dispersion withincreasing grid size.
Possible HP1 operator-splitting errors were quantified in terms of: (i) total L2 = (∑L22)0.5, (ii) wet
L2, being the same as total L2 but limited to concentration profiles between 1.8 and 2.2 years (arelatively wet period with mainly infiltration), and (iii) max L2. Wet L2 was used to evaluate theaccuracy of the sequential non-iterative approach when errors related to grid discretization are
expected to be low.

Fig. 3. Time series of L2 values for the transport of an inert solute obtained with HYDRUS-1D for three grid sizes and amaximum time step of 1 day (reference solution is for Δx = 0.5 cm and a maximum time step of 0.1 day).

206

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

3.1.3. Transport subject to equilibrium adsorption
Next, the three contaminants were considered to undergo only equilibrium adsorption (linear
for A1, nonlinear for A2 and A3), but no kinetic degradation. No contaminants were assumed to be
initially present and possible density effects after their application were neglected. Since nodegradation took place (no degradation products from A1), the concentration of each contaminant
in the inflowing rain water was assumed to be simply 1 mol l− 1 for all precipitation events.
Simulations were carried out with different grid sizes, maximum time steps and varying
performance indices ωs (see Table 2). Fig. 4 shows total L2, wet L2, and max L2 for A1 and A2 as
a function of simulation time (A3 concentrations were very similar to those of A2, and are notfurther shown here). Results for the HYDRUS-1D simulations with different spatial
discretizations and a maximum time step of 1 day are shown as straight lines. Note that while
Table 2
Overview of test simulations of the transport of the first-order solute decay chain
Run

Δx
(cm)

Maximumtime step
(day)

ωs

Simulation time (min)a
Only equilibrium
reactions

Only kinetic
reactions

Mixed equilibrium
kinetic reactions

Reference case HYDRUS-1D
0.5
0.1

2

4.6

2.4

4.1

HYDRUS-1D runs
a
2
b
1
c
0.5

1
1
1

2
2
2

0.8
1.5
3.1

0.5
1.0
2.1

0.7
1.4
1.3

HP1 runs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

1
1
1
0.25
1
1
1
0.75
0.50
0.25
0.125
1
0.75
0.50
0.25
0.125
1
1
1
1
0.25
1

1
0.5
0.25
0.25
0.1
1
0.5
0.5
0.5
0.5
0.5
0.25
0.25
0.25
0.25
0.25
0.1
1
0.5
0.25
0.25
0.1

a

2
2
2
2
2
1
1
1
1
1
1
1
1
1
1
1
1
0.5
0.5
0.5
0.5
0.5

16
22
19
20
28
33
34
35
35
38
50
38
38
38
40
51
60
73
73
80
86
116

82
84
94
104
131
164
172
172
179
195
229
187
189
190
202
235
275
377
378
417
391
526

76
78
85
108b
117
152
157
162
163b
205b
264b
171b
172
173
187
237b
269
382
385
422b
507
598b

On a Pentium(R) 4, 3.06 GHz.
These runs did not finish in one time, but were restarted from the last saved information before failure. All runs finished using
this two-step approach. Hypothetical time the simulation would need if completed in one run; value was calculated based on thenumber of time steps multiplied by the average computational time needed for one time step during the first part of the simulation.
b

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

207

Fig. 4. L2 measures for simulations of the transport of adsorbing contaminants (only equilibrium reactions) for A1 (a–c)
and A2 (d–e) (square: Δx = 0.5 cm; open star: Δx = 1 cm, ωs = 1; solid star: Δx = 1 cm, ωs = 0.5; open diamond: Δx = 1 cm,
ωs = 0.25; solid diamond: Δx = 1 cm, ωs = 0.1; open circle: Δx = 2 cm). Also shown are HYDRUS-1D results for
Δx = 0.5 cm (solid lines), Δx = 1 cm (short dashed lines), and Δx = 2 cm (long dashed black lines). The long curved dashedgrey lines in (a) and (d) are indicative of Pareto fronts.

we did not identify every simulation separately in Fig. 4, the smaller ωs or the maximum time
step, the smaller generally the accuracy measures. Using the computational times listed in Table 2,
it is possible to identify the individual simulation runs as needed.
Contaminant A1 was found to be the most mobile, and A3 the least (results again not shown for thelatter). Errors associated with the HYDRUS-1D grid discretization were found to be positively
correlated with the mobility of the contaminants: largest for A1 and smallest for A3. There are several
reasons for this: (1) since the more mobile solutes travelled deeper in the soil profile compared to theless mobile solutes, more non-zero terms are included in L2, (2) due to the higher sorption,
concentrations in the liquid phase are smaller for the less mobile solutes, thus making L2 smaller, and(3) errors in the fluxes for larger grid sizes have less impact on the transport of the less mobile solutes.
By comparison, the less mobile contaminants produced larger errors for the operator-splitting
scheme relative to the grid discretization errors of HYDRUS-1D. HP1 simulations with Δx = 1 or
2 cm reached the total L2 value for the HYDRUS-1D simulation for A1, but not with Δx = 1 cm forA2. The HP1 simulations with Δx = 0.5 always produced larger total L2 values as compared to
HYDRUS-1D (run c). These results suggest that operator-splitting errors during transient flow may
increase with decreasing grid size (as opposed to grid discretization errors as evaluated with
HYDRUS-1D) and decreasing mobility (or increasing nonlinearity for small grid sizes).

208

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

We now focus on the HP1 simulations (symbols in Fig. 4) in order to provide insight in theoperator-splitting errors, and to identify possible strategies for maximizing accuracy with only alimited increase in computational time. Values of total L2 for A1 for Δx = 1 and Δx = 0.5 cm werequite similar. This suggests that the same level of accuracy can be obtained with a coarser grid(Δx = 1 cm with adjusted ωs value and maximum time step) and with less computational time ascompared to the finer grid (Δx = 0.5 cm). Even coarser grids could be used for the less mobile
nonlinearly adsorbing contaminants (A2 and A3), which would lead to significantly smaller
computational times (Table 2).
A comparison of wet L2 and max L2 reveals several other aspects of the accuracy of HP1 fortransient flow problems. As a reference, lines (with arrows) connect simulations with the same
maximum time step (1 day) and ωs (but for different Δx): ωs = 0.25 (circle, open diamond andsquare for Δx = 2, 1, and 0.5 cm, respectively) and ωs = 0.1 (circle, solid diamond and square for
Δx = 2, 1, and 0.5 cm, respectively). Simulations with Δx = 0.5 cm performed better than thecoarser grid sizes during the main infiltration periods for a given combination of ωs andmaximum time step (wet L2 in Fig. 4). On the other hand, the fine discretization produced a larger
max L2, especially for larger ωs values. As was the case with the HYDRUS-1D runs, the max L2
occurred during high Ep boundary conditions (Fig. 3 for the HYDRUS-1D simulations). The
increases in max L2 for A1 when Δx decreased from 1 to 0.5 cm were 40 (i.e., from 0.5488 to
0.7687) and 46% (i.e., from 0.2076 to 0.3039) for ωs 0.25 and 0.1, respectively. However, sinceerrors are smaller for the finer grid during wet boundary conditions, total L2 was only 25% larger
for the finer grid when ωs = 0.25, whereas the total L2 is similar when ωs = 0.1. The max L2
decreased significantly when we used smaller ωs values (Δx = 0.5 cm) or a combination of smaller
maximum time steps and ωs values (Δx = 1 cm). For the less mobile nonlinearly adsorbing
contaminants (A2 and A3), smaller ωs values for the simulations with Δx = 0.5 cm did not produce
smaller max L2-values than those obtained with the coarsest discretization, as was the case for A1.
The simulations with Δx = 1 cm allow us to identify how operator-splitting errors can be reduced
in a computationally efficient way by adjusting the maximum time step and/or ωs. Wet L2 can be
reduced most efficiently by reducing ωs. Decreasing ωs from 1 to 0.25 (with a maximum time step of
1 day) produced significantly lower wet L2 values without a significant increase in computational
effort. On the other hand, decreasing the maximum time step resulted in a significant decrease in max
L2, with a negligible increase in simulation time (maximum time steps of 1, 0.75, 0.5 and 0.25 forωs = 0.5 and 0.25). This was to be expected since wet L2 was calculated during a period with highinfiltration rates, a wet soil profile, and high water fluxes, and thus likely high Courant numbers. The
performance index ωs determines the maximum allowed Cr (for a fixed Pe) and thus the time step.
Reducing ωs hence leads to smaller time steps when water fluxes are high. As discussed earlier, max
L2 occurs during periods with high Ep (and low precipitation), and hence when the water fluxes in thesoil are relatively small. The time step is then not controlled by ωs, but by the maximum time step.
Decreasing the latter thus should have its largest effect on L2 during dry periods.
The plots in Fig. 4 can be interpreted also in terms of a multi-objective optimization with accuracy(as estimated with the different L2 measures) and computational time as the two objectives. The grey
dashed lines in Fig. 4a and d show so-called Pareto fronts (Stadler, 1988). Sharp Pareto fronts for allcontaminants indicate that the accuracy can be significantly improved with minor increases incomputational time. However, further improvement in accuracy beyond a given level may requireconsiderably more computational time. Several non-Pareto optimal runs appear to exist, including 4of the 5 runs for Δx = 0.5, some runs for Δx = 1 cm for A1, and many runs for Δx = 1 cm for A2.
Relatively good accuracy for all contaminants may be obtained by using a grid discretizationof 1 cm, ωs = 0.25, and maximum time steps of 0.25 day resulting in a computational time of

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

209

Fig. 5. Concentration profiles of A1 (a) and A2 (b) for the transport of adsorbing contaminants (only equilibrium reactions)for the reference case (black lines) and run 15 (grey lines with dots) at t = 0.53 year and t = 1.81 years. These times
correspond to the time with max L2 and the largest L2-value during the wet period (between 1.8 and 2.2 years).

40 min. The accuracy could be further improved slightly by reducing the maximum time step to0.125 day, which would increase the computational time by 10 min. While this increase may seem
relatively small, the computational time would increase to 60 min when the equilibrium reactions
are combined with the kinetic degradation reactions (to be discussed later, and shown in Table 2).
Fig. 5 shows concentration profiles for A1 and A2 for the max L2 value and the highest L2
value during the ‘wet’ period. Differences between the reference case and run 15 (Table 2) werevery small. The largest difference occurred at the soil surface when Ep was large.
3.1.4. Transport subject to kinetic degradation
The same runs as in the previous example were carried out for a problem involving only kinetic
degradation reactions (using parameters as defined in Table 1). Computational times were much

Fig. 6. Total L2 for simulations of the transport of a first-order decay chain of three contaminants (only kinetic reactions)for A1. Symbols are explained in Fig. 4.

210

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

Fig. 7. Water content and concentration profiles at selected times for the transport of a first-order decay chain of adsorbing
contaminants obtained with the reference run (black lines) and run 15 (grey lines with dots).

larger than those needed for the runs with only equilibrium reactions. Since kinetic reactions in
PHREEQC are solved using the Runga–Kutta technique, the computational time for solving thegeochemical part of the problem was much larger than when only equilibrium reactions arepresent. Optimizing the grid size and time steps is hence much more important when kinetic
reactions are involved.
Table 3
Soil hydraulic properties and cation exchange capacities of the five soil horizons of the hypothetical sandy soil used for test
example 1
HorizonLayer thickness (cm)θrθsα (cm− 1)nKs (cm day− 1)lCation exchange capacity(eq/1000 cm3 soil)

A
E
Bh1
Bh2
Bh/C

13
10
5
5
17

0.065
0.035
0.042
0.044
0.039

0.476
0.416
0.472
0.455
0.464

0.016
0.015
0.016
0.028
0.023

1.94
3.21
1.52
2.01
2.99

93
311
39
860
1198

0.5
0.5
0.5
0.5
0.5

0.01830.01140.06640.05420.0116

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

211

Table 4
Hypothetical pH and solution concentrations of the initial soil solution and the infiltrating solution (μmol l− 1) for example 2
Solution

pH

Na⁎

K

Ca

Mg

Br

Cd

Pb

Zn

0–28 cm depth
28–50 cm depth
Applied water

8.5
8.5
3.5

401.9
454.0
127.5

120
120
120

98
98
98

5
5
5

780
780
780

0.8
0.0
0.0

2.5
0
0

50
0
0

⁎ Concentration of Na was adjusted to obtain the desired pH.

Fig. 6 shows the total L2 for A1 for different runs. Results for the other L2 measures and the othertwo contaminants were very similar and are not shown here. Notice that the L2 values are larger thanthose for the simulations with only equilibrium adsorption reactions (Fig. 4). When only kinetic
degradation occurs and no sorption, the contaminants are more mobile, leading to larger errors for thecoarser spatial discretizations (see discussion in Section 3.1.3.), with the errors becoming of the same
magnitude as for the inert tracer simulations (Section 3.1.2, Fig. 3).
HP1 errors were close to those obtained with HYDRUS-1D for different grid discretizations,indicating that operator-splitting errors are small. This is due to the fact that the first-orderdegradation constants are small relative to the time steps. Carrayrou et al. (2004) showed thatwhen the dimensionless time NOS1 (= kΔt with k the first-order degradation constant) is smaller
than 0.1, mass balance errors will be less than 1%, also for the non-iterative operator-splitting
approach. For a maximum time step of 1 day, NOS1 in our example was less than 0.1, whichshould lead to relatively small mass balance errors. Also, the time steps were generally smaller
than the maximum time step for the transient flow simulations, thus leading to even smaller errors.
3.1.5. Transport subject to combined equilibrium adsorption and degradation
This last step of our analysis considers a combination of equilibrium adsorption and kinetic
degradation using the full set of parameters defined in Table 1. Not all runs finished the 3 years
simulation in one time, notably those for Δx = 0.5 cm, and those with a maximum time step of0.125 day or ωs = 0.1 for Δx = 1 cm (the reason for this failure was not clear since no convergenceerror in PHREEQC occurred). Each run, however, could be finished when it was restarted from thelatest save of the profile information. The various accuracy measures were limited to only the first
1.8 year (i.e., shortest run), while wet L2 was calculated between 0.74 and 1.11 years. Table 2 lists

Table 5
Overview of aqueous equilibrium reactions and corresponding equilibrium constants for example 2 (data from the phreeqc.
dat database, Parkhurst and Appelo, 1999)
Nr

Aqueous speciation reaction
−

Log_K
− 14
− 14.18− 14.46− 12.78− 11.44

+

(1)
(2)
(3)
(4)
(5)

H2O = OH + H
Na+ + H2O = NaOH + H +K+ + H2O = KOH + H+
Ca2+ + H2O = CaOH + + H +Mg2+ + H2O = MgOH + + H +

(6)
(7)
(8)
(9)

Me2+ + H2O = MeOH + + H +Me2+ + 2H2O = Me(OH)2 + 2 H +Me2+ + 3H2O = Me(OH)−3 + 3 H ++Me2+ + 4H2O = Me(OH)2−
4 +4 H

Cd

Pb

Zn− 10.08− 20.35− 33.30− 47.35− 7.71− 17.12− 28.06− 39.70− 8.96− 16.90− 28.40− 41.20

212

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

Table 6
Log K parameters for the multi-site exchange complex used in example 2
Y-exchanger a

NaY

KY

MgY2

CaY2

CdY2

PbY2

ZnY2

HYb

− 1.0
HYa
1.65

−0.3
HYb
3.3

− 0.4
HYc
4.95

− 0.2
HYd
6.85

−0.2
HYe
9.6

0.05
HYf
12.35

− 0.2

a

The value for NaY was taken from Appelo et al. (1998). Values for the other complexes were taken from the phreeqc.
dat database (Parkhurst and Appelo, 1999) and adapted relative to the K for NaY.
b
Values taken from Appelo et al. (1998).

computational times for those runs that finished normally, as well as extrapolated computationaltimes for those that were terminated early. The computational times were similar to those obtained
for the kinetic reactions alone.
Although the L2 measures were calculated for a different time period, similar Pareto fronts wereobserved as for simulations with only equilibrium reactions (results not shown). Again, a goodcompromise between required computational time and accuracy was obtained with Δx = 1 cm,
ωs = 0.25, and a maximum time step of 0.25 day. It is important to note that the wet L2 for Δx = 0.5 cm
and ωs = 0.1 also resulted in the lowest value for all runs.
The water content and concentration profiles at selected times for run 15 were found to be ingood agreement with the reference HYDRUS-1D simulations (Fig. 7). In general, errors
associated with both operator-splitting and grid discretization were relatively small, leading toaccurate results for HP1. Notice that the effects of root water uptake are clearly visible in Fig. 7 by
the lower water contents in the top 20 cm of the profile at t = 225 day.
3.2. Heavy metal transport with a pH-dependent cation exchange complex
In order to describe the competition between protons, major cations and heavy metals foradsorption sites as a function of pH, Appelo et al. (1998) previously proposed a multi-site cationexchange complex consisting of several sites, each having a different selectivity coefficient for theexchange of protons. Their model captures the typical behaviour and increased capacity of thesurfaces to adsorb major cations and heavy metals at higher pH values. The model for cationexchange in the presence of organic matter and Fe-oxyhydroxides consists of six sites. The multi-site
cation exchange complex model was used in the second test example to simulate the transport ofseveral major cations (Na, K, Ca, and Mg) and three heavy metals (Cd, Zn, and Pb) in a heavy-metal

Table 7
Overview of test simulations of the transport of heavy metals with a pH-dependent cation exchange complex
Δx = 0.5 cm

Δx = 1 cm

Δx = 2 cm

Maximum time step (day)ωsMaximum time step (day)ωsMaximum time step (day)ωs

0.01450.01080.00730.0036

0.1
0.075
0.05
0.025

0.02900.217
0.01460.00720.00430.0029

0.2
0.15
0.1
0.05
0.03
0.02

0.05790.04340.02920.01450.00870.00580.0029

0.4
0.3
0.2
0.1
0.06
0.04
0.02

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

213

Fig. 8. Total L2 summed over 7 cations for simulations of the transport of major cations and heavy metals in a soil with apH-dependent cation exchange complex during steady-state flow. The grey dashed line indicates the Pareto front. Numbers
indicate Pareto optimal runs shown in Fig. 9.

contaminated soil that initially had a high pH (8.5), but was subjected to infiltration with an acidheavy-metal-free solution (pH 3).
The 60-cm deep soil profile was assumed to contain five distinct soil horizons with different soil
hydraulic properties and exchange capacities (Table 3). The reported CEC is the sum of the CEC of
the six sites, with each site having the same CEC. A constant water flux of 0.05 m day− 1 and freedrainage were taken as the upper and lower boundary conditions, respectively. This situation leads tounsaturated flow with the water content varying with depth. We used a dispersivity of 0.05 m and adiffusion coefficient of 9.2 · 10− 10 m2 s− 1. Table 4 lists the concentrations of the initial soil solution
and the infiltrating solution. Note that only the top 30 cm of the soil profile was initially contaminatedwith the three heavy metals. The considered aqueous equilibrium reactions (with their corresponding
equilibrium constants) and the parameters for the multi-site exchange complex are given in Tables 5
and 6, respectively. The six sites are identified in Table 6 as Ya, Yb, …, Yf, with equilibrium constantsfor the same element (Na, K, Mg, Ca, Cd, Pb and Zn) being identical (first two rows in Table 6, where
Y is understood to refer to Ya, …, Yf). The sites all have different equilibrium constants for theexchange of the protons, as shown by the third and fourth rows in Table 6.
HP1 simulations in terms of L2 will be compared with results obtained with CRUNCH (Steefel, 2000),
which solves the problem using a global implicit approach and therefore is considered to be the reference
run. The reference run used a spatial discretization Δx=0.5 cm and a maximum time step of 0.055 day.
HP1 simulations were performed for three spatial discretizations (Δx=0.5, 1, and 2 cm) and for different
maximum time steps (d) with corresponding ωs-values (Table 7). L2 values were calculated for theconcentration profiles of each component for every 5 days between day 0 and day 365.
Fig. 8 shows a plot of the total L2 (= (∑ L22)0.5) for all 7 cations combined (Na, K, Ca, Mg, Cd,
Pb, Zn) versus the required computer time (for a Pentium 4 3.06 GHz computer). Different
simulations for a given Δx can be identified by the longer computational time for smaller ωsvalues. For steady-state flow conditions, a finer spatial discretization produced a smaller L2 norm.
This is in contrast to the transient flow simulations which showed slightly larger L2 norms withfiner grids (see Section 3.1.2), but similar to the wet L2 norm defined in Section 3.1.2 (onlyequilibrium reactions). Neglecting periods with high evapotranspiration (as was done for the wet
L2 norm), the finer discretizations produce more accurate results. For the steady-state flow
situations considered here, this effect was very clear. The coarser grid sizes never reached theaccuracy of the smaller grid sizes, even after reducing the maximum time step for the coarser

214

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

grids. However, the difference between Δx = 2 cm and Δx = 1 cm was much larger than the
difference between Δx = 1 and Δx = 0.5 cm.
A Pareto front drawn through the data in Fig. 8 shows the non-optimal runs. Cd outflow curves
and profiles of Ca, Cd and pH for three Pareto front runs (Δx = 1 cm, ωs = 0.05; Δx = 0.5 cm,
ωs = 0.1 and 0.05) are shown in Fig. 9 together with the reference run using CRUNCH. Infiltration
with the low-pH solution causes a decrease in the capacity of the solid surfaces to adsorb heavy
metals in favour of protons. This will lead to desorption of heavy metals and subsequent leachingfrom the soil profile (Fig. 9a, showing leaching of Cd).
3.3. Strategy for reducing operator-splitting errors
The spatial and time discretizations suggested earlier hold only for the specific combination ofsoil, time series of ETp and P, and the particular geochemical properties of the contaminants used inour examples. A more elaborate analysis covering a range of contaminant properties, soil types, andinfiltration and evapotranspiration rates is required to provide more general guidelines. Still, ouranalyses indicates that (i) operator-splitting errors are small for kinetic reactions when the kinetic

Fig. 9. Calculated Cd outflow curves (a) and profiles of pH (b), Ca concentrations (c) and Cd concentrations (d) for thetransport of major cations and heavy metals in a soil with a pH-dependent cation exchange complex during steady-state
flow. HP1 (Δx = 0.5 cm, ωs = 0.05), HP1 (Δx = 0.5 cm, ωs = 0.1), and HP1 (Δx = 1.0 cm, ωs = 0.05) correspond to thePareto optimal runs 1, 2, and 3 indicated in Fig. 8.

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

215

rates are small relative to the time steps. Given that time steps during transient flow simulations are
small, operator-splitting errors can be easily controlled for these type of problems; (ii) in case of
equilibrium heterogeneous reactions (linear and nonlinear sorption and cation exchange studiedhere), numerical errors are reduced by reducing the grid size for steady-state flow conditions or
during wet periods with mainly infiltration; and (iii) during transient flow (especially during periods
with mainly evaporation), operator-splitting errors depend in a complex manner on the griddiscretization and the time stepping. In that case, the finest grid discretization with the smallest
ωs value and maximum time step will lead to the smallest wet L2 value. Our simulation results hence
may serve as a preliminary guide for selecting a reasonable spatial and temporal discretization.
In view of our simulations we suggest the following procedure for selecting time step
constraints for other problems: (1) select a short period from a time series of P and ETp containing
both wet and dry periods with high ETp (e.g., 2 years out a 30-year simulation period) or one
simulation out of a set of many (e.g., in uncertainty analyses), (2) run the model with a relatively
fine grid, low ωs and a small maximum time step to obtain a reference solution (which may
require a long computational time), (3) run the model for several cases with coarser grids and alarger ωs, (4) calculate wet L2 using the reference run to analyze the effect of grid discretizationand ωs, and determine an optimal Δx, (5) run the model for several cases with larger maximum
time steps for a fixed (optimum) grid discretization and possibly also a fixed ωs, and (6) calculate
max L2 using the case with the smallest maximum time step as reference for a fixed griddiscretization to analyze the effect of the maximum time step on the error.
4. Conclusions
A non-iterative sequential coupling approach may be used to combine separate water flow andsolute transport codes with geochemical speciation codes. Although such coupling has many
advantages from a programmer's perspective, its most important disadvantage is the possible
occurrence of operator-splitting errors. This paper investigates several ways to reduce these typeof errors in transient flow problems by optimizing grid discretization and time stepping, the latterin terms of two controlling parameters: the maximum time step and the performance index ωs.
The most critical system we investigated involved heterogeneous equilibrium reactions duringtransient flow conditions (example 1). Contrary to the global-implicit approach, we found that finer
spatial discretizations do not always lead to smaller numerical errors, especially for less mobile
solutes with fine grids during dry periods with a high evapotranspiration demand. This result is incontrast with steady-state flow problems where finer grids did produce smaller operator-splitting
errors. Also, for (relatively slow) kinetic reactions, coarser grids introduced larger numerical errors
for transient flow situations, but results were independent of the operator-splitting errors due to thesmall values of the dimensionless rate constant NOS1 (smaller than 0.1, Carrayrou et al., 2004). We
note here that time steps during transient variably-saturated flow simulations often are alreadyrelatively small to enable accuracy in solutions of the flow equation, thus frequently ensuring small
NOS1 values. The observation that numerical errors are in some cases larger during transient flow
simulations when finer grids are used indicates that operator-splitting errors are more significant
during periods with high Ep boundary conditions with finer grids.
Adequate time step control is required for transient flow simulations when using a reactive
coupled code. An optimal combination of the two criteria determining the time steps (maximum
time step and ωs) is critical. To reduce numerical errors during periods of infiltration, ωs is themost appropriate and efficient parameter. The maximum time step is the more important factor
determining numerical errors during time periods with high evapotranspiration rates.

216

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

Since our analysis was carried out for only one soil type and a limited number of contaminantswith different adsorption properties, results should not be used to provide general guidelines for
optimal choice of the performance index and the maximum time step. Following the procedures usedin this paper, a fine spatial grid discretization (e.g., Δx = 0.5 cm) could be used in combination with asmall ωs value (e.g., 0.125) and a small maximum time step (e.g., 0.125 day depending on the total
simulation time and water fluxes) to first generate an accurate numerical solution during infiltration
events at high computational cost. This reference solution could then be used to test the accuracy ofalternative spatial and time discretizations for a particular problem at a lower computational cost, butwith increased accuracy during periods of high evapotranspiration.
Acknowledgements
This work was partly supported by SAHRA (Sustainability of semi-Arid Hydrology andRiparian Areas) under the STC Program of the National Science Foundation, Agreement No.
EAR-9876800 and the Terrestrial Sciences Program of the U.S. Army Research Office (Terrestrial
Processes and Landscape Dynamics and Terrestrial System Modeling and Model Integration).
Additional support was obtained through the bilateral agreement project “Development andevaluation of a coupled geochemical transport model between SCK•CEN and USDA-ARS(Agreement No. 58-5310-0-F105), between the SCK•CEN and the University of California,
Riverside (Agreement No. C0-90001412.01), and by the SCK•CEN R and D program on“Environmental Remediation” (project E022031 and CO91002). We acknowledge the insightfulcomments and suggestions provided by two anonymous reviewers and the editor, which
considerably improved the quality of the paper.

References
Adler, M., 2001. Interaction of claystone and hyperalkaline solutions at 30 °C: A combined experimental and modelingstudy. Ph. D., Universitaet Bern, 120 p.
Appelo, C.A.J., Verweij, E., Schäfer, H., 1998. A hydrogeochemical transport model for an oxidation experiment with
pyrite/calcite/exchangers/organic matter containing sand. Appl. Geochem. 13, 257–268.
Ayora, C., Taberner, C., Saaltink, M.W., Carrera, J., 1998. The genesis of dedolomites: a discussion based on reactive
transport modeling. J. Hydrol. 209, 346–365.
Barry, D.A., Miller, C.T., Culligan-Hensley, P.J., 1996. Temporal discretisation errors in non-iterative split-operator
approaches to solving chemical reaction/groundwater transport models. J. Contam. Hydrol. 22, 1–17.
Barry, D.A., Miller, C.T., Culligan, P.J., Bajracharya, K., 1997. Analysis of split operator methods for nonlinear andmultispecies groundwater chemical transport models. Math. Comput. Simul. 43, 331–341.
Carrayrou, J., Mosé, R., Behra, P., 2004. Operator-splitting procedures for reactive transport and comparison of massbalance errors. J. Contam. Hydrol. 68, 239–268.
Davis, J.A., Meece, D.E., Kohler, M., Curtis, G.P., 2004. Approaches to surface complexation modeling of Uranium(VI)adsorption on aquifer sediments. Geochim. Cosmochim. Acta 68, 3621–3641.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of field water use and crop yield. John Wiley and Sons, New
York, NY.
Hundsdorfer, W., Verwer, J.G., 1995. A note on splitting errors for advection-reaction equations. Appl. Numer. Math. 18,
191–199.
Hupet, F., Lambot, S., Javaux, M., Vanclooster, M., 2002. On the identification of macroscopic root water uptake
parameters from soil water content observations. Water Resour. Res. 38, 1300, doi:10.1029/2002WR001556.
Jacques, D., Šimůnek, J., 2005. User manual of the multicomponent variably-saturated flow and transport model HP1:
Description, Verification, and Examples (Version 1.0). SCK•CEN, Mol, Belgium, BLG-998, 79 p.
Kaluarachchi, J.J., Morshed, J., 1995. Critical assessment of the operator-splitting technique in solving the advection–
dispersion-reaction equation: 1. First-order reaction. Adv. Water Resour. 18, 89–100.

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

217

Kent, D.B., Abrams, R.H., Davis, J.A., Coston, J.A., LeBlanc, D.R., 2000. Modelling the influence of variable pH on the transportof zinc in a contaminated aquifer using semi-empirical surface complexation models. Water Resour. Res. 36, 3411–3425.
Kohler, M., Curtis, G.P., Kent, D.B., Davis, J.A., 1996. Experimental investigation and modeling of uranium(VI) transport
under variable chemical conditions. Water Resour. Res. 32, 3539–3551.
Langergraber, G., Šimůnek, J., 2005. Modeling variably-saturated water flow and multicomponent reactive transport inconstructed wetlands. Vadose Zone J. 4, 924–938.
Lichtner, P.C., 1996. Continuum formulation of multicomponent–multiphase reactive transport. In: Lichtner, P.C., Steefel,
C.I., Oelkers, E.H. (Eds.), Reactive transport in porous media. Reviews in Mineralogy, vol. 34. Mineralogical Societyof America, Washington, DC, pp. 1–81.
Mayer, K.U., 1999. A numerical model for multicomponent reactive transport in variably saturated porous media. Ph.D.
thesis, Department of Earth Sciences, University of Waterloo, 285 p.
Mayer, K.U., Frind, E.O., Blowes, D.W., 2002. Multicomponent reactive transport modeling in variably saturated mediausing a generalized formulation for kinetically controlled reactions. Water Resour. Res., doi:10.1029/2001WR000682.
Meeussen, J.C., Scheidegger, A., Hiemstra, T., van Riemsdijk, W.H., Borkovec, M., 1998. Predicting multicomponent
adsorption and transport of fluoride at variable pH in a goethite–silica sand system. Environ. Sci. Technol. 30, 481–488.
Morel, F., Hering, J., 1993. Principles and applications of aquatic chemistry. New York, John Wiley and Sons Inc.
Morshed, J., Kaluarachchi, J.J., 1995. Critical assessment of the operator-splitting technique in solving the advection–
dispersion-reaction equation: 2. Monod kinetics and coupled transport. Adv. Water Resour. 18, 101–110.
Öztürk, H.S., Özkan, İ., 2004. Effects of evaporation and different flow regimes on solute distribution in soil. Transp.
Porous Media 56, 245–255.
Parkhurst, D.L., Appelo, C.A.J, 1999. User's guide to PHREEQC (Version 2) — A computer program for speciation,
batch-reaction, one-dimensional transport and inverse geochemical calculations. Water-Resources Investigations,
Report 99-4259, Denver, Co, USA, 312 pp.
Parkhurst, D.L., Kipp, K.L., Engesgaard, P., Charlton, S.C., 2004. PHAST — A program for simulating ground-water flow,
solute transport and multicomponent geochemical reactions. U.S.G.S. Techniques and Methods, vol. 6-A8. 154 pp.
Perrochet, P., Berod, D., 1993. Stability and the standard Crank–Nicolson–Galerkin scheme applied to the diffusion–
convection equation: some new insights. Water Resour. Res. 29, 3291–3297.
Prommer, H., Barry, D.A., Zheng, C., 2003. MODFLOW/MT3DMS based reactive multi-component transport modelling.
Ground Water 41, 247–257.
Saaltink, M.W., Battle, F., Ayora, C., Carrera, J., Olivella, S., 2004. RETRASO, a code for modeling reactive transport insaturated and unsaturated porous media. Geol. Acta 2, 235–251.
Samper, J., Juncosa, R., Delgado, J., Montenegro, L., 2000. CORE2D: A code for nonisothermal water flow and reactive
solute transport. User manual version 2. ENRESA Technical Publication 06/2000, 131 pp.
Schäfer, D., Schäfer, W., Kinzelbach, W., 1998. Simulation of reactive processes related to biodegradation in aquifers. 1.
Structure of the three-dimensional reactive transport model. J. Contam. Hydrol. 31, 167 186.
Šimůnek, J., Valocchi, A.J., 2002. Geochemical transport, In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis,
Part 1: Physical Methods, Chapter 6.9, Third ed. SSSA, Madison, WI, pp. 1511–1536.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1998. The HYDRUS-1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variably-saturated media. Version 2.0, IGWMC-TPS-70,
International Ground Water Modeling Center. Colorado School of Mines, Golden, CO. 202 pp.
Stadler, W., 1988. Fundamentals of multicriteria optimization. In: Stadler, W. (Ed.), Multicriteria optimization in
Engineering and the Sciences. Plenum Press, New York, pp. 1–25.
Steefel, C.I., 2000. New directions in hydrogeochemical transport modeling: incorporating multiple kinetic andequilibrium pathways. In: Bentley, L.R., Sykes, J.F., Brebbia, C.A., Gray, W.G., Pinder, G.F. (Eds.), Computational
Methods in Water Resources XIII. A.A. Balkema, Rotterdam, pp. 331–338.
Steefel, C.I., MacQuarrie, K.T.B., 1996. Approaches to modeling of reactive transport in porous media. In: Lichtner, P.C.,
Steefel, C.I., Oelkers, E.H. (Eds.), Reactive transport in porous media. Reviews in Mineralogy, vol. 34. MineralogicalSociety of America, Washington, DC, pp. 83–129.
Steefel, C.I., Yabusaki, S.B., 1996. OS3D/GIMRT, Software for multicomponent–multidimensional reactive transport,
user manual and programmer's guide. PNL-11166, Pacific North West Laboratory, Richland, Washington.
Van Dam, J.C., Feddes, R.A., 2000. Numerical simulation of infiltration, evaporation and shallow groundwater levels with
the Richards' equation. J. Hydrol. 233, 72–85.
van der Lee, J., De Windt, L., 2001. Present state and future directions of modeling of geochemistry in hydrogeologicalsystems. J. Contam. Hydrol. 47, 265–282.
van der Lee, J., De Windt, L., Lagneau, V., Goblet, P., 2003. Module-oriented modeling of reactive transport with HYTEC.
Comput. Geosci. 29, 265–275.

218

D. Jacques et al. / Journal of Contaminant Hydrology 88 (2006) 197–218

Vanclooster, M., Viaene, P., Diels, J., Christiaens, K., 1994. WAVE: a mathematical model for simulating water andagrochemicals in the soil and vadose environment, reference and user's manual (release 2.0). Inst. For Land and Water
Manage. Katholieke Univ. Leuven, Leuven, Belgium.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Sci. Soc. Am. J. 44, 892–898.
Walter, A.L., Frind, E.O., Blowes, D.W., Ptacek, C.J., Molson, J.W., 1994. Modeling of multicomponent reactive transport
in groundwater. 1. Model development and evaluation. Water Resour. Res. 30, 3137–3148.
Yeh, G.-T., Cheng, H.-P., 1999. 3DHYDROGEOCHEM: A 3-dimensional model of density-dependent subsurface flow
and thermal multispecies–multicomponent HYDROGEOCHEMical transport. EPA/600/R-98/159, 150 pp.