Inverse Mobile–Immobile Modeling of Transport During Transient Flow:
Effects of Between-Domain Transfer and Initial Water Content
J. Maximilian Köhne,* Sigrid Köhne, Binayak P. Mohanty, and Jirka Šimůnek

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTJ.M. Köhne and B.P. Mohanty, Department of Biological & Agricultural Engineering, Texas A&M University, Scoates Hall, College Station, TX 77843-2117; S. Köhne, 1912 Vinewood, Bryan, TX 77802; J.
Šimůnek, University of California, Riverside, Dep. of Environmental
Sciences, 900 University Avenue, A135 Bourns Hall, Riverside, CA
92521. Received 16 Sept. 2003. Original Research Paper. *Corresponding author (mkoehne@cora.tamu.edu).

solute concentrations, as was demonstrated using laboratory experiments with a macroporous soil column(Köhne and Mohanty, unpublished data, 2004).
The initial or antecedent water content (i.e., the moisture status of a soil before solute and rain application)is likely to affect the infiltration of rain and of surfaceapplied chemicals. However, investigations on the effectof the initial water content on solute transport haveyielded ambiguous results. For example, White et al.
(1986) and Shipitalo et al. (1990) found that lower initialwater contents favor physical nonequilibrium during
fast solute leaching, while other studies indicated thathigh initial water contents enhance leaching (Tillman
et al., 1992; Clothier and Green, 1994; Granovsky etal., 1993). In this study we analyzed the effect of theantecedent water content on observations and simulation results of transient water and Br⫺ movement in an
aggregated Ap soil.
Two- and multiregion transport models can be used
for the risk assessment of physical nonequilibrium transport of contaminants into groundwater. The classicalMIM comprises a mobile region subject to steady-stateflow and convective–dispersive solute transport, alongwith diffusive solute exchange with a stagnant (immobile) region (e.g., van Genuchten and Wierenga, 1976).
The original MIM is not easily applicable to field transport situations involving transient flow. More complex
dual-permeability models (e.g., Jarvis and Leeds-Harrison, 1987; Jarvis et al., 1991; Chen and Wagenet, 1992;
Gerke and van Genuchten, 1993; Kohler et al., 2001)
and multiple permeability models (e.g., Gwo et al., 1995;
Hutson and Wagenet, 1995) may be used to describetransient, variably saturated flow and transport in two
or more mobile regions, coupled with water and solute
mass transfer terms. The application of dual- or multipermeability models is hampered by their large numbersof model parameters. Only few studies have been presented recently where model parameters were obtainedby independent measurement (e.g., Köhne et al., 2002a,
2002b; Logsdon, 2002; Castiglione et al., 2003) or by
inverse estimation (e.g., Schwartz et al., 2000; Kätterer
et al., 2001; Dubus and Brown, 2002; Roulier and Jarvis, 2003).
The MIM approach in several studies has been extended to describe transport under variably saturated
flow conditions (e.g., Russo et al., 1989; Zurmühl andDurner, 1998; Clothier et al., 1998; Larsson et al., 1999).
However, advective transfer between regions was considered only in the study of Larsson et al. (1999). The
most comprehensive MIM accounting for water andsolute transfer between regions at variable saturations

Published in Vadose Zone Journal 3:1309–1321 (2004).
© Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: BTC, breakthrough curve; MIM, mobile–immobile
models.

Mobile–immobile models (MIM) have rarely been used for inverse
simulation of measurements of variably saturated water flow andcontaminant transport. We evaluated two MIM approaches with watertransfer across the mobile and immobile regions either based onrelative saturation (Se) differences, MIM(Se), or pressure head (h )differences, MIM(h ), for inverse simulation of transient water flow
and Br⫺ transport in aggregated soil. Six undisturbed Ap soil columns(14.7-cm length and diameter) at wet, medium, and dry initial watercontents were subjected to application of 0.005 L Br⫺ solution of1000 mg L⫺1 and subsequent irrigation of 1 cm h⫺1 for 3 h. Two similar
irrigations were applied after 7 and 14 d. Measurements comprisedpressure heads at depths of 2.8 and 12.8 cm, average soil column
water contents, outflow, and effluent solute concentrations. This experimental information was used for simultaneous optimization of
van Genuchten soil hydraulic parameters for mobile and immobile
regions, the dispersivity, and the first-order rate coefficients for waterand solute transfer. In total, eight parameters were estimated for
MIM(Se) and 10 parameters for MIM(h ). The inverse MIM approaches described the various hydraulic and transport data adequately. Physical nonequilibrium was more pronounced for wet anddry than for intermediate initial moisture, while initial moisture had no
obvious effect on the total Br⫺ lost. For wet and dry initial conditions,parameter estimates seemed fairly reliable, with the exception of thehighly correlated saturated water contents in mobile and immobile
regions. MIM(h ) yielded parameters that appeared physically moreconsistent with observations, but required smaller time steps thanMIM(Se) to overcome oscillations of pressure heads. Both MIM approaches were found to be suitable for inverse simulation of physicalnonequilibrium transport during variably saturated flow.

During the past several decades, increasing concentrations of surface-applied agrochemicals inground and surface waters have given rise to environmental concerns and have intensified research efforts
to understand underlying transport processes (Flury,1996). Various studies showed that in fine-textured soils,
chemicals may be rapidly channeled to groundwater
through preferential pathways caused by aggregation,
cracking, decaying plant roots, and earthworm burrows(e.g., Flury et al., 1994; Mohanty et al., 1997; Lennartz
et al., 1999). One manifestation of preferential flow insoil is physical nonequilibrium, that is, distinct local(centimeter to decimeter scale) differences with respect
to flow velocities, pressure heads, water contents, and1309

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1310

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

was recently presented by Šimůnek et al. (2003). Herewe will extend their approach by utilizing inverse parameter estimation to analyze experimental data for physical nonequilibrium solute transport under transient
flow conditions.
An important aspect of this study is the effect of theadvective component in the MIM solute transfer termon the inverse transport simulation results. Solute transfer has been identified as crucial for accurate predictionof nonequilibrium transport (e.g., Flühler et al., 1996;
Wilson et al., 1998; Gerke and Köhne, 2004). Advective
transfer directly depends on water transfer that in its
turn can be either assumed to be proportional to thedifferences in saturation (e.g., Jarvis et al., 1991; Kohler
et al., 2001) or the pressure head (e.g., Gerke and van
Genuchten, 1993; Gwo et al., 1995) across regions. However, the effect of having two different advective terms
on MIM simulation results for solute transport has notbeen evaluated to date.
The objectives of this work were (i) to study the effectof the initial soil moisture on physical nonequilibrium
water flow and solute transport as observed and simulated for six aggregated Ap soil columns and (ii) toevaluate two MIM approaches in terms of their suitability for inverse simulation of measurements of variably
saturated flow and physical nonequilibrium Br⫺ transport. The approaches assume water and advective transfer either based on differences between the immobile
and mobile regions in pressure heads, approach MIM(h),or saturations, approach MIM(Se).
MATERIALS AND METHODSSix undisturbed soil columns (14.7-cm diameter and 14.7cm length) were taken from the Ap horizon of a loamy soil
in a field near Kiel (Northern Germany). The Ap soil consistedof 2.3, 43.8, and 53.9% (w/w) clay, silt, and sand, respectively,
and showed subangular and crumbed peds of 0.5- to 1-cm
diameter. Total C and CO3⫺–C of the soil materials were determined by dry-ashing at 1200⬚C and measuring the volume ofCO2 produced after reaction with HCl, respectively. Organic
C (1.83%) was calculated as the difference between total andCO3⫺–C. A soil pH of 6.7 was measured in 0.01 M CaCl2
solution at a soil/solution ratio of 1:2.5. The average bulk
density, ␳ (1.37 g cm⫺3), and porosity, P (0.47), were determined from 20 100-cm3 samples. Assuming a mean particle
density (␳p) of 2.6 g cm⫺3, porosity was calculated as P ⫽1
1 ⫺ ␳␳⫺p .
Tracer displacement experiments were conducted to study
the effect of different initial water contents on transport behavior during transient flow. We used two replicates for eachof three different initial water saturations, henceforth called“wet” (Columns Ap1, Ap2), “medium” (Ap3, Ap4), and “dry”(Ap5, Ap6). Conditions of selected experiments are summarized in Table 1. Rainfall at the top of the column was simulated by means of an automated sprinkler, while constantsuctions in a range from 15 to 60 hPa were applied to thecolumns through the lower boundaries. Pressure heads at two
depths (2.8 and 12.8 cm) were recorded using tensiometers.
Pressure heads were assumed to represent the mobile region
since the tensiometers’ ceramic tip of 5-cm length was intersecting several interaggregate sections. Columns Ap1, Ap3,
Ap4, and Ap6 were placed on automatically weight registering
scales to monitor depth-averaged water contents. At the end

Table 1. Experimental conditions of the transport experiments.
Rainfall depth¶
in event
Column

␪ i†

␪max‡
3

Ap1
Ap2
Ap3
Ap4
Ap5
Ap6

