Journal of Contaminant Hydrology 142–143 (2012) 109–125

Contents lists available at SciVerse ScienceDirect

Journal of Contaminant Hydrology
journal homepage: www.elsevier.com/locate/jconhyd

Inverse optimization of hydraulic, solute transport, and cation exchange
parameters using HP1 and UCODE to simulate cation exchange
Diederik Jacques a,⁎, Chris Smith b, Jiří Šimůnek c, David Smiles babcInstitute for Environment, Health, and Safety, Belgian Nuclear Research Centre (SCK•CEN), Mol, Belgium
CSIRO Land and Water, Canberra, ACT, Australia
University of California Riverside, Riverside, CA, USAa r t i c l e

i n f oArticle history:
Received 21 November 2011
Received in revised form 15 March 2012
Accepted 30 March 2012
Available online 5 April 2012
Keywords:
Fertilization
Coupled reactive transport
Unsaturated water flow
Calibrationa b s t r a c t
Reactive transport modeling is a powerful tool to evaluate systems with complex geochemical
relations. However, parameters are not always directly measurable. This study represents oneof the first attempts to obtain hydrologic, transport and geochemical parameters from anexperimental dataset involving transient unsaturated water flow and solute transport, using an
automatic inverse optimization (or calibration) algorithm. The data come from previously
published, controlled laboratory experiments on the transport of major cations (Na, K, Mg, Ca)during water absorption into horizontal soil columns that were terminated at different times.
Experimental data consisted of the depth profiles of water contents (θ), Cl concentrations, andtotal aqueous and sorbed concentrations of major cations. The dataset was used to optimize
several parameters using the reactive transport model, HP1 and the generic optimization code,
UCODE. Although the soil hydraulic and solute transport parameters were also optimized, thestudy focused mainly on the geochemical parameters because the soil columns wereconstructed from disturbed soil. The cation exchange capacity and the cation exchange
coefficients for two exchange models (Gapon and Rothmund–Kornfeld) were optimized. The
results suggest that both calibrated models satisfactorily described the experimental data,
although the Rothmund–Kornfeld model fit was slightly better. However, information content
and surface response analyses indicated that parameters of the Gapon model are wellidentifiable, whereas those of the Rothmund–Kornfeld model were strongly correlated. The
calibrated geochemical parameters were validated using an independent dataset. In agreementwith the identifiability analysis, the Gapon approach was better than the Rothmund–Kornfeld
model at calculating the observed concentrations of major cations in the soil solution and on theexchange sites.
© 2012 Elsevier B.V. All rights reserved.

1. Introduction
Various human activities, such as irrigation of agricultural lands, may alter the geochemistry of soil systems(e.g. Gonçalves et al., 2006; Ramos et al., 2011). In this study,
we evaluate data collected for the analysis of the wet flush
system for intensive piggeries. Such systems produce large⁎ Corresponding author.
E-mail address: DJACQUES@SCKCEN.BE (D. Jacques).
0169-7722/$ – see front matter © 2012 Elsevier B.V. All rights reserved.
doi:10.1016/j.jconhyd.2012.03.008

volumes of effluent from the urine, feces and waste feed. In
Australia, irrigation with treated effluent is the most productive way to utilize this waste product. In addition to highnutrient content, the effluents also contain high concentrations of K and Na (Smiles and Smith, 2004a,b). Potassium is
becoming the dominant monovalent ion released into theenvironment, as piggery operators, concerned about theeffects of Na on soil properties, have replaced Na with K
salts in the feed. The use of the K-rich effluent for irrigation
has increased K concentrations in the soil solution and on the
cation exchange sites (Smiles and Smith, 2004a).

110

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

In principle, the fate of cations in solution and sorbed bysoils under natural conditions can be predicted using models
that simulate variably-saturated water flow and solutetransport, while considering reaction chemistry and/or chromatographic theory. Examples are LEACHM (Hutson andWagenet, 1992) and/or HP1 (Jacques et al., 2006), which
couples HYDRUS-1D (Šimůnek et al., 2008a) with PHREEQC(Parkhurst and Appelo, 1999). However, the use of suchcomplex numerical models requires specification of a large
number of soil hydraulic, solute transport, and biogeochemical reaction input parameters. In the application involving
irrigation with piggery effluent, cation exchange parameters
are of particular interest. Several approaches exist to obtainthese parameters. For example, the competitive exchange ofNa, K, and Ca during absorption by initially homoionic soils
was studied by Bond and Phillips (1990a,b), Bond (1997), andBond and Verburg (1997). In other studies, cation exchange
parameters were estimated or calibrated using detailed information on binary sorption isotherms and then applied toternary (or multication) transport experiments (e.g., Voegelin
et al., 2000, 2001; Vulava et al., 2000). Alternately, Appelo andPostma (2005; pages 252–260) outlined an approach inwhich cation exchange can be modeled without the need ofprior knowledge of the adsorption isotherms for each pair ofcompeting cations in the presence of all the others.
In this study, we propose to obtain flow, transport, andgeochemical parameters by calibrating an appropriate model
against experimental data involving reactive transport underunsaturated transient conditions in natural soils. Suitabledatasets for such analysis are rarely available (e.g., Mansell etal., 1993). The experimental data of Smiles and Smith(2004b), involving water absorption with piggery effluent
into initially unsaturated soil columns, provided detailed
profiles of water contents and aqueous and sorbed concentrations. This dataset gives us the opportunity to investigate
the use of inverse optimization to estimate geochemicalparameters.
The use of parameter estimation techniques is wellestablished for determining soil hydraulic properties (e.g.,
Hopmans et al., 2002; Vrugt et al., 2008) or basic solutetransport and reaction parameters (e.g., Šimůnek et al., 2002,
2008b; and references given there). This approach has been
widely used for various laboratory and field experiments.
Estimation of soil hydraulic properties typically requires
transient, variably-saturated flow experiments and numericalmodels such as HYDRUS-1D. On the other hand, solutetransport parameters are often derived from column experiments under steady-state water flow conditions (e.g., NkediKizza et al., 1984). Analytical solutions of the transport
equations are then fitted to experimental breakthrough curves
by using parameter estimation codes such as CXTFIT (Toride etal., 1995) or STANMOD (Šimůnek et al., 2008b). Obtaining
solute transport parameters for conditions for which no
analytical solutions exist, such as for transient flow conditions
and/or nonlinear biogeochemical reactions such as adsorption
or cation exchange, can be accomplished only using numericalsolutions (Dai and Samper, 2004; Šimůnek et al., 2002, 2008b).
The current study represents one of the first attempts toobtain soil hydrologic, transport, and geochemical parameters with an automatic inverse optimization algorithm andexperimental data obtained under transient water flow

