Journal of Hydrology 259 (2002) 15±31

www.elsevier.com/locate/jhydrol

Calibration of Richards' and convection±dispersion equations to®eld-scale water ¯ow and solute transport under rainfall conditions
Diederik Jacques a,*, Jirka SÏimuÊnek b, Anthony Timmerman a, Jan Feyen aaInstitute for Land and Water Management, Katholieke Universiteit Leuven, Vital Decosterstraat 102, 3000 Leuven, BelgiumbGeorge E. Brown Jr Salinity Laboratory, USDA-ARS, 450 W. Big Springs Road, Riverside, CA 92507, USA
Received 21 July 2000; revised 30 October 2001; accepted 12 November 2001

Abstract
In this paper, the applicability of Richards' equation for water ¯ow and the convection±dispersion equation for solutetransport is evaluated to model ®eld-scale ¯ow and transport under natural boundary conditions by using detailed experimental
data and inverse optimization. The data consisted of depth-averaged time series of water content, pressure head and resident
solute concentration data measured several times a day during 384 d. In a ®rst approach, effective parameters are estimatedusing the time series for one depth and assuming a homogeneous soil pro®le. In a second approach, all time series were usedsimultaneously to estimate the parameters of a multi-layered soil pro®le. Water ¯ow was described by the Richards' equationand solute transport either by the equilibrium convection±dispersion (CDE) or the physical non-equilibrium convection±
dispersion (MIM) equation. To represent the dynamics of the water content and pressure head data, the multi-layered soil
pro®le approach gave better results. Fitted soil hydraulic parameters were comparable with parameters obtained with othermethods on the same soil. At larger depths, both the CDE- and MIM-models gave acceptable descriptions of the observed
breakthrough data, although the MIM performed somewhat better in the tailing part. Both models underestimated signi®cantly
the fast breakthrough. To describe the breakthrough curves at the ®rst depth, only the MIM with a mixing layer close to the soil
surface gave acceptable results. Starting from an initial value problem with solutes homogeneously distributed over the mobile
and immobile water phase was preferable compared to the incorporation of a small layer with only mobile water near the soil
surface. q 2002 Elsevier Science B.V. All rights reserved.
Keywords: Parameter optimization; Numerical modeling; Unsaturated zone; Water and solute transport

1. Introduction
The Richards' equation for water ¯ow and theconvection±dispersion equation (CDE) for solutetransport are frequently used to describe water andsolute transport experiments in undisturbed soilcolumns and ®eld soils (e.g. Schulin et al., 1987;* Corresponding author. Address: SCK.CEN, Studiecentrum voor
Kernenergie, Waste and Disposal, Boeretang 200, 2400 Mol,
Belgium. Tel.: 132-14-333-209; fax: 132-14-323-553.
E-mail address: djacques@sckcen.be (D. Jacques).

Butters and Jury, 1989; Roth et al., 1991; Flury,
1993; Radcliffe et al., 1996). Sometimes, however,signi®cant differences are observed between observations and the results of model calculations such as
higher measured concentrations at the beginning of
the solute breakthrough or tailing (e.g. Snow et al.,
1994; Jacques et al., 1998) which is especially thecase in structured or macroporous soils. An alternativeis the non-equilibrium convection±dispersion
equation to describe solute transport in such structured soils (e.g. Jacobsen et al., 1992; Mallants,
1996; Schoen et al., 1999; Vanderborght et al., 2000).

0022-1694/02/$ - see front matter q 2002 Elsevier Science B.V. All rights reserved.
PII: S 0022-169 4(01)00591-1

16

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