cm
0.316
0.288
0.254
0.250
0.222
0.208

cm⫺3
0.359
0.334
0.330
0.323
0.312
0.315

Suction§
hPa
⫺30
⫺22
⫺60
⫺60
⫺60
⫺60

1

2

2.93
2.96
2.95
2.94
2.87
2.96

cm
2.97
2.95
3.04
2.89
2.93
2.97

3

Applied Br⫺

2.92
2.99
3.01
2.87
2.91
2.91

mg L⫺1
1010
1001
1001
1010
1001
1001

† ␪i, initial water content.
‡ ␪max, maximum water content reached during a rainfall event.
§ Suction applied to the columns’ bottom of the columns.
¶ Rainfall intensity was 1 cm h⫺1 (i.e., rainfall depth in cm was equal torainfall duration in hours).

of the transport experiment, the column was weighed and thesoil was dried at 105⬚C to determine the gravimetric watercontent. The final volumetric water content was calculated bymultiplying the gravimetric value with the bulk density, which
was obtained by dividing the dry weight of the soil by thecolumn volume. This final volumetric water content and thebalance registrations were used to back-calculate depth-average volumetric water contents during the experiments. Since
only four scales were available for automated measurements,
for the Ap2 and Ap5 columns only the initial and the maximum
water contents could be determined from the initial weights
and the weights measured after the first irrigation.
The transport experiment was initiated by applying 5 cm3(0.03 cm) of a KBr solution with a concentration of 1000
mg L⫺1 Br on the soil surface within 10 min by means of achromatography syringe. Starting 20 min later, 3 cm of rain
was sprinkled onto the columns at a rate of 1 cm h⫺1. Twomore similar irrigation events were applied in 7-d intervals.
The effluent of each column was collected in 16- to 18-cm3
volume fractions and analyzed for Br⫺ (flux type) concentrations. For a complete mass balance, soil aliquots were takenfrom 3-cm-thick layers of the Ap columns to determine watercontent and the mass of Br⫺ residing in these layers. The latter
was determined by extracting 50 g dry soil with 50 cm3 H2Oby overhead shaking for 1 h. Subsequently, the batch samples
were rotated in a centrifuge and the supernatant solution
was analyzed for Br⫺ (resident type) concentration. Bromide
concentrations were determined using an ion-chromatographysystem with a conductivity detector (GAT Wescan, Kontron
Instruments, Exchingen, Germany). The solid phase was a 100mm long anion cartridge (Methrom, Super-Sep 6.1009.000,
Deutsche Methrom GmbH & Co. KG, Filderstadt, Germany)
together with a 20-mm guard-cartridge (Bischoff chromatography, part 6302137, Metrohm AG, Herisau, Switzerland). The
mobile phase was a 2.5-mM phthalic acid eluant (pH 4.2) thatwas degassed via an automatic degasser (Degasys, Chrom
Tech, Apple Valley, MN). The detection limit of the system
was 0.2 mg L⫺1 Br. Additional information about the experiments can be found in Meyer-Windel (1998). Bromide andherbicide transport in the wet Ap1 and Ap 2 soil columns
have been discussed in detail by Meyer-Windel et al. (1999).

MODELS
The MIM approach assumes partitioning of the soil watercontent, ␪, into a mobile region water content, ␪m, and an
immobile region water content, ␪im (van Genuchten and Wierenga, 1976):␪ ⫽ ␪m ⫹ ␪im[1]

The Richards’ equation is used to describe one-dimensional

1311

www.vadosezonejournal.org

transient water flow in the mobile region, while a source–sinkterm accounts for advective water exchange with the variably
saturated immobile region (slightly modified after Šimůnek
et al., 2003):

冤

冢

冣冥

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

⳵
⳵␪m
⳵h
⫽
K(hm) m ⫺ 1 ⫺ ⌫w
⳵t
⳵z
⳵z
⳵␪im
⫽ ⌫w
⳵t

[2]

where t (T) is time, z (L) is vertical distance positive downward,
K (L T⫺1) is the hydraulic conductivity as a function of pressure
head in the mobile region, hm (L), and ⌫w (T⫺1) is the watertransfer rate between the mobile and the immobile region.
We used two alternative approaches for calculating thewater transfer term ⌫w. First, ⌫w was assumed to be proportional to the difference in effective saturation, Se (0 ⱕ Se ⱕ1), between the mobile and immobile regions, MIM(Se)⌫w ⫽ ␻(Sem ⫺ Seim)[3]

where ␻ (T⫺1

) is a first-order water transfer rate coefficientand where Sem (Seim) is the effective saturation of the mobile(immobile) region defined as␪m ⫺ ␪m,r
␪m,s ⫺ ␪m,r
␪ ⫺ ␪im,rSeim ⫽ im␪im,s ⫺ ␪im,r
Sem ⫽

[4b]

[5]

where ␪m and nm are parameters defining the shape of theretention function for the mobile region.
Alternatively, ⌫w was assumed to be proportional to thedifference in pressure heads between the mobile and immobile
regions, MIM(h ).

⌫w ⫽ ␻*(hm ⫺ him)
where ␻* (L⫺1

⫺1

[6]

T ) is a first-order water transfer rate coefficient. The MIM(h ) approach requires a retention function of
the immobile region analogous to Eq. [5]. Water transfer terms
based on pressure head differences (Eq. [6]) are commonly
perceived as being more physically realistic than water transfer
terms based on saturation differences (Eq. [3]), since the difference in pressure heads is the actual driving force (Šimůnek
et al., 2003). However, note that the first-order term (Eq. [6])approximates the real, unknown pressure head profile as adifference between only two values (i.e., hm and him). Due tothis first-order simplification, no substantial conceptual differences exist between first-order terms based on saturation andpressured head. In fact, Eq. [6] and [3] would become identical
if in Eq. [3] the effective saturations in both mobile and immobile regions were calculated using van Genuchten (1980) relations (Eq. [5]) with the same parameter values for ␣ and n.
Hence, the main difference between these two approaches isthat while MIM(Se) assumes similarity in the retention curvesof the mobile and immobile regions, the MIM(h ) approachallows different retention curves for each region.
To assess the potential effect of initial moisture on themobile water content, the fraction of mobile water at saturation, ␤s, was calculated as␪m,s␪m,s ⫹ ␪im,s[7]

The MIM formulation for nonreactive solute transport comprises a convection–dispersion equation for the mobile region
coupled via a sink–source term to the immobile region (Šimůnek et al., 2003):

冢

冣

⳵
⳵␪mcm
⳵q c
⳵c
⫽
␪mDm m ⫺ m m ⫺ ⌫s
⳵t
⳵z
⳵z
⳵z

[8a]

⳵␪imcim
⫽ ⌫s
⳵t

[8b]

⌫s ⫽ ␣s(cm ⫺ cim) ⫹冦⌫⌫ ccw mw im⌫w ⬎ 0⌫w ⬍ 0[8c]where cm and cim (M L⫺3) are solute concentrations in themobile and immobile domains, respectively, qm (L T⫺1) is theflux density in the mobile region, ⌫s (M L⫺3 T⫺1) is the solutetransfer rate between mobile and immobile region, with ⌫w
calculated either using Eq. [3] for MIM(Se) or Eq. [6] forMIM(h ), ␣s (T⫺1) is the first-order solute transfer rate coefficient, and Dm (L2 T⫺1) is the dispersion coefficient in themobile region defined as[4a]where ␪m (␪im) (L3 L⫺3) denotes the mobile (immobile) watercontent with residual and saturated constraints of ␪m,r (␪im,r)and ␪m,s (␪im,s), respectively. To describe ␪m in Eq. [4a], thefollowing form of van Genuchten’s (1980) equation was
assumed␪m ⫽ ⎣1 ⫹ (␣m|h|)nm⎦⫺(1⫹1/nm)

␤s ⫽

Dm ⫽

␭ m qm
⫹ D0 ␶m
␪m[9a]

in which
7⁄

␶m ⫽

␪m3
␪2m,s

[9b]

where ␭m (L) is the dispersivity of the mobile region, D0 (L2T⫺1) is the molecular diffusion coefficient, ␶m is the tortuosity
coefficient of the mobile domain evaluated according to Millington and Quirk (1961), and ␪m,s is the saturated water content of the mobile domain.
The above MIM equations cannot be solved analytically.
They were implemented in the numerical HYDRUS-1D computer model (Šimůnek et al., 2003). To ensure an accurate
numerical solution, the numerical grid consisted of 200 elements of only 0.07-cm length, while the maximum time-step
was restricted to 0.001 d in a time-step adaptive Picard iteration solution scheme. The effect of using different spatial andtemporal resolutions on numerical simulation results was analyzed (see Results section).
The initial condition for water flow was represented by apressure head profile linearly interpolated between measurements of tensiometers installed at the 2.8- and 12.8-cm depth.
We assumed initial equilibrium with identical pressure heads
in the mobile and immobile regions. The upper flow boundary
condition was a variable flux-type condition equal to the timevariable rates of irrigation (22 cm d⫺1) and evaporation (0.3cm d⫺1) set according to the experimental conditions. Upon
ponding this boundary condition may evolve into a programcontrolled head-type condition. For the lower boundary, an
unsaturated seepage face at a user-specified pressure head(constant head imposed at the lower boundary of soil column)
was implemented in HYDRUS-1D. The initial condition forsolute transport was set to zero. The upper and lower boundaries for solute transport were represented by specified fluxand (zero) gradient conditions, respectively.

Parameter Estimation
For describing mobile–immobile water flow using Eq. [2],[3], [4], and [5], MIM(Se) requires eight parameters (Ks, ␪m,r,Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1312

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

␪m,s, ␣m, nm, ␪im,r, ␪im,s, and ␻). The MIM(h ) approach (Eq. [2],[4], [6]) additionally requires ␣im and nim. To reduce the numberof unknowns, ␪m,r and ␪im,r were fixed to values of zero consistent with ␪r values of zero fitted to water retention and hydraulic conductivity data for all Ap soil columns (Meyer-Windel,
1998). This limited the total number of unknown water flow
parameters to six [MIM(Se)] and eight [MIM(h )], respectively. To calculate solute transport using Eq. [8] and [9], bothinverse MIM approaches additionally require estimation of ␣sand ␭, whereas D0 (1.797 cm2 d⫺1) was assumed to be a known
value for Br (Atkins, 1990). Hence, full descriptions of mobile–immobile water flow and solute transport assuming zero residual water contents require 8 and 10 parameters for MIM(Se)and MIM(h ), respectively.
The inverse parameter estimation was performed by minimization of the objective function ⌽ (Šimůnek et al., 1998):mn⌽(b) ⫽ 兺 vj 兺 wi, j [Oj (z,ti ) ⫺ Ej (z,ti,b)]2
j⫽1

[10]

i⫽1

where m is the number of different sets of measurements, nrepresents the number of observations in a particular measurement set, Oj(z,ti) are observations such as pressure heads inthe mobile region at the 2.8- and 12.8-cm depths, average soilcolumn water contents (averaged over the column depth andover mobile and immobile regions), cumulative water fluxes
across the lower boundary, and effluent concentrations at time
ti at location z, Ej(z,ti,b) are the corresponding estimated spacetime variables for the vector b of optimized parameters (␪im,r,␪im,s, ␪m,r, ␪m,s, ␣m, nm, Ks, ␻, ␭, ␣s, and additionally ␣im, nim forMIM(h ), and vj and wi,j are weighting factors associated witha particular measurement set or point, respectively. We assumed in this study that the weighting coefficients wi,j in Eq.
[10] are equal to one; that is, the variances of the errors withina particular measurement set are the same. Only data thatwere measured at larger time intervals (and hence were underrepresented as compared with the more frequent measurements) were assigned larger weights wi,j, calculated as theratio of their measurement time interval to the average timeinterval. The weighting factors vj were obtained in two steps.
An initial value for vj was calculated for each inverse simulation as (Clausnitzer and Hopmans, 1995)

vj ⫽1
nj ␴ 2j[11]

which assumes that vj is inversely related to the variance ␴ 2j
within the jth measurement set and to the number of nj measurements of the set. However, using Eq. [11] was not sufficientto yield equal contributions of the four measurement sets
(i.e., water contents, pressure heads, water outflow, and Br⫺
concentrations in the outflow) to the optimized objective function (Eq. [10]) and resulted in poor simulations of those measurement sets that were underrepresented in the objectivefunction. Therefore, the weight vj of every measurement set
was fine-tuned by repeated manual calibration in consecutive
inverse simulation runs to achieve equal individual contributions of the four measurement sets to the overall sum ofsquared residuals.
The Levenberg–Marquardt algorithm was used to minimize
the objective function (Eq. [10]) and to calculate confidence
intervals and correlation matrices for the parameters (Šimůnek et al., 1998). Our inverse modeling strategy relied on
simultaneous estimation of hydraulic and transport parameters from combined water flow and solute transport information. This approach takes advantage of all information contained in the data, since concentrations are a function of waterflow also (Medina and Carrera, 1996). For coupled single-

domain flow and transport, simultaneous estimation of hydraulic and transport parameters was found to yield smaller
estimation errors and confidence intervals for model parameters than sequential inversion of hydraulic properties fromwater content and pressure head data followed by inversionof transport properties from concentration data (Šimůnek etal., 2002).
To facilitate direct comparisons of the goodness-of-fit fordifferent data sets, the RMSE between observations (Oi) andmodel estimates (Ei) for a particular data set was normalizedwith the arithmetic mean of observations, O, to yield thecoefficient of variation (CV, %) as冢冣100 兺in⫽1(Oi ⫺ Ei)2 ⁄2 100
CV ⫽
⫽
RMSE
O
n
O
1

[12]

An average coefficient of variation, CV, over the four measurement sets of a particular Ap column was calculated as astatistical measure of the overall agreement between modeland observations as follows:CV ⫽兺km⫽1CVkm[13]

where CVk is the coefficient of variation (Eq. [12]) for measurement set k, and m denotes the total number of CV values
or data sets (in this study m ⫽ 4).
To enable a more physically based interpretation of thesimulated MIM transport processes, we analyzed the degree ofphysical (non-)equilibrium for different initial water contents.
For steady-state flow, an index frequently used for defining
physical nonequilibrium is the Damköhler number (Damköhler, 1936; Bahr and Rubin, 1987; Vanderborght et al., 1997;
Michalak and Kitanidis, 2000). However, the Damköhler number cannot be used for transient flow conditions. In this study,
we introduce an empirical physical equilibrium (PE) index toquantify the degree of physical equilibrium during transient
flow by describing the combined effects of advective and diffusive first-order mass transfers. The PE index (%) is defined asPE ⫽ 100

␻ ⫹ ␣s␻max ⫹ ␣ s,max[14]

where ␻ (T⫺1) (Eq. [3], [6]) is a first-order water transfer rate
coefficient, ␣s (T⫺1) (Eq. [8]) is a first-order solute transfer
rate coefficient (Table 1), and ␻max and ␣s,max are maximum
parameter constraints for ␻ and ␣s, respectively, indicative ofwater and solute transfer equilibrium. Theoretically, instantaneous equilibrium exists only if ␻ and ␣s are infinite. Practically, parameter constraints ␻max and ␣s,max could be operationally defined, for example, by using values that force watercontents and solute concentrations in mobile and immobile
regions to reach 99% equilibrium within the experimental timescale. However, such a definition would require individual
simulations for every set of parameters, and initial and boundary conditions to identify ␻max and ␣s,max. Using a more practical
approach, ␻max and ␣s,max in this study were both simply set
equal to a large value of 10 d⫺1; hence, ␻max⫹ ␣s,max equaled
20 d⫺1. This was found to effectively ensure equilibration ofwater contents or pressure heads and concentrations within
the experimental time scale. For the parameter combinations
tested in this study, increasing values for ␻ or ␣s to above 10d⫺1 did not substantially change the results (not further shownhere). By definition, the PE index is inversely related to thedegree of physical nonequilibrium and varies between 0%(perfect nonequilibrium, no mass transfer) to 100% (interregion equilibrium due to very fast mass transfer).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 1. Water flow data and simulation results obtained with MIM(Se)and MIM(h ) for the initially wet columns Ap1 (left) and Ap2 (right).
(a, b) Water contents, (c, d) pressure heads in the mobile region
at depths of 2.8 cm and (e, f) 12.8 cm, and (g, h) the cumulative
outflow.

RESULTS AND DISCUSSION
Water Flow Results
Initially Wet Ap Soil Columns (Ap1, Ap2)
Experimental water flow data and model simulations
for the initially wet soil columns are shown in Fig. 1.
Water contents, pressure heads in the mobile region at
depths of 2.8 and 12.8 cm, and the cumulative outflow
exhibited a recurrent pattern of increasing values during
infiltration and decreasing or constant (cumulative outflow) values during redistribution (Fig. 1a–1h). For theAp2 column, only the initial water content and a value
after the first irrigation event were known (Fig. 1b).
Observations were comparable for Ap1 and Ap2, reflecting similar water flow dynamics in the initially wet
replicate soil columns. According to Fig. 1, the MIM(Se)and MIM(h) approaches approximated observed watercontents, pressure heads, and cumulative outflow equally
well. Since infiltration occurs on a smaller time scalethan drainage, Fig. 2 shows a detailed view of the shortterm responses to one irrigation event. The third irrigation event is presented here to illustrate the maximum
deviations between model results and data, as opposedto the first event, where known initial conditions resulted in a closer match with the data (not shown). The
detailed view still shows reasonable agreement between
hydraulic data measured during the third irrigation
event and both MIM approaches (Fig. 2).

1313

Fig. 2. Same as Fig. 1, but with detail for the third irrigation event only.

Medium Initial Water Content in Ap Soil Columns(Ap3, Ap4)
For the replicate Ap3 and Ap4 soil columns of medium initial water content, hydraulic dynamics as reflected by water contents, mobile region pressure heads
at the 2.8-cm depth, and cumulative outflow were again
comparable (Fig. 3). All observations were matchedsatisfactorily by MIM(Se) and MIM(h), both on averageduring the experimental duration of 20 d (Fig. 3) andduring single irrigation events (not shown).
Initially Dry Ap Soil Columns (Ap5, Ap6)The Ap5 and Ap6 soil columns were characterized
by low initial water contents (see Table 1). Hydraulicdynamics were comparable for Ap5 and Ap6 (Fig. 4),
except for differences in water content variations during
the first 7 d of the experiment (no water contents could
be measured for Ap5 after 7 d) (Fig. 4a). However, it
is likely that these apparent differences were causedby less reliable measurements of Ap5 water contents.
Except for the water contents for Ap5, hydraulic data
were reasonably well described with both MIM(Se) andMIM(h) during the 20-d experiment (Fig. 4a–4f) andduring the irrigation events (detailed view not shown).
Comparison of Hydraulic Data, Simulations,and Parameters
The maximum water contents that were reached after
the irrigations declined in the order: Ap1,2 (wet); Ap3,4(medium); Ap5,6 (dry) (see Fig. 1a, 1b, 3a, 3b, 4a, andReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1314

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

Fig. 3. Water flow data and simulation results obtained with MIM(Se)and MIM(h ) for the initially medium wet columns Ap3 (left) andAp4 (right). (a, b) Water contents, (c, d) pressure heads in the mobile region at the 2.8-cm depth, and (e, f) the cumulative outflow.

Fig. 4. Water flow data and simulation results obtained with MIM(Se)and MIM(h ) for the initially dry columns Ap5 (left) and Ap6 (right).
(a, b) Water contents, (c, d) pressure heads in the mobile region
at the 2.8-cm depth, and (e, f) the cumulative outflow.

4b), indicating incomplete saturation during irrigationof the initially drier soil columns. The MIM(h) simulation results for the saturated water content, ␪s, summed
for the mobile and immobile regions (␪s,m ⫹ ␪s,im), reflected the decreasing order of maximum water contents
with decreasing initial water content, whereas respective
␪s values obtained with MIM(Se) did not consistently
decrease with drier initial conditions (Table 2). The ␪svalues obtained with MIM(h) were hence more consistent with the observations, even if statistical CV values
indicated that MIM(Se) and MIM(h) on average described hydraulic observations of water contents, pressure heads, and outflow equally well (Table 3).
For the initially wet columns Ap1 and Ap2, independently measured and inverse estimates of retention functions, ␪(h), and hydraulic conductivity functions, K(h),are displayed in Fig. 5. Measurements of ␪(h) and K(h)were obtained at static and steady-state conditions, respectively (Meyer-Windel, 1998). The estimated bimodal functions represent the summation (Durner,
1993) of van Genuchten subcurves based on parameters
obtained by the inverse numerical solution of MIM(h)and MIM(Se) using all available transient flow andtransport data. Good agreement between independentmeasurements and predictions of ␪(h) and K(h) was
achieved particularly for MIM(h), while the inverse
MIM(Se) approach overestimated ␪s for Ap1 (Fig. 5).
Incomplete saturation during infiltration may indicate
the presence of immobile regions. Hence, one may assume that the initially dry columns that reached smaller
maximum saturations also had a smaller fraction of mobile water at saturation, ␤s (Eq. [7]). However, ␤s assumed values between 11 and 41%, independent of theinitial water content (Table 4). Moreover, the actual
fraction of mobile water at partial saturation may havebeen variable during transient flow. Overall, the saturation differences were apparently too small (comparedwith other causes for variation) to exert a notable effecton ␤s.
Among the fitted hydraulic mobile–immobile model
parameters, only the MIM(h) parameters ␣im, nim, and␪s (␪s,m ⫹ ␪s,im) showed a tendency to decrease for drier
initial conditions (Table 2), thus indicating an effect of
initial water content on the water flow process. The
remaining parameters (␣m, nm, Ks, ␪s,m, ␪s,im) neither depended on the invoked model [MIM(h) vs. MIM(Se)]nor on the initial water content. The water transfer
coefficients were found to strongly affect the solutetransport simulation, and therefore will be included in
the discussion below.

Solute Transport Results
For all Ap columns, the breakthrough curves (BTC)
in terms of the fraction of applied Br⫺ concentration
vs. cumulative outflow are shown in Fig. 6. The BTCs
exhibited an asymmetric shape with a long tail typicalof physical nonequilibrium. The low peaks of only 0.004
to 0.008 (4–8 mg L⫺1) of the applied concentration(⬇1000 mg L⫺1, Table 1) indicate that only mild physicalnonequilibrium existed for all columns, as opposed topreferential flow, where much higher concentration
peaks would be expected. The initially wet column Ap1
and the initially dry Ap5 and Ap6 columns exhibitedan early rise, whereas the Ap2 column and the initiallymedium wet Ap3 and Ap4 columns displayed a somewhat delayed initial rise. Note that the Ap2 column
had in fact an initial water content somewhere between
those of the initially wet Ap1 and medium wet Ap3
and Ap4 columns (Table 1). Overall, the MIM(Se) andMIM(h) described the BTC reasonably well. According

1315

www.vadosezonejournal.org

Table 2. Estimated van Genuchten (1980) model parameters for the mobile region (vG-Mobile) and the immobile region (vG-immobile),and saturated water content ␪s(␪m,s ⫹ ␪im,s).
vG-mobile (⫾95% confidence interval)
␪m,sColumn␣m

vG-immobile (⫾95% confidence interval)

nm

Ks

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cm⫺1
Ap1
MIM(Se)MIM(h )Ap2
MIM(Se)MIM(h )Ap3
MIM(Se)MIM(h )Ap4
MIM(Se)MIM(h )Ap5
MIM(Se)MIM(h )Ap6
MIM(Se)MIM(h )

␪s

1.25 ⫾ 0.01
1.25 ⫾ 0.04

22.0 ⫾ 0.2
22.5 ⫾ 0.9

0.310 ⫾ 0.009
0.214 ⫾ 0.014–
0.016 ⫾ 0.002–
2.00 ⫾ 0.31

0.43
0.36

0.064 ⫾ 0.005
0.061 ⫾ 0.002

0.047 ⫾ 0.003
0.056 ⫾ 0.003

1.30 ⫾ 0.06
1.43 ⫾ 0.01

20.3 ⫾ 1.3
20.9 ⫾ 0.6
Medium0.317 ⫾ 0.009
0.300 ⫾ 0.003–
0.016 ⫾ 0.001–
1.61 ⫾ 0.15

0.36
0.36

0.044 ⫾ 0.008
0.042 ⫾ 0.002

0.029 ⫾ 0.003
0.039 ⫾ 0.002

1.45 ⫾ 0.05
1.55 ⫾ 0.08

20.5 ⫾ 1.7
20.5 ⫾ 1.1

0.343 ⫾ 0.015
0.295 ⫾ 0.006–
0.015 ⫾ 0.001–
1.51 ⫾ 0.05

0.39
0.34

0.123 ⫾ 0.010
0.111 ⫾ 0.019

0.048 ⫾ 0.001
0.057 ⫾ 0.000

1.26 ⫾ 0.02
1.35 ⫾ 0.07

21.6 ⫾ 0.1
21.6 ⫾ 2.0
Dry0.248 ⫾ 0.015
0.223 ⫾ 0.018–
0.013 ⫾ 0.001–
1.40 ⫾ 0.09

0.37
0.33

0.130 ⫾ 0.005
0.090 ⫾ 0.010

0.043 ⫾ 0.001
0.021 ⫾ 0.005

1.33 ⫾ 0.02
1.23 ⫾ 0.10

24.9 ⫾ 0.9
22.1 ⫾ 1.8

0.296 ⫾ 0.011
0.234 ⫾ 0.014–
0.011 ⫾ 0.002–
1.53 ⫾ 0.12

0.43
0.32

0.086 ⫾ 0.006
0.055 ⫾ 0.006

0.055 ⫾ 0.004
0.070 ⫾ 0.017

1.20 ⫾ 0.01
1.36 ⫾ 0.08

53.6 ⫾ 2.9
41.0 ⫾ 24.6

0.255 ⫾ 0.009
0.263 ⫾ 0.007–
0.013 ⫾ 0.002–
1.25 ⫾ 0.01

0.34
0.32

CV
␪(t )

hm(z ,t )qcum(t )cout(t )

CV†

%
Wet

Ap5
MIM(Se)MIM(h)Ap6
MIM(Se)MIM(h)

nim

0.068 ⫾ 0.004
0.034 ⫾ 0.002

Table 3. Coefficient of variation between measurement sets andmodel estimates for cumulative outflow, qcum(t ), Br⫺ concentration in outflow, cout(t ), pressure heads in the mobile region at
depths of 2.8 and 12.8 cm, hm(z,t), and water contents averaged
over the column depth, ␪(t ).

Ap3
MIM(Se)MIM(h)Ap4
MIM(Se)MIM(h)

␪imcm⫺10.115 ⫾ 0.005
0.147 ⫾ 0.014

to CV values (Table 3), Br⫺ concentrations of Ap1,
Ap2, Ap3, and Ap4 were better matched with MIM(Se),whereas the MIM(h) approach better described the
BTC of the initially dry Ap5 and Ap6 columns.
The fraction of mobile water (discussed above), thedispersivity for the mobile region of the fitted tworegion models, ␭m, the rate coefficients for transfer ofwater, ␻ (␻*), and solute, ␣s, and the PE index are given
in Table 4. The ␭m values varied between 1 mm and3.75 cm, which is within the range of ␭m values reported

Ap1
MIM(Se)MIM(h )Ap2
MIM(Se)MIM(h)

␪im,s

cm d⫺1
Wet

2.2
2.7

20.5
46.0

2.7
2.8

37.4
33.2

15.7
21.3

–
–

64.8
54.9
Medium

8.3
10.6

91.3
101.2

54.8
55.6

5.2
2.1

40.2
45.6

6.9
3.7

77.5
98.3

34.1
37.5

2.8
3.6

43.8
80.8
Dry

3.5
7.5

36.9
43.1

21.8
33.8

3.7
5.6

33.7
37.8

3.5
10.0

100.8
62.6

35.4
29.0

3.4
1.3

43.0
49.5

5.1
1.8

60.5
45.6

28.0
24.6

† Average CV value (Eq. [13]).

by others for column scale studies (e.g., Schwartz et al.,2000). Estimated values for ␻, ␻*, and ␣s for MIM(h)were smaller for Ap1, Ap5, and Ap6 than for Ap2, Ap3,and Ap4 (Table 4). This indicates restricted advective and diffusive transfer (nonequilibrium) for initially
wet and dry as compared with initially intermediatemoisture conditions (equilibrium). For MIM(h), the ␣svalue actually reached its maximum constraint (␣s ⫽ 10)
for Ap2, Ap3, and Ap4 (Table 4), indicating enhanced
diffusive solute transfer for concentration equilibration
across regions. For Ap2, caution must be taken to interpret ␣s with respect to solute transport processes since
the BTC of Ap2 was only described adequately for the
first 2 cm of outflow, whereas the low tailing concentrations could not be described (Fig. 6). For MIM(Se),optimized values for water and solute transfer rate coefficients were not as consistent as for MIM(h). For
example, MIM(Se) gave a large ␻ value for Ap6 and asmall ␻ for Ap5 for similar initial water contents (Table 4).
The PE values obtained with MIM(h) were small forthe high (Ap1) and low (Ap5 and Ap6) initial moisture
conditions, while PE values were large for medium

Fig. 5. (a) Measured water retention and (b) hydraulic conductivity
data and bimodal Durner (1993) functions estimated inversely using
MIM(Se) and MIM(h ) for the Ap1 and Ap2 columns.

1316

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

Table 4. Estimated fraction of mobile water at saturation, ␤s, dispersivity of the mobile region, transfer rate coefficients, and physicalequilibrium index (PE).
Transfer rate coefficients
Column␤s

Dispersivity ␭m

␻†

␣s

PE

cm

d ⫺1

d⫺1

%

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Wet
Ap1
MIM(Se)MIM(h)Ap2
MIM(Se)MIM(h )

0.27
0.41

1.08 ⫾ 0.19
2.26 ⫾ 0.89

0.40 ⫾ 0.04
0.19 ⫾ 0.04

2.20 ⫾ 0.14
1.34 ⫾ 0.27

12.5
7.6

0.17
0.17

1.41 ⫾ 1.08
3.75 ⫾ 1.15

4.51 ⫾ 0.71
3.56 ⫾ 0.36

4.65 ⫾ 1.56
10 #

45.8
67.8

MediumAp3
IM(Se)MIM(h )Ap4
MIM(Se)MIM(h)

0.11
0.13

0.13 ⫾ 0.04
0.02 ⫾ 0.002

4.91 ⫾ 0.76
5.35 ⫾ 0.20

3.53 ⫾ 0.45
10‡

42.2
76.8

0.33
0.33

1.91 ⫾ 0.53
3.59 ⫾ 1.19

1.32 ⫾ 0.19
2.70 ⫾ 0.27

10‡
10‡

56.6
63.5

Dry
Ap5
MIM(Se)MIM(h)Ap6
MIM(Se)MIM(h )

0.31
0.28

1.70 ⫾ 1.36
2.72 ⫾ 1.96

0.73 ⫾ 0.04
0.22 ⫾ 0.02

3.78 ⫾ 1.70
0.72 ⫾ 0.29

22.6
4.7

0.25
0.17

0.97 ⫾ 0.70
0.30 ⫾ 0.38

5.58 ⫾ 0.33
0.11 ⫾ 0.03

2.09 ⫾ 0.57
1.58 ⫾ 0.27

38.4
8.5

† For the MIM(h) approach it is ␻* (cm⫺1 d⫺1).
‡ Fixed at upper constraint.

(Ap3, Ap4) initial moisture, indicating more physicalnonequilibrium for initially wet and dry than for intermediate conditions. For MIM(Se), PE showed sametendencies, but differences were not as pronounced
among PE values for the various water contents (Table4). The simulations with low PE values (Ap1, Ap5, Ap6)
displayed secondary concentration peaks typical of physical nonequilibrium; these were caused by diffusive backFig. 6. Relative Br⫺ concentration vs. cumulative outflow: data for
all Ap columns and simulation results obtained with MIM(Se) andMIM(h ).

transfer from the immobile into the mobile region during no-flow periods. Secondary peaks were also observed for the experimental BTC of Ap2. Considered
together, the experimental data, the model results, and
the PE values indicate a higher degree of physical nonequilibrium for initially wet and dry conditions thanfor medium initial soil moisture. Note, however, thatdifferent initial moisture contents did not systematicallyaffect the total mass loss of Br for the Ap columns understudy, which amounted to 62, 47, 51, 73, 41, and 45%of the applied Br mass for Column Ap1, 2, 3, 4, 5, and6, respectively.
We now discuss differences between MIM(h) andMIM(Se) on the basis of simulation results for Ap1,
Ap4, and Ap6 after 70 min (i.e., 60 min after solution
application and 50 min after initiating the first rain
event), and we suggest a mechanistic interpretation for
the transport at different initial water contents. Effective
saturations (Se) as simulated with MIM(Se) (Fig. 7a)
indicate that the wetting front (Se ⫽ 1) in the mobile
region at 70 min almost reached the bottom of the Ap1column (14.7 cm), followed by Ap4 and Ap6. Behind
the respective wetting fronts, saturation differences, assumed as the driving force for water transfer in MIM(Se), were smaller for the initially wet Ap1 than for
initially intermediate wet Ap4 and initially dry Ap6(Fig. 7a). Additionally, for MIM(Se) the water transfer
coefficient, ␻, increased in the order: Ap1, Ap4, andAp6 (Table 4). Accordingly, water transfer into the immobile region was smallest for Ap1, larger for Ap4, andlargest for Ap6 (Fig. 7b). For the initially dry Ap6, water
laterally imbibed into the immobile region, thereby
slowing down vertical water flow in the mobile region.
By contrast, mobile water in the wet Ap1 could percolate faster due to less horizontal transfer into the rela-1317

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 9. Bromide transfer rate (⌫s) into the immobile (im) region in
the 0- to 2-cm depth layer at the end of the Br⫺ application (10 min)for the Ap1, Ap4, and Ap6 columns. Simulation results obtainedwith (a) MIM(Se) and (b) MIM(h ).

Fig. 7. Depth profiles at 70 min characterizing water flow for the Ap1,
Ap4, and Ap6 columns as simulated with MIM(Se) (top): (a) water
saturation and (b) water transfer rate (⌫w), and as simulated with
MIM(h) (bottom): (c) pressure heads for Ap6 and (d) for Ap1 andAp4, and (e) water transfer rate. Symbols: mo, mobile region; im,immobile region.

tively wet immobile region. The differences in Se values
between the mobile and immobile regions of the Ap6column (Fig. 7a) were two orders of magnitude smaller
than similar differences in h values (Fig. 7c). However,water transfer calculated using both MIM approaches
was comparable for the Ap6 column (Fig. 7b, 7e), since␻ fitted with MIM(Se) was 50 times that of ␻* fittedwith MIM(h) (Table 4). For Ap1 and Ap4, estimatedvalues for ␻ and ␻* were of the same order of magnitude(Table 4), but water transfer in the lower profile of
Ap1 and Ap4 was larger for MIM(h) (Fig. 7e) than for
MIM(Se) (Fig. 7b).
Figure 8 shows simulated concentration profiles andsolute transfer rates at 70 min. For MIM(Se), differences
between Br⫺ in the mobile and immobile regions were
somewhat more pronounced for Ap1 than for Ap4 (Fig.
8a, 8b), corresponding to a larger overall (advective–
diffusive) Br⫺ transfer rate ⌫s for Ap4 (Fig. 8c). For
MIM(h), the Ap4 concentration profiles (Fig. 8d, 8e)
and Br⫺ transfer rate (Fig. 8e) indicate an advanced
state of equilibration throughout the profile. The above
results support the interpretation of having physicalnonequilibrium for high and low initial moisture andtransport close to equilibrium at medium initial moisture contents.
The following explanation can be given for solutetransfer at 70 min being directed into the immobile
region, except in the upper 2 cm (Fig. 8c, 8f). The advective transfer component was always directed into theimmobile region. At the end of the tracer application(at 10 min), diffusive mass transfer in the uppermost
2 cm of the soil profile was also directed into the immobile region, thus leading to a large overall solute transferrate, ⌫s (Fig. 9a, 9b). The first irrigation subsequently
flushed the Br⫺ downward through the mobile region,
thus leaving higher concentrations in the immobile region, which reversed the concentration gradient in theupper 2 cm and directed diffusive transfer back into themobile region after 70 min (Fig. 8c, 8f).

Parameter Uniqueness and Overall
Model Evaluation

Fig. 8. Depth profiles at 70 min characterizing Br⫺ transport for theAp1, Ap4, and Ap6 columns as simulated with MIM(Se) (top): (a)Br⫺ concentration in mobile (mo) and (b) immobile (im) region,(c) Br⫺ transfer rate (⌫s), and as simulated with MIM(h ) (bottom),(d) Br⫺ concentration in mobile and (e) immobile region, (f) Br⫺
transfer rate (⌫s).

Having a large number of flow and transport parameters—8 MIM(Se) and 10 MIM(h) parameters—that
needed to be identified at the same time evokes thequestion of whether or not the estimated parametervalues were unique. Unfortunately, the typical analysis
of parameter uniqueness using two-dimensional response surfaces of the objective function as a functionof parameter pairs (e.g., Šimůnek and van Genuchten,
1996) is not readily possible for such a large numberof parameters. The number of parameter pairs to beanalyzed using response surfaces in our study would be
28 [MIM(Se)] and 45 [MIM(h)], respectively. Moreover, these response surfaces as two-dimensional crosssections through the 8- or 10-dimensional parameter
space may not accurately reflect the true behavior of
the objective function. To assess the reliability of thepresented parametric values, we (i) evaluated confidence intervals for the estimated parameter values, (ii)1318

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Table 5. Correlation, r, between MIM(Se) parameters for inverse simulation of the Ap1, Ap4, and Ap6 transport studies.†
Ap1
␪m,s
␣m
nm
Ks
␪im,s
␻
␭m
␣s
Ap4
␪m,s
␣m
nm
Ks
␪im,s
␻
␭m
␣s
Ap6
␪m,s
␣m
nm
Ks
␪im,s
␻
␭m
␣s

␪m,s

␣m

nm

Ks

␪im,s

␻

␭m

␣s

1
⫺0.01
⫺0.51
⫺0.38
(⫺0.89)
⫺0.16
0.60
⫺0.11

1
⫺0.46
0.08
0.14
⫺0.36
0.01
⫺0.02

1
0.19
0.65
0.07
⫺0.36
0.09

1
0.31
0.28
⫺0.33
0.13

1
⫺0.06
⫺0.56
0.10

1
⫺0.1
0.00

1
0.09

1

1
⫺0.11
⫺0.44
0.41
(⫺0.78)(⫺0.76)
⫺0.12
⫺0.22

1
⫺0.40
0.16
⫺0.07
0.38
0.04
0.00

1
⫺0.56
0.64
⫺0.06
0.04
0.10

1
⫺0.36
⫺0.34
⫺0.21
⫺0.09

1
0.40
0.31
0.39

1
0.19
0.17

1
(0.80)

1

1(⫺0.76)⫺0.71⫺0.46(⫺0.95)
⫺0.16
0.12
0.01

1
⫺0.39
⫺0.19
0.77
(⫺0.38)
⫺0.11
⫺0.03

1
0.56
0.70
0.42
0.00
⫺0.06

1
0.40
0.74
⫺0.03
⫺0.06

1
0.10
⫺0.10
⫺0.02

1
0.07
0.01

1
0.56

1

† Italics denote |r| ⬎ 0.5; parentheses denote |r| ⬎ 0.75.

analyzed correlation matrices, and (iii) studied how different initial parameter values affected final parameter results.
Confidence intervals of the parameters were in general relatively small. Typically, an approach with small
CV values like MIM(Se) for Ap1 (Table 3) yielded
smaller confidence intervals (Table 2). For both MIM(Se)and MIM(h), smaller confidence intervals for the watertransfer rate coefficients were obtained for initially wet(Ap1) and dry (Ap5,6) conditions than for initiallymedium (Ap3,4) and medium wet (Ap2) conditions (Table 4). This indicates that for initially wet and dry conditions, where physical nonequilibrium was more pronounced, parameter estimation using the MIM approaches
was more reliable than for medium initial moisture conditions.
Correlation matrices for Ap1, Ap4, and Ap6 are presented in Table 5 for the MIM(Se) approach and in
Table 6 for the MIM(h) approach. In Tables 5 and 6,
moderately large correlation coefficients (italics denote|r| ⬎ 0.5) and large coefficients are identified (parentheses denote |r| ⬎ 0.75). Correlations between parameters
obtained with MIM(Se) in general appear relatively
small: for Ap1, 96% are below 0.75 and 82% are smaller
than 0.5 (Ap4: 89% ⬍ 0.75, 82% ⬍ 0.5; Ap6: 93% ⬍
0.75, 71% ⬍ 0.5) (Table 5). Particularly for Ap1 andAp6, small correlations for water and solute transfer
coefficients indicate reliable estimates of ␻ and ␣s (Tables 5 and 6). By contrast, for Ap4 two large correlations(|r| ⬎ 0.75) obtained for ␻ and ␣s indicate that for intermediate initial moisture, ␻ and ␣s may have been lesswell defined. For both MIM(Se) (Table 5) and MIM(h)(Table 6), correlations between ␪m,s and ␪im,s were always
large, indicating that their simultaneous estimation did
not yield reliable and unique parameter values. This inturn explains that no consistent effect of initial watercontent on ␪m,s and ␪im,s could be found.

Experimental methods for independent, approximate
experimental determination of ␪m and ␪im have previously been proposed (e.g., Clothier et al., 1995). However, we did not utilize these methods here since weintended to elucidate the full capacity of the inverse
parameter identification method by analyzing the extent
to which parameters can be identified at the same time
without resorting to additional, time-consuming experimental methods. Apart from the correlation between
␪m and ␪im, correlations between parameters obtainedwith MIM(h) were generally moderate for Ap1 andAp6: 96% were below 0.75, and 71% (Ap1) and 82%(Ap6) were less than 0.5. For Ap4, high correlation
values for nim, ␣im, nm, and ␣m indicate that these parameters may not have been uniquely identified for intermediate initial water contents. Correlations for ␻* and ␣sdisplayed mostly moderate values, indicating proper parameter identification as required for calculations of the
PE index.
Simulations using MIM(Se) and MIM(h) were repeated for Ap1, Ap4, and Ap6 by increasing each initial
parameter by 30% of its previous value. In general, final
optimized parameters were within ⫾10% of their prior
optimized values, thus indicating reliable parameter
identification for Ap1 and Ap6. As an example, thesimulation for Ap1 using MIM(Se) with 30% higher
initial parameter values gave the following final estimates (previous optimization result are given in parentheses): ␪m,s ⫽ 0.128 (0.115); ␣m ⫽ 0.09 (0.07) cm⫺1; nm ⫽1.22 (1.25); Ks ⫽ 25.8 (22.0) cm d⫺1; ␪im,s ⫽ 0.312 (0.310);
␻ ⫽ 0.41 (0.40) d⫺1; ␭m ⫽ 1.17 (1.08) cm; ␣s ⫽ 2.29 (2.20)
d⫺1, while for MIM(h) the following final estimates wereobtained for Ap1: ␪m,s ⫽ 0.150 (0.147); ␣m ⫽ 0.059 (0.034)cm⫺1; nm ⫽ 1.25 (1.25); Ks ⫽ 22.5 (22.5) cm d⫺1; ␪im,s ⫽0.242 (0.214); ␣im ⫽ 0.021 (0.016) cm⫺1; nim ⫽ 1.68 (2.00);
␻ ⫽ 0.06 (0.19) d⫺1; ␭m ⫽ 2.04 (2.26) cm; ␣s ⫽ 1.68 (1.34)

1319

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Table 6. Correlation, r, between MIM(h ) parameters for inverse simulation of the Ap1, Ap4, and Ap6 transport studies.†
Ap1
␪m,s
␣m
nm
Ks
␪im,s
␣im
nim
␻*
␭m
␣s
Ap4
␪m,s
␣m
nm
Ks
␪im,s
␣im
nim
␻*
␭m
␣s
Ap6
␪m,s
␣m
nm
Ks
␪im,s
␣im
nim
␻*
␭m
␣s

␪m,s

␣m

nm

Ks

␪im,s

␪im

nim

␻*

␭m

␣s

1
⫺0.66
⫺0.55
0.17
(⫺0.96)
0.61
⫺0.16
0.61
0.52
⫺0.32

1
0.14
⫺0.57
0.56
⫺0.58
0.12
⫺0.38
⫺0.12
0.21

1
⫺0.15
0.64
⫺0.41
⫺0.23
⫺0.05
⫺0.14
0.17

1
⫺0.13
0.25
0.38
⫺0.32
⫺0.06
⫺0.12

1
⫺0.49
0.05
⫺0.54
⫺0.46
0.25

1
⫺0.39
0.54
0.38
⫺0.42

1
(⫺0.83)
⫺0.10
0.05

1
0.36
⫺0.16

1
⫺0.25

1

1
⫺0.55
⫺0.73
(0.78)(⫺0.92)
⫺0.60
0.74
⫺0.05
⫺0.34
⫺0.21

1
0.42
(⫺0.79)
0.48
0.32
⫺0.53
0.25
0.22
0.09

1
⫺0.40
(0.83)(0.93)(⫺0.88)
⫺0.52
0.42
0.32

1
⫺0.66
⫺0.27
0.50
⫺0.46
⫺0.38
⫺0.07

1
0.74
(⫺0.86)
⫺0.12
0.43
0.27

1
(⫺0.90)
⫺0.53
0.39
0.29

1
0.21
⫺0.41
⫺0.38

1
⫺0.02
0.02

1
0.39

1

1
⫺0.13
⫺0.68
0.20
(⫺0.85)
0.01
0.51
⫺0.48
0.41
⫺0.22

1
⫺0.05
(0.85)
0.27
⫺0.19
0.04
0.14
⫺0.40
0.09

1
⫺0.26
0.74
⫺0.06
⫺0.65
0.43
⫺0.48
0.20

1
0.04
⫺0.22
0.40
0.24
⫺0.22
⫺0.01

1
0.18
⫺0.59
0.69
⫺0.44
0.21

1
⫺0.46
0.36
0.26
⫺0.09

1
⫺0.24
0.26
⫺0.15

1
⫺0.14
0.05

1
0.11

1

† Italics denote |r| ⬎ 0.5; parentheses denote |r| ⬎ 0.75.

d⫺1. For Ap4, parameter identification using MIM(Se)with 30% higher initial parameter values was somewhat
less accurate than for Ap1 and Ap6: ␪m,s ⫽ 0.178 (0.128);␣m ⫽ 0.049 (0.048) cm⫺1; nm ⫽ 1.29 (1.26); Ks ⫽ 34.9
(21.6) cm d⫺1; ␪im,s ⫽ 0.216 (0.248); ␻ ⫽ 0.50 (1.32) d⫺1;
␭m ⫽ 2.35 (1.91) cm; ␣s ⫽ 10 (10) d⫺1. The inverse
solution of MIM(h) for input parameters increased by
30% was not completed within 12 d, at which time weaborted the calculation.
The inverse MIM approaches overall gave acceptable
representations of the data. Most of the 8 [MIM(Se)]and 10 [MIM(h)] flow and transport parameters could
be reasonably well identified at the same time for initially wet and dry conditions. The parameters ␪m,s and␪im,s could not be reliably identified at the same timeand hence should be determined independently. For
medium initial water contents, large correlations between several parameters and dependence of the final
on the initial parameter values demonstrated that simultaneous fitting did not represent a well-posed problem for both MIM approaches. This is because medium
initial moisture conditions produced near-equilibriumtransport where MIM parameters are less well identifiable due to overparameterization of the model. It has
previously been found that for soils with spherical aggregates, MIM parameters are not always identifiable frommiscible displacement experiments (Vanderborght etal., 1997). The problem of parameter nonuniqueness isoften caused not only by the large number of optimized
parameters, but also by the lack of sufficient experimental information. Parameter identifiability can therefore
be increased either by reducing the number of optimized

parameters, or by including appropriate additional experimental information (Hopmans and Šimůnek, 1999),
such by using simultaneously information about multiple variables. Note that we defined the objective function using four different sets of variables (i.e., pressure
heads, water contents, cumulative water flow, and soluteconcentrations). This combination of variables yieldeda relatively good identification of most parameters inour study.

Numerical Discretization Effects on
Simulation Results
For the above comparison of inverse MIM(h) andMIM(Se) approaches, we always used the same smallvalues for the maximum time step dtmax (0.001 d) andthe space increment dz (0.074 cm). Numerical pressure
head oscillations for MIM(h) occurred when the time
step adaptive iteration scheme was not constrained bya small value for dtmax and calculated larger time steps
during drainage. A forward simulation for the Ap5 column shows that numerical oscillations in the pressure
head, h, were more effectively reduced by decreasing
dtmax than by decreasing dz (Fig. 10). Keeping dz at 0.5
cm (29 finite elements), the oscillations in h could be
eliminated by reducing dtmax from 0.03 to 0.01 d (Fig.
10a, 10b, 10c). Constraining dtmax to 0.03 d and reducing
dz to 0.037 cm (400 elements) did not eliminate oscillations (Fig. 10d, 10e, 10f). Even for the largest selecteddz (0.5 cm) and dtmax (0.3 d) the MIM(Se) approach did
not show oscillations in h (Fig. 10a). The need forsmaller time steps to eliminate oscillations makes theReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1320

VADOSE ZONE J., VOL. 3, NOVEMBER 2004

large values for the related PE index. The respective
parameters for MIM(Se) displayed similar tendencies
but were less consistent. Overall, the MIM(Se) inverse
approach was found to be well suited for simulating
physical nonequilibrium solute transport during variably saturated flow. Results indicate that the MIM(h)approach can yield higher accuracy along with moreconsistent parameter values if the observations containsufficient information for successful fitting. A future
theoretical analysis using synthetic datasets could complement this study to further clarify such inverse optimization issues as parameter sensitivity and uniquenessfor both MIM(Se) and MIM(h).
Fig. 10. Pressure heads for the Ap5 column as simulated with MIM(Se)and MIM(h ) assuming different maximum time steps (dtmax) while
keeping the space increment (dz ) constant at (a, b, c) 0.5 cm and(d, e, f) assuming different space increments while keeping time
steps constant at 0.03 d.

MIM(h) approach more time-consuming than the MIM(Se) approach. Numerical oscillations for MIM(h) wereclearly caused by the mass transfer term, which for larger
time steps could lead to overshooting of the equilibrium
state between the two pore regions. Numerical oscillations could have been eliminated if the time stepping
scheme had been adjusted to prevent this overshooting.
This was done only for the MIM(Se) approach for which
the adjustment was trivial due to the similarity assumption between the two pore regions, and for which accordingly no oscillations were registered.

CONCLUSIONS
Including experimental data of the pressure head,water content, cumulative outflow and effluent concentration into the objective functions of the MIM(Se) andMIM(h) approaches allowed us to identify model parameters, most of which appeared to be reliable as reflected by relatively small confidence intervals, generally moderate linear correlation between parameters,
and only minor dependence of the final optimized parameter values on initial parameter estimates. As an
exception, strong parameter correlations were foundbetween ␪m,s and ␪im,s, indicating that these parameters
should be determined independently. In spite of having
two more adjustable parameters, MIM(h) (10 parameters) produced results that were superior to those of
MIM(Se) (8 parameters) only for dry initial conditions,which indicates that at least for other than dry initialconditions, the inverse problem for MIM(h) was not
well-posed. Numerical oscillations in the pressure heads
for MIM(h) could be eliminated by reducing the maximum time step. The results further show that for aggregated soils, physical nonequilibrium can be more pronounced for wet and dry than for intermediate initialwater contents. Even if the initial water content influences the transport process, no distinct effect was foundon total nonreactive solute loss. For MIM(h), physicalnonequilibrium was reflected by smaller values for thewater and solute transfer coefficients and relatively

REFERENCES
Atkins, P.W. 1990. Physical chemistry. 4th ed. W.H. Freeman, New York.
Bahr, J.M., and J. Rubin. 1987. Direct comparison of kinetic and local
equilibrium formulations for solute transport affected by surface
reactions. Water Resour. Res. 23:438–452.
Castiglione, P., B.P. Mohanty, P.J. Shouse, J. Šimůnek, M.Th. van
Genuchten, and A. Santini. 2003. Lateral water diffusion in an
artificial macroporous system: Modeling and experimental evidence. Available at www.vadosezonejournal.org. Vadose Zone J.
2:212–221.
Chen, C., and R.J. Wagenet. 1992. Simulation of water and chemicals
in macropore soils. Part I. Representation of the equivalent macropore influence and its effect on soilwater flow. J. Hydrol. (Amsterdam) 130:105–126.
Clausnitzer, V., and J.W. Hopmans. 1995. Non-linear parameter estimation: LM_OPT. General-purpose optimization code based on
the Levenberg-Marquardt algorithm. Land, Air and Water resources Paper 100032. University of California, Davis.
Clothier, B.E., and S.R. Green. 1994. Rootzone processes and efficientuse of irrigation water. Agric. Water Manage. 25:1–12.
Clothier, B.E., L. Heng, G.N. Magesan, and I. Vogeler. 1995. The
measured mobile-water content of an unsaturated soil as a functionof hydraulic regime. Aust. J. Soil Res. 33:397–414.
Clothier, B.E., I. Vogeler, S.R. Green, and D.R. Scotter. 1998. Transport in unsaturated soil: Aggregates, macropores, and exchange.
p. 273–295. In H.M. Selim and L. Ma (ed.) Physical nonequilibrium
in soils. Ann Arbor Press, Chelsea, MI.
Damköhler, G. 1936. Einflüsse der Strömung, Diffusion und des Wärmeüberganges auf die Leistung von Reaktionsöfen. Z. Elektrochem. 42:846–862.
Dubus, I.G., and C.D. Brown. 2002. Sensitivity and first-step uncertainty analyses for preferential flow model MACRO. J. Environ.
Qual. 31:227–240.
Durner, W. 1993. Predicting the unsaturated hydraulic conductivity
using multi-porosity water retention curves. In M.Th. van Genuchten et al. (ed.) SHYPFIT, user’s manual. Res. Rep. 93.11, Soil
Physics. ITOE, ETH Zürich.
Flühler, H., W. Durner, and M. Flury. 1996. Lateral solute mixing
processes–A key for understanding field-scale transport of waterand solutes. Geoderma 70:165–183.
Flury, M., H. Flühler, W.A. Jury, and J. Leuenberger. 1994. Susceptibility of soils to preferential flow of water: A field study. Water
Resour. Res. 30:1945–1954.
Flury, M. 1996. Experimental evidence of transport of pesticides
through field soils—A review. J. Environ. Qual. 25:25–45.
Gerke, H.H., and J.M. Köhne. 2004. Dual-permeability modeling ofpreferential bromide leaching from a tile drained glacial till agricultural field. J. Hydrol., accepted.
Gerke, H.H., and M.Th. van Genuchten. 1993. A dual-porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29:305–319.
Granovsky, A.V., E.L. McCoy, W.A. Dick, M.J. Shipitalo, and W.M.
Edwards. 1993. Water and chemical transport through long-term
no-till and plowed soils. Soil Sci. Soc. Am. J. 57:1560–1567.
Gwo, J.P., P.M. Jardine, G.V. Wilson, and G.T. Yeh. 1995. A multiple-

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

pore-region concept to modeling mass transfer in subsurface media.
J. Hydrol. (Amsterdam) 164:217–237.
Hopmans, J.W., and J. Šimůnek. 1999. Review of inverse estimationof soil hydraulic properties. p. 643–659. In M.Th. van Genuchten
et al. (ed.) Characterization and measurement of the hydraulicproperties of unsaturated porous media. University of California, Riverside.
Jarvis, N., P.-E. Jansson, P.E. Dik, and I. Messing. 1991. Modelling
water and solute transport in macroporous soils. I. Model description and sensitivity analysis. J. Soil Sci. 42:59–70.
Jarvis, N.J., and P.B. Leeds-Harrison. 1987. Modelling of water movement in drained clay soil. I. Description of the model, sample
output and sensitivity analysis. J. Soil Sci. 38:487–498.
Hutson, J.L., and J.R. Wagenet. 1995. A multiregion model describing
water flow and solute transport in heterogeneous soils. Soil Sci.
Soc. Am. J. 59:743–751.
Kätterer, T., B. Schmied, K.C. Abbaspour, and R. Schulin. 2001.
Single- and dual-porosity modeling of multiple tracer transport
through soil columns: Effects of initial moisture and mode of application. Eur. J. Soil Sci. 52:25–36.
Kohler, A., K.C. Abbaspour, M. Fritsch, M.Th. van Genuchten, and
R. Schulin. 2001. Simulating unsaturated flow and transport in amacroporous soil to tile drains subject to an entrance head: Modeldevelopment and preliminary evaluation. J. Hydrol. (Amsterdam)
254:67–81.
Köhne, J.M., H.H. Gerke, and S. Köhne. 2002a. Effective diffusion
coefficients of soil aggregates with surface skins. Soil Sci. Soc. Am.
J. 66:1430–1438.
Köhne, J.M., S. Köhne, and H.H. Gerke. 2002b. Estimating soil hydraulic functions of dual-permeability models from bulk soil data.
Water Resour. Res. 38(7). DOI:10.1029/2001WR000492.
Larsson, M.H., N.J. Jarvis, G. Torstensson, and R. Kasteel. 1999.
Quantifying the impact of preferential flow on solute transport
to tile-drains in a sandy field soil. J. Hydrol. (Amsterdam) 215:
116–134.
Logsdon, S. 2002. Determination of preferential flow model parameters. Soil Sci. Soc. Am. J. 66:1095–1103.
Lennartz, B., J. Michaelsen, W. Wichtmann, and P. Widmoser. 1999.
Time variance analysis of preferential solute movement at a tiledrained field site. Soil Sci. Soc. Am. J. 63:39–47.
Medina, A., and J. Carrera. 1996. Coupled estimation of flow andsolute transport parameters. Water Resour. Res. 32:3063–3076.
Meyer-Windel, S. 1998. Transport processes of bromide and herbicides
under steady-state and transient flow conditions. Ph.D. diss. Schriftenreihe Institut für Wasserwirtschaft und Landschaftsökologie der
Christian-Albrechts-Universität zu Kiel.
Meyer-Windel, S., B. Lennartz, and P. Widmoser. 1999. Bromide andherbicide transport under steady-state and transient flow conditions. J. Eur. Soil Sci. 50:23–33.
Michalak, A.M., and P.K. Kitanidis. 2000. Macroscopic behavior andrandom-walk particle tracking of kinetically sorbing solutes. Water
Resour. Res. 36:2133–2146.
Millington, R.J., and J.P. Quirk. 1961. Permeability of porous soils.
Trans. Faraday Soc. 57:1200–1207.
Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van Genuchten. 1997. New piecewise-continuous hydraulic functions for1321

modeling preferential flow in an intermittent-flood-irrigated field.
Water Resour. Res. 33:2049–2063.
Roulier, S., and N. Jarvis. 2003. Analysis of inverse procedures for
estimating parameters controlling macropore flow and solute transport in the dual-permeability model MACRO. Available atwww.vadosezonejournal.org. Vadose Zone J. 2:349–357.
Russo, D., W.A. Jury, and G.L. Butters. 1989. Numerical analysis ofsolute transport during transient irrigation. 2. The effect of immobile water. Water Resour. Res. 25:1785–1794.
Schwartz, R.C., A.S.R. Juo, and K.J. McInnes. 2000. Estimating parameters for a dual-porosity model to describe non-equilibrium,
reactive transport in a fine-textured soil. J. Hydrol. (Amsterdam)
229:149–167.
Shipitalo, M.J., W.M. Edwards, W.A. Dick, and L.B. Owens. 1990.
Initial storm effects on macropore transport of surface appliedchemicals in no-till soil. Soil Sci. Soc. Am. J. 54:1530–1536.
Šimůnek, J., D. Jacques, J.W. Hopmans, M. Inoue, M. Flury, andM.Th. van Genuchten. 2002. Solute transport during variably-saturated flow—Inverse methods. p. 1435–1449. In J.H. Dane and G.C.
Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5.
SSSA, Madison, WI.
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs.
2003. Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadose zone. J.
Hydrol. (Amsterdam) 272:14–35.
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D software package for simulating the one-dimensional
movement of water, heat, and multiple solutes in variably-saturated
media. Version 2.0. IGWMC- TPS- 70. International Ground Water
Modeling Center, Colorado School of Mines, Golden.
Šimůnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated
soil hydraulic properties from tension disc infiltrometer data bynumerical inversion. Water Resour. Res. 32:2683–2696.
Tillman, R.W., D.R. Scotter, B.E. Clothier, and R.E. White. 1992.
Solute transport during intermittent water flow in a field soil andsome implications for irrigation and fertilizer application. Agric.
Water Manage. 20:119–133.
Vanderborght, J., D. Mallants, M. Vanclooster, and J. Feyen. 1997.
Parameter uncertainty in the mobile–immobile solute transportmodel. J. Hydrol. (Amsterdam) 190:75–101.
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., and P.J. Wierenga. 1976. Mass transfer studies
in sorbing porous media. I. Analytical solutions. Soil Sci. Soc.
Am. J. 40:473–481.
White, R.E., J.S. Dyson, Z. Gerstl, and B. Yaron. 1986. Leaching of
herbicides through undisturbed cores of a structured clay soil. Soil
Sci. Soc. Am. J. 50:277–283.
Wilson, G.V., J.P. Gwo, P.M. Jardine, and R.J. Luxmoore. 1998.
Hydraulic and Physical Nonequilibrium Effects on Multiregion
Flow. p. 37–61. In H.M. Selim and L. Ma (ed.) Physical nonequilibrium in soils. Ann Arbor Press, Chelsea, MI.
Zurmühl, T., and W. Durner. 1998. Determination of parameters forbimodal hydraulic functions by inverse modeling. Soil Sci. Soc.
Am. J. 62:874–880.