Journal of Hydrology 264 (2002) 185–200
www.elsevier.com/locate/jhydrol

Sensitivity analysis of physical and chemical properties affecting
field-scale cadmium transport in a heterogeneous soil profile
P. Seuntjens*, D. Mallantsb, J. Šimůnekc, J. Patyna, D. JacquesbaVlaamse Instelling voor Technologisch Onderzoek, Flemish Institute for Technological Research, Boeretang 200, B-2400 Mol, Belgiumb
SCK-CEN, Belgian Nuclear Research Centre, Boeretang 200, 2400 Mol, BelgiumcGeorge E. Brown, Jr., Salinity Laboratory, USDA, ARS, 450W. Big Springs Road, Riverside, CA 92507-4617, USA
Received 27 August 2001; revised 12 March 2002; accepted 28 March 2002

Abstract
Field-scale transport of reactive solutes depends on spatially variable physical and chemical soil properties. The quantitative
importance of physical and chemical parameters required for the prediction of the field-scale solute flux is generally unknown.
A sensitivity analysis is presented that ranks the importance of spatially variable water flow and solute transport parameters
affecting field-scale cadmium flux in a layered sandy soil. In a Monte-Carlo simulation approach, partial rank correlation
coefficients were calculated between model parameters and cadmium flux concentrations at various time steps. Data on theheterogeneity of flow and transport parameters were obtained from a 180 m-long and 1 m-deep Spodosol transect. Each soil
layer was described in terms of probability density functions of five model parameters: two shape parameters of van
Genuchten’s water retention curve, saturated hydraulic conductivity, dispersivity and soil – water distribution coefficient. The
results showed that the cadmium flux concentrations at the bottom of the soil profile were most sensitive to the cadmium
deposition rate and the soil– water distribution coefficient of all soil horizons. The maximum cadmium flux concentrations werealso affected by variations in hydraulic conductivity of the humic topsoil horizons. Variations in shape parameters of the water
retention curve did not significantly affect the field-scale cadmium flux. Variations in the dispersivity of the subsoil significantly
influenced the early time cadmium concentrations. Monte-Carlo simulations involving non-linear sorption showed that
cadmium flux concentrations were dominated by variations in the sorption constant and in the exponent of the Freundlich
isotherm. q 2002 Elsevier Science B.V. All rights reserved.
Keywords: Cadmium; Solute transport; Monte-Carlo; Sensitivity analysis; Uncertainty analysis; Stream-tube models

1. Introduction
Transport of reactive chemicals in soils is affectedby complex interactions involving biological, chemical and physical processes (Brusseau, 1989). The
sorption and transport processes of cadmium (Cd) in* Corresponding author.
E-mail address: piet.seuntjens@vito.be (P. Seuntjens).

soils were intensively studied (e.g. Chardon, 1984;
Christensen, 1984; Selim et al., 1992; Boekhold,
1992; Temminghoff et al., 1995 among many others).
For acid sandy soils like those focused on in this work,it was shown that sorption of Cd can be describedempirically using a Freundlich isotherm with parameters that account for the effects of pH, calcium
concentration, ionic strength and inorganic complexation. These ‘scaled’ isotherms were successfully used

0022-1694/02/$ - see front matter q 2002 Elsevier Science B.V. All rights reserved.
PII: S 0 0 2 2 - 1 6 9 4 ( 0 2 ) 0 0 0 7 1 - 9

186

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

by Temminghoff et al. (1995) to predict sorbed Cd
concentrations in samples from a sandy soil differingin pH, ionic strength and soil water composition.
Soils are heterogeneous, thus causing soil properties to vary from one place to another (Jury, 1985;
Jacques et al., 1999). Soil heterogeneity has a
profound effect on the transport of contaminantsfrom the soil surface to the groundwater table (van der
Zee and van Riemsdijk, 1987). This physical andchemical heterogeneity is often represented byvariations in hydraulic conductivity (Jury, 1985) andsoil – water distribution coefficient (Robin et al.,1991). Macro-scale heterogeneity can be accounted
for by stream-tube models that conceptualize a field asconsisting of a number of parallel non-interacting soilcolumns (Dagan and Bresler, 1979; Toride and Leij,
1996). Each stream-tube is envisaged as a single
realization of a multivariate probability density
function of stochastic variables. Van der Zee andvan Riemsdijk (1987) used a stream-tube model topredict Cd transport at the field-scale. A stochastic
distribution coefficient was calculated from variations
in pH and organic matter content, while scaling theory
was used to derive the probability density function of
the water flux. Boekhold and van der Zee (1991) intheir sensitivity analysis showed that field-scale Cd
transport is largely influenced by deposition rate,water flow rate and sorption constant. A second typeof models, stochastic-continuum models, use random
space functions (RSF) which are characterized by amean, variance and a specified autocorrelation(Chrysikopoulos et al., 1990; Destouni and Cvetkovic,
1991). Using numerical Monte-Carlo simulations andmoment analysis, Bosma (1994) found that theposition of a non-linearly interacting contaminant
plume is sensitive to the physical and chemical
heterogeneity of the flow domain and to non-linear
sorption, whereas the shape of the plume is dominatedby the degree of non-linear sorption.
Studies that quantitatively estimate the sensitivityof a reactive transport model output to variations inmodel parameters, are scarce. In the present work,
field-scale Cd transport is calculated using MonteCarlo simulations. The field is viewed as a series ofparallel non-interacting vertical stream-tubes. Localscale transport in each stream-tube is described
deterministically by an equilibrium convection –dispersion model. Soil samples from a Cd-contami-

nated field were taken to determine the statistical
moments (mean and variance) of the model parameters. Few transport studies thus far considered thetemporal variations of model sensitivity. As part ofthis study, this time-dependency was addressed by theuse of partial rank correlation coefficients between
model parameter and model output for each calculation time step.