However, a throughout evaluation of models basedon Richards' equation and the (non-) equilibrium
convection±dispersion equation to describe water¯ow and solute transport under transient ®eld conditions is still limited, mainly because detailed experimental data at that scale are scarce. A strict evaluationof the processes described by a model should be doneusing independent experimentally obtained parameters. For some parameters, however, no techniques
are available to measure them at the ®eld-scale. Alternatively, models are calibrated using measured
variables (such as water content, solute concentration,¯uxes) to obtain a set of model parameters able toreproduce the measurements, also called inverse optimization. The calibrated model is then evaluated by
judging its ability to predict the dynamics of the states
variables for different initial and boundary conditions
or for other state variables (e.g. resident versus ¯ux
concentrations, Vanderborght et al., 2000). Whendetailed (in time and/or space) measurements of(several) variables are available, it is possible toperform an evaluation of a conceptual or mathematical model based on how good the model reproduces
the dynamics of the measurements and a physicalinterpretation of the ®tted parameters.
In general, a relatively large amount of model parameters has to be estimated since most soil pro®lesconsist of soil layers having different properties. For
each soil layer, the interaction between water content,pressure head, and unsaturated hydraulic conductivity
has to be determined. Instead of using a `subjective'
trial and error approach, more objective and/or automated algorithms are used to estimate soil hydraulicparameters in multi-layered soils (Lehman andAckerer, 1997; Kasteel, 1997; Abbaspour et al.,1999, 2000), solute transport parameters under natural
boundary conditions (e.g. Ventrella et al., 2000) andcoupled estimation of ¯ow and solute transport parameters for ground water ¯ow problems (e.g. Mishra
and Parker, 1989; Sun and Yeh, 1990; Medina andCarrera, 1996; Mayer and Huang, 1999). However,studies simultaneously estimating water ¯ow andsolute transport parameters for transient variablysaturated media are less common (e.g. Inoue et al.,2000), especially in a layered soil-pro®le and for®eld conditions.
The main objectives of the paper were to evaluate
the applicability of Richards' equation and the equi-

librium or non-equilibrium convection dispersion
equation to describe ®eld-scale water ¯ow and solutetransport in a layered soil pro®le under natural rainfall
boundary conditions and to determine simultaneously
soil hydraulic and solute transport parameters fordifferent horizons. The evaluation is based on the
combination detailed measurements of time series ofwater content, pressure heads, and solute concentration, an semi-automatic inverse optimization technique, and a numerical model.
2. Theory
2.1. Water ¯ow model
In this study, one-dimensional vertical water ¯ow
in variably-saturated porous medium is describedusing Richards' equation:


2u
2
2c
2K c
K c
1

2
2z
2z
2t
2z
where u is the volumetric water content [L 3L 23], t thetime [T], z the vertical coordinate (positive downwards) [L], c the pressure head [L], and K thehydraulic conductivity [LT 21]. The upper boundary
condition is de®ned here as an atmospheric condition.
Under this condition, water in®ltration corresponds tohourly rainfall rates as long as there is no ponding ofwater at the soil surface. Surplus of water is assumed
to be removed by run-off. Also, note that we neglect
evaporation since it was assumed that it is negligible
during the experiment due to the presence of a thinlayer of gravel (see Section 3). The lower boundary isde®ned as a zero pressure head gradient, i.e. freedrainage. The initial condition is given as a linear
interpolation of the measured pressure heads between
the ®ve depths at day 0.
We used the empirical model of van Genuchten(1980) to describe the water retention characteristic,u c :uu s 2 ur 
1 urn
1 1 uacu m

c # 0

2

where u r and u s [L 3L 23] are the residual and thesaturated water content, respectively, a [L 21], n,and m [±] are shape parameters. Combined withthe statistical pore-size distribution model for theD. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

unsaturated hydraulic conductivity of Mualem (1976)
and the restriction m  1 2 1=n; an analytical modelfor the hydraulic conductivity curve, K c; becomes(van Genuchten, 1980):
K c  K s

u1 2 aucu mn 1 1 aucu n 2m u2
1 1 aucu n ml3

c # 0
where Ks [LT 21] is the saturated hydraulicconductivity, and l the pore connectivity andtortuosity factor [±]. A preliminary analysis showedthat u r did not in¯uence the simulated time series of u
and c for the experimental conditions and was ®xed at
0.078 which is the average value for a loamy soil asgiven by Carsel and Parrish (1988). Therefore, thevector of unknown parameters of the soil hydraulicfunctions contains ®ve unknowns and is denoted asb h  {us ; a; n; Ks ; l}; in this paper.
2.2. Solute transport model
In this paper, the equilibrium convection±
dispersion equation (denoted as CDE throughout the
paper) represents a ®rst model, while the physicalnon-equilibrium convection±dispersion equation,
also called the mobile±immobile model (MIM),
represents a second model. The mathematicalrepresentation of the CDE for a non-sorbing, conservative solute is:


2uC
2
2C
2vuC

4
uD
2
2t
2z
2z
2z
where C is the solute concentration [ML 23], D thedispersion coef®cient [L 2T 21], and v the averagepore water velocity [LT 21]. The dispersion coef®cientis de®ned as (ignoring molecular diffusion):D  luvu

5

where l is the dispersivity [L]. The dispersivity isviewed as a material constant independent of the¯ow rate. Since v is obtained from the numerical solution of the water ¯ow model (the water ¯ux q [LT 21]
divided by u ), the vector of unknown parameters for
the CDE-model, bse, contains only l .
In the MIM-model, it is assumed that a part of theliquid phase is immobile and is not available forconvective transport. However, solutes can be

17

exchanged between the mobile and immobile regions
by means of a ®rst-order exchange process. The
mathematical formulation of the MIM for transientwater ¯ow for a non-sorbing, conservative trancer is(van Genuchten and Wierenga, 1976):
2um Cm
2u im Cim
2
2C
2vum Cm1um D m 2
2z
2t
2t
2z
2z
6
2uim Cim
 v Cm 2 Cim 
2t7

where u m and u im are the mobile and immobile watercontents, respectively, Cm and Cim the resident soluteconcentrations in the mobile and immobile regions,
respectively, and v the ®rst-order exchange coef®cient [T 21]. Note that v is the average pore water¯ux in the mobile region, thus v  q=um : Three
unknown parameters should be estimated for the
MIM transport model, i.e. bsne  {l; uim ; v}: The
upper solute transport boundary condition was de®ned
as a Cauchy boundary condition (SÏimuÊnek et al.,1998). Details are given below (see Section 3). A
zero solute concentration gradient was assumed as
the lower solute transport boundary condition. No
solute was present at the beginning of the simulation.
2.3. Formulation of the inverse problem
During the inverse procedure, the unknown parameter vectors are optimized by minimizing a given
objective function. The objective function is de®nedhere as (SÏimuÊnek et al., 1998):
( nj
)
m
i2
X
X h p
F g; buc 
8
vj
wij gj ti  2 gj ti ; buc
j1

i1

where nj is the number of observations for the jthmeasurement set, g the vector containing themeasured values (see below), buc the vector b containing the unknown parameters estimated given the ®xed
and known parameters in vector c, m the number ofparameter sets, gpj ti  the measured variable at time ti
for the jth measurement set, gj ti ; buc the corresponding model prediction, and vj and wij are the weights forthe jth measurement set and the individual observations, respectively. In this study, each wij is assumed

18

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

Table 1
Average soil textures measured at different depths in the trench
Depth (cm)

No of samples. 50 mm (%)

50±20 mm (%)

20±10 mm (%)

10±2 mm (%), 2 mm (%)

15
35
55
75
95

7
8
5
3
4

58.6
56.7
57.3
49.6
43.8

19.3
18.9
17.6
21.2
30.03

6.4
7.8
6.6
7.9
7.4

4.5
3.2
3.7
4.4
4.5

11.1
13.3
14.8
17.4
14

to be equal to 1 and vj is given by:
vj 

1
nj s j2

9

where s j2 is the measurement variances for the jthmeasurement set. This formulation of the inverse
problem is similar to those of SÏimuÊnek et al. (1999,
2000) and Inoue et al. (2000). The Hydrus-1D code(SÏimuÊnek et al., 1998) was used for the numericalsolution of the governing equations and the non-linear
parameter optimization. The algorithm proposed by
Marquardt (1963) is used here to minimize theobjective function (Eq. (8)).
3. Materials and methods
3.1. Site description and experimental design
The experimental ®eld is located at Bekkevoort,
Belgium. It is situated on a gentle slope and is covered
with a meadow. The soil is classi®ed as a Eutric
Regosol. Typically, the top 100 cm show three soil
horizons: an Ap horizon between 0 and 25 cm, a C1
horizon between 25 and 55 cm and a C2 horizon

between 55 and 100 cm (Mallants, 1996). Dye studies
revealed the occurrence of macropores throughout
the soil pro®le, mainly from biological origin(Vanderborght et al., 2000). Previous studies characterized the soil hydraulic and solute transport parameters on different scales and under different ¯owconditions (Mallants et al., 1996; Vanderborght,
1997; Jacques et al., 1997, 1998).
A trench, 1.2 m deep and 12 m long, was dug at the®eld site. Textural analyses were done at different
depths in the trench (Table 1). Next to the trench,
the vegetation was removed from the soil surface
over an area of 8 £ 2 m 2 and the surface was leveled.
An overview of the experimental set-up and the datacollected is given in Table 2. One side of the trench
was covered using a plastic sheet as a barrier tohydraulically isolate the disturbed trench zone.
Along a transect of 5.5 m, two-rod TDR-probes(25 cm long, 0.5 cm rod diameter, 2.5 cm rod spacing)were installed each 50 cm (12 locations) at 5 depths(15, 35, 55, 75, and 95 cm deep). Tension cups (6 mmdiameter, 25 mm long) were installed at a horizontal
distance of 10 cm from each TDR-probe. Temperature probes (6 cm long) were installed at ®ve depths
on three locations at 1.3, 2.7, and 4.2 m from the ®rst

Table 2
Overview of measured variables, used experimental devices, and frequency of measurements (A: automated, M: manual, d: number of depths, p:number of positions)
Quantity

DeviceFrequency

MethodTotal positions

Water content, uResident solute concentration,
Cr
Pressure head, cTemperature, T
Water ¯uxes
Solute ¯uxes
Rainfall

TDR-probes
TDR-probes

Every 2 hEvery 2 hAA

5d £ 12p
5d £ 12p

Tension cups
Temperature probes
PCAPS
PCAPS
PluviographEvery hEvery hEvery 2/3 dEvery 2/3 dContinuously

A
A
M
M
A

5d £ 12p
5d £ 3p
2d £ 3p
2d £ 3p
1

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

probe. The travel time of an electromagnetic wave
along the lines of the TDR-probe and the re¯ectioncoef®cient were automatically measured every 2 husing the system described by Heimovaara and de
Water (1991). Each probe was calibrated using thetechniques of Heimovaara et al. (1995) to obtain theapparent dielectric constant, Ka, and the bulk soil
electrical conductivity, ECa (dS m 21), from themeasurements. The water content, u , was calculated
from Ka using a site-speci®c calibration curve(Jacques et al., 1999). The resident solute concentration was determined from converting ECa to theelectrical conductivity of the soil water, ECw, by theempirical model of Vogeler et al. (1996) (see Jacques,2000 for details). Finally, ECw [dS m 21] wasconverted to C r using a linear relationship with aslope of 0.667 and an intercept of 20.52 (Jacques etal., 1998). Pressure heads were automaticallyrecorded using pressure transducers in the set-up of
Bertuzzi et al. (1996). Rainfall was continuously
measured with a rainfall recorder (200 cm 2) with a¯oated pen system on a paper (0.1 mm intervals,
rotation speed of 1 cm h 21).
After all devices were installed, the trench wasre®lled and a thin layer of gravel (between 1 and2 cm) was put on the experimental area to decrease
the erosive effect of raindrops on the bare soil surface,minimize evaporation and decrease the weed growthwhich were regularly removed. Measurements werestarted on March 11, 1998, de®ned as day 0. On an
area of 16 m 2, 6000 g of a tracer (CaCl2´2H2O) was
manually applied on August 28, 1998 (day 168.5)
over a period of 3.5 h using a water depth of 0.5 cm
with a total concentration of 75 g l 21. The experiment
ended on March 31, 1999 (day 384). A throughout
discussion of the experimental results and thequalitative interpretation was given in Jacques (2000).
3.2. Overview of measurement sets used in the inverse
optimization
We used four different measurement sets in the
inverse optimization: gu;d ; gc;d ; gq ; and gc;d ; where(1) gu;d is a vector containing the water content data
averaged by depth for the dth depth (d  1 for the
15 cm depth, d  2 for the 35 cm depth,¼, d  5
for the 95 cm depth). Only four daily observations(at 0, 6, 12, and 18 h) were included in the objective19

function. (2) gc;d is a vector containing the pressure
head data averaged by depth for the dth depth and four
observations a day. (3) gq is a vector containing thehourly rainfall rates that are smaller then 0.5 cm h 21.
The decision to include gq in F was based on preliminary optimizations (results not shown). The
upper boundary was de®ned as hourly rainfall rates
as de®ned above. In preliminary simulations with nosurface run-off (thus with surface ponding), a build-up
of water of about 10 cm height during some heavy
rainfall periods was simulated. However, ponding ofwater was not observed at the experimental ®eld. On
the other hand, only a small part of the precipitationwater was actually in®ltrating into the soil in simulations with surface run-off (and no g q in F). The occurrence of excessive run-off was not observed at the®eld. By including g q in F, the optimization will
search for parameters for which nearly all waterin®ltrates during moderate rainfall (q , 0.5 cm h 21).
(4) g c;d is a vector containing the time-integralnormalized
resident concentrations averaged bypdepth, Cdrt t [T 21] de®ned as (Vanderborght et al.,
1996; Jacques et al., 1998):

u tCdr t
Cdrtp t  Z1 d
ud uCdr udu100where u d and Cdr are the water content and the resident
solute concentration averaged by depth, respectively.
For the transient ¯ow problem in this study, a sensitivity analysis has shown that the effect of the calibration coef®cients for ECw on Cdrtp t is very small andsmaller than when Cdr t is used.
3.3. Inverse procedures
Two different approaches were used for the estimation of the water ¯ow and solute transport parameters.
In the ®rst approach (Section 4.1), parameters are
estimated using time series of observations for a particular depth assuming a homogeneous soil pro®le.
Parameters are interpreted as effective parameters
representative for an equivalent homogeneous porous
medium between the soil surface and a particulardepth. Thus, we seek to minimize, for example, theobjective function F g u;d ; gc;d ; bhde ; i.e. to optimize
the effective (subscript e) water ¯ow parameters of thedth depth using time series of the averaged u and c

20

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

for the dth depth. In the second approach (Section4.2), a layered soil pro®le with different parameters
for each layer is assumed. The layers are de®ned
between two observation depths. For example, theobjective function F g u;d ; gc;d ; gu;d21 ; gc;d21 ; gu;d22 ;gc;d22 ; bhd ubhd21 ; bhd22  leads to the estimation of
the water ¯ow parameters for the layer between the(d 2 1)th and dth depth. The time series of averaged u
and c for the (d 2 2)th, (d 2 1)th and dth depth areused in this inversion procedure and parameters of theupper two layers are assumed to be previously
estimated.
In both approaches, parameters are estimated in two
steps. In approach 1, bhde is ®rst estimated using gu ,d,
gc ,d and gq, followed by the estimation of the solutetransport parameters (bse or bsne) using gc,d and bhde. A
similar sequential approach is used in the second
approach in which the different layers are subsequently added to a soil pro®le of 150 cm deep fromthe top to the bottom layer. Parameters for particularlayers, bhd are estimated in a sequence starting fromthe top layer. When all sets of bh,d are estimated, somefurther tuning is done by re-estimating combinationsof bh,d from different layers. When all hydraulic parameters are estimated, the same stepwise procedure isused to estimate the solute transport parameters.
Although Mishra and Parker (1989) recommended
to simultaneously estimate the soil hydraulic andsolute transport parameters (see Inoue et al., 2000),
we used the sequential approach to minimize

problems of non-uniqueness, especially for approach
2 having up to 36 unknown parameters. To test
uniquenessoftheoptimized
parameters,
F g u;d ; gc;d ; bhde  for d  1; 2, and 3 were minimized
for different initial parameter sets. Final estimates
were close to each other and their 95% con®dence
interval overlapped (results not shown). Since solutetransport is in¯uenced by water ¯ow and thus by thesoil hydraulic properties (Medina and Carrera, 1996),
hydraulic and solute transport parameters werealso optimized simultaneously by minimizing
F g u;d ; gc;d ; gc;d ; bhde ; bsede  with d  2 and 3 using
the previously optimized parameters from the sequential approach as initial parameters. No signi®cantdifferences in parameter estimates were foundbetween the simultaneous and sequential approach(results not shown). Based on these ®ndings, weused the sequential approach and only one initial
parameter set for all optimization problems.

4. Results
4.1. Approach 1: equivalent homogeneous soil pro®leFirst, we discuss the effective hydraulic parameters
estimated by minimizing F g u;d ; gc;d ; gq ; bhde : Results
are shown in Figs. 1±4 and Tables 3±5. Since watercontents and pressure heads in the ®fth depth showed
little temporal variability, we used the time series ofFig. 1. Measured cumulative rainfall and simulated cumulative actual in®ltration for the different objective functions as de®ned in Table 3.

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±3121Fig. 2. Measured water contents, and water contents simulated using the equivalent homogeneous approach (approach 1, solid line), and multilayered soil approach (approach 2, Fall in Table 3, dashed line) for (a) 15 (F1 in Table 3 for approach 1), (b) 35 (F2 in Table 3 for approach 1),(c) 55 (F3 in Table 3 for approach 1), (d) 75 (F4 in Table 3 for approach 1), and (e) 95 (F5 in Table 3 for approach 1) cm depth.

the fourth and ®fth depth simultaneously to estimate
the parameters for the ®fth depth (Table 3). There is
only a small difference between the parameters of thedifferent depths. A detailed discussion of the
optimized parameters will be given in Section 4.2.

For the ®ve optimizations, the total simulated
cumulative in®ltration is only between 50.2 and
83.8% of the total amount of rain (Table 3). Fig. 1
shows that the largest difference occurred during theheavy rainfall on day 185, and smaller differences

22

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31Fig. 3. Measured pressure heads, and pressure heads simulated using the equivalent homogeneous approach (approach 1, solid line), and multilayered soil approach (approach 2, Fall in Table 3, dashed line) for (a) 15 (F1 in Table 3 for approach 1), (b) 35 (F2 in Table 3 for approach 1),(c) 55 (F3 in Table 3 for approach 1), (d) 75 (F4 in Table 3 for approach 1), and (e) 95 (F5 in Table 3 for approach 1) cm depth.

during moderate rainfall conditions (e.g. between
days 100 and 150, and 250 and 350, except for
F g u;1 ; gc;1 gq ; bh1e : These differences may be due
to either surface run-off across the slope, either in®ltration in the soil in macropores, or a combination of

both. It was not possible to discriminate between bothprocesses because run-off was not measured (although
no large run-off ¯ow was observed during the experiment), and measurements of water content andpressure head generally do not indicate preferential

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31¯ow. However, macropore in®ltration might be an
important process. During the same experiment,water and solute ¯uxes were measured at two depths
and three locations with passive capillary wicksamplers. Qualitative interpretation of these measurements indicated the occurrence of signi®cant macropore water ¯ow as was discussed in Jacques (2000).
The dynamics of water content and pressure head timeseries gradually decreased with depth, but the effect ofheavy rainfalls was still clearly observed at deeper
depths by a sudden increase in u after heavy rain
events. A low conductivity is needed to model theslow water ¯ow dynamics at deeper depths, whereasa high conductivity is required to model the fast
responses after heavy rainstorms. The model has
solved this dilemma by generating a large amount of
run-off, even during moderate rainfall conditions.
Simulations very well reproduced the dynamics of
the water content with time, especially between days
80 and 140 and during the last 100 days of the experiment (Fig. 2). At the fourth and ®fth depths, the calculated water contents overpredicted the dynamics ofthe observed water contents, especially between
days 25 and 150. Furthermore, the dynamics of thecalculated water contents are somewhat retarded
compared with the observations during the last
150 d of the experiment. Visual inspection of Fig. 3
shows that pressure head simulations are less goodthan the ones for the water content. However, pressure
heads that were measured in smaller soil volumes may
be associated with larger uncertainties than watercontents that are measured in larger soil volumes
using TDR-probes. In addition, external factors suchas variation in temperature and the long capillarytubes between suction cup and pressure transducer
may result in larger uncertainties of the pressure
head measurements. Calculated pressure heads
decreased less compared to observations during thedry period in the summer (day 140±160). Althougha thin layer of gravel was put on the soil surface, somewater loss due to evaporation may have caused thisdiscrepancy between measurements and observations.
Note that the difference between measured and simulated water contents during this period was lesspronounced. At the second depth, pressure heads are
signi®cantly overpredicted between days 150 and 200.
This may be rather due to measurement errors anduncertainty than due to an incorrect model. Water

23

contents during this period were as large as aroundday 100 or during the last 100 d of the experiment.
However, during these periods, the pressure heads
were larger than between days 150 and 200.
Subsequently, the transport parameters were determined using the solute breakthrough data. We will
®rst discuss the results for the CDE- and MIM-models
for the second to ®fth depth (Fig. 4(b)±(e)). For thesecond, third and fourth depths, the MIM-model
describes the observed breakthrough curves better
than the CDE-model, especially in the tailing part.
In contrast to the simulations presented here, one
usually expects that the tailing of the MIM-model
will be more pronounced than that of the CDEmodel. The optimization algorithm may have resulted
in an overestimation of the dispersivity values for the
CDE-model in order to describe the fast breakthrough.
Consequently, the tailing of the breakthrough curve isalso overestimated. At the ®fth depth, the MIM-model
simulated better the early breakthrough but overpredicted the tailing compared to the CDE-model.
Fitted effective solute transport parameters are
given in Table 5. Dispersivity values increase signi®cantly with depth, which was also observed forsteady-state lysimeter (Vanderborght et al., 2000)
and in situ (Jacques et al., 1998), and transientlysimeter (Vanderborght et al., 2000) solute transport
experiments on the same soil. Immobile watercontents were high at the second and ®fth depth
compared to small values at the other two depths.
Extremely low v values were obtained for these two
depths compared to values reported in the literature(see compiled values in Vanderborght et al., 1997).
We could not describe the solute breakthroughcurve at the ®rst depth using the CDE or the MIMmodel. The `best' ®t of the MIM-model showed a verylarge peak just after solute application (with a maximum C rtp of 0.195) followed by an extreme long tailing (with l  643 cm; F g C;1 ; bsne1e ubh1e   0:821;
results not shown). A major reason for this failure
may be the difference between the actual in®ltration
processes at the soil surface during solute applicationand the processes described in the MIM-model. The
soil consists of different ¯ow domains that can be
subdivided into two groups: (i) preferential ¯owpaths, mainly macropores in this soil, and (ii) the
soil matrix consisting of domains of slowly moving
water and immobile water (see Jacques, 2000). In our

24

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31Fig. 4. Time-integral-normalized resident concentrations measured and simulated for (a) 15, (b) 35, (c) 55, (d) 75, and (e) 95 cm depth for theobjective functions as de®ned in Table 5.

®eld experiment, solutes were applied after a long dry
period, and it is reasonable to assume that water
entered both the soil matrix due to sorptivity effects
and in the preferential ¯ow domain consisting mainlyof macropores (water may ¯ow into macropores underno ponding conditions as shown by Godrathi et al.
(1999)). Water in the preferential ¯ow domainreached the 15 cm depth immediately after solute
application. A second peak in the breakthroughcurve indicated the solute arrival in the slowly ¯owing

water in the soil matrix, while tailing may indicate
diffusive solute exchange between mobile andimmobile water (Jacques, 2000). All solute appliedat the soil surface will move only into the mobile
water phase for the MIM-model. Only a small amountof solute will be exchanged to the immobile watercontent since the transport distance is small and theexchange time short. A dual- or multi-porosity transport model (e.g. Jarvis et al., 1991; Gerke and van
Genuchten, 1993; Steenhuis et al., 1990; Skopp andD. Jacques et al. / Journal of Hydrology 259 (2002) 15±3125hydraulic properties as the rest of the soil pro®le,but with only mobile water (i.e. CDE). The dispersivity of this ®rst mixing layer and parameters of
the MIM-model for the other layer were estimatedusing the observed breakthrough data at the 15 cm
depth. A value of 2.2 was obtained for l of the mixing
layer (see Table 5 for other parameters). Although thevalue of the objective function for the two-layered
system of the ®rst depth is the largest, an acceptable
description of the observed breakthrough curve isobtained (Fig. 4(a)) given the complex transportmechanisms near the soil surface.
In conclusion for the estimated effective parameters: (i) time series of water contents, pressure
heads and resident solute concentrations are reasonably well described, (ii) the MIM-model described thepeak and tailing parts of the observed solute breakthrough curves better than the CDE-model, and (iii)
goodness-of-®t of the BTCs showed a decreasing
trend with depth (except for the breakthrough curveat the ®rst depth).

Fig. 4. (continued)

Gardner, 1992) may be more appropriate to represent
the transport processes near the soil surface.
Alternatively, a mixing zone near the soil surface,in which solutes are distributed over the different ¯owregions, may mimic these processes. We de®ned a
5 cm thick layer at the soil surface with similar

4.2. Approach 2: multi-layered soil pro®leBased on the results of the ®rst approach, four
layers for hydraulic properties (fourth and ®fth

Table 3
Summary of statistical measures of goodness-of-®t for time series of water contents, pressure heads and precipitation rates (MSEÐmean
squared error). F is the objective function as de®ned in Eq. (8) for the three measurement sets (water content, pressure head, and rainfall rates).
Fu , Fc , and Fq are the objective functions as de®ned in Eq. (8) using only one measurement set j  1; (i.e. water content, pressure head, andrainfall rates, respectively) calculated for the parameters obtained by minimizing FF
Approach 1: equivalent homogeneous soil pro®leF g u;1 ; gc;1 ; gq ; bhle   F 1
1.287
F g u;2 ; gc;2 ; gq ; bh2e   F 2
1.984
F g u;3 ; gc;3 ; gq ; bh3e   F 3
2.415
2.543
F g u;4 ; gc;4 ; gq ; bh4e   F 4
F g u;4 ; gc;4 ; gu;5 ; gc;5 ; gq ; bh5e   F 5
2.925
Approach 2: layered soil pro®leF all b F all
15 cm depth
35 cm depth
55 cm depth
75 cm depth
95 cm depth
a
b

Fu

Fc

Fq

R2

%TF a

MSEu
( £ 10 25)

MSEc

0.656
0.607
0.812
1.05
1.385

0.542
0.717
0.725
0.992
1.073

0.089
0.66
0.877
0.5
0.467

0.86
0.854
0.813
0.913
0.748

83.8
60.1
50.2
66.4
67.4

12.5
4.1
3.98
6.21
3.79

317.9
205.7
45.4
31.1
79.2

11.6
3.85
4.1
4.21
3.46

350.6
243.3
42.4
27.9
120.1

78.6

The ratio of the total simulated cumulative ¯ux at the top boundary over the cumulative amount of rain (109 cm) multiplied by 100.
The total objective function of the second approach includes all three-time series. All hydraulic parameters are estimated simultaneously.

26

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

Table 4
Estimated hydraulic parameters based on time series of water contents, pressure heads and precipitation ratesa (cm 21)n

Ks (cm d 21)

Approach 1: equivalent homogeneous soil pro®le0.390
F g u;1 ; gc;1 ; gq ; bh1e 
F g u;2 ; gc;2 ; gq ; bh2e 
0.369
F g u;3 ; gc;3 ; gq ; bh3e 
0.378
F g u;4 ; gc;4 ; gq ; bh4e 
0.399
0.399
F g u;4 ; gc;4 ; gu;5 ; gc;5 ; gq ; bh5e 

0.01730.01370.02510.02480.0355

1.43
1.38
1.41
1.435
1.435

5.18
1.64
1.31
2.31
2.55

4.70
0.13
5.62
5.25
5.25

Approach 2: multi-layered soil pro®leF(all) a
15 cm depth
35 cm depth
55 cm depth
75 1 95 cm depth

0.01370.01370.02460.0219

1.56
1.55
3.35
1.57

3.89
1.55
1.22
2.35

0.395
0.12
0.001
15.28

us

a

0.378
0.371
0.37
0.391

l

The total objective function of the second approach includes all three-time series. All hydraulic parameters are estimated simultaneously.

depth were combined in one layer) and ®ve layers fortransport properties were de®ned. In addition, we did
not use the CDE-model, and de®ned a 5 cm thickmixing layer containing only mobile water near thesoil surface. Thus, 20 soil hydraulic and 16 solutetransport parameters were estimated (approach 2a).
Alternatively, the mixing layer was de®ned as an
initial value problem (approach 2b). Solutes were
evenly distributed over the immobile and mobile
water content to a depth of 1.44 cm (  0.5 cm applied
water divided by the average water content in that
depth interval, 0.347 cm 3 cm 23) and only 15 transportparameters were estimated. Results are summarized inTables 3±5 and Figs. 1±4.
The total actual in®ltration at the soil surface isclose to the value obtained by minimizing
F g u;1 ; gc;1 ; gq ; bh1e  and is 85.7 cm. This value issubstantially larger than the total actual in®ltration
obtained for the other four objective functions used
for estimating effective parameters for deeper depths.
As discussed above, in order to represent the decreasein the water ¯ow dynamics with depth in the equivalent homogeneous approach, a large part of the rainfall did not in®ltrate in the soil. The multi-layered soil
pro®le had a dampening effect on the water ¯owreducing the dynamics deeper in the soil.
Visual inspection of Figs. 2 and 3 shows thatapproach 1 described the observed time series somewhat better than approach 2 (see also MSE in Table3). There are only small differences between the two
approaches in terms of the MSE. MSEc are larger forthe ®rst, second and the last depth for approach 2
compared to approach 1. Except for the third depth,
approach 2 performed somewhat better than approach
1 for water contents based on MSEu .
A quantitative comparison between the threemethods to estimate the MIM-model parameters(approach 1, 2a, and 2b) is based on the MSE forthe ®ve individual observation depths (Table 5). The
results from approach 1 are worse than those obtained
for approaches 2a and 2b. Approach 2b performed
better at all depths compared to approach 2a.
However, the difference between the 2a and 2b
decreased with depth: the MSE of 2b is 38.5, 55.1,
77.1, 82.7, and 98.6% of that of 2a. The two
approaches used to mimic the complex ¯ow patterns
in the topsoil have the largest impact on the simulations of solute breakthrough in the upper horizons
with an obvious preference for approach 2b.
However, the better simulation results obtainedwith approach 2 are not consistent with depth, as
can be seen in Fig. 4. For example, approach 2b
obviously resulted in a better description of the tailingpart of the breakthrough curve at the ®rst depth or of
the solute peak between days 220 and 240 for thesecond depth. On the other hand, approach 2a and2b failed to describe C rtp between days 185 and 210
in contrast to approach 1 (MIM-model). At the otherthree depths, there are only small differences between
approaches 2a and 2b, although 2b described themeasured concentrations between days 185 and 210
somewhat better, especially at the third depth.

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

27

Table 5
Summary of statistical measures of goodness-of-®t and estimated solute transport parameters for time series of resident solute concentrations(MSEÐmean squared error, naÐnot available)
R2

MSEC ( £ 10 26)

l (cm)

u im (cm 3 cm 23)v (d 21)

Approach 1: equivalent homogeneous soil pro®leCDE-modelna
F g C;1 ; bse1e ubh1e   F 1 CDE0.146
F g C;2 ; bse2e ubh2e   F 2 CDEF g C;3 ; bse3e ubh3e   F 3 CDE0.215
F g C;4 ; bse4e ubh4e   F 4 CDE
0.376
F g C;5 ; bse5e ubh5e   F 5 CDE
0.314

na
0.854
0.785
0.648
0.794

na
3.22
1.11
1.83
0.84

na
8.45
19.96
15.84
21.72

±
±
±
±
±

±
±
±
±
±

MIM-model
F g C;1 ; bsne1e ubh1e  a F 1 MIMF g C;2 ; bsne2e ubh2e   F 2 MIMF g C;3 ; bsne3e ubh3e   F 3 MIMF g C;3 ; bsne3e ubh3e   F 3 MIMF g C;4 ; bsne4e ubh4e   F 4 MIM

0.697
0.958
0.878
0.806
0.755

19.7
1.1
0.82
1.35
0.69

3.227
1.98
11.05
10.6
22.65

0.32
0.32
0.053
0.048
0.325

0.019
0.059
3.0 10 27
4.7 10 27
0.038

12.1
3.01
0.51
1.22
0.64

0.067
8.74
7.65
53.9
3.28

0.32
0.32
0.32
0.32
0.32

0.05
0.089
0.045
0.028
0.027

4.66
1.66
0.39
1.01
0.62

0.433
5.66
6.88
35.45
2.99

0.327
0.322
0.3033
0.32
0.32

0.12
0.086
0.027
0.027
0.02

F

0.306
0.05
0.158
0.279
0.260

Approach 2: multi-layered soil pro®leMIM-model; approach 2a0.1514F(´)  Fall(MIM, approach 2a)
15 cm depth
35 cm depth
55 cm depth
75 cm depth
95 cm depth
MIM-model: approach 2bF(´)  Fall(MIM, approach 2b)
15 cm depth
35 cm depth
55 cm depth
75 cm depth
95 cm depth
a

0.078

0.849

0.922

Estimated with a 5 cm thick surface layer consisting of only mobile water (i.e. CDE).

Compared to approach 1, approaches 2a and 2b
predicted earlier breakthrough, although there is stilla signi®cant difference with the observed early breakthrough. The tailing part is well described by all
methods using the MIM-model. The opposite isobserved for the ®fth depth: approach 1 described
better the early breakthrough, whereas the approaches
2a and 2b ®tted better the peak and tailing of the
breakthrough.
4.3. Comparison of estimated parameters withprevious parameter values
The a -coef®cient of Eqs. (2) and (3) showed anincreasing trend with depth, while n is constant withdepth, except at the 55 cm depth for the layered soil
pro®le (Table 4). At that depth, the high n value issomewhat compensated by the low l-parameter. The
estimated values, reported in Table 4, are comparedwith parameter values determined in earlier studies on
the same soil (Table 6). The equilibrium desorption
u 2 c measurements resulted in smaller a -values andsimilar n-values. The contrary was observed for thetransient drainage experiment: similar a -values(except between 5 and 15 cm), but smaller n-values.
Furthermore, a and n parameter values were different
for K c based on the drainage experiment. The
results of Mallants et al. (1997a) indicate using thesame parameters in Eqs. (2) and (3) may be too
restrictive for modeling transient water ¯ow in soils.

28

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

Table 6
Overview of parameters of the water retention characteristic and the unsaturated hydraulic conductivity described with Eqs. (2) and (3),respectively, determined in previous experiments on the same soil (upper and lower positions refers to depths in the pro®le for which parameters
were estimated)
Experimental information
Methodu (c ) Desorption technique

Parameter values

Reference

Diameter Height Upper Lower Parameter Mean Minimum Maximum(cm)(cm) position position
(cm)
(cm)
5

5

10

90

Transient drainage experiment 30

100

15
60
90
±
90
±
±

K(c ) Instantaneous pro®le method

30

100

5
30
5
0

Constant head method

5
20
30

5
20
100

10
0
0

The Ks-values measured on the lysimeters of 100 cm
height are comparable with the ®tted ones from our®eld data. Note the different boundary conditions used
to estimate Ks: a laboratory method with ponding
conditions for the former in contrast to an in situ
method with natural boundary for the latter.
The estimated immobile water content for
approach 2 is high in each layer (Table 5). Vanderborght et al. (2000) also found small mobile watercontents in their transient ¯ow experiments on thesame soil. In contrast to our study, they noticed adecrease in u m with depth. For approach 2b, transfer
coef®cients, v , show a decreasing trend with depth
and are in the same order as values reported for thesame soil for steady-state saturated solute transport(0.0384±0.744 d 21, Mallants et al., 1996). These aresmaller than the value (0.89 d 21) obtained for thetransient ¯ow experiments of Vanderborght et al.
(2000). The dispersivity of the fourth depth is large
compared to values derived for the other depths. Thisis probably due to an artifact, since the optimization
algorithm, based on minimizing the sum of squared
errors, tried to describe both the fast breakthroughand the long tailing of the observed breakthroughcurve. The small dispersivity value of the ®fth
depth is due to a combination of the high dispersivity
value of the fourth depth and the incomplete
measurement of the tailing part of the breakthroughcurve.

a
n
a
a
n
a
n
Ks
Ks
Ks

0.0038
1.311
0.103
0.017
1.05
0.038
2.63

333.6
126.7
11.8 0.96

0.013
1.535
0.147
0.031
1.14
0.065
2.85
12.7

Mallants et al., 1997a
Mallants et al., 1997a
Mallants et al., 1997a
Mallants et al., 1997a
Mallants et al., 1997a
Mallants et al., 1997a
Mallants et al., 1997a
Mallants et al., 1997b
Mallants et al., 1997b
Mallants et al., 1997b

5. Final comments and conclusions
Although modeling based on inverse optimization
could describe some of the main features of theexperimental data, and although estimated parameters
are within the range of earlier determined ¯ow andtransport parameters, this study also shown that the
combination of Richards' equation and the MIMmodel failed to describe important features such as`internal run-off' and early solute breakthroughunder natural ¯ow conditions in a structured soil.
Some model improvements are proposed here,
which are to be considered in further studies. A
unimodal description (e.g. Eqs. (2) and (3)) of thehydraulic characteristics may not describe the fast
increase in hydraulic conductivity near saturation
when macropores or structural elements are present
in the soil. When a medium with multiple ¯owdomains (e.g. matrix and macropores) is represented
by an equivalent continuum (van Genuchten et al.,1999), the hydraulic characteristics can be describedby multi-modal functions (e.g. Durner, 1994;
Mohanty et al., 1997, 1998). The hydraulic conductivity then increases drastically near saturation
without a large increase in pressure head or watercontent resulting in a larger amount of actual in®ltration under wet condition or during heavy rainfall as
observed in our study.
Transient ¯ow processes may signi®cantly

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

complicate the parametric description of the solutetransport processes and result in ¯ow-dependency ofsolute transport parameters. Vanderborght et al.
(2001) compiled the solute transport parameters fordifferent tracer experiments in key Belgian soil types.
They found that longitudinal dispersivity increased
with increasing water ¯ux Jw, implying that dispersivity cannot be viewed entirely as a material
constant, but is also dependent on the ¯ow conditions.
For the soil type used in this study, Vanderborght et al.
(2000) found that for conditions close to saturation
this relation is l  3:44J w0:83 : Based on the relations
between dispersivity, depth and ¯ow rates at one hand
and the length of the pores at the other hand (Skopp
and Gardner, 1992), Vanderborght et al. (2000)
concluded that long continuous pores were activatedwhen ¯ow rates increased. The experimental studies
of Reedy et al. (1996) and Bajracharya and Barry(1997) indicated that v may also increase with ¯owvelocity. A modeling study of ZurmuÈhl and Durner(1996) showed that treating u im as an dynamic
variable may signi®cantly in¯uence the simulationof solute transport through structured soils.
The MIM-model did not represent the complex¯ow and transport mechanisms, especially near thesoil surface (Fig. 4). Snow et al. (1994) also observedthat solute transport near the soil surface was different
from transport processes deeper in the soil. They
stated that applied solutes move quickly into theimmobile water and then slowly diffuse back to themobile water. This process was mimicked here byapproach 2b. An alternative is to use dual-porosity
models assuming advective transport in both ¯owdomains (e.g. Jarvis et al., 1991; Gerke and van
Genuchten, 1993). A second problem is that theMIM-model assumes that solutes are uniformly
distributed in each water phase. As was experimentally demonstrated by Ghodrati et al. (1999), only asmall part of the matrix domain interacts with thepreferential ¯ow domain, also called the active matrix
layer (Wallach and Steenhuis, 1998). Furthermore,lateral mixing processes between different ¯owdomains change with depth (FluÈhler et al., 1996;
Flury et al., 1994; de Rooij, 1996). Multi-region transport models (e.g. Hutson and Wagenet, 1995) may beneeded to describe the complex lateral mixing
processes between macropores, active matrix layer,
and mobile and immobile regions.

29

As a general conclusion, observed time series ofwater content, pressure head and solute concentration
could reasonable be described by applying the
Richards' equation for water ¯ow and the nonequilibrium convection±disperison equation forsolute transport given a layered soil pro®le andde®ning a kind of mixing layer for solute near thesoil surface. We applied successfully the sequential
parameter estimation approach to estimate the large
number of hydraulic and solute transport parameters
needed to describe the ®eld-scale water ¯ow andsolute transport under transient ¯ow conditions.
However, the applied models failed to describe
some speci®c features such as macropore water ¯ow(failure of Richards' equation) and mixing of solutes
between different ¯ow domains (failure of theconvection±dispersion equation).
Acknowledgements
This study was done with ®nancial support of a
Research Grant of the Fund for Scienti®c ResearchFlanders (Belgium) (F.W.O.- Vlaanderen). The ®rst
author acknowledges the ®nancial support of ascholarship of the Flemish Institute for the
EncouragementofScienti®c±TechnologicalResearch in the Industry (IWT). The inverse optimization was performed at the US Salinity
Laboratory, USDA-ARS, made possible through a
travel grant of the Fund for Scienti®c ResearchÐ
Flanders (Belgium) (F.W.O.-Vlaanderen) and a
grant of the National Science FoundationÐUSA.
We acknowledge the family Deckers, owners of theexperimental ®eld, for the permission to perform theseexperiments on their property. The assistance of J.
Schaerlaekens, F. Serneels, and J. Vanderborght inthe ®eld and the laboratory is highly appreciated.
The comments of the anonymous reviewer to improve
the overall quality of the paper are also highly
appreciated.
References
Abbaspour, K., Sonnleitner, M., Schulin, R., 1999. Uncertainty inestimation of soil hydraulic parameter by inverse modeling:example lysimeter experiments. Soil Sci. Soc. Am. J. 63,
501±509.

30

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31

Abbaspour, K., Kasteel, R., Schulin, R., 2000. Inverse parameter
estimation in a layered unsaturated ®eld soil. Soil Sci. 165, 109±
123.
Bajracharya, K., Barry, D.A., 1997. Nonequilibrium solute transportparameters and their physical signi®cance: numerical andexperimental results. J. Cont. Hydrol. 24, 185±204.
Bertuzzi, P., Gauda, J.C., Bourlet, M., Mohrath, D., Chanzy, A.,
1996. Methode de Wind: Notice Technique. Institute Nationalde Recherche Agronomique, UniteÂ du Sol d'Avignon, France,
54 pp.
Butters, G.L., Jury, W.A., 1989. Field scale transport of bromide in
an unsaturated soil, 2. Dispersion modeling. Water Resour. Res.
25, 1583±1589.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability
distributions of soil water retention characteristics. Water
Resour. Res. 24, 755±769.
de Rooij, G.H., 1996. Preferential ¯ow in water-repellent sandy
soils: model development and lysimeter experiments. PhD
Dissertation, University of Wageningen, The Netherlands, 229
pp.
Durner, W., 1994. Hydraulic conductivity estimation for soils withheterogeneous pore structure. Water Resour. Res. 30, 211±223.
FluÈhler, H., Durner, W., Flury, M., 1996. Lateral solute mixing
processesÐa key for understanding ®eld-scale transport ofwater and solutes. Geoderma 70, 165±183.
Flury, M., 1993. Transport of bromide and chloride in a sandy and aloamy ®eld soil. PhD Dissertation ETH, NO. 10185, 136 pp.
Flury, M., FluÈhler, H., Jury, W.A., Leuenberger, J., 1994. Susceptibility of soils to preferential ¯ow of water: a ®eld study. Water
Resour. Res. 30, 1945±1954.
Gerke, H.H., van Genuchten, M.Th., 1993. A dual-porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29, 305±319.
Ghodrati, M., Chendorain, M., Chang, Y., 1999. Characterization ofmacropore ¯ow mechanisms in soil by means of a split macropore column. Soil Sci. Soc. Am. J. 63, 1093±1101.
Heimovaara, T., de Water, E., 1991. A computer controlled TDR
system for measuring water content and bulk electrical conductivity of soils. Laboratory of Physical Geography and Soil
Science, University of Amsterdam, Report no. 41, 27 pp.
Heimovaara, T., Focke, A.C., Bouten, W., Verstraten, J.M., 1995.
Assessing temporal variations in soil water composition withtime domain re¯ectometry. Soil Sci. Soc. Am. J. 59, 689±698.
Hutson, J.L., Wagenet, R.J., 1995. A multiregion model describing
water ¯ow and solute transport in heterogeneous soils. Soil Sci.
Soc. Am. J. 1995, 743±751.
Inoue, M., SÏimuÊnek, J., Shiozawa, S., Hopmans, J.W., 2000. Simultaneous estimations of soil hydraulic and solute transportparameters from transient in®ltration experiments. Adv. Water
Resour. 23, 677±688.
Jacobsen, O.H., Leij, F.J., van Genuchten, M.Th., 1992. Parameter
determination for chloride and tritium transport in undisturbed
lysimeters during steady ¯ow. Nordic Hydrol. 23, 89±104.
Jacques, D., 2000. Analysis of water ¯ow and solute transport at the®eld-scale. PhD, no 454, Faculteit Land bouwkundige enToegepaste Biologische Wetenschappen, K.U. Leuven,
Belgium, 255 pp.

Jacques, D., Vanderborght, J., Mallants, D., Kim, D.-J., Vereecken,
H., Feyen, J., 1997. Comparison of three stream tube models
predicting ®eld-scale solute transport. Hydrol. Earth Syst. Sci. 4,
873±893.
Jacques, D., Kim, D.-J., Diels, J., Vanderborght, J., Vereecken, H.,
Feyen, J., 1998. Analysis of steady-state chloride transport
through two heterogeneous ®eld soils. Water Resour. Res. 34,
2539±2550.
Jacques, D., Timmerman, A., Feyen, J., 1999. Experimental study ofwater ¯ow and solute transport in a macroporous soil undernatural boundary conditions: Bekkevoort-site description,
experimental layout and calibration procedures. Institute forLand and Water Management, K.U. Leuven, Internal Report
No. 56, 99 pp.
Jarvis, N., Jansson, P.-E., Dik, P., Messing, I., 1991. Modelling
water and solute transport in macroporous soil. I. Modeldescription and sensitivity analysis. J. Soil Sci. 42, 305±319.
Kasteel, R., 1997. Solute transport in an unsaturated ®eld soil:
describing heterogeneous ¯ow ®elds using spatial distributionof hydraulic properties. Dissertation ETH. No. 12,477, 108p.
Lehmann, F., Ackerer, P., 1997. Determining soil hydraulicproperties by inverse method in one-dimensional unsaturated¯ow. J. Environ. Qual. 26, 76±81.
Mallants, D., 1996. Water ¯ow and solute transport in a heterogeneous soil pro®le. PhD, no 309, Faculteit Landbouwkundige
en Toegepaste Biologische Wetenschappen, K.U. Leuven,
Belgium, 312 pp.
Mallants, D., Vanclooster, M., Feyen, J., 1996. Transect study on
solute transport in a macroporous soil. Hydrol. Proc. 10, 55±70.
Mallants, D., Jacques, D., Tseng, P.-H., van Genuchten, M.Th.,
Feyen, J., 1997a. Comparison of three hydraulic property
measurement methods. J. Hydrol. 199, 295±318.
Mallants, D., Mohanty, B.P., Vervoort, A., Feyen, J., 1997b. Spatial
analysis of saturated hydraulic conductivity in a soil withmacropores. Soil Technol. 10, 115±131.
Marquardt, D.W., 1963. An algorithm for least-squares estimationof nonlinear parameters. SIAM J. Appl. Math. 11, 431±441.
Mayer, A.S., Huang, C., 1999. Development and application of acoupled-process parameter inversion model based on the maximum likelihood estimation method. Adv. Water Resour. 22,841±853.
Medina, A., Carrera, J., 1996. Coupled estimation of ¯ow and solutetransport parameters. Water Resour. Res. 32, 3063±3076.
Mishra, S., Parker, J.C., 1989. Parameter estimation for coupledunsaturated ¯ow and transport. Water Resour. Res. 25, 385±
396.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten,
M.Th., 1997. New piecewise-continuous hydraulic functions formodeling preferential ¯ow in an intermittent-¯ood-irrigated®eld. Water Resour. Res. 33, 2049±2063.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., SÏimuÊnek, J., vanGenuchten, M.Th., 1998. Preferential transport of nitrate to atile drain in an intermittent-¯ood-irrigated ®eld: model development and experimental evaluation. Water Resour. Res. 34,
1061±1076.
Mualem, Y., 1976. A new model for predicting the hydraulic

D. Jacques et al. / Journal of Hydrology 259 (2002) 15±31
conductivity of unsaturated porous media. Water Resour. Res.
12, 513±522.
Radcliffe, D.E., Tillotson, P.M., Hendrix, P.F., West, L.T., Box,
J.E., Tollner, E.W., 1996. Anion transport in a piedmont ultisol:
I. Field-scale parameters. Soil Sci. Soc. Am. J. 60, 755±761.
Reedy, O.C., Jardine, P.M., Wilson, G.V., Selim, H.M., 1996.
Quantifying the diffusive mass transfer of nonreactive solutes
in columns of fractured saprolite using ¯ow interruption. Soil
Sci. Soc. Am. J. 60, 1376±1384.
Roth, K., Jury, W.A., FluÈhler, H., Attinger, W., 1991. Transport ofchloride through an unsaturated ®eld soil. Water Resour. Res.
27, 2533±2541.
Schoen, R., Gaudet, J.P., Elrick, D.E., 1999. Modelling of solutetransport in a large undisturbed lysimeter during steady-statewater ¯ux. J. Hydrol. 215, 82±93.
Schulin, R., van Genuchten, M.Th., FluÈhler, H., Ferlin, P., 1987. An
experimental study of solute transport in a stony ®eld soil. Water
Resour. Res. 23, 1785±1794.
SÏimuÊnek, J., SÏejna, M., van Genuchten, M.Th., (1998). The
HYDRUS-1D software package for simulating the onedimensional movement of water, heat, and multiple solutes invariably-saturated media. Version 2.0, IGWMC-TPS-70, International Ground Water Modeling Center, Colorado School ofMines, Golden, Colorado, 202 pp.
SÏimuÊnek, J., Wendroth, O., van Genuchten, M.Th., 1999. Estimating unsaturated soil hydraulic properties from laboratory tension
disc in®ltrometer experiments. Water Resour. Res. 35, 2965±
2979.
SÏimuÊnek, J., Jacques, D., Hopmans, J.W., Inoue, M., Flury, M., vanGenuchten, M.Th., 2001. Solute transport during variablysaturated ¯owÐinverse methods. In: Dane, J.H., Topp, G.C.
(Eds.). Methods of Soil Analysis, Part 1, Physical Methods.
3rd ed. SSSA, Madison, WI in press.
Skopp, J., Gardner, W.R., 1992. Miscible displacement: an interacting ¯ow region model. Soil Sci. Soc. Am. J. 56, 1680±1686.
Snow, V.O., Clothier, B.E., Scotter, D.R., White, R.E., 1994. Solutetransport in a layered soil: experiments and modelling using theconvection±disperion approach. J. Cont. Hydrol. 16, 339±358.
Steenhuis, T.S., Parlange, J.-Y., Andreini, M.S., 1990. A numericalmodel for preferential solute movement in structured soils.
Geoderma 46, 193±208.
Sun, N.-Z., Yeh, W.W.-G., 1990. Coupled inverse problems ingroundwater modeling. 1. Senstivity analysis and parameter
identi®cation. Water Resour. Res. 26, 2507±2525.

31

Vanderborght, J., 1997. Experimental and numerical study on nonreactive solute transport in soils. PhD no. 349, Faculteit
Landbouwkundige en Toegepaste Biologische Wetenschappen,
K.U. Leuven, Belgium, 310 pp.
Vanderborght, J., Vanclooster, M., Mallants, D., Diels, J., Feyen, J.,
1996. Determining convective lognormal solute transport parameters from resident concentrations. Soil Sci. Soc. Am. J. 60,
1306±1317.
Vanderborght, J., Mallants, D., Vanclooster, M., Feyen, J., 1997.
Parameter uncertainty in the mobile±immobile solute transportmodel. J. of Hydrol. 190, 75±101.
Vanderborght, J., Timmerman, A., Feyen, J., 2000. Solute transport
for steady-state and transient ¯ow in soils with and without
macropores. Soil Sci. Soc. Am. J. 64, 1305±1317.
Vanderborght, J., Vanclooster, M., Timmerman, A., Mallants, D.,
Jacques, D., Hubrechts, L., Gonzalez, C., Feyen, J., Diels, J.,
Deckers, J., 2001. Overview of inert tracer experiments in key
Belgian soil types: relation between transport, morphologicaland hydraulic soil properties. Water Resour. Res. in press.
van Genuchten, M.Th., 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc.
Am. J. 44, 892±898.
van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer studies
in sorbing porous media I. Analytical solutions. Soil Sci. Soc.
Am. J. 40, 473±480.
van Genuchten, M.Th., Schaap, M.G., Mohanty, B.P., SÏimuÊnek, J.,
Leij, F.J., 1999. Modeling ¯ow and transport processes at the
local scale. In: Feyen, J., Wiyo, K. (Eds.). Modelling ofTransport Processes in Soils at Various Scales in Time andSpace, pp. 23±45.
Ventrella, D., Mohanty, B.P., SÏimuÊnek, J., Losavio, N., vanGenuchten, M.Th., 2000. Water and chloride transport in a®ne-textured soil: ®eld experiments and modeling. Soil Sci.
165, 624±631.
Vogeler, I., Clothier, B.E., Green, S.R., Scotter, D.R., Tillman,
R.W., 1996. Characterizing water and solute movement bytime domain re¯ectometry and disk permeametry. Soil Sci.
Soc. Am. J. 60, 5±12.
Wallach, R., Steenhuis, T.S., 1998. Model for nonreactive transport
in structured soils with continuous preferential ¯ow paths. Soil
Sci. Soc. Am. J. 62, 881±886.
ZurmuÈhl, T., Durner, W., 1996. Modeling transient water and solutetransport in a biporous soil. Water Resour. Res. 32, 819±829.