conditions. This is accomplished using the universal optimization code UCODE (Poeter et al., 2005) with the numericalmodel HP1. The previously acquired experimental dataset ofSmiles and Smith (2004b) is used to calibrate: i) the soil
hydraulic parameters of the van Genuchten–Mualem model(van Genuchten, 1980); ii) the solute transport parameters(dispersivity); and, most importantly, iii) the geochemicalparameters (major cation exchange coefficients). The identifiability of the cation exchange parameters is investigated
for the presented experimental setup by analyzing parameter
sensitivities and response surfaces.
2. Material and methods
2.1. Experimental setup
Two laboratory experimental datasets were used in thisstudy. In the first set of experiments, soil columns wereplaced in contact with an artificial piggery effluent (Smiles
and Smith, 2004b). Because of the destructive nature ofutilized methods, three soil columns were prepared tomeasure concentration profiles at three different times.
Data collected at these three columns were used to calibrate
flow, transport, and exchange parameters. The second
dataset was obtained from a series of experiments, in which
two soil columns were placed in contact with gypsumsaturated solutions (Smiles and Smith, 2008). This dataset
was used to validate optimized parameters.
Details of the method used to produce the first dataset are
given in Smiles and Smith (2004b). In brief, the topsoil(0–100 mm) from a Deep Red Chromosol [site 5 of Smiles
and Smith, 2004a] was used. The soil had not previously been
irrigated with piggery effluent. The soil contains 15–20% clay(by mass) and has a cation exchange capacity (CEC) of
46 mmolc/kg. This is equal to a CEC value of 0.077 ±0.002molc/dm 3, providing the average soil bulk density was
1.68 ± 0.05 kg/dm 3. More information on the soil is given in
Smiles and Smith (2004a).
In the first set of laboratory experiments, sieved soil
(b2-mm) with a water content (θi) of ≈0.05 g/g was packed
into acrylic columns with a 20-mm internal diameter and20 cm long. The columns were made from sections 4, 6, and
10 mm long with the shortest sections close to the inflow
end. The soil was added in increments of 2–3 g and packed
using a small drop hammer to ensure uniformity. The
resulting soil bulk density was 1.68 ± 0.05 kg/dm 3 throughout the column. A solution with a cation composition similar
to piggery effluent, but with Cl − as the sole anion, wasapplied at zero water potential to one end of the horizontal
column. The artificial effluent had the following composition:
Na=10 mmolc/L; K = 20.1 mmolc/L; Ca= 4.69 mmolc/L, and
Cl=35 mmolc/L (Smiles and Smith, 2004b). The pH of theartificial effluent was 7.7. Initial pore water concentrationsin the columns were: Na = 5 mmolc/L, K= 2 mmolc/L,
Ca= 40 mmolc/L, Mg =8 mmolc/L, and Cl =4 mmolc/L.
The experiments were terminated by sectioning thecolumns after 36, 106, and 144 min. Each moist section was
weighed and a water sample was removed by centrifuging in
the presence of 5 cm 3 of 1,1,2-trichloro-1,2,2-trifluoro ethanethat had a specific gravity of 1.57 (Phillips and Bond, 1989).
The soil sample was washed with glycerated alcohol toD. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

remove excess water-soluble cations, and then exchangeable
cations were extracted with 1 M NH4Cl. After washing withalcohol, the CEC was estimated by measuring the remaining
NH4+ (Rayment and Higginson, 1992).
Water-soluble cations were measured by inductively
coupled plasma spectroscopy (ICP) of the water sample
after diluting a small volume of soil solution (0.2–0.5 cm 3) to
10 cm 3. Exchangeable cations were measured in the NH4Cl
extract by atomic absorption spectroscopy (AAS). The NH4+
representing the exchange capacity was measured using an
autoanalyser (Rayment and Higginson, 1992). Dilutions weretracked gravimetrically and all quantities were referred tothe oven-dry weight of soil in each section, which wasmeasured by oven drying following all extractions.
The second dataset is from Smiles and Smith (2008).
Briefly, the acrylic columns used in the first series ofexperiment were packed with potassic soil to a bulk densityof 1.76 ± 0.2 kg/dm 3. This soil was the same as the one used
to develop the optimized geochemical parameters, exceptthat it had been irrigated with K-rich piggery effluent andthus had a high exchangeable K ratio. A gypsum-saturated
solution was then applied at one end of the horizontal
columns from a Mariotte bottle set at a small suction(~5 mm), and the experiments were terminated after 200
and 400 min. The gypsum solution consisted of 24 mmolc/Lof Ca and 24 mmolc/L of S(6). The initial soil solution
concentrations were: 28 mmolc/L of Na, 30 mmolc/L of K,
40 mmolc/L of Ca, 9 mmolc/L, and 13 mmolc/L of Cl. The
experimental methods were the same as in the first set ofexperiments.

2.2. Model description
Modeling the transport of major cations through the soilcores involves modeling the following processes: a) water
flow; b) solute transport (advection and dispersion); andc) geochemical reactions (aqueous complexation and cationexchange).
Water absorption is modeled using a numerical solutionof the Richards equation, which requires knowledge of thewater retention curve (the relation between the watercontent θ [cm 3 cm − 3] and the pressure head h [cm]), andthe unsaturated hydraulic conductivity function (the relation
between the hydraulic conductivity K [cm min − 1] and thepressure head). These nonlinear relations are parameterizedusing the closed-form system of equations proposed by van
Genuchten (1980):

θðhÞ ¼ θr þ

θs −θrð1 þ jαhjn Þm

h

 i
l
1=m m 2:
K ðhÞ ¼ K s Se 1− 1−Se

ð1Þ

ð2Þ

These equations have six adjustable parameters: thesaturated water content θs [cm 3 cm− 3], the residual watercontent θr [cm 3 cm− 3], two shape parameters α [cm − 1] andn [−], the saturated hydraulic conductivity Ks [cm min− 1], anda pore-connectivity and tortuosity factor l [−]. The degree of111

saturation Se and the parameter m are obtained from the otherparameters: Se = (θ − θr)/(θs − θr) and m = 1 − 1/n.
Transport of the aqueous components (that is, the total
aqueous concentration of an element) is described with theadvection–dispersion equation:


∂θC i
∂
∂qC i
w ∂C i
¼
−
þ Ri
θD
∂t
∂x
∂x
∂x

ð3Þ

where t is time [min], x is spatial distance [cm], i (= 1,…,Nm)is the aqueous component number (Nm is the total number ofaqueous components), Ci is the total aqueous concentrationof the ith component [mol cm − 3], q is the volumetric flux
density [cm min − 1], D w is the dispersion coefficient in theliquid phase [cm 2 min − 1], and Ri is the general source/sink
term due to geochemical reactions [cm cm − 3 min − 1]. In
general, this sink/source term (Ri) contains heterogeneous(between different phases) equilibrium reactions and homogeneous (within one phase) and heterogeneous kinetic
reactions. In the current study, only heterogeneous equilibrium reactions are accounted for. The dispersion coefficient isdefined as:w

θD ¼ λjqj þ θD

w0

τw

ð4Þ

w0

where D is the molecular diffusion of aqueous components
in free water [cm 2min − 1], λ is the longitudinal dispersivity[cm], and τw is a tortuosity factor in the liquid phase [−],
calculated using the model of Millington and Quirk (1961):
τw ¼

θ7=3
:
θ2s

ð5Þ

In the multicomponent transport model described above(Eqs. (3)–(5)), only one unknown transport parameter(i.e., the dispersivity, λ) is needed to characterize the flow
field.
Two geochemical processes were considered in thecurrent study: aqueous speciation reactions and cationexchange. The former is described by a set of aqueous
complexation reactions and corresponding mass action
constants. The phreeqc.dat database (Parkhurst and Appelo,1999) is used for the aqueous complexation reactions, themass action constants, and activity correction models.
The relation between the concentrations of ions in thepore water and on the surface exchange sites may bedescribed using the Gapon convention (Appelo and Postma,2005). The equation of the exchange reaction on surface site
X, involving cations N and M, with charges of n and mrespectively, is:mþ

N1=n X þ 1=mMnþ

¼ M1=m X þ 1=nN :ð6Þ

The selectivity coefficient KGMN is given:hinþ 1=nM1=m X Ni
K GMN ¼ h
N1=n X ½M mþ 1=m

ð7Þwhere the square brackets [] represent activities. The activity ofeach exchangeable species is equal to its equivalent fraction.
Here the subscript G refers to the Gapon formulation.

112

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

An alternative exchange convention is the Rothmund–Kornfeld equation (Bloom and Mansell, 2001; Bond, 1995).
For exchange reactions written as:m

nmNXn þ nM ¼ nMXm þ mN

ð8Þ

the selectivity coefficient, KRKMN, is defined as:½MX m n
K RKMN ¼½NX n m

 n m !NMNN
½M m n(Marquardt, 1963) embedded within HYDRUS-1D is used forthis purpose. In the next step, the geochemical parameters are
estimated by HP1 using the model-independent optimization
tool, UCODE (Poeter et al., 2005). The objective function is
obtained as a weighted sum of squared differences (SSD)
between observations and experimental data. Weights aredefined to minimize differences in weighting between different data types (Clausnitzer and Hopmans, 1995):

ð9Þ

where NMN is an empirical exponent for the activity ratio of
the aqueous species. This exponent enables the model todescribe sigmoidal shapes of binary isotherms, which reflect
the sorption of species to strong-binding sites when concentrations of one of the species is low (Appelo and Postma,2005).
2.3. Simulation tool
The HP1 code (Jacques et al., 2006, 2008a,b; Šimůnek etal., 2006, 2009) was used as the simulator to solve thecoupled flow, transport, and geochemical problem. HP1(version 2.3) couples the water flow, multiple solutetransport and heat transport code HYDRUS-1D (Šimůnek etal., 2008a) with the geochemical equilibrium and kinetic
code PHREEQC (version 2.17.5) (Parkhurst and Appelo,1999). The transport simulator solves the water flow andsolute transport equations, whereas the geochemical simulator calculates chemical speciation in the aqueous phase andon the exchange complexes.
The Gapon exchange reactions are defined as half reactions (X − + 1/m M m + = M1/mX with a thermodynamic
constant KGM) relative to a reference half reaction (inPHREEQC; the half reaction with Na has log KGNa = 0). Notethat each Gapon exchange reaction (Eq. (6)) can be written
as a linear combination of the defined half reactions (given in
Appendix 1). As such, the thermodynamic constants of thehalf reactions for the other cations are related to the Gapon
selectivity coefficient as log KGMN = log KGM − logKGN. The
Rothmund–Kornfeld equation is implemented following
example 6.6 in Appelo and Postma (2005) (see also
Appendix 1).
2.4. Inverse optimization
Unknown parameters of the described model are: θs, θr, α,n, Ks, l, λ, KGK, KGCa, KGMg, and CEC for the Gapon approach;and θs, θr, α, n, Ks, l, λ, KRKKNa, KRKCaNa, KRKMgNa, NKNa, NCaNa,
NMgNa and CEC for the Rothmund–Kornfeld approach. These
parameters are estimated using the measured profiles ofwater contents, Cl concentrations, aqueous Na, K, Ca, and Mg
concentrations, and sorbed Na, K, Ca, and Mg concentrations,by minimizing the sum of squared differences between
observations and corresponding simulated values.
A sequential parameter estimation procedure is used here.
HYDRUS-1D alone is used first to estimate the soil hydraulicand solute transport parameters because the flow andtransport simulator of HP1 is essentially identical to the
HYDRUS-1D code. The inverse optimization algorithm basedon the Levenberg–Marquardt nonlinear minimization method

wj ¼

1
nj σ jð10Þ

where wj is the weight for the jth data point, nj is the number ofobservations for a given data type and σj is the measurement
variance within a given data type. Twenty four different datatypes are used when optimizing cation exchange parameters:aqueous and sorbed concentration profiles of four cations atthree times.
The initial and boundary conditions used in the modelingexercise were as follows: the initial volumetric water content
was 0.084 cm3 cm− 3, the volumetric water content at the inletof the column was set to 0.305 cm3 cm− 3 (inferred from theexperimental data close to the inlet end of the horizontal
column), and the outlet boundary was set to free drainage. The
initial estimate of CEC of the soil was 0.085 mol/dm 3. In thisparticular example, the initial estimate of CEC is not criticalsince this parameter is well-defined, as discussed below in thesection “Response surface analysis of the objective function”.
Cation exchange capacity was optimized since the bulk density(and thus CEC) varied slightly within and between soilcolumns. One additional optimization was also carried out
without fitting CEC to evaluate if the use of the measured CEC
has an effect on fitted selectivity coefficients. The initial
estimates of the Gapon selective coefficients for Ca/Na, Ca/K,
and Ca/Mg exchange were KGCaNa = 2.9, KGCaK = 0.2, andKGCaMg = 1.2 respectively, as the middle values of the ranges
given in MDH (2003). Because of the lack of additional
information, the same initial values were also used for the
Rothmund–Kornfeld selectivity coefficients and n was initially
set at 0.7, 0.9 and 0.9 for nKNa, nCaNa, and nMgNa, respectively.
It is important to emphasize that the hydraulic and solutetransport parameters are fitted only to reproduce the waterflow and solute transport conditions so that the reactive
constants could be determined from these transient experiments. Because the experiments were done for disturbed
soils, these parameters may not be representative of
undisturbed soil conditions. As such, uniqueness and correlation between the hydraulic and transport parameters (seeparameter identifiability) are not further assessed here. The
reader is referred to other studies on identifiability of
hydraulic and transport parameters (e.g., Šimůnek and van
Genuchten, 1996; Toorman et al., 1992).
2.5. Parameter identiﬁability
Identification of parameters using an inverse optimization
approach requires that there is enough information in theexperimental data to obtain unique and/or well-definedparameters. Two types of analysis were performed in order
to investigate the identifiability of the cation exchange
parameters in this particular experimental set up. The firstD. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

113

analysis uses dimensionless and composite scaled sensitivities. The dimensionless scaled sensitivity of the ith observation and jth parameter, dssji, is defined as (Hill and Tiedeman,2007):
dssji ¼!
∂yj 
1=2
 jb jw∂bi b i j

ð11Þ

where yj is the simulated value corresponding to a given
observation, bi is the ith estimated parameter, b is the vector,
which contains the parameter values at which the sensitivities are evaluated, and wj is the weight of the jth observation
point as defined above. The sensitivity is the derivative of thesimulated value with respect to the parameter bi. It provides
information on the relative importance of an observation to
the estimation of a single parameter, compared to the otherobservations. Larger values of dssji indicate observations
containing more information. Alternatively, the dimensionless scaled sensitivity also provides a measure that compares
the relative importance of a given parameter for a given
simulated value. The information content of different sets ofdata types is quantified by taking the sum of the absolutevalues of the dimensionless scaled sensitivities. The averagevalue is used to compare the different data types (e.g., profiles
at a given time of the aqueous or sorbed concentrations, orthe sum of the three profiles for a given aqueous or sorbed
species).
The composite-scaled sensitivity for the ith parameter(cssi) is calculated as (Hill and Tiedeman, 2007):
cssi ¼"
N obs 
X
j¼1

dssji


#1=2
2 

b =N obs


ð12Þ

where Nobs is the number of observations.
The composite-scaled sensitivity provides a means toevaluate the amount of information in all experimental data
to estimate a given parameter. For example, Friedel (2005)
interpreted css for different combinations of data types in acoupled water–heat–solute transport model.
The second analysis is based on the evaluation of response
surfaces (e.g., Šimůnek and van Genuchten, 1996; Toorman
et al., 1992). A response surface is a two-dimensional plot of
the objective function obtained by varying two parameters,
usually on a regular grid (80 × 80 in this study), while
keeping the other parameters at their fitted values. Since
only the cation exchange parameters are included in thisanalysis, there are 6 contour plots for the Gapon approachand 21 plots for the Rothmund–Kornfeld approach.
3. Result and discussion
Experimental data from Smiles and Smith (2004b) arepresented in Fig. 1 (water contents and chloride concentrations) and in Fig. 2 (aqueous and sorbed cation concentrations). Note that due to the large scatter in the 36-minute K
and Cl concentrations (attributed to the contamination of soil
samples), these data were discarded from the fitting
procedure. The data are plotted as a function of the
Boltzmann variable (distance/time 0.5) because, for the given
boundary conditions, water content and ion concentration

Fig. 1. Simulated [lines] profiles of water contents (a) and chloride
concentrations (b) compared to experimental data [symbols] from water
absorption and solute transport experiments with piggery effluent (Smiles
and Smith, 2004b).

profiles should coalesce to a single curve when graphed interms of the Boltzmann variable (Smiles and Smith, 2004b).
For a full interpretation of the experimental dataset, thereader is referred to that study. Only the key findings arelisted below in order to facilitate interpretation of modeled
experimental results. The piston front (the location in thecolumn where the infiltrating solution completely displaces
the initial water) is between 0.75 and 0.85 cm/min 0.5 (seeSmiles and Smith, 2004b). Non-reactive solutes (such as Cl)
should have the center of their diffusive front around thislocation (Fig. 1b). A change in aqueous concentrations of thereactive cations (Mg and Ca in Fig. 2c and d) around thepiston front occurs without any noticeable change in sorbed
concentrations. This is because of the concentration differences between the initial solution (the total sum of cations isabout 35 mmolc/L) and the infiltrating solution (the totalsum of cations is about 60 mmolc/L). Across the piston front,
the cation ratios are buffered by the exchange processes sothat the total concentration adjusts while the sorbed andaqueous cation ratios appear to stay constant. Concentration
changes around the piston front are thus associated with the

114

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125
Table 1
Estimated parameter values and linear 95% confidence intervals (CI) for the
Gapon approach.
Parameter

Upper
95% CI

Lower
95% CI

Standard
deviation

Optimization with estimated CEC
Log KGK
0.91
Log KGCa
0.56
Log KGMg
0.16
CEC (mol/dm³)
0.072

0.96
0.59
0.19
0.073

0.85
0.53
0.12
0.071

0.027
0.015
0.018
5.4 × 10− 4

Optimization with fixed CECa
Log KGK
0.99
Log KGCa
0.56
Log KGMg
0.19

1.05
0.60
0.22

0.93
0.53
0.15

0.029
0.016
0.018

a

Final
value

115

Table 2Estimated parameter values and linear 95% confidence intervals (CI) for the
Rothmund–Kornfeld approach.
Parameter

Final
value

Upper
95% CI

Lower
95% CI

Standard
deviation

Log KRKKNa
Log KRKCaNa
Log KRKMgNa
NKNa
NCaNa
NMgNa
CEC (mol/dm³)

1.14
1.85
1.27
0.307
0.717
0.843
0.079

1.18
1.98
1.42
0.383
0.784
0.929
0.081

1.10
1.72
1.12
0.230
0.647
0.756
0.079

0.010
0.069
0.076
0.039
0.035
0.044
7.9 × 10− 4

The value of CEC was fixed at 0.077 mol/dm3.

displacement of the initial solution. In the region where theBoltzmann variable is smaller than 0.25 cm/min 0.5, typicalcation exchange reactions occur, during which Na and K
replaces Ca and Mg. Since all experimental curves coalesce inthis region after applying the Boltzmann transformation, theexchange processes are in chemical equilibrium (i.e., timescales for chemical reactions are short relative to time-scalesof solute movement).
3.1. Optimization results
Fig. 1 compares the water content (θ) and chloride
concentration profiles measured during the experimentswith the optimized model predictions. Both the Richards
equation and the advection–dispersion equation were able toreproduce very well the experimental water content andchloride profiles, respectively. Note that the simulated Cl
profiles are a slightly ahead of the experimental data aroundthe piston front. Fitted parameters are θs = 0.307 cm 3 cm − 3,
α = 0.230 cm − 1, n = 1.45, Ks = 0.157 cm min − 1, l = 0.99, and
λ = 0.0694 cm.
The modeled profiles of the major cations in the aqueousand sorbed phases for the three times are compared to themeasured data in Fig. 2. In general, both approaches are able
to describe the experimental data reasonably well, althougheach model missed some particularities of the experimental
data.
Fitted parameters are given in Tables 1 and 2 for the
Gapon and the Rothmund–Kornfeld approaches, respectively.
The weighted sum of squared differences (SSD) is about 20%
lower for the Rothmund–Kornfeld approach than for the
Gapon approach. The lower SSD for the Rothmund–Kornfeld
approach is to be expected, because of the higher number ofadjustable parameters in the geochemical model (7 versus 4
in the Gapon approach). Nevertheless, both the Akaike's
information criterion (AIC) and the Bayesian information
criterion (BIC) (Hill and Tiedeman, 2007) are also lower for
the Rothmund–Kornfeld approach (− 959 and −926 for AIC
and BIC, respectively) than for the Gapon approach (−875

and −855 for AIC and BIC, respectively). The difference
between the two approaches is sufficiently large enough(Burnhum and Anderson, 2004) to give preference to theRothmund–Kornfeld model. Based on the SSD, AIC, and BIC,
the Rothmund–Kornfeld model fits the experimental data
better than the Gapon model.
The Gapon approach was additionally used to optimize
the exchange constants while setting CEC equal to themeasured value (0.077 mol/dm 3). Fitted values are again
given in Table 1. The aqueous concentrations were almost
identical to those obtained when the CEC was optimized. The
largest differences were in the initial sorbed Ca and Mg
concentrations, which were about 6.5 and 13% higher,
respectively, when the CEC was not optimized. However,the 95% confidence intervals of the three exchange parameters obtained using the two approaches (optimized andfixed CEC) overlap.
The fitted CEC was 0.072 and 0.079 mol/dm 3 for the
Gapon and the Rothmund–Kornfeld approaches, respectively.
These values are, respectively, slightly lower and greater thanthe measured CEC of 0.077 ± 0.002 mol/dm 3. The linear 95%
confidence interval (CI) of the CEC is small, which indicates awell-defined value. Note that the measured CEC falls within
the 95% confidence interval of the CEC estimated using theRothmund–Kornfeld approach, but not the Gapon approach.
Correlation coefficients of the CEC with the other parameters
are all less than 0.1 in both approaches (Table 3).
The correlation coefficients for other parameters in the
Gapon approach were also relatively small, with only one
being larger than 0.75 (between log KGCa and log KGMg). Eachcorrelation coefficient between two log KGM was positive(Table 3). The parameters of the Gapon exchange equations,calculated from the fitted half-reaction constants, were:
KGCaNa = 3.55, KGCaK = 0.46, and KGCaMg = 2.51. Uncertaintyof the estimated parameters was slightly larger in theRothmund–Kornfeld approach, as evidenced by a greaternumber of correlation coefficients exceeding 0.75 (Table 4).
Further analysis of the identifiability and uniqueness of thefitted parameters is given below using the response surface
analysis of the objective function.
The Gapon approach was able to describe Na in the
aqueous and exchange phase within the scatter of theFig. 2. Simulated profiles of aqueous (a–d) and sorbed (e–h) cation concentrations obtained with the Gapon approach [solid lines] or the Rothmund–Kornfeld
approach [dashed lines] compared to experimental data [symbols] from water absorption and solute transport experiments with piggery effluent (Smiles andSmith, 2004b).

116

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

Table 3
Correlation coefficients between parameters used in the Gapon approach.

Log KGK
Log KGCa
Log KGMg
CEC (mol/dm³)

Log KGK

Log KGCa

Log KGMg

CEC (mol/dm³)

1
0.383
0.327
0.206

1
0.811
6.41 × 10− 2

1
0.241

1

experimental data. However, the exchangeable K (sorbed)
concentration was under-predicted. The description of
exchangeable Ca was generally good, with the exception of
the final predicted exchangeable Ca when the Boltzmann
variable is less than 0.2 cm min 0.5. The description of
exchangeable Mg was generally within the variability in theexperimental data, although close inspection of the aqueous
Mg data suggests that the simulated concentration front isslightly ahead of the measured data between the piston front
and the region of cation exchange (i.e., between 0.25 and0.75 cm min − 0.5; Fig. 2d). The Gapon approach seemed to beunable to describe sorbed K concentrations, especially when
the K concentrations are low and the concentrations of thebivalent cations (in particular Ca) are high (initial condition).
Thus, at low K concentrations, sorption is apparently stronger
than can be described by the Gapon approach. Although notinvestigated in detail here, a possible reason for high K sorption
at low concentrations may have been the presence of illite in
the soil sample (20% of the clay, Smiles and Smith, 2004b). Illite
typically has some high-preference sorption sites for K, that are
related to the formation of K inner-sphere complexes between
the clay surfaces (Evangelou and Philipps, 2005).
As suggested earlier by the weighted sum of squared
differences, the Akaike's information criterion, and theBayesian information criterion, the Rothmund–Kornfeld
approach appears to describe the experimental data better
than the Gapon approach. While the overall description of
the cations' behavior was similar between the Rothmund–Kornfeld and Gapon conventions, the simulated initial andfinal exchangeable K and Ca more closely matched measured
values when using the Rothmund–Kornfeld approach.
Another way of representing the data is to show thesorbed concentrations against the aqueous concentrations(Fig. 3). The redistribution of aqueous Ca and Mg without any
visible changes in the cation exchange composition aroundthe piston front is clearly illustrated in this figure (vertical

Table 4
Correlation coefficients between parameters used in the Rothmund–Kornfeld approach.
Log
KGK

Log
KGCa

Log
KGMg

NKNa

NCaNa

NMgNa

CEC

1
Log
KGK
0.268
1
Log
KGCa
0.221
0.722
1
Log
KGMg
− 0.422
0.467
0.354
1
NKNa
NCaNa
− 0.03− 0.926 − 0.676 − 0.556
1
NMgNa − 0.01
− 0.739 − 0.922 − 0.478
0.825
1
CEC
0.580
0.223
0.229 − 0.490 − 0.178 − 0.174 1

Fig. 3. Simulated relations between sorbed and aqueous concentrations
obtained with the Gapon approach [solid lines] or the Rothmund–Kornfeld
approach [dashed lines] compared to experimental data [symbols] fromwater absorption and solute transport experiments with piggery effluent(Smiles and Smith, 2004b).

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

lines in Fig. 3c, d) and nicely represented by both model
approaches. Given the large scatter in the experimental
data, it is impossible to make any conclusions about theadequacy of the two models when comparing their resultswith the Na, Mg, and Ca data. However, the Rothmund–Kornfeld approach performs much better in reproducing
the K isotherm than the Gapon approach. At high sorbed(and aqueous) K concentrations, both approaches perform
equally well. However, at low K concentrations, sorption in
the Rothmund–Kornfeld approach is larger than in the
Gapon approach, indicating some preferential sorption of K
at low concentrations.
3.2. Information content in experimental data
The composite-scaled sensitivity coefficients (css) for the
Gapon parameters are given in Table 5. The experimental datacontain the most information for the estimation of the cationexchange capacity (CEC), although the css for KGCa is also high.
Less information is present in the experimental data to estimate KGK and KGMg. The sorbed data provide more information
than the aqueous data, as indicated by higher values of theaverage dimensionless scaled sensitivity coefficients (dss)(Table 5), except for K. From the four cations, Ca has thelargest average dimensionless scaled sensitivity.
In many cases, identifiability of a particular parameter
benefits from a combination of different data types, because asingle parameter influences multiple output variables. In
addition, the dimensionless scaled sensitivities change withobservations as a function of the Boltzmann variable. Let us
provide interpretation of dss for KGCa, which has high values
for aqueous Na and K concentrations around their moving
fronts and for sorbed Ca concentrations (Fig. 4). While the Na
front with increasing aqueous Na concentrations passes closebehind the piston front, the K front, where K replaces Ca andMg on the exchange sites, is nearer to the inlet boundary.
Higher KGCa will result in less K on the exchange complex,
which explains the high dss for aqueous K between 0.25 and0.75 cm/min 0.5. Similarly, a higher KGCa will result in more Na
in the aqueous phase, which explains the high dss foraqueous Na concentrations. Also, the value of KGCa defines
the amount of the sorbed Ca, which explains the high css for117

the sorbed Ca for the whole profile. All other dss profiles canbe interpreted in a similar way.
This information content is similar for all aqueous data.
Aqueous concentrations do not contribute to the information
content ahead of the piston flow (i.e., for the Boltzmann
variable larger than 0.75 cm/min 0.5) where dss is close to
zero. Aqueous Na concentrations provide useful information
between 0.25 and 0.75 cm/min 0.5 where the Na front islocated. The three other cations provide the most information
around the K front where sorbed K replaces sorbed Ca andMg (between 0 and 0.25 cm/min 0.5). On the other hand, theinformation content of the sorbed concentrations is quite
independent of the Boltzmann variable, at least when it islarger than 0.25 cm/min 0.5, since in this region there are nochanges in sorbed concentrations of K, Ca, and Mg. However,
the initial sorbed concentrations (especially Ca and Mg) haveimportant information for several parameters because they
constrain the initial distribution of cations between aqueousand sorbed concentrations. Additional information is contained in sorbed concentrations at the inlet where theexchange between K, Ca, and Mg is complete (i.e., at small
Boltzmann variables). Additionally, the concentration profiles measured at 144 min provide more information than theother two datasets since they contain more data points withhigher sensitivity.
A similar analysis was also done for the Rothmund–Kornfeld approach (Table 6 and Fig. 5). The css values for the
Rothmund–Kornfeld approach are larger than those for the
Gapon approach. In fact, css of the parameters involving K
exchange are smaller compared to those of the otherparameters. The profiles of the dimensionless scaled sensitivities of the observations for KRKKNa, KRKCaNa, KRKMgNa, andthe CEC (results not shown) are very similar to those for the
Gapon approach, except that they have higher absolutevalues. For the N coefficients in the Rothmund–Kornfeld approach, sorbed Ca and Mg concentrations are the most informative experimental data ahead of the piston front. For
Boltzmann variables smaller than 0.25 cm/min 0.5, the aqueous Ca and Mg concentrations provide similar amount of
information compared to the sorbed Ca and Mg concentrations. As expected from the css of NKNa, the dimensionless
scaled sensitivity for this parameter is small.

Table 5Composite-scaled sensitivities (css, top value in the second column), percentage of css compared to maximum values (bottom value in the second column) andthe averaged dimensionless scaled sensitivities (dss) for each measured cation for both the aqueous data (top value of columns 3–6) and sorbed data (bottom
value of columns 3–6) for the Gapon approach.
Parameter

css

Na

K

Ca

Mg

Top value(Absolute value)(Aqueous)(Aqueous)(Aqueous)(Aqueous)

Bottom value(% of max)

(Sorbed)(Sorbed)(Sorbed)(Sorbed)Log KGK

0.65
0.27
1.11
0.46
0.25
0.10
2.4
1.00

0.222
0.300
0.357
0.709
0.023
0.124
0.374
0.768

0.504
0.425
0.266
0.213
0.011
0.012
0.423
0.289

0.241
0.679
0.119
1.822
0.030
0.387
0.226
5.916

0.095
0.187
0.222
1.814
0.039
0.538
0.197
1.870

Log KGca
Log KGMg
CEC (mol/dm³)

Bold values indicate the experimental dataset which provides the largest information for a given parameter.

118

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

Fig. 4. Profiles of dimensionless-scaled sensitivities (dss) for the parameters in the Gapon approach for the aqueous (left) and sorbed (right) data.

The above discussed information content can be used
when designing the experimental set up and procedures. The
highest information content is in, and thus it should bedesirable to measure, initial sorbed concentrations andaqueous concentrations of different cations at different
critical positions behind the piston front (i.e., around the Na

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

119

Table 6Composite-scaled sensitivities (css, top value in the second column), percentage of css compared to maximum values (bottom value in the second column) andthe averaged dimensionless scaled sensitivities for each measured cation for both the aqueous data (top value of columns 3–6) and sorbed data (bottom value ofcolumns 3–6) for the Rothmund–Kornfeld approach.
Parameter

css

Na

K

Ca

Mg

Top value(Absolute value)(Aqueous)(Aqueous)(Aqueous)(Aqueous)

Bottom value(% of max)

(Sorbed)(Sorbed)(Sorbed)(Sorbed)Log KRKKNa

1.44
0.41
3.48
1.00
2.04
0.58
0.22
0.06
3.46
0.99
3.08
0.88
2.47
0.71

0.315
0.403
0.482
1.251
0.061
0.282
0.045
0.172
0.216
1.094
0.072
0.318
0.320
0.863

0.692
1.449
0.637
0.999
0.051
0.122
0.200
0.177
0.301
0.849
0.127
0.202
0.486
0.662

0.396
2.273
0.269
6.386
0.174
3.064
0.094
0.285
0.470
6.414
0.545
4.514
0.247
5.900

0.173
0.722
0.600
5.838
0.348
4.349
0.052
0.093
1.128
5.925
1.201
6.393
0.219
1.879

Log KRKCaNa
Log KRKMgNa
NKNa
NCaNa
NMgNa
CEC (mol/dm³)

Bold values indicate the experimental dataset which provides the largest information for a given parameter.

front and around the K front when K displaces Ca and Mg on
the exchange sites). In general, the measurement of a single
profile at a large time should be sufficient because suchmeasurement provides more informative data points.
3.3. Response surface analysis of the objective function
Two-dimensional contour plots of the objective function
can be used to investigate the identifiability of optimized
parameters and their mutual correlations. In the soilhydrological literature, contour plots are often used toinvestigate the identifiability, uncertainty, and correlation ofsoil hydraulic parameters in an inverse optimization problem(e.g., Šimůnek and van Genuchten, 1996; Toorman et al.,1992). Such analysis often indicated that many soil hydraulicparameters are not well defined or strongly correlated. For
example, for the one-step outflow experiment, Toorman et al.
(1992) identified close correlations between parameters α-nand α-Ks. Similarly, Šimůnek and van Genuchten (1996)
found close correlations between parameters α-n, α-Ks, andn-Ks for the tension disk infiltration experiment. In both ofthese experiments, only water flux was measured, and theparameter identifiability dramatically improved when additional information (such as pressure heads) was provided.
The severity of unidentifiability or correlation thus depends
strongly on the boundary conditions applied and/or themeasured variables (e.g., water contents, water fluxes,pressure heads) used to constrain the inverse optimization.
Fig. 6 shows three selected contour plots (out of six) forthe parameters in the Gapon approach. Well identifiable
parameters show a well-defined single minimum (i.e., asmall circular region). For example, parameter uncertainty
and correlations are small in the log KGK–log KGCa (Fig. 6a)and log KGK–CEC (Fig. 6c) planes. Less-identifiable parameters show contour plots with an extended region around theminimum or multiple minima, indicating that the same valueof the objective function can be obtained for different

combinations of optimized parameters (e.g., the log KGC–log
KGMg contour plot, Fig. 6b). Extended regions of small valuesof the objective function indicate an increased uncertainty of
optimized parameters. Overall, Fig. 6 shows that the experimental set-up used in this study permits the estimation ofwell-defined reaction parameters, which is in contrast to
many experimental setups used for determining soil hydraulic parameters.
The same analysis was also performed also for theparameters of the Rothmund–Kornfeld approach. The analysis
for seven parameters results in 21 contour plots, nine of which
are shown in Fig. 7. Contour plots not involving NKNa have awell-defined minimum, although a strong correlation exists
between two parameters (e.g., between KRKCaNa and KRKMgNa or
between KRKCaNa and NCaNa, as previously discussed). The
contour plots with NKNa show contours parallel with the NKNa
axis, which indicates that NKNa is a rather insensitive parameter(as was also clear from the analysis of css and the dimensionless scaled sensitivity values). Nevertheless, although NKNa
itself is a rather insensitive parameter, there is a significant
improvement in applying the Rothmund–Kornfeld approachrather than the Gapon approach to the K isotherm (Fig. 3).
Although the analysis of the sensitivities (discussed above)
indicated that the experimental data contained significant
information to estimate the Rothmund–Kornfeld parameters,
the two-dimensional contour analysis shows that someparameters are ill-defined.
3.4. Validation of derived parameters
In the present experiment, the main interest was inestimating the cation exchange parameters. Soil hydraulicand solute transport parameters were fitted in order to havea proper description of the water flow and hydrodynamic
dispersion processes, thus removing confounding effects
from the results. The fitted hydraulic and solute transportparameters, which are often affected by the soil structure, are120

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

Fig. 5. Profiles of dimensionless-scaled sensitivities (dss) for selected parameters in the Rothmund–Kornfeld approach for the aqueous (left) and sorbed (right) data.

of minimal relevance to similar but undisturbed field soil
because the laboratory experiments were carried out using ahighly disturbed soil material. In contrast, cation exchange
parameters estimated using the disturbed soil are assumed tobe applicable to other disturbed and undisturbed soils of thesame type and chemical composition, since these parameters
are affected mainly by the soil texture and mineralogy. Thishypothesis was tested by applying the cation exchange
parameters derived from experimental data of Smiles andSmith (2004b) against a new experimental dataset (Smiles
and Smith, 2008) that uses the same soil, with the exceptionthat it had been irrigated with K-rich effluent. In thisexperiment, water equilibrated with gypsum was absorbed
into a soil column that initially had a high exchangeable K
ratio (Smiles and Smith, 2008).
The performance of the model and its optimized parameters were evaluated by comparing calculated and measured

concentration profiles for major cations in the aqueous phaseand on the exchange sites (Fig. 8). The larger differences
between the initial sorbed measured and simulated concentration compared to the differences between the initial
aqueous measured and simulated concentrations are because
the initial sorbed concentrations are calculated with thegeochemical model, whereas the measured initial aqueousconcentrations are direct input in the model. The calculatedaqueous concentrations of Na, K, Ca, and Mg agree well withthe measured data when the Gapon coefficients are used. In
contrast, there were larger variations between the measured
and calculated data using the Rothmund–Kornfeld approach.
In general, the Gapon exchange convention provides better
agreement between the simulated and measured datacompared to the Rothmund–Kornfeld convention. Although
the Rothmund–Kornfeld approach fitted the first set ofexperiments better (considering the Akaike criterion), a

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

121

they were calibrated. The Gapon convention, which uses
fewer parameters, appears to be more robust.
Three different analyses discussed above (i.e., information
content analysis, response surfaces analysis, and validation)
showed that under conditions involving the transport ofmultiple cations (i.e., Na, K, Ca, and Mg), the aqueous andsorbed concentrations are not sensitive to the NKNa parameter of the monovalent exchange reaction. Predicted concentrations are more sensitive to parameters describing theheterovalent exchange. However, PHREEQC calculations ofbinary monovalent exchange between K and Na showed that
the K-isotherm in particular is quite sensitive to NKNa (results
not shown). It may be worthwhile to investigate if theinverse problem for the Rothmund–Kornfeld approach would
be better defined when independent measurements of binary
monovalent exchange for the K-isotherm (covering theentire range of the sorbed mole fraction) are included.
4. Summary and conclusions

Fig. 6. Contour plots of the objective function defined using the aqueous andsorbed concentrations for the Gapon approach. (a) logKGK–log KGCa, (b) log
KGCa–log KGMg, and (c) log KGK–CEC.

larger number of optimized parameters led to higher
uncertainty and lower identifiability, as was shown using
the analysis of the dimensionless scaled sensitivity coefficients and response surfaces. That is, the Rothmund–Kornfeld
model is slightly over–parameterized given the experimental
information content. This means that some optimized
parameters may not have physical meaning and may not be
reliable when used in conditions other than those for which

Previously published laboratory experiments by Smiles
and Smith (2004b) that examined the transport of major ions
during water absorption in horizontal soil columns wereused here to estimate soil hydraulic, solute transport, andsolute reaction parameters. Soil hydraulic parameters were
estimated first using the measured water content profiles atthree different times with the inverse option of HYDRUS-1D.
Solute transport parameters were then similarly identified
using the measured Cl profiles. Finally, the geochemicalparameters of the cation exchange were estimated using
liquid-phase and exchange-phase concentrations of four
major cations (Na, K, Mg, and Ca), the biogeochemical HP1
code, and the universal optimization code, UCODE. The cationexchange parameters were identified for the Gapon andRothmund–Kornfeld models. While both approaches describedthe experimental data reasonably well, the Rothmund–Kornfeld model described the experimental data better thanthe Gapon model during the calibration stage. Both the Akaike's
information criterion and the Bayesian information criterionindicated that, despite having more parameters and consequently providing more flexibility in fitting experimental data
than the Gapon model, the Rothmund–Kornfeld model wassuperior.
The information content and response surface analyses
indicated that while all parameters of the Gapon model were
well identifiable, there was a close correlation between
several parameters of the Rothmund–Kornfeld model (Figs. 6
and 7). Thus, using the Rothmund–Kornfeld convention, theseparameters cannot be identified simultaneously from theexperimental data and some need to be measured, orotherwise estimated, independently.
These findings were confirmed by a subsequent validation
exercise, in which geochemical reaction parameters estimated from the initial set of column experiments were used to
predict liquid-phase and exchange-phase concentrations offour major cations from another set of experiments involving
absorption of gypsum-saturated water into horizontal soilcolumns. While the parameters of the Gapon model allowed
for a good fit between the modeled and empirical aqueous
data, the parameters of the Rothmund–Kornfeld model
resulted in a poorer fit, attributed to the fact that not all

122

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

Fig. 7. Contour plots of the objective function defined using the aqueous and sorbed concentrations for the Rothmund–Kornfeld approach. First column: (top) log
KK–log KCa, (middle) log KCa–log KMg, and (bottom) log KK–CEC; Second column: log KK–nK, log KCa–nCa, and log KK–nCa; Third column: log KK–nCa, log KCa–nK, andlog KK–nMg.

parameters were identifiable, as shown by the analysis of theinformation content and response surfaces.
The combined use of the flow and transport code
HYDRUS-1D, the biogeochemical code HP1, and the universal
optimization code UCODE proved to be a powerful tool for
analyzing experimental data involving transient variablysaturated water flow, solute transport, and geochemical
reactions. A flexible code, such as HP1, allows one to analyze

different geochemical conceptual models (in our case, theGapon and Rothmund–Kornfeld models) using one code. A
key element however, is that experimental data should
contain enough information to estimate all parameters. For
example, the sensitivity analysis indicated that the current
experimental data (Smiles and Smith, 2004b) does notcontain enough information to identify all parameters in theRothmund–Kornfeld model.

Fig. 8. Experimental [symbols] and simulated [solid lines for the Gapon approach and dashed lines for the Rothmund–Kornfeld approach] profiles of water solubleand exchangeable cations during the absorption of gypsum-saturated solution into a soil column.

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

123

124

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125

Appendix 1
1. A section of the HP1 input file where the cation exchange
is defined using the Gapon convention:
EXCHANGE_MASTER_SPECIESG GEXCHANGE_SPECIES
G- = G-; log_k 0
G- + Na + = NaG; log_k 0
G- + K + = KG; log_k 1.16
G- + 0.5 Ca+2 = Ca0.5G; log_k 0.462
G- + 0.5 Mg+2 = Mg0.5G; log_k 0.383.
Where log_k is log K of the reaction (values given here are
the initial estimates). The reference half-reaction for eachexchanger is Na; the half reaction is Na + + X − = NaX andit has a log K of 0.0.
2. A section of the HP1 input file where the cation exchange
is defined using the Rothmund–Kornfeld equations:
EXCHANGE_MASTER_SPECIESRREXCHANGE_SPECIES
R- = R-; log_k 0
Na+ + R- = NaR; log_k 0
0.318 K+ + NaR = KR + .318 Na+; log_k 1.12;-no_check; -mole_balance KR
0.736
Ca+2 + 2NaR = CaR2 + 1.472Na+;log_k1.804; -no_check; -mole_balance CaR2
0.868 Mg+2 + 2NaR = MgR2 + 1.712 Na+;log_k 1.224;-no_check; -mole_balance MgR2.

References
Appelo, C.A.J., Postma, D., 2005. Geochemistry, Groundwater and Pollution,
2nd ed. AA Balkema Publishers, Leiden. The Netherlands a member of
Taylor & Francis Group plc.
Bloom, S.A., Mansell, R.S., 2001. An algorithm for generating cation exchange
isotherms from binary selectivity coefficients. Soil Science Society ofAmerica Journal 65, 1426–1429.
Bond, W.J., 1995. On the Rothmund–Kornfeld description of cation exchange.
Soil Science Society of America Journal 59, 436–443.
Bond, W.J., 1997. Competitive exchange of K+, Na+ and Ca2 + duringtransport through soil. Australian Journal of Soil Research 35, 739–752.
Bond, W.J., Phillips, I.R., 1990a. Cation exchange isotherms obtained withbatch and miscible-displacement techniques. Soil Science Society ofAmerica Journal 54, 722–728.
Bond, W.J., Phillips, I.R., 1990b. Approximate solutions for cation transport
during unsteady, unsaturated soil water flow. Water Resources Research
26, 2195–2205.
Bond, W.J., Verburg, K., 1997. Comparison of methods for predicting ternaryexchange from binary isotherms. Soil Science Society of America Journal
61, 444–455.
Burnhum, K.P., Anderson, D.R., 2004. Multimodal inference: understanding
AIC and BIC in model selection. Sociological Methods & Research 33,
261–304.
Clausnitzer, V., Hopmans, J.W., 1995. Non-linear parameter estimation:
LM_OPT. General purpose optimization code based on the Levenberg–
Marquardt algorithm. Land, Air and Water Resources Paper No. 100032.
University of California Davis, CA.
Dai, Z., Samper, J., 2004. Inverse problem of multicomponent reactive
chemical transport in porous media: formulation and applications.
Water Resources Research 40, W07407. doi:10.1029/2004WR003248.
Evangelou, V.P., Philipps, R.E., 2005. Cation exchange in soils. In: Tabatabai,
M.A., Sparks, D.L. (Eds.), Chemical Processes in Soil. Soil Science Societyof America, pp. 343–410.

Friedel, M.J., 2005. Coupled inverse modeling of vadose zone water, heat, andsolute transport: calibration constraints, parameter nonuniqueness, andpredictive uncertainty. Journal of Hydrology 312, 148–175.
Gonçalves, M.C., Šimůnek, J., Ramos, T.B., Martins, J.C., Neves, M.J., Pires, F.P.,
2006. Multicomponent solute transport in soil lysimeters irrigated withwaters of different quality. Water Resources Research 42, W08401.
doi:10.1029/2006WR004802, 17 pp.
Hill, M.C., Tiedeman, C.R., 2007. Effective Groundwater Model Calibration
with Analysis of Data, Sensitivities, Predictions and Uncertainty. John
Wiley & Sons, NJ, USA.
Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2002. Inverse modelingof transient water flow, In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil
Analysis, Part 1, Physical Methods, Chapter 3.6.2, Third ed. SSSA,
Madison, WI, pp. 963–1008.
Hutson, J.L., Wagenet, R.J., 1992. LEACHM. A process based model of water andsolute movement, transformations, plant uptake and chemical reactions in
the unsaturated zone. Version 3. Cornell University Department of Soil,
Crop and Atmospheric Sciences Research Series No. 92-3.
Jacques, D., Šimůnek, J., van Genuchten, M.Th., Mallants, D., 2006. Operatorsplitting errors in coupled reactive transport codes for transient variably
saturated flow and contaminant transport in layered soil profiles.
Journal of Contaminant Hydrology 88, 197–218.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.Th., 2008a. Couplinghydrological and chemical processes in the vadose zone: a case study onlong term uranium migration following mineral P-fertilization. Vadose
Zone Journal 7, 698–711.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.Th., 2008b. Modelling
coupled water flow, solute transport and geochemical reactions affecting
heavy metal migration in a podzol soil. Geoderma 145, 449–461.
Mansell, R.S., Bloom, S.A., Bond, W.J., 1993. Simulating cation transport during
water flow in soil: two approaches. Soil Science Society of America Journal
57, 3–9.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear
parameters. SIAM Journal on Applied Mathematics 11, 431–441.
MDH Engineered Solutions Corp., 2003. Evaluation of computer models forpredicting the fate and transport of salt in soil and groundwater. Phase II
Report.
Millington, R.J., Quirk, J.M., 1961. Permeability of porous solids. Transactionsof the Faraday Society 57, 1200–1207.
Nkedi-Kizza, P., Biggar, J.W., Selim, H.M., van Genuchten, M.Th., Wierenga, P.J.,
Davidson, J.M., Nielsen, D.R., 1984. On the equivalence of two conceptualmodels for describing ion exchange during transport through anaggregated oxisol. Water Resources Research 20, 1123–1130.
Parkhurst, D.L., Appelo, C.A.J., 1999. User's guide to PHREEQC (version 2)—acomputer program for speciation, batch-reaction, one-dimensionaltransport and inverse geochemical calculations. Water Resources
Investigation, Report 99-4259, Denver, CO, USA.
Phillips, I.R., Bond, W.J., 1989. Extraction procedure for determining solution
and exchange ions in the same soil sample. Soil Science Society ofAmerica Journal 53, 1294–1297.
Poeter, E.P., Hill, M.C., Banta, E.R., Mehl, S., Steen, C., 2005. UCODE_2005 andsix other computer codes for universal sensitivity analysis, calibration
and uncertainty evaluation. U.S. Geological Survey Techniques andMethods 6-A11.
Ramos, T.B., Šimůnek, J., Gonçalves, M.C., Martins, J.C., Prazeres, A.,
Castanheira, N.L., Pereira, L.S., 2011. Field evaluation of a multicomponent solute transport model in soils irrigated with saline waters. Journalof Hydrology 407 (1–4), 129–144 2011.
Rayment, G.E., Higginson, F.R., 1992. Australian Laboratory Handbook of Soiland Water Chemical Methods. Inkata Press, Melbourne.
Šimůnek, J., van Genuchten, M.Th., 1996. Estimating unsaturated soil
hydraulic properties from tension disc infiltrometer data by numericalinversion. Water Resources Research 32, 2683–2696.
Šimůnek, J., Jacques, D., Hopmans, J.W., Inoue, M., Flury, M., van Genuchten,
M.Th., 2002. Solute transport during variably-saturated flow—inverse
methods, In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis, Part
1, Physical Methods, Chapter 6.6, Third ed. SSSA, Madison, WI, pp.
1435–1449.
Šimůnek, J., Jacques, D., van Genuchten, M.Th., Mallants, D., 2006. Multicomponent geochemical transport modeling using HYDRUS-1D and HP1.
Journal of the American Water Research Association 46, 1537–1547.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th., 2008a. The
Hydrus-1D Software Package for Simulating the One-dimensional
Movement of Water, Heat, and Multiple Solutes in Variably-saturated
Media. Department of Environmental Sciences, University of California
Riverside, Riverside, California.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2008b. Development andapplications of the HYDRUS and STANMOD software packages, andrelated codes. Vadose Zone Journal. doi:10.2136/VZJ2007.0077, Special
Issue “Vadose Zone Modeling”.7, 587–600.

D. Jacques et al. / Journal of Contaminant Hydrology 142–143 (2012) 109–125Šimůnek, J., Jacques, D., Twarakavi, N.K.C., van Genuchten, M.Th., 2009.
Selected HYDRUS modules for modeling subsurface flow and contaminant transport as influenced by biological processes at various scales.
Biologia 64, 465–469.
Smiles, D.E., Smith, C.J., 2004a. A survey of the cation content of piggery
effluents and some chemical consequences of their use to irrigate soils.
Australian Journal of Soil Research 42, 231–246.
Smiles, D.E., Smith, C.J., 2004b. Absorption of artificial piggery effluent by
soil: a laboratory study. Australian Journal of Soil Research 42, 96–975.
Smiles, D.E., Smith, C.J., 2008. Absorption of gypsum solution by a potassic
soil: a data set. Australian Journal of Soil Research 46, 67–75.
Toorman, A.F., Wierenga, P.J., Hills, R.G., 1992. Parameter estimation of
hydraulic properties from one-step outflow data. Water Resources
Research 28, 3021–3028.
Toride, N., Leij, F.J., van Genuchten, M.Th., 1995. The CXTFIT code for
estimating transport parameters from laboratory or field tracer experiments. Version 2.0, Research Report No. 137. U. S. Salinity Laboratory,
USDA, ARS, Riverside, CA.

125

van Genuchten, M.Th., 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal 44, 892–898.
Voegelin, A., Vulava, V.M., Kuhnen, F., Kretzschmar, R., 2000. Multicomponent transport of major cations predicted from binary adsorption
experiments. Journal of Contaminant Hydrology 46, 319–338.
Voegelin, A., Vulava, V.M., Kretzschmar, R., 2001. Reaction-based model
describing competitive sorption and transport of Cd, Zn and Ni in anacidic soil. Environmental Science & Technology 35, 1651–1657.
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V.V., 2008.
Inverse modeling of subsurface flow and transport properties: a review
with new developments. Vadose Zone Journal 7, 843–864.
Vulava, V.M., Kretzschmar, R., Rusch, U., Grolimund, D., Westall, J.C.,
Borkovec, M., 2000. Cation competition in a natural subsurface material:modelling of sorption equilibria. Environmental Science & Technology
34, 2149–2155.