2. Materials and methods
2.1. Sampling and laboratory measurements
Soil samples were taken from a 180 m-long and
1 m-deep transect at the ‘Kattenbos’ experimental
field (Lommel, Belgium). The transect comprises a
Spodosol toposequence in an area that has been
contaminated for more than 100 years by atmosphericdeposition of heavy metals emitted by non-ferrous
industry. Soils are classified as sandy Spodosols, i.e. atypic Placohumod at the dry end of the toposequence
and a typic Aquod at the wet end of the transect (SoilSurvey Staff, 1998).
To optimize the sampling and data collection
effort, a sequential sampling strategy was adopted.
Purpose of the first sampling was to determine thevertical variability in both the water flow and transportparameters, using samples from each soil horizon. For
this purpose 1 m-long cores were taken in triplicate at
three distinct locations along the transect (at the upper
part, in the middle and at the lower end of the transect,
see Fig. 1). Organic matter content, pH and metal
content were determined on 60 samples collectedfrom each soil horizon in the three cores (i.e. six soil
horizons in the wet Spodosol, seven in the middle andseven in the dry Spodosol).
Close to the position where the soil cores weretaken, triplicate in situ measurements of fieldsaturated hydraulic conductivity Kfs were conductedin each soil horizon using a pressure infiltrometer(Elrick and Reynolds, 1992). Following each conductivity measurement, small 0.0001 m3-soil cores were
driven in the soil volume inside the infiltrometer forthe determination of the drainage branch of the water
retention curve. The water retention data wereobtained with a sand-box apparatus for soil water
potential h values of 2 0.01, 2 0.03, 2 0.10, andP. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

187

Fig. 1. Overview of the field sampling scheme. Arrows indicate locations where the soil cores were taken for dispersivity measurements.
Numbers at the right side of the transect indicate the sequence of layers used in the Monte-Carlo Analysis.

2 0.32 m, and in low and high pressure chambers forh values of 2 2.00 and 2 6.31 m and h values of
2 25.12 and 2 158.48 m, respectively. Water retention parameters a, n, ur, us were determined using thenon-linear least-squares optimization code RETC
(van Genuchten et al., 1991). This was done foreach of the 60 retention curves, of which 18 were
attributed to Layer 1, 30 to Layer 2 and 12 to Layer 3.
The data set on physical properties was completed bydetermining dispersivities in large soil columns taken
at the wet and dry end side of the transect.
Displacement experiments with a non-sorbing tracer(chloride) in two undisturbed 1 m-long and 0.8 mdiameter soil columns were conducted to determine
the dispersivity of each soil horizon. Details of thedisplacement experiments can be found Seuntjens
et al. (2001a).
A second sampling was aimed at determining thehorizontal variation in metal content and chemical
properties (i.e. organic matter content and pH) ofindividual soil layers along the transect. Soil cores
were taken to a depth of 1 m at equally spaced
intervals of 5 m along the transect. At six locations
short-distance-samples were taken in nested samplings spaced 0.5, 1 and 2 m apart. In this way, 66
1 m-deep soil cores were sampled (Fig. 1). From thesecores the thickness of each soil horizon was noted. To
reduce the number of samples, individual horizons
were grouped into functional layers according toprocedures described in Seuntjens et al. (1999). A

functional property, i.e. the Cd travel time, was used
in a statistical analysis to derive the functional layers
from the soil horizons. Based on the non-parametricMann – Whitney U-test (Mann and Whitney, 1947),
horizons having statistically the same functional
property were grouped, whereas horizons showing
significantly different Cd travel times were separated
into different functional layers. Three functional
layers were defined: Layer 1, containing the humic
A and Bh horizons of both the Aquod soil (wet sandysoil) and the Humod soil (dry sandy soil); Layer 2,
comprising the E, Bh2 and Bh þ C horizons andLayer 3, made up of the parent material (C1 and C2
horizon of Humod soil). As a result of the soil horizon
grouping procedure, the entire field is conceptualized
as a five-layered soil with mean layer thickness andwithout distinction between wet and dry Spodosols.
As shown in Fig. 1, the boundary between the wet andthe dry Spodosol is situated at approximately 15 mfrom the beginning of the 180 m-long transect,indicating that the wet Spodosol only accounts for aminor part of the transect. Therefore, we assumed thefunctional layer sequence of the dry Spodosol, as
shown on the right side of the transect in Fig. 1, to bevalid for the entire field.
2.2. Transport model
Solute transport calculations were made using the
HYDRUS1D numerical finite-element code (Šimůnek

188

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

et al., 1998). A five-layer profile is considered
consisting of a sequence of three functional layers:
Layer 1 (0 – 8 cm) –Layer 2 (8 –18 cm) – Layer 1 (18 –25 cm) –Layer 2 (25 –48 cm) – Layer 3 (48 – 100 cm).
To eliminate possible numerical oscillations in thecalculations due to spatial and temporal discretization,
the stability criterion developed by Perrochet andBerod (1993) was used (i.e. Peclet number , 2 andCourant number Cr , 1), resulting in Dz ¼ 1 cm,
while the time-step Dt was adjusted to satisfy thecriterion. Solute transport was simulated in an initiallycadmium-free soil profile through a simultaneous
solution of the Richards’s (1931) equation (Eq. (1))
and the convection– dispersion equation (Eq. (2))(Lapidus and Amundson, 1952), given by:
 

›u
›
›h
¼
K
21
ð1Þ
›z
›t
›z

r

›S
›C
›2 C
›C
þu
¼ uD 2 2 J w
›t
›t
›z
›z

ð2Þ

where K is the unsaturated hydraulic conductivity(L T21), h the soil water pressure head (L), S theamount of Cd adsorbed (M M21), C the concentrationor activity of Cd in the soil solution (M L23), t thetime (T), z the distance from the soil surface (L), r thesoil bulk density (M L23), u the soil water content(L3 L23), D the hydrodynamic dispersion coefficient(L2 T21) and Jw is the volumetric water flux density(L T21). All flow and transport parameters wereassumed to be constant in time. The boundary
conditions for water flow are determined by a constantwater flux at the soil surface and free drainage (zeropressure gradient) at the bottom boundary:›h2 1 ¼ Jw
hðz; tÞ ¼ hi ðzÞ;2K
z ¼ 0;
›z

›h
¼0
›z

z¼L
ð3Þ

where hi ðzÞ is the initial steady state pressure head
distribution, determined by a stationary flux of
0.064 cm d21, which is the long-term average precipitation surplus (total rainfall minus actual evapotranspiration) for the region (Patyn, 1997).
A zero initial Cd concentration distribution wasused together with a third-type flux boundary
condition at the top and a second-type zero concen-

tration gradient at the bottom:
cðz; 0Þ ¼ 0;

›c
¼0
›z

2 uD

›c
þ Jw c ¼ Jw co ðtÞ›z

z ¼ 0;

z¼L
ð4Þ

in which Jwco is the actual Cd deposition flux. This
flux was estimated from annually averaged mean
solute deposition data as obtained from a nearby
station (Lommel, station MN301) from 1980 to 1998
(VMM, 1999). The overall historical deposition rate
was estimated by integrating the observed solutedepth profiles up to 1 m following procedures
proposed by Cernik et al. (1994). The time-dependent
historic deposition rate was calculated from a
production history pattern for the regional non-ferrous
industry (LISEC, 1989). For the purpose of thesensitivity analysis, the industrial emissions andconsequently the deposition rate were assumed to be
proportional to the overall production rate. The
production itself was assumed to increase exponentially for 85 years from the time the production started(i.e. 1895), after which a change in the production
process was introduced and additional actions weretaken to reduce heavy metal emissions. The change inthe production process caused a sudden decrease in
the deposition flux, as shown by the deposition data
between 85 and 100 years in Fig. 2. Since theestimation of the historical deposition flux is uncertain, the calculated input concentration was multipliedwith an uncertainty factor I. We assumed I to betriangularly distributed between 0.5 and 2, with its
probability distribution function pdf reaching amaximum at I ¼ 1.5.
Two deposition scenarios were considered (bothstarting at 85 years after production started): (i) the
deposition rate will remain at a constant level equal tothe average deposition flux of the last 10 years (i.e.
0.02 mg l21 between 90 and 100 years), (ii) after
T ¼ 85 years, the deposition rate will drop to zero.
The total simulation time was 400 years.
2.3. Calculation of the soil – water distribution
coefficient
The soil – water distribution coefficient Kd was

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

189

Fig. 2. Estimated former Cd deposition (before T ¼ 85 years after the contamination started) and two alternative future Cd deposition scenarios
expressed in flux concentrations in the infiltrating rain water.

estimated from soil properties. The organic carbon
content and pH of each soil layer and a predefined Ca
activity were used to calculate the sorption constant kfrom Eq. (5):k ¼ kp ocðHþ Þ20:69 ðCaÞ20:34

ð5Þ

where k p (1023.77 mol(12n ) ln kg21) is a universal
sorption constant (n in the dimension of k is theFreundlich exponent) determined in batch experiments (Seuntjens et al., 2001b), oc is the organiccarbon content (%), (Hþ) is the proton activity in soil
solution (mol l21) and (Ca) is the Ca activity in soil
solution (mol l21). Eq. (5) is derived from the ThreeSpecies-Freundlich model of Temminghoff et al.
(1995) which was successfully used to describe Cd
sorption in acid sandy soils. In our study, two different
Ca activities were used to calculate k: a high Caactivity of 0.01 M and a low Ca-activity of 0.0003 M.
A low Ca concentration generates a high sorption
constant k, whereas a high Ca concentration results ina low k. The distribution coefficient kd was calculated
by linearizing the non-linear Freundlich isotherm (i.e.
s ¼ kcnCd ):kd ¼ kcnCd 21

ð6Þ

where c is an average Cd concentration for the entire
field, i.e. c ¼ 0.05 mg l21, whereas the exponent nCdof the Freundlich isotherm was estimated from batch
experiments conducted on samples from selectedhorizons of the same field soil, i.e. nCd ¼ 0.86
(Seuntjens et al., 2001b).

2.4. Monte-Carlo sensitivity analysis
In the Monte-Carlo analysis, each functional layer
was characterized by a multivariate probability
density function (pdf) of five model parameters (Ks,a, n, Kd, l ). Variations in layer thickness were notincluded in the analysis, but may be important. The
Shapiro – Wilk W-test (Shapiro and Wilk, 1965) wasused to determine whether the variables werenormally or lognormally distributed. A W-value
close to one indicates that the probability density
function (pdf) closely resembles a normal distribution. If the W-statistic is significant ( p , 0.05), thehypothesis that the respective distribution is normal, isrejected.
A latin-hypercube-sampling (LHS) was used togenerate the multivariate parameter distributions for
the Monte-Carlo simulations. Details of the inductionof a correlation matrix on a multivariate input randomvariable are given by Iman and Conover (1982). The
matrix operations are performed by the PREP utilityof the LISA-code (Saltelli, 1987). Iman and Helton(1985) stated that a sample size between 2m and 5m
provides good results, where m is the number ofvariables. Since 16 variables were considered in theMonte-Carlo analysis (five flow and transport parameters for three layers and one parameter describing
the variability for the input function), a sample size of
N ¼ 100 was assumed to be sufficiently large. Thismeans that the field is conceptualized as a set of 100
non-interacting stream-tubes.

190

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

Table 1
Overview of Monte-Carlo simulation transport scenarios
Scenario

Deposition

Kd

nCda

1
2
3
4
5

Continuous sourceb
Continuous source
Finite sourcee
Finite source
Finite source

Highc
Lowd
High
High
High

Constant, 1
Constant, 1
Constant, 1
Constant, 0.86
Variable, m ¼ 0.86; s ¼ 0.05anCd: Freundlich isotherm exponent.
Future deposition remains at actual deposition.
c
A high Kd corresponds to a low Ca activity.
d
A low Kd corresponds to a high Ca activity.
e
Future deposition becomes zero.
b

To estimate the sensitivity of the model output toan individual input variable, a partial rank correlation
coefficient (PRCC) was calculated (Helton et al.,1991). Rank regression is preferred over usual
regression when poor linear fits are caused by nonlinear data. The rank regression (PRCC) analysis isconducted with the SPOP utility of the LISA code
(Saltelli, 1987). The advantage of the PRCC over
several other sensitivity indicators is that it iscalculated at each time step, which allows to introduce
a time dependency in the sensitivity analysis. This
type of sensitivity analysis offered the possibility toshow that the model result could be sensitive tovariations of a certain parameter at one time but notnecessarily at other times.
Monte-Carlo simulations were conducted to determine the sensitivity of the Cd flux concentrations at a
depth of 1 m to the model parameters of eachfunctional layer. Five cases were considered (Table1). The simulations with a high or low mean Kd(scenarios 1 and 2) were conducted to assess the effectof the mean Kd on the sensitivity of the model output
to variations in the remaining parameters (i.e. whether
or not the same sensitivity is observed when the mean
sorption is low or high). Two scenarios wereconsidered that corresponded to two different management options: a continuous source (referred to as
scenario 1 in Table 1: the actual deposition rate exists
for another 315 years) and a finite source (referred toas scenario 3 in Table 1: deposition rate drops to zeroat T ¼ 85 years after production started). Whereas
scenarios 1– 3 dealed with linear sorption, scenarios 4
and 5 were conducted to assess the effect of non-linear
sorption. In case of Freundlich type sorption iso-therms, the degree of non-linearity was represented by
the value of the Freundlich exponent nCd. The resultsof scenario 4 (fixed nCd ¼ 0.86) were compared toscenario 3 (nCd ¼ 1) for a finite pollutant source.
Scenario 5 assumed the stochastic Freundlich exponent, similar for all three layers, to be normally
distributed with a mean of 0.85 and a standard
deviation of 0.05.

3. Results and discussion
3.1. Variation in the estimated soil –water distribution
coefficient Kd
The variation of Kd along the transect in the four
sampled soil horizons is shown in Fig. 3. The highest
Kd values were found for the A and spodic Bhhorizons, due to a higher organic carbon content. The
eluvial E horizon exhibited a similar Kd as the C
horizon. The effect of a higher organic carbon content
of the E horizon as compared to the C horizon wascompensated by a lower pH. Note from Fig. 3 that theKd of the C horizon varied less in comparison withthose of the topsoil horizons. Because the wet
Spodosol has no C horizon, data for the first 15 m of
the transect were not shown in Fig. 3.
3.2. Multivariate probability density functions
Mean and standard deviation of the untransformed
and loge-transformed model parameters for the threefunctional layers and the corresponding results of theW-test are given in Table 2. Table 2 shows that, based

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

Fig. 3. Horizontal variation in the estimated distribution coefficient
Kd of four sampled horizons of the experimental field.

on the significance of the W-statistic, the water
retention parameters were equally well described bya normal and a lognormal pdf, except for shapeparameter a of Layer 2. The saturated hydraulicconductivity Ks was lognormally distributed in Layers
1 and 2, contrary to Layer 3 which showed a normally
distributed Ks. The W-value closer to unity for the
loge-transformed dispersivity l in Layer 2 indicatedthat the pdf of l in this layer closely resembles alognormal distribution.
Layer 3 showed the highest saturated hydraulicconductivity Ks and the lowest distribution coefficient
Kd of the three functional layers. The coefficients of

191

variation (CV) ofpthetransformed variables were
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
calculated using expðs2 Þ 2 1; where s 2 is thevariance of the loge-transformed variable. Variability,
as expressed by the CV, was larger for Layers 1 (the A
and Bh horizons) and 2 (the E and Bh2 horizon), as
compared to Layer 3. Layer 3 comprised the parent
soil material which was much more homogeneousthan the topsoil horizons. The saturated hydraulicconductivity displayed a CV ranging from 32 to
175%, which is in the range of field CVs reported by
Jury (1985). We found a relatively high variability indispersivity l, having a CV of 46 –83%. Mallants
et al., 1996 found a CV of 61% for the logtransformed
dispersion coefficient. The CV of 27– 40% for the
distribution coefficient Kd lies well within the range
reported by Jury (1985). The lowest variability was
found for the saturated water content us. This is inagreement with observations in other studies compiled by Jury (1985). Because of the small variability
in us and since ur has little influence on water flow ascompared to other hydraulic parameters (Durner,
1994), both parameters were assumed to be constantin the subsequent Monte-Carlo simulations.
Correlations among the five variables (Ks, a, n, Kd,l ) were calculated for each layer (Table 3). For thispurpose, data from the first sampling (N ¼ 60) wereused to calculate the correlations because only here all
five parameters were determined at the same location.
Significant ( p , 0.05) positive correlations were
found between Ks and n of Layer 1, between Ks anda of Layer 2, and between n and l of Layer 3.
Significant negative correlations occurred between n
and l of Layer 1, and n and Kd of Layer 2.
Insignificant correlations were set to zero. The LHS
technique generated a spurious correlation between Ks
and l of Layer 1 (2 0.63). This correlation was notdefined in Table 3, but is the result of the LH samplingtechnique. The correlation between the shape parameter n of the water retention curve and l on the one
hand and Ks on the other, most likely induced thecorrelation between Ks and l.
3.3. Sensitivity of the Cd flux concentrations
In the following discussion, the term ‘concentration’ refers to ‘flux concentrations at a depth of
1 m’. Fig. 4 shows the field-average Cd breakthroughcurve for scenario 1 (Table 1). Also shown is the192

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

Table 2
Statistics of the untransformed and loge-transformed model parameters. Also given are the results of the Shapiro–Wilk W-statistic (N is numberof samples, m is mean, CV is coefficient of variation, W is Shapiro– Wilk’s W-value and p is the p-level of the W-statistic). The statistical
distribution used in the subsequent Monte-Carlo analysis is displayed in boldloge-Transformed

UntransformedNm

CV (%)Wpm

CV (%)Wp

Layer 1
Ksa, cm d21
urusa, cm21n
Kdb, l kg21
Kdc, l kg21
l, cm

78
18
18
18
18
132
132
15

382
0.097
0.488
0.015
2.07
9.93
38.1
1.84

119.1
46.4
7.2
26.7
25.1
41.5
43.3
65.8

0.73
0.92
0.97
0.96
0.93
0.81
0.83
0.91

0.00
0.16
0.85
0.51
0.19
0.00
0.00
0.14

408
0.098
0.488
0.015
2.08
10.0
38.1
1.48

174.9
54.9
7.3
27.4
24.8
39.8
41.9
79.9

0.98
0.97
0.93
0.97
0.95
0.92
0.94
0.97

0.39
0.85
0.21
0.82
0.38
0.00
0.00
0.83

Layer 2Ks, cm d21
urusa, cm21n
Kdb, l kg21
Kdc, l kg21
l, cm

31
30
30
30
30
76
76
40

563
0.059
0.43
0.020
2.67
4.34
17.2
2.45

66.6
45.8
8.4
25.0
17.0
46.3
47.0
134.7

0.83
0.84
0.93
0.81
0.98
0.91
0.90
0.48

0.00
0.00
0.05
0.00
0.80
0.00
0.00
0.00

564
0.059
0.43
0.020
2.68
4.35
17.2
2.21

66.7
45.5
8.7
23.7
18.0
47.3
48.3
83.3

0.95
0.89
0.92
0.93
0.96
0.97
0.95
0.87

0.20
0.01
0.03
0.06
0.42
0.19
0.02
0.00

Layer 3Ks, cm d21
urusa, cm21n
Kdb, l kg21
Kdc, l kg21
l, cm

13
12
12
12
12
56
56
16

1221
0.03
0.41
0.021
4.34
2.73
15.6
1.80

30.4
20.0
5.1
9.5
12.2
27.5
27.4
57.2

0.95
0.95
0.97
0.94
0.96
0.96
0.95
0.71

0.52
0.55
0.85
0.42
0.69
0.09
0.06
0.00

1235
0.03
0.41
0.021
4.34
2.74
15.6
1.80

35.4
22.3
5.1
7.7
12.6
27.6
27.5
46.2

0.89
0.93
0.97
0.95
0.94
0.98
0.98
0.92

0.09
0.37
0.85
0.58
0.46
0.82
0.73
0.15

a

Saturated hydraulic conductivity from both laboratory (i.e. N ¼ 66) and field measurements (i.e. N ¼ 12).
High Ca activity.
c
Low Ca activity.

btime-dependency of the PRCC for the three functional
layers. The shaded area represents the range of PRCC
values, which do not significantly differ from zeroaccording to a two-sided t-test. The PRCC value
between Cd concentration and the deposition rateincreased gradually with time and reached a constantvalue close to unity. After 400 years, steady state wasreached in the soil profile and evidently, Cd
concentrations were perfectly correlated with the Cd
deposition rate.
The PRCC between Kd and Cd concentration
increased rapidly with time at the start of thesimulation and became relatively constant andnegative between 20 and 140 years. The PRCC
subsequently changed from negative to positive,
reached a maximum positive value of 0.4 (Layer 1)
and 0.7 (Layers 2 and 3) 240 years after thecontamination started, and decreased towards theend of the simulation. A negative correlation between
Kd and Cd concentration was expected since higher
sorption results in a lower Cd concentration. The
negative correlation between Kd and Cd concentration
was maximum at about 140 years, i.e. at the mean
breakthrough time of the adsorption front of the Cd
plume. When the Cd flux concentration decreased as aresult of the decrease in deposition rate, the PRCC

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200
Table 3
Cross-correlations among flow and transport parameters. Statistically significant correlations ( p , 0.05) are displayed in bold. N is
the number of samples
Ks
Layer 1 (N ¼ 30)
1
Ks
a
20.05
n
0.81
Kd
0.39
l
20.40
Layer 2 (N ¼ 18)
Ks
1
a
0.63
n
20.06
Kd
0.04
l
0.17
Layer 3 (N ¼ 12)
Ks
1
a
20.33
n
0.01
Kd
20.34
l
0.30

a

1
20.13
20.07
20.09

1
20.35
20.13
0.22

1
20.28
0.11
20.57

n

1
0.34
20.74

1
20.41
20.18

1
20.11
0.65

Kd

1
20.27

1
0.10

1
20.10

Table 4
Effect of mean distribution coefficient Kd on the maximum Cd flux
for a continuous source. Statistically significant correlations( p , 0.05) are displayed in bold

l

1

1

1

between Kd and Cd concentration became positive.
Stream-tubes with higher Kd have accumulated more
Cd and released more Cd when Cd deposition was
reduced. Finally, the PRCC between Kd and Cd
concentration was insignificant when steady state wasreached. The temporal changes in the PRCC between
individual model parameters and Cd concentration
were further shown in Fig. 5. Keeping all parametersconstant, the effects on Cd concentration of a varying
Kd were shown in Fig. 5(a). Before 120 years, a low
Kd value caused higher solute concentrations than ahigh Kd. After 120 years the concentration of Cd
displaying a higher Kd exceeded the concentration
corresponding to a lower Kd, which illustrates the
positive correlation between Cd concentration and Kd.
The PRCC calculated between dispersivity l andCd concentration showed some interesting features. In
Layer 1, the Cd concentrations were uncorrelated withl, even at early times. The most likely explanation forthis insensitivity was the small contribution of Layer 1
to the overall longitudinal dispersion. In contrast to
Layer 1, Layers 2 and 3 showed rapidly increasing
PRCCs between l and Cd concentration at earlytimes. The PRCC reached a maximum after 20 years,
and then decreased gradually. A higher dispersivity

193

High Kd

Low KdLayer 1
a
n
Ks
l
Kd

0.12
20.17
0.26
20.11
20.43

0.15
20.16
0.25
20.05
20.47

Layer 2
0.06
n
Ks
l
Kd

0.17
0.08
0.13
20.04
20.53

0.08
0.11
20.11
20.45

Layer 3
a
n
Ks
l
Kd

0.05
20.05
20.06
0.08
20.35

0.14
20.09
20.04
0.12
20.20

Jwco

0.79

0.95

caused Cd to reach greater depths at early times, thus
implying a positive correlation between Cd concentration and l. Fig. 5(b) showed that at early times(before 120 years), higher concentrations werepredicted when a larger dispersivity was assumed.
This lead to a positive correlation between dispersivity and Cd flux. Fig. 5(b) also showed that negative
PRCCs are expected when the maximum Cd concentration breaks through (i.e. the higher l produces thelower peak concentration).
Fig. 4 also showed that variation in hydraulicproperties (Ks, a, n ) of all layers did not significantly
affect Cd concentrations. Only variations in thehydraulic conductivity of the humic topsoil layers
affected the maximum Cd concentration. Boekhold
and van der Zee (1991) found a significant effect ofpore water velocity on Cd leaching. A possible
explanation may be the larger variation in porewater velocity v (CV ¼ 100%) for their experimental
field as compared to the CV of v in our field(CV ¼ 19%). Secondly, they considered only leaching from a topsoil, ignoring the properties of deeper
soil layers.

194
P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200
Fig. 4. Field-scale Cd flux concentration at 1 m depth for a continuous future Cd source. Evolution with time of the Partial Rank Correlation Coefficient (PRCC) between model
parameters and Cd flux for the three functional soil layers. The PRCC values outside the shaded area are significant at the 95% level.

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

195

Fig. 6. Sensitivity of the field-scale Cd flux to the saturated
hydraulic conductivity Ks, the dispersivity l and the distribution
coefficient Kd for a finite deposition rate and assuming linear
sorption.

3.4. Effect of the mean distribution coefficient Kd

Fig. 5. Illustration of the effects of: (a) the soil–water distribution
coefficient Kd (high Kd corresponds to a low Ca activity, low Kdcorresponds to high Ca activity) (b) the dispersivity l (high l is theupper 95% confidence limit on the mean, low l is the lower 95%
confidence limit on the mean) and (c) the Freundlich exponent nCd
on calculated Cd breakthrough flux concentrations. The vertical line
indicates the transition from positive (negative) to negative(positive) PRCCs.

Results of the sensitivity analysis thus far werebased on the high Kd-values (scenario 1). Additional
simulations were performed using low Kd-values(scenario 2) that were on average four times smaller.
Table 4 illustrates the effect of the mean Kd on themodel sensitivity, expressed in terms of the PRCC.
The maximum Cd concentration (i.e. at 100 and 176
years for the low and the high mean Kd, respectively)
was used as model output. The PRCC between
deposition rates and the maximum Cd concentrations
was larger for a low mean Kd. A decrease in the mean

196

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

sorption in cultivated and natural sandy soils. Bosma
et al. (1993), using a mean Kd varying over a factor of50, showed that the longitudinal spreading of alinearly reacting solute in a physical and chemical
heterogeneous porous medium was significantly
affected by the mean Kd. For equal degrees ofphysical and chemical heterogeneity it was shownthat chemical heterogeneity dominates the spreading
process for large Kd values, whereas physicalheterogeneity was dominating in case of a smallmean Kd.
3.5. Effect of deposition scenario

Fig. 7. Sensitivity of the field-scale Cd flux to the saturated
hydraulic conductivity Ks, the dispersivity l and the sorption
constant k for a finite deposition rate and assuming non-linear
sorption.

sorption increased the sensitivity of the maximum Cd
concentrations to variations in the deposition rate. The
maximum Cd flux was affected by the saturated
hydraulic conductivity of Layer 1 when the mean Kd
was both high and low. The variation in the model
variables of the C horizon (Layer 3, see Table 1) had
no significant effect on the maximum Cd
concentration.
In general, changing the mean Kd only minimallyaffected the sensitivity for other parameters. Of notehere is that the mean Kd values were purposely
restricted to values that were representative for Cd

The effect of two different future deposition rates(i.e. a continuous deposition and a finite deposition)on the sensitivity of the model results was evaluated.
Fig. 6 shows the evolution with time of the PRCC
between the Cd concentration and the relevant model
parameters (l, Kd and Ks) for a finite pollution source(scenario 4). The PRCC between the deposition rate
and Cd concentration increased almost linearly withtime, reached a maximum at the same time when themean flux reached a maximum, and decreased
towards the end of the simulation. Judging from the
PRCC in all three layers, Cd concentration was moresensitive to variations in the Kd than to variations in
the deposition rate, except when Cd reached its
maximum concentration.
Similar as for the continuous pollution source, the
PRCC between l and Cd concentration in Layers 2
and 3 (not shown here) increased rapidly towards a
maximum during the early stages of contamination,
reached a minimum when the breakthrough concentration was maximum and increased again towards theend of the simulation. A high l results in increased
spreading and tailing of the breakthrough curve. Thismeans that at later times higher concentrations will beobserved when l is larger. Conversely, low l values
imply less tailing and hence lower concentrations atlater times.
The PRCC between Kd and Cd concentration
showed more or less the same behavior for all threelayers. At the beginning, the PRCC quickly reachednegative values of 2 0.6 for Layer 1 and 2 0.9 forLayers 2 and 3 (not shown). For Layer 1, the PRCC
further decreased until a minimum value of 2 0.8 was
obtained after 145 years. Until T ¼ 225, the PRCC

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

197

Fig. 8. Effect of the degree of non-linearity expressed by the Freundlich exponent nCd on Cd concentration and on the sensitivity of the model toother transport parameters (dispersivity, Kd) and the deposition rate.

followed the same pattern for cases 1 (continuous
deposition) and 4 (finite deposition). At later times,
the PRCC between Kd and Jsf for the finite source
remained high (or decreased only slightly for Layer3), while decreasing to nearly zero for the continuous
source (Fig. 4).
Summarizing the results for the finite source weshowed that: (i) the dispersivity l and the distribution
coefficient Kd affected the Cd flux before the breakthrough of the adsorption front and after the passageof the desorption front (ii) Kd affects the mean
breakthrough of the adsorption and desorption frontsand (iii) the deposition rate affected mostly the
maximum solute flux. Similar to the continuous Cdsource, the hydraulic conductivity Ks of Layer 1 only
affected the maximum Cd flux.
3.6. Effect of non-linear sorption
The effect of non-linear sorption on field-scale Cd
transport was studied using the Monte-Carlo simulation approach, but with the linear isotherm(nCd ¼ 1) being replaced by a non-linear Freundlichtype isotherm (nCd ¼ 0.86). The effect of non-linear

sorption on the shape of the breakthrough curve isshown in Fig. 7. The breakthrough curve became
asymmetric by showing considerable tailing. This wasdue to increased retardation resulting from lower
concentrations in the back-end of the breakthroughcurve. The solute front moved somewhat slower forthe non-linear case compared to linear sorption, while
the maximum concentration was smaller. These
results are consistent with those by Bosma (1994)
and van Genuchten and Cleary (1982) who showedthat non-linear sorption causes more solute spreading.
Fig. 7 further shows the behavior of the PRCC withtime. The PRCC values for all layers remained zero
until T ¼ 70 years. Due to front sharpening, Cd
appeared at the bottom of the soil profile only after 70
years. From then on, the PRCC between Kd and Cd
concentration increased slower for non-linear sorption(Fig. 7) as compared to linear sorption (Fig. 6). This
indicated that non-linear sorption had a dominating
effect on Cd breakthrough at early times. Also theeffect of variations in l at early times diminished ascompared to the linear sorption case. The effect ofvariations in l on the early time Cd flux concentrations is counteracted by non-linear sorption. On the

198

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

other hand, non-linear sorption had a minor effect on
the PRCC in the tail of the breakthrough curve.
To quantify the impact of non-linear sorption on
the field-scale Cd breakthrough, a sensitivity analysis
was conducted treating the Freundlich exponent nCd
as a stochastic variable. The calculated field-scale Cd
flux is given in Fig. 8. Results suggested that aspatially variable Freundlich exponent nCd caused Cd
to break through more rapidly (at T ¼ 50 years) and toreach a higher mean concentration (kcl ¼ 0.11
mg l21) as compared to a spatially uniform nCd. The
breakthrough curve displayed a narrow and nearly
symmetrical shape in contrast to the deterministic
exponent case. The predicted Cd concentration
variance was larger in the ascending part of the
breakthrough curve than in the descending part.
The PRCC between nCd and Cd concentration
increased during the early stages of contamination to a
relatively constant and positive value. Its value laterchanged from positive to negative when the soluteconcentration reached a maximum, attaining its
lowest values (2 0.6 in Layer 1 and 2 0.8 in Layers
2 and 3) at T ¼ 200 years. These values did notchange anymore in the remaining part of thesimulation. The (positive or negative) effects of nCd
on the solute flux were most significant at relatively
early and late breakthrough times, and insignificant
when the center of Cd mass passed through the soil.
The positive effect of nCd on the solute concentration at early times can be deduced from the shape ofthe solute front, shown in Fig. 5(c) by the twosimulation runs with different exponents. When nCd is
low (e.g. nCd ¼ 0.75), the sorption front is selfsharpening and early time concentrations increase
rapidly with time. When nCd is larger (i.e.
nCd ¼ 0.95), the sorption front is wide and the soluteconcentration increases gradually. This means that at
early times the flux concentration is higher when nCd
is higher. On the other hand, the desorption part of abreakthrough curve is less sharp for the lower nCd as
compared to the higher exponent. Consequently, low
nCd values will lead to higher concentrations at theend of the contamination process, resulting in anegative correlation between model parameter nCd
and the Cd concentrations.
The effect of variations in nCd was larger for
Layer 3 than for Layers 1 and 2. This is most
likely due to the moderate variability of the waterflow and solute transport model parameters ofLayer 3. The effect of the other model parameters
and the deposition rate were much less when nCd
was treated as a stochastic variable. The effects ofKd and l were especially much less during theearly and final stages of contamination. Similar
results were obtained by Bosma (1994), whofound that variations in the degree of non-linear
sorption had a dominating effect on the dimensions of the contaminant plume, whereas physicaland chemical heterogeneity (i.e. sorption constant)hardly affected the plume geometry. The present
sensitivity analysis confirmed the important effectof non-linear sorption on the breakthrough ofreactive solutes.

4. Conclusions
Variations in field-scale Cd flux was found to be
dominated by the variations in deposition rate and theparameters of the Freundlich sorption isotherm. In
case of a pulse-type contaminant source, it was shown
that: (i) variations in the deposition rate mostly
affected the maximum Cd concentrations leachingfrom the soil profile, (ii) variations in dispersivity
affected the Cd concentrations before the mean Cd
adsorption front reached the bottom of the soil profile
and after the mean desorption front traveled throughthe medium, and (iii) variations in soil – water
distribution coefficient affected the position of themean adsorption and desorption front.
The effect of variations in the water retention
curve parameters on the field-scale Cd flux was
insignificant in nearly all scenarios considered.
Only variations in the hydraulic conductivity of
the humus topsoil layers (i.e. the A and Bhhorizon) affected the maximum Cd flux. The effectof variations in dispersivity on Cd leaching wasrestricted to deeper soil layers.
Non-linear sorption was found to have a dominating effect on the Cd flux concentrations. Non-linear
sorption lead to a compression of the solute plume at
the adsorption front and increased solute spreading at
the back-end of the breakthrough curve. It was also
shown that the sensitivity to other model parameters
was significantly reduced when non-linear sorption
was treated stochastically.

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

References
Boekhold, A.E., 1992. Field scale behaviour of cadmium in soil.
PhD dissertation, Agricultural University Wageningen, The
Netherlands, 181 pp.
Boekhold, A.E., van der Zee, S.E.A.T.M., 1991. Long term effects
of soil heterogeneity on cadmium behaviour in soil. J. Contam.
Hydrol. 7, 371–390.
Bosma, W.J.P., 1994. Transport of reactive solutes in heterogeneous
porous formations. PhD dissertation, Agricultural University
Wageningen, The Netherlands, 229 pp.
Bosma, W.J.P., Bellin, A., van der Zee, S.E.A.T.M., Rinaldo, A.,
1993. Linear equilibrium adsorbing solute transport in physically and chemically heterogeneous mporous formations. 2.
Numerical results. Water Resour. Res. 29, 4031–4043.
Brusseau, M.L., 1989. Modeling the transport of solutes influenced
by multiprocess nonequilibrium. Water Resour. Res. 9,
1971–1988.
Cernik, M., Federer, P., Borkovec, M., Sticher, H., 1994. Modelingof heavy metal transport in a contaminated soil. J. Environ.
Qual. 23, 1239–1248.
Chardon, W.J., 1984. Mobiliteit van cadmium in de bodem. PhD
dissertation, Agricultural University Wageningen, The Netherlands.
Christensen, T.H., 1984. Cadmium soil sorption at low concentrations: I. Effect of time, cadmium load, pH and calcium.
Water, Air, Soil Pollut. 21, 105–114.
Chrysikopoulos, C.V., Kitanidis, P.K., Roberts, P.V., 1990.
Analysis of one-dimensional solute transport through porous
media with spatially variable retardation factor. Water Resour.
Res 26 (3), 437 –446.
Dagan, G., Bresler, E., 1979. Solute dispersion in unsaturated
heterogeneous soil at field-scale: I. Theory. Soil Sci. Soc. Am. J.
43, 461 –467.
Destouni, G., Cvetkovic, V., 1991. Field scale mass arrival ofsorptive solute into the groundwater. Water Resour. Res. 27,
1315–1325.
Durner, W., 1994. Hydraulic conductivity estimation for soils withheterogeneous pore structure. Water Resour. Res. 27, 967–981.
Elrick, D.E., Reynolds, W.D., 1992. Infiltration from constant headwell permeameters and infiltrometers. Advances in measurement of soil physical properties: bringing theory into practice,
SSSA Special Publication No. 30.
Helton, J.C., Garner, J.W., McCurley, R.D., Rudeen, D.K., 1991.
Sensitivity analysis techniques and results for performanceassessment at the waste isolation pilot plant. US Department
Energy, SAND90-SA7103.
Iman, R.L., Conover, W.J., 1982. A distribution-free approach toinducing rank correlation among input variables. Commun.
Stat.—Simulat. Comput. 11 (3), 311 –334.
Iman, R.L., Helton, J.C., 1985. A comparison of uncertainty andsensitivity analysis techniques for computer models. US
Nuclear Regulatory Commission, NUREG/CR-3904.
Jacques, D., Kim, D.-J., Diels, J., Vanderborght, J., Vereecken, H.,
Feyen, J., 1999. Analysis of steady state chloride transport
through two heterogeneous field soils. Water Resour. Res 34(10), 2539–2550.

199

Jury, W.A., 1985. Spatial Variability of soil physical parameters insolute migration: a critical literature review. Research Report
EA-4228, Electrical Power Research Institute, Palo Alto,
California.
Lapidus, L., Amundson, N.R., 1952. Mathematics of adsorption in
beds. VI. The effect of longitudinal diffusion in ion exchange
and chromatographic columns. J. Phys. Chem. 56, 984–988.
LISEC, 1989. Inventory of existing measurements and production
data. Report Subgroup Air-Heavy metals in Northern Limburg(in Dutch).
Mallants, D., Vanclooster, M., Feyen, J., 1996. Transect study on
solute transport in a macroporous soil. Hydrol. Proc. 10, 55–70.
Mann, H.B., Whitney, D.R., 1947. On a test of whether one of two
random variables is stochastically larger than the other. Ann.
Math. Stat. 18, 50–60.
Patyn, J., 1997. Hydrodynamic model SCR Sibelco (in Dutch).
Perrochet, B., Berod, D., 1993. Stability of the standard Crank–
Nicholson Galerkin scheme applied to the diffusion– convection
equation: some new insights. Water Resour. Res. 29,
3291–3297.
Richards, L.A., 1931. Capillary conduction of liquids in porous
mediums. Physics 1, 318.
Robin, M.J.L., Sudicky, E.A., Gilham, R.W., Kachanoski, R.G.,
1991. Spatial variability of strontium distribution coefficientsand their correlation with hydraulic conductivity in the
Canadian Forces Base Borden Aquifer. Water Resour. Res.
27, 2619–2632.
Saltelli, A., 1987. PREP and SPOP utilities. Two Fortran programsfor sample preparation, uncertainty analysis and sensitivity
analysis in Monte Carlo simulation, Programs Description andUser’s Guide, Commission of the European Committee, JRC,
Ispra, EUR 11034 EN.
Selim, H.M., Büchter, B., Hinz, C., Ma, L., 1992. Modeling thetransport and retention of cadmium in soils: multireaction andmulticomponent approaches. Soil Sci. Soc. Am. J. 56,
1004–1015.
Seuntjens, P., Cornelis, C., De Brucker, N., Geuzens, P., 1999.
Derivation of functional layers in a podzol toposequence forsimulating cadmium transport. Phys. Chem. Earth 24 (7),869 –873.
Seuntjens, P., Mallants, D., Toride, N., Cornelis, C., Geuzens, P.,
2001a. Grid lysimeter study of steady state chloride transport in
two Spodosol types using TDR and wick samplers. J. Contam.
Hydrol. 51, 13– 39.
Seuntjens, P., Tirez, K., Šimůnek, M., van Genuchten, Th.,
Cornelis, C., Geuzens, P., 2001b. Aging effects on cadmiumtransport in undisturbed contaminated sandy soil columns.
J. Environ. Qual. 30, 1040–1050.
Shapiro, S.S., Wilk, M.B., 1965. An analysis of variance test fornormality (complete samples). Biometrika 52 (3 and 4),
591 –611.
Šimůnek, J., Huang K., van Genuchten, M.Th., 1998. The
HYDRUS code for simulating the one-dimensional movementof water, heat and multiple solutes in variably saturated media.
Research Report No. 144, US Salinity Laboratory, USDA,
Riverside, California, 142 pp.

200

P. Seuntjens et al. / Journal of Hydrology 264 (2002) 185–200

Soil Survey Staff, 1998. Keys to Soil Taxonomy. USDA, Natural
Resources Conservation Service, Seventh Edition.
Temminghoff, E., van der Zee, S.E.A.T.M., de Haan, F.A.M., 1995.
Speciation and calcium competition effects on cadmium
sorption by sandy soils at various pH levels. Eur. J. Soil Sci.
46, 649 –655.
Toride, N., Leij, F.J., 1996. Convective–dispersive stream tube
model for field-scale solute transport; II. Examples andcalibration. Soil Sci. Soc. Am. J. 60, 325– 361.
van der Zee, S.E.A.T.M., Van Riemsdijk, W.H., 1987. Transport ofreactive solutes in spatially variable soil systems. Water Resour.
Res 23 (11), 2059– 2069.
van Genuchten, M.Th., Cleary, R.W., 1982. Movement of solutes insoil: computer-simulated and laboratory results. In: Bolt, G.H.,(Ed.), Soil Chemistry, B. Physico-Chemical Models, Elsevier,
Amsterdam, pp. 349–386.
van Genuchten, M. Th., Leij, F.J., Yate, S.R., 1991. The RETC codefor quantifying the hydraulic functions of unsaturated soils. US
Salinity Laboratory, Riverside, California.
VMM, 1999. Air quality in Flanders 1998 (in Dutch), 178 pp.