Journal of Hydrology 403 (2011) 141–156

Contents lists available at ScienceDirect

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

Coupled simulation of surface runoff and soil water ﬂow using multi-objective
parameter estimation
John Maximilian Köhne a,⇑, Thomas Wöhling b,c, Valérie Pot d, Pierre Benoit d, Sophie Leguédois e,f,
Yves Le Bissonnais e,g, Jirka Šimůnek haHelmholtz Centre for Environmental Research – UFZ, Theodor-Lieser-Straße 4, D-06120 Halle (Saale), Germany
Lincoln Ventures Ltd., Lincoln Environmental Research Division, Private Bag 3062, Hamilton, New Zealand
Water and Earth Science Systems Research Centre, University of Tübingen, Institute for Geosciences, Sigwartstrasse 10, 72076 Tübingen, GermanydINRA, AgroParisTech, UMR 1091 EGC, F-78850 Thiverval-Grignon, FranceeINRA, UR0272, Unité de Science du Sol, Centre de Recherche Orléans, F-45075 Orléans Cedex 2, FrancefINRA, UMR Laboratoire Sols et Environnement, Nancy-Université, 2 Avenue de La Forêt de Haye, BP 172, 54505 Vandoeuvre-lès-Nancy, FrancegINRA, UMR LISAH, 2 Place Pierre Viala, F-34060 Montpellier cedex 1, FrancehDepartment of Environmental Sciences, University of California Riverside, Riverside, CA, USAbca r t i c l e

i n f oArticle history:
Received 28 September 2010
Received in revised form 7 March 2011
Accepted 1 April 2011
Available online 7 April 2011
This manuscript was handled by Dr.
A. Bardossy, Editor-in-Chief, with theassistance of Vazken Andréassian, Associate
Editor
Keywords:
Overland ﬂow
Surface runoff
Multi-objective global parameter
optimization
Mobile–immobile model
Simulation
Flow channel experiments u m m a r y
A comprehensive description of water ﬂow in environmental and agricultural systems requires anaccount of both surface and subsurface pathways. We present a new model which combines a 1D overland ﬂow model and the 2D subsurface ﬂow HYDRUS-2D model, and uses the multi-objective global
search method AMALGAM for inverse parameter estimation. Furthermore, we present data frombench-scale ﬂow experiments which were conducted with two 5-m long replicate soil channels. Whilerainfall was applied, surface runoff was recorded at the downstream end of the soil channel, subsurface
drainage waters were sampled at three positions equally spaced along the channels, and pressure heads
were recorded at ﬁve depths. The experimental observations were used to evaluate the performance ofour modeling system. The complexity of the modeling approach was increased in three steps. First, only
runoff and total drainage were simulated, then drainage ﬂows from individual compartments were additionally evaluated, and ﬁnally a surface crust and immobile soil water were also considered. The results
showed that a good match between measured and observed surface runoff and total drainage does notguarantee accurate representation of the ﬂow process. An inspection of the Pareto results of different
multiobjective calibration runs revealed a signiﬁcant trade-off between individual objectives, showingthat no single solution existed to match spatial variability in the ﬂow. In spite of the observed crust formation, its consideration in the more complex model structure did not signiﬁcantly improve the ﬁt
between the model and measurements. Accounting for immobile water regions only slightly improved
the ﬁt for one of the two replicate soil channels. Discrepancies between relatively complex model simulations and seemingly simple soil channel experiments suggest the presence of additional unknowns,
such as heterogeneity of the soil hydraulic properties. Nevertheless, with its versatile subsurface options
and powerful inverse method, the model system shows promise for studying hillslope ﬂow problems
involving both surface runoff and subsurface ﬂow.
Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction
Overland ﬂow during storm events may affect the water andmass balance of environmental and agricultural systems. However,current vadose zone ﬂow models include at most a simpliﬁed
description of this process, such as instant removal of saturation/
inﬁltration excess water as surface runoff (e.g., Ahuja et al.,⇑ Corresponding author. Tel.: +49 345 5585 5406.
E-mail address: max.koehne@ufz.de (J.M. Köhne).
0022-1694/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2011.04.001

2000). As an alternative, combined models of overland and subsurface ﬂow were developed.
The diffusion wave or kinematic wave approximations of theSaint-Venant equation are frequently used to describe overlandﬂow, while the Green-Ampt or Richards equations are used tomodel inﬁltration or subsurface ﬂow. Earlier approaches usuallyrelied on analytical solutions, limiting their applications torelatively simple boundary conditions and homogeneous soil
(Govindaraju et al., 1990; Wallach et al., 1997). Numerical solutions
then allowed various generalizations, such as time-dependent
boundary conditions (e.g., Govindaraju and Kavvas, 1991; Motha142

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

and Wigham, 1995; Woolhiser et al., 1990) or spatial variability ofsoil hydraulic properties (e.g., Corradini et al., 1998; Merz andPlate, 1997), and allowed consideration of runon in addition torainfall (e.g., Nahar et al., 2004). While many approaches reliedon 1D overland ﬂow and pseudo-2D (staggered independent 1Dvertical columns) subsurface ﬂow, extensions to two dimensions
were also used for both overland ﬂow (e.g., Gandolﬁ and Savi,
2000; Merz and Plate, 1997) and subsurface ﬂow (e.g., Govindaraju
and Kavvas, 1991; Singh and Bhallamudi, 1998). A 1D surface andquasi-3D subsurface (2D planes) model was introduced by
Wöhling and Schmitz (2007) and successfully applied by Wöhlingand Mailhol (2007). A comprehensive coupled 2D surface and 3D
subsurface watershed model was introduced by Panday andHuyakorn (2004). A few models also considered the effect ofinﬁltration and vertical gravity ﬂow in macropores (e.g., Léonard
et al., 2001; Merz and Plate, 1997), also in combination withheterogeneity in soil hydraulic properties (Herbst et al., 2006). A
review of various surface–subsurface ﬂow approaches has recently
been presented by Furman (2008).
None of the above studies focused on inverse parameter identiﬁcation, and most approaches do not even have a provision for inverse estimation. However, while prediction of vadose zone
dynamics is already a challenge, adding surface ﬂow to the model
analysis further increases the complexity. The main problem is tocharacterize the spatial variability of soil hydraulic properties
caused by variations in soil texture, structure, and layer interfaces(Kodesová et al., 2009; Kulli et al., 2003; Vogel et al., 2005), among
other effects. The difﬁculty of independent model parameterization was the driving motivation behind exploring the beneﬁt ofusing inverse parameter estimation in this study.
For parameter estimation of vadose zone models, local searchmethodologies are often employed that begin their iterative searchfrom a single initial point in the parameter space. A prominent
example is the nonlinear, gradient-based Levenberg–Marquardt
search algorithm (Marquardt, 1963). However, such local derivative-based search methods often fail to evolve towards the global
optimum in situations where the response surface exhibits multiple local optima in the parameter domain (e.g., Schwefel, 1993). In
such cases, the ﬁnal solution of the local search method is nonunique and essentially dependent on the starting point in theparameter space. While multi-start local search methods may help
to overcome this problem under some circumstances, the convergence to the global solution cannot be guaranteed. As a responseto this problem, global search algorithms were developed thatuse multiple concurrent searches from different starting points
to efﬁciently reduce the chance of getting stuck in a single areaof attraction (e.g., Deb et al., 2002; Madsen, 2003; Vrugt andRobinson, 2007).
Since any model is a simpliﬁcation of reality and the underlying
reality is never known exactly, model simulation results are inherently uncertain. One expression of this uncertainty is equiﬁnality,
characterized by multiple model structures and parameter setsyielding acceptable ﬁts to observed data (Beven, 2006). To better
understand the uncertainty and limitations of the model structure,
the optimization problem can be formulated in a multi-objective
context, where information contained in several data-sets is
exploited simultaneously (Gupta et al., 1998; Vrugt et al., 2003;
Wöhling and Vrugt, 2008). Among different multi-objective global
parameter optimization procedures, A MultiALgorithm Genetically
Adaptive search Method (AMALGAM) developed by Vrugt andRobinson (2007) was found to be the most efﬁcient algorithm forﬁnding best attainable parameter sets with the least number ofmodel evaluations (Wöhling et al., 2008).
The overall objective is to evaluate the usefulness of linking asurface ﬂow model with a comprehensive vadose zone model(Simunek et al., 2003, 2008) and advanced multi-objective optimi-

zation techniques for studying the interacting ﬂow-runoff behavior
in a sloping soil. In this context, data from the controlled benchscale ﬂow experiment were utilized.

2. Material and methods
2.1. Overland ﬂow model
The overland ﬂow model is based on unpublished work by Šimůnek (2003). Hortonian overland ﬂow is typically described using
the kinematic wave equation, which is a simpliﬁcation of theSaint-Venant equations and an excellent approximation of mostoverland ﬂow conditions (Morris and Woolhiser, 1980; Singhet al., 2005; Woolhiser et al., 1990):@ho @Q o
þ
¼RI
@t
@x

ð1Þ

where ho(x, t) is the unit storage of water (or mean depth forsmooth surfaces) in the overland ﬂow (L), Qo(x, t) is the overlanddischarge per unit width (L2 T1), t is time (T), x is the distance coordinate over the soil surface (L), R(x, t) is the local effective rainfall
rate (rain minus evaporation) (L T1), and I(x, t) is the local inﬁltration (or exﬁltration) rate (L T1). The discharge Qo per unit width
can be calculated as follows:m

Q o ¼ aho

ð2Þ

where a (L2m T1) and m (–) are parameters related to slope, surface roughness, and ﬂow conditions (laminar or turbulent ﬂow).
Overland ﬂow is often turbulent, with a large Reynolds number
Re = Qo/v > 500, where v is the kinematic viscosity of water(L2 T1). There is some evidence that turbulent ﬂow conditions
may prevail during rainfall even at lower Reynolds numbers (Denget al., 2005). Flow for our experimental conditions was assumed tobe turbulent and the coefﬁcients in Eq. (2) were thus evaluatedusing the empirical Manning hydraulic resistance law and assuming
m = 5/3 (Wallach et al., 2001):

a¼k

S1=2
nM

ð3Þ

where S (–) is the slope, nM (–) is the Manning’s roughness coefﬁcient for overland ﬂow, and k is a conversion constant which in SI
units has the value of 1 and units of (L1/3 T1). From Eqs. (2) and(3) it can be seen that the nM coefﬁcient depends not only on surface
roughness, but also on the slope and the overland discharge. For
bare ﬁeld soils, nM typically ranges between 0.01 for a rathersmooth surface (ﬂat slope and slow overland discharge), and
0.2 for a rough soil surface (steeper slope and faster discharge)(Li and Zhang, 2001; Mohamoud, 1992; Sepaskhah and Bondar,
2002).
By substituting (2) into (1), and substituting R  I by the net
rate of local input, q(x, t) (L T1), the following equivalent expressions describing overland ﬂow are obtained:m@ho @ ahoþ
¼ qðx; tÞ
@t
@x

ð4Þ

@ho
m1 @hoþ amho¼ qðx; tÞ
@t
@x

ð5Þ

A numerically stable, fully implicit four-point ﬁnite difference
method, similar to the one used in the Kineros model (Woolhiser
et al., 1990), was used to solve Eq. (1). The following discretization
scheme was applied (subscript o is omitted for convenience).

143

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156
jþ1

j

!

jþ1
j
Q jþ1  Q jþ1
Q j  Q ji
1 hi  hi hiþ1  hiþ1i
þ e iþ1þþ ð1  eÞ iþ1
2
Dt
Dt
Dx
Dx


1
1
j
j
¼ e qjþ1þ qjþ1iiþ1 þ ð1  eÞðqi þ qiþ1 Þ
2
2

ð6Þ

Subscripts i are nodal numbers, superscripts j and j + 1 represent previous and current times, Dt and Dx are temporal and spatial discretization steps, and e is the temporal weighting
coefﬁcient; e = 1 for the implicit scheme (the Crank–Nicholson
scheme with e = 0.5 gave similar results). The resulting nonlinear
system of equations is solved using the Picard iterative approach,similar to the one used for the solution of the Richards equationin HYDRUS-2D (Šimůnek et al., 1999). The numerical model was
validated by comparing its results with the corresponding analytical solution of Eqs. (4), (5) (Šimůnek, 2003).
The one-dimensional equation for overland water ﬂow requires
initial conditions, plus up to two boundary conditions; one for runon (if any) at the upstream entrance and another for runoff at thedownstream exit. The initial conditions for overland ﬂow are represented by the initial ponding depths, ho0(x), assumed to be equal
to the initial soil water pressure heads, h0(x, zsurf), along the soil
surface with the z-coordinate zsurf (L):ho ðx; 0Þ ¼ ho0 ðxÞ ¼ h0 ðx; zsurf Þ for x; zsurf eCO

ð7Þ

when the pressure head is positive (the soil surface is ponded), or
equal to zero when the pressure head is negative (the soil surface
is unsaturated).

Fig. 1. Conceptual schematic of the coupled runoff-subsurface ﬂow model based on
HYDRUS-2D. An immobile region was assumed in the ‘MIM + Crust’ scenario,
whereas in the remaining scenarios water in the entire soil domain was assumed tobe mobile.

The boundary conditions for water ﬂow represent water runon
at the entrance side (a Dirichlet boundary condition) with a prescribed ponding head, and a zero gradient for the ponding depth
is assumed at the outlet side of the plot:ho ð0; tÞ ¼ ho0

for x ¼ x0 ;

t>0

ð8Þ

@ho ðn; tÞ¼ 0 for x ¼ xn ;
@x

t>0

ð9Þ

If there is no runon, ho0(t) in (8) becomes zero. Since the lateralrunoff discharge, Qo, is calculated as a function of ho according to(2), the runon is fully deﬁned by Eq. (8), which is therefore equivalent with a ﬂux-type (Neumann) condition.
The numerical solution of the overland ﬂow equation was iteratively coupled with HYDRUS-2D (Šimůnek et al., 1999, 2008).
Interactions with the HYDRUS-2D subsurface ﬂow model are represented by the sink-source term I in Eq. (1). Coupling to overland
ﬂow is accomplished via dynamic atmospheric boundary conditions in the computation module of HYDRUS-2D as follows. During
unsaturated conditions, there is no surface runoff and water inﬁltration is similar to the applied rainfall rate, which is describedusing a speciﬁed ﬂux (Neumann) upper boundary condition. If
the pressure at the surface exceeds zero, i.e. ponding occurs, thenthe upper boundary switches to a speciﬁed head (Dirichlet) condition (Šimůnek et al., 1999). Inﬁltration, I, is then calculated in thesubsurface model using the Richards equation, where the pressure
head in the subsurface domain at the upper boundary, h0(x, zsurf), is
assumed to be equal to the ponding water depth, ho0(x), which inturn is calculated as overland ﬂow depth according to (1). The
same Neumann–Dirichlet switch in boundary conditions was
implemented for runon at any location.
Finally, each drainage water outlet is represented by a single
node at the bottom of the transport domain with a seepage faceboundary condition. During each iteration, a saturated node of
the potential seepage face is treated as a prescribed pressure headboundary with h = 0, while an unsaturated node is treated as a prescribed ﬂux boundary with Q = 0. While in the former case the unknown boundary (drainage) ﬂux Q is calculated, in the latter casethe unknown pressure head is calculated.
The soil structure generally invokes zones of different water
mobility in the soil. Often, preferential ﬂow paths form highly mobile regions, contrasting with relatively immobile zones in the soil
matrix. We coupled the model for overland ﬂow with the dualporosity mobile-immobile water model (assuming the saturation-based water transfer) in HYDRUS-2D. Then the pressure headin the mobile region at the surface is assumed to be equal to the
ponding water depth, ho0(x). A schematic of the conceptual modelis shown in Fig. 1. The model was set-up to represent the experimental system (Fig. 2).

Tensiometers
Drain 1a

1m
1m

0.3 m

geotextile

Drain 2a

BoxR

BoxL
5m
r
ato

Upstream separ
1.4 m
5% slope

Drain 1

Tensiometers

t or
ara

Midstream sep
1.6 m

Drain 2

Downstream
2.0 m

Drain 3Fig. 2. Schematic overview of the soil channel. The grey lines denoted ‘separator’ represent metal sheets of a few cm height preventing lateral ﬂow along the bottom between
the up-, mid-, and down-stream segments.

144

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156
Table 1
Bulk densities of the Fontaine silt loam soil ﬁlled in the left (BoxL) and right (BoxR)
channels.
BoxL

aFig. 3. Photo of the soil ﬂow channel at INRA Orléans.

2.2. Runoff experiment2.2.1. Experimental set-up
A runoff experiment was conducted in the Soil Science Laboratory of the Institut National de la Recherche Agronomique (INRA)
at Orléans (France). The experimental set-up was made of two replicate soil-ﬁlled boxes integrated in a channel of 5 m length, 2 m
total width (1 m per replicate), and 0.3 m height, with a slope of
5% (Figs. 2 and 3). The left and right replicate boxes (Fig. 2) are subsequently referred to as BoxL and BoxR. They were ﬁlled withsieved soil obtained from a silt loam and silty clay loam, respectively. The experiment was designed to study the transport of pesticides with different sorbing properties as affected by interrill
erosion (Leguédois et al., 2006; Leguédois and Bissonnais, 2004).
In this study, however, we utilized only data from the silt loam soil,
which carries the local designation ‘Fontaine soil’ and is found in
the western part of the Paris Basin, Normandy (France). Fontaine
soil has an aggregated structure, but aggregates are unstable andthe soil is susceptible to erosion (Fig. 4) (Le Bissonnais et al.,
1998; Leguédois and Bissonnais, 2004). The soil material was sampled from the Ap horizon of a plot on cultivated arable land near
Fontaine (Bourville, Seine-Maritime), from a depth of 0–0.3 m.
Filling the soil boxes was accomplished in two steps, as follows.
First, a 25-cm thick layer was built up by pouring dry-sieved(<5 cm aggregates) soil in 5-cm increments. Natural compactionof this 25-cm layer was realized by alternately applying about 10
recurring cycles of rainfall (until saturation) and drainage. Carewas taken to avoid any ponding, erosion, or crust formation. Whenthe soil stopped settling, the remaining 5 cm to the top were ﬁlledwith sieved (<2 cm) soil in two 2.5-cm increments. This ﬁlling procedure invoked slight differences in the soil bulk density and watercontent with depth, as was measured on soil cores taken after theexperiment (Table 1). The soil surface was reshaped into a furrow-

BoxR

Depth (cm)qb (g/cm3)

haDepth (cm)qb (g/cm3)

ha

0–5
7.3–12.3
13.5–18.5
19.5–24.5
25–30

1.31
1.35
1.02
1.36
1.47

0.33
0.33
0.30
0.34
0.47

0–5
6.5–11.5
13–17
20–25
25–30

1.29
1.37
1.16
1.44
1.38

0.33
0.28
0.20
0.37
0.44

Water content measured a few hours after the end of the experiment.

ridge micro-topography (Leguédois and Bissonnais, 2004) parallelto the slope with a 2-cm height-amplitude and a 10-cm wavelength in a cross-section perpendicular to the slope (Fig. 4). The
channel box was equipped with ﬁve tensiometers at ﬁve soil
depths between 2 and 29 cm, centered at the mid position of
x = 2.5 m.
The bottom of the box was covered with highly permeable geotextile. At upstream (1.4 m), midstream (3 m), and downstream(5 m) locations, metal sheets a few cm high at the bottom of thebox impeded lateral ﬂow along the bottom. They served as gutterschanneling drainage water to the outlets (drains) representing different compartments, i.e., drain 1 – upstream (0–1.4 m), drain 2 –midstream (1.4–3 m), and drain 3 – downstream (3–5 m). Surface
runoff was collected at the downstream end.
Soil hydraulic properties were independently measured in thelab on undisturbed ﬁeld-sampled soil cores (1 dm3) using theWind’s (1966) method (Tamari et al., 1993). The Wind method
yields uncertain results close to saturation, since the derived saturated hydraulic conductivity, Ks, is an extrapolated value. Additionally, Ks was measured using a permeameter method on samples of1 dm3 in size (Table 2). The custom-made permeameter measured
the saturated water ﬂux for an imposed constant head as a function of time. The ﬁnal steady-state ﬂux was assumed to be equal
to Ks.

Table 2
Soil hydraulic van Genuchten–Mualem parameters for measurements on the Fontaine
silt loam soil (wind and permeameter measurements on undisturbed ﬁeld samples:two replicates; permeameter measurements on box soil samples: mean ± standard
deviation for 10 measurements).
Wind’s method (evaporation experiment)

hr
0.000
0.000
a
b

hs
0.483
0.464

n
1.36
1.56

Permeameter
Field
sampleBox sampleaqb(1/cm)(g/cm3)Ks(cm/d)

Ks(cm/d)

Ks(cm/d)

0.01120.0105

1.37
1.34

12.2
12.8

16.3
16.3

146.0 ± 43.6a
9.1 ± 5.5b

Soil samples below crust.
Soil with crust.

Fig. 4. The soil ﬂow channel with aggregates and furrow-ridge surface micro-topography, prior to rainfall (left) and thereafter (mid). Detailed view of a cross-section througha rill revealing the surface crust (right).

145

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156
Table 3Mass balance information for the left (BoxL) and right (BoxR) soil boxes at the end of the ﬁrst and second irrigation, respectively.
Time (min)Rain (l)Runoff (l)Drain 1 (l)Drain 2 (l)Drain 3 (l)DrainR (l)Storage (l)DrainR (%)

Runoff (%)

Storage (%)

First rain event
BoxL
32.5
BoxR
32.5

79.9
74.5

9.1
16.3

0.0
0.0

0.0
0.0

0.0
0.0

0.0
0.0

70.8
58.2

0.0
0.0

11.4
21.9

88.6
78.1

Second rain event
BoxL
100
BoxR
100

245.9
229.2

175.9
197.9

0.0
3.2

0.7
14.2

11.4
16.4

12.1
33.9

57.9
2.5

4.9
14.8

71.5
86.3

23.5
1.1

However, since hydraulic properties measured for undisturbed
soil may not represent the sieved soil material used for ﬁlling thesoil box, additional permeameter measurements were also conducted on soil samples (5 cm height and diam.) collected fromthe box after the experiments. Samples were taken from minifurrows with surface crust (11 measurements) and from small
ridges (10 measurements) (Table 2). Samples taken from ridges
had an order-of-magnitude higher Ks than samples taken fromfurrows. The surface crust thus strongly reduced Ks. However, Ksof samples from ridges were also much larger than Ks measured
on ﬁeld samples. A possible explanation is that the ﬁlling procedure may have caused development of inter-aggregate regionsthat had increased Ks. Moreover, a relatively small sample size
may have caused an overestimation of mean Ks and its variability(Mallants et al., 1997).
2.2.2. Experiment
Two rainfall events were applied to each of the replicate soilboxes. The rainfall simulator was designed to simulate rain eventsof spatially uniform intensity. It was made of ﬁve oscillating nozzles at 1.1 m horizontal spacing and was placed at 6.4 m above
the soil surface. The device produced raindrops with 1.3 mm median diameter and a kinetic energy of 16 J mm1 m2. The rainfall
intensity received by the box was veriﬁed by placing a plastic sheeton the soil surface and measuring the runoff. The spatial variabilityof the rain intensity was measured within the experimental area
using rainfall gauges; the coefﬁcient of variation was 15%. The
mass balance for both rainfall events and the two soil boxes is given in Table 3. The rainfall rate differed slightly between 29 mm/h
(147.6 l/h) for BoxL (the total surface area of 5.09 m2), and28.4 mm/h (137.5 l/h) for BoxR (the total surface area of
4.84 m2). The difference between left and right surface areas isdue to an imprecision when dividing the box into two parts byplacing a metal sheet of 5 m length in the middle. The width of
BoxL is 1.02 m (upstream) and 1.015 m (downstream), while for
BoxR it is 0.97 m (upstream) and 0.975 m (downstream),
respectively.
The ﬁrst rainfall was applied to pre-wet the soil and did not produce any drainage. Runoff started after 20 min and 21 min for theleft and right boxes, respectively. Rainfall was stopped after
32.5 min when the surface runoff rate reached 10% of the applied
rainfall rate. Rapid formation of a seal crust in the mini-furrows
was visually observed within 10 min of the ﬁrst simulated rainfall.
In the crust, boundaries between aggregates disappeared and thesurface appeared smoother and of lighter color (Fig. 4). After a
24-h interruption, the second rainfall of 100 min was applied,
using the same rate as for the ﬁrst event.
2.3. Multi-objective optimization
2.3.1. Theory
The aim of the inverse modeling procedure is to ﬁnd values ofmodel parameters that provide the best attainable ﬁt between
model simulation and corresponding observations. The multi-

objective framework used herein comprised three different
criteria:2

F 1 ðuÞ

3

6
7
min FðuÞ ¼ 4 F 2 ðuÞ 5
F 3 ðuÞ

ð10Þ

where F1 to F3 are deﬁned as the normalized Root Mean Square
Error (nRMSE) objective functions of the differences between
observed and simulated drainage ﬂuxes at the midstream (F1) andlower (F2) drains, and the surface runoff (F3) collected at the lower
end of the sand box, and u is the vector of model parameters to be
estimated. Here, u includes both the soil hydraulic van Genuchten–
Mualem (VGM) parameters (van Genuchten, 1980) and the Manning’s roughness coefﬁcient, nM, usually a not well-known overland
ﬂow parameter.
A normalization is required to compensate for differences inmagnitude and variability of the objective functions, i.e., to balance
their weights in the global search. The objective function values aretransformed so that the minima of the objectives in the initialsample of parameter sets exhibit a similar distance to the origin(Madsen, 2000; Mertens et al., 2004; Wöhling and Vrugt, in press):

Fi ¼

fi

ri

þ /i

i ¼ 1; 2; 3

ð11Þ

where fj is the non-transformed (RMSE) objective function, rj is thestandard deviation associated with fj of the initial sample, and /j is atransformation constant given by


 
fj
;
/i ¼ max min

rj


 
fi
j ¼ 1; 2; 3  min
i ¼ 1; 2; 3

ri

ð12Þ

The AMALGAM method (A MultiALgorithm Genetically Adaptive search Method), recently developed by Vrugt and Robinson(2007), was used for parameter estimation of the coupled surface–subsurface model. AMALGAM is a population-based global
search algorithm that combines two concepts of a simultaneous
multi-algorithm search method with genetically adaptive offspringcreation for estimating the Pareto set of solutions of a given multiobjective inverse problem. The Pareto optimal or Pareto efﬁcientsolutions represent trade-offs among different objectives, having
the property of moving from one solution to another results in
the improvement of one objective while causing deterioration in
one or more of the others (Gupta et al., 1998; Vrugt et al., 2003).
An illustrative summary of the basic features of AMALGAM appears
Table 4
Optimization scenarios. MIM – mobile and immobile regions. The + sign indicates anobjective.
Scenario

Uniform-2

Uniform-3

MIM + Crust

Objectives
RunoffTotal drainage
Midstream drainage
Downstream drainage

2
+
+



3
+

+
+

3
+

+
+

146

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

in terms of an equilibrium pressure proﬁle, which was estimatedby interpolating measured tensiometer values.
Scenarios with different complexities should reveal the effects
of data available for model optimization, and of particular soil
properties. Among the six scenarios originally studied, three distinct scenarios were selected for this paper (Table 4). The ‘Uniform’scenarios were built upon the assumption of homogeneous soil
hydraulic properties in the entire domain, represented by a single
set of VGM parameters. The ‘Uniform-2’ scenario includes twoin Wöhling et al. (2008), while a more detailed description can befound in Vrugt and Robinson (2007), and is therefore not repeatedhere.

2.3.2. Optimization procedure and scenarios
A reliable model simulation requires the realistic deﬁnition of
initial and boundary conditions. Boundary conditions were set
equal to the applied rainfall rate. Initial conditions were deﬁned

Table 5
Range of the optimized parameters.
hsna (1/cm)

Ks (cm/d)

sb

nM

hs,ma

hs,imaxa (1/d)

Min
Max

0.15
0.70

1.001
9 (4a)0.001
2 (1a)

0.000144
14,4003
30.0001
0.5

0.05
0.5

0.05
0.5

0
1

Rain, runoff (dm³/min)

Only in the MIM + Crust scenario.
Exponent in the van Genuchten Mualem model.

0.5

BoxL

3
2.5

0.4

2

0.3

0.5

BoxR

0.4

2

0.3

1.5

1.5
0.2

1

0.1

0.5

0
140

0
0

Cumulative flow (dm³)

3
2.5

20

40

60

80

100

120

0.2

1

0.1

0.5
0

1000

1000

100

100

10

10

1

1

0.1

0.1

0.01
0

20

40

60

80

100

120

140

0.01

0

20

40

20

40

Time (min)
Rain

60

0
80 100 120 14060

80 100 120 140

0Drainage (dm³/min)ab

Parameters

Time (min)RunoffDrain 2

Drain 1

Drain 3Pressure head (cm)Fig. 5. Experimental data for the BoxL (left) and BoxR (right) replicates as measured during the second rainfall event. Applied rainfall, surface runoff, and drainage at threedrains are given in terms of rates (top) and cumulative ﬂows (bottom), respectively.

0

0

-20

-20

-40

-40

-60

-60

-80

-80

3 cm

-100

7 cm
15 cm
23 cm

-100
BoxL

28 cm

BoxR

-120

-120
0

50

100

Time (min)

150

0

50

100

150

Time (min)Fig. 6. Pressure head response to rainfall as observed at different depths using tensiometers installed in the middle of the replicates BoxL and BoxR.

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

objective functions for runoff and total drainage, i.e., the sum of thethree individual drainage ﬂows. This reduces Eq. (10) to two components, i.e., F1 for drainage and F2 for runoff. The ‘Uniform-3’ scenario includes three objective functions as deﬁned in Eq. (10). Bothscenarios were compared to assess the potential importance ofusing spatially resolved subsurface ﬂow information.
The ‘MIM + Crust’ scenario considers mobile and immobile regions in a soil with a thin crust layer placed on top of the soil. Crustformation was observed during the experiment. The presence ofmobile and immobile regions was assumed to be an effective
way of accounting for any type of ﬂow heterogeneity, such as ﬁngering ﬂow in the layered soil system, inhomogeneous inﬁltration
due to the micro-funneled surface, or inter-aggregate percolation
through the aggregated soil. Furthermore, the ‘MIM + Crust’ optimization includes 3 objective functions. Additionally, a ‘Crust’(without MIM) scenario was simulated, but since its results did
not yield much additional information, we did not report them.
From the resulting Pareto surface, we selected four Pareto
points. The ﬁrst three Pareto points were the best individual solu-

147

tions with respect to the individual objectives, which are henceforth denoted as Pareto extremes, F^i . There is one Pareto extreme
for midstream drainage (F^1 ), one for downstream drainage (F^2 ),and one for runoff (F^3 ). The extremes are a subset of the Pareto
optimal or efﬁcient solutions deﬁned above. The fourth Pareto
point is subsequently referred to as the compromise solution,
where the sum of all three normalized root mean square error
objectives, Fi, is at its minimum (F^0 ):F^0 ¼ min3X

F i ðuÞ

ð13Þ

i¼1

Only data for the second rain event for both BoxL and BoxR
were used, since the ﬁrst rain event produced only little runoff
and drainage. For the same reason, drainage (drain 1) at the upstream position was not considered. Flux data for the recession
phase and tensiometer recordings were not included in the parameter estimation, but were utilized for independent model performance evaluation.

Fig. 7. Hydraulic functions (left) and corresponding model simulation results for surface runoff and subsurface drainage (right) based on different parameters a (a–d), hs (e–h), Ks (i–l), n (m–p), and Manning’s nM (q and r).

148

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

F1

The different AMALGAM scenarios described herein were terminated after the algorithm converged on a stable Pareto surface,which occurred within 20,000–50,000 model evaluations, depending on the number of converging runs. Each scenario required 1–
2 weeks of run time in sequential mode using the Matlab R2008a(b) BOXR(a) BOXL

F2
Fig. 8. Uniform-2 scenario (two objective functions): Pareto trade-off planes
between objective functions F1 (runoff) and F2 (total drainage rate) for (a) BoxLand (b) BoxR. Small grey dots – Pareto solutions, Big (red) dots – Pareto optimal
solutions, a cross – compromise solution. (For interpretation of the references tocolour in this ﬁgure legend, the reader is referred to the web version of this article.)-1

Drainage (dm³ min )-1

Runoff (dm³ min )

^
F1

^
F0

data

Pressure head (cm)

(64 bit) Microsoft Windows™ XP Professional (64 bit) modelingenvironment on a Dell™ Precision 390 workstation with a QuadCore IntelÒ Core™2 Extreme processor QX6700 (2.67 GHz) and4 GB of RAM. The run time could be reduced by at least a factorof three by running AMALGAM in parallel mode on all four coresof the workstation. The goodness-of-ﬁt was evaluated by comparing (i) simulated and observed runoff and drainage, and (ii) F^i values of Pareto extreme and compromise solutions in the objectivefunction space as visualized in Pareto plots.
The total number of optimized parameters depended on the
scenario (Table 5). In the VGM model, the relation m = 1  1/n
was assumed, along with an air entry value of 2 cm (Ippischet al., 2006). The residual water content was set to zero, hr = 0. With
these assumptions, six parameters remained to be optimized in theuniform scenario. These include Manning’s roughness coefﬁcient,nM, and ﬁve VGM parameters: the saturated water content, hs (–),
the curve shape parameters n (–) and a (L1), the saturated hydraulic conductivity, Ks (L T1), and the exponent in the hydraulic conductivity mode, s (–). The ‘MIM + Crust’ scenario considers mobile
and immobile soil water regions. Here, a version of the MIM model
was chosen in which the mass transfer between domains is basedon differences in saturation instead of pressure head (Köhne et al.,
2004; Simunek et al., 2003). With the residual water content assumed as zero, the only additional parameters are the mobile
and immobile saturated water contents, hs,m and hs,im, and thetransfer coefﬁcient, x. Thus, the ‘MIM + Crust’ scenario required
14 parameters, i.e., nM and x plus the set of six parameters (hs,m,hs,im, n, a, Ks,s) for both the soil and crust. A large range was assigned to each parameter to permit AMALGAM to search the entire
possible space (Table 5).

^
F2

0.6

0.6

BoxR

BoxL
0.4

0.4

0.2

0.2

0

0
0

50

100

150

2

2

1.5

1.5

100

150

1

1
fit

0.5

0.5

predict

0
0

50

100

150

fit
0

50

predict

100

150

50
100
Time (min)

150

40

40
0

0

-40

-40

-80

-80

-120

50

2.5

2.5

0

0

0

50
100
Time (min)

3 cm

7 cm

150

-120

15 cm

0

23 cm

^
F0

Fig. 9. Uniform-2 scenario (two objective functions): observations and simulations based on parameters at the Pareto extreme solutions for drainage (F^1 ) and surface runoff(F^2 ), and the compromise solution (F^0 ). The objective function comprised only runoff and drainage until 100 min. Later observations and pressure heads are predictions. The
pressure head prediction is shown based on the compromise solution.

149

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

3. Results
3.1. Experimental results
Only the second rainfall event was studied, since the ﬁrst event
did not produce any drainage and pressure head values indicatedvery inhomogeneous moisture distribution prior to the ﬁrst event.
Fig. 5 shows the experimental data for the second rainfall event,including runoff and drainage responses in both replicates. Surface
runoff and drainage began approximately 5 min after the start of
rainfall. Final steady ﬂow rates of about 2 dm3/min were 10 times
larger for surface runoff than for the largest drainage rates, whichwere observed at the downstream drains (Fig. 5). The sum of
cumulative drainage and runoff from BoxR was approximately
equal to the rainfall amount of the second experiment after
100 min. For BoxL, a gap of 57 dm3 (23%) in the water balance suggests a water content increase of 0.038. This hypothesis was qualitatively conﬁrmed by the increase of pressure heads during a
rainfall event (Fig. 6). For BoxL, initial pressure heads ranged between 60 cm and 120 cm (Fig. 6a), suggesting more unsaturated
moisture conditions than in BoxR, where values of 50 to 60 cm
were measured (Fig. 6b). Moreover, in BoxL, the initial response
times for the pressure heads increased with soil depth and showeda slow increase at depths below 7 cm. By contrast, the tensiometers in BoxR showed a rapid succession of their initial response,
and the increase at the 23- and 28-cm depths was faster than at
the 15-cm depth (Fig. 6b), suggesting preferential ﬂow in more
permeable soil regions embedded in less permeable (‘immobile’)
regions. The 15-cm depth was characterized by a low bulk density
and relatively dry moisture content recorded a few hours after theexperiment (Table 1). Apart from preferential ﬂow caused by surface topography or layer interfaces, ﬁlling the box with sieved(<5 cm) soil aggregates might have created inter-aggregate ﬂow
paths. Pressure heads could then show more rapid or slower response, depending on whether the tensiometer is located in or between soil clods. Such immobile soil regions could also explain theincongruity between the apparently closed water balance for BoxRand the observed increase in pressure heads suggesting an increasein water storage. Since immobile regions saturate and drain more
slowly than the remaining soil, they may not have reached equilibrium during the experiment.
After an initial fast increase, the surface runoff rate slowly continued to increase. This gradual increase in the runoff rate implies a
related decrease in the inﬁltration rate, which in turn suggests a

corresponding decrease of surface hydraulic conductivity overtime. Indeed, the erosion of soil particles and the formation of acrust on the soil surface were observed during the ﬁrst rainfall
event, particularly for BoxL. During the experiment this crust
may have further settled and compacted, with a concurrent decrease in the soil surface hydraulic conductivity. For BoxR, a noticeable increase of the runoff rate at approximately 50 min coincides
with the lowest pressure head measurement (at a 15-cm depth)
reaching saturation. While replicate experiments will never show
identical results, the differences between the hydraulic responsesof the two boxes appear to be larger than would be expected froma random variation. Both received identical pre-treatment with 10
rainfall events. So what is the explanation that the two boxes do
not show more similar results? The differences in bulk densities(Table 1) suggest either some spatial variability in the soil reconstruction or some heterogeneity in the soil material.
Limited data were also collected to characterize the recessing
limb of the drainage and runoff hydrographs after the end of rainfall at 100 min (Fig. 5). For BoxR, drainage ceased within 2, 20, and30 min, at drains 1–3, respectively (Fig. 5c), which is the opposite
sequence of the onset of drainage at the beginning (Fig. 5d). For
BoxL, drainage was essentially zero for drain 1, ceased within
2 min for drain 2, and lingered for 30 min for drain 3.
3.2. Model sensitivity
A qualitative sensitivity study was conducted in order to assess
the effect of changing parameter values on the simulation results
for runoff and drainage ﬂuxes. The VGM parameters a, hs, Ks, n, s,
and Manning’s nM of the base case scenario were obtained by manually calibrating the model to observations from BoxL (Fig. 7). The
calibrated parameters are mostly within the range of experimental
values (Table 2). For sensitivity analysis, all parameters were varied by +/50%. As an exception, the lower threshold for n was set
to 25%, since n has to be larger than unity. Using the same relative
variations for all parameters has the advantage of giving a quickﬁrst impression of the sensitivity. However, the results need to
be seen in the light of differing range of parameter values fordescribing different soils. This range is often smallest for hs, larger
for n, a, s, nM, and largest for Ks. Fig. 7 shows the changes in the soil
hydraulic functions and model responses to a change in the input
parameter values. Hydraulic functions are plotted in terms ofh(h) and K(h), since h(h)  K(h) and Se(h)  K(Se) do not respondto a variation of some of the parameters.

Table 6
Optimized parameters for the uniform scenario with two objectives optimization. RMSE – normalized root mean square error, F0 – compromise solution, Fr – Pareto extreme forsimulated runoff, Fd – Pareto extreme solution for the simulation of total drainage.
Pareto resultshsna (1/cm)

Ks (cm/d)

s

nM

F1

BoxL
Fd
Fr
F0

0.451
0.405
0.467

1.80
1.51
1.83

0.003240.00285
0.00322

28.2
90.6
41.3

2.99
1.34
1.74

0.00110.00190.0008

0.222
0.717
0.281

0.542
0.194
0.421

0.764
0.911
0.702

1.66
1.37
1.37

0.00506
0.004890.00489

76.5
63.5
77.6

0.41
0.88
0.91

0.114
0.00010.0001

0.226
0.316
0.240

2.462
0.256
0.281

2.688
0.574
0.521

1.5

0.003

31.7

1

0.007

hs

n

a

Ks

0.483
0.464

1.36
1.56

0.01120.0105

BoxR
Fd
0.467
Fr
0.499
F0
0.500
BoxL – manual calibration (Fig. 7)
0.4

F2

FR

MeasuredabSoil samples with crust.
Soil below crust.

Wind

Permeameter (ﬁeld)

Permeameter (lab)

12.2
12.8

16.3
16.3

146.0 ± 43.6a
9.1 ± 5.5b

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

3.2.1. Uniform-2 scenario
The Uniform-2 scenario considers a homogeneous soil and two
objective functions, i.e., F1 (F2) based on the normalized RMSE between simulated and observed total drain ﬂow (runoff). First, inverse model results are studied in terms of the F1–F2 Pareto
plane (Fig. 8). Among all F1–F2 pairs (grey dots), the Pareto efﬁcientparameter solutions (red solid circles) form Pareto fronts, whichlook different for the two replicates. For BoxL (left panel ofFig. 8), the front is diagonal. For BoxR, the Pareto front shows a distinct corner located at F1 = 0.24 and F2 = 0.27, representing thecompromise solution, F^0 . These comparable values for F1 and F2
suggest a successful normalization. For BoxR, a parameter set canbe identiﬁed that ﬁts both normalized objective functions equally
well. This is not the case for BoxL, where the Pareto front exhibitsa more distinguished trade-off between F1 and F2.
Observations and the corresponding model simulations arecompared in Fig. 9. The trade-off between Pareto extreme solutions
for BoxL is clearly illustrated (Fig. 9, left panel): a reduction indrainage goes along with an increase in surface runoff, and viceversa. For both Pareto extremes, the respective complementary
ﬂuxes were overestimated. Accordingly, the compromise solution,F^0 , slightly overestimates both drainage and surface runoff bysome 5%. This reveals that for scenario 1 and BoxL, no parameterset exists for which the ﬂow model will accurately represent theexperiment. Simulation results could be improved if experimental
data were not considered as error-free, but to some extent uncertain. For example, if a 5% lower effective rainfall rate was assumed

in the model input, to account for lateral spray losses as suggestedby the slight gap in the water balance, model optimization would
certainly better match runoff and drain ﬂow observations for BoxL.
For BoxR (Fig. 9, right panel), the Pareto extremes invoked very
different ﬂow simulations. The model with parameters at Pareto
extreme F^1 incorrectly predicted prolonged drainage after theend of rainfall (where data were not included in the optimization),along with an entirely wrong runoff pattern. The F^0 compromise
solution ﬁtted both runoff and drainage well, and also predicteddrainage recession.
The prediction of pressure heads, which were not included in
the model calibration, is shown in Fig. 9 (bottom panels). The simulations match the initial tensiometer response quite well. However, some observed pressure heads lag considerably behind
others. As discussed previously, some soil regions appear to be lessaccessible to inﬁltrating water, which cannot be represented in themodel that assumes uniform homogeneous conditions. Moreover,

BoxL

BoxR

2

F2

1.5

1

0.5

0

0

0.5

1

1.5

2 0

0.5

F1

1

1.5

2

1.5

2

1.5

2

F1

2

1.5

F3

Decreasing a shifts the point of the steepest slope of h(h) towards more negative (larger absolute) values of h (van Genuchten,1980). Hence, within the pressure head range of the experiment,
decreasing a increases h(h). While changing a does not affect theshape of the K(h) curve, decreasing a indirectly increases K(h) bylimiting it to a more saturated range (Fig. 7a and b). Accordingly,
the soil saturates faster, leading to an earlier start of drainage(Fig. 7c) and runoff (Fig. 7d).
A variation of hs by +/50% results in seemingly very different
hydraulic functions (Fig. 7e and f), which invoke almost the same
ﬂow responses as corresponding changes of a (Fig. 7g and h). The
explanation is that when a or hs are reduced, the slope of h(h) is lesssteep in both cases, so that K(h) stays closer to its saturated endpoint, Ks.
Doubling or bisecting the value of Ks itself (Fig. 7i and j) does notaffect the onset and cessation of ﬂows more than a or hs. However,Ks is the single parameter that controls the partitioning between
runoff and drainage at saturated steady-state conditions (Fig. 7k
and i). When additionally considering its potentially largest rangeof variation, Ks emerged as the most sensitive parameter.
The variation of n shows relatively complex results. As opposed
to a or hs, a decrease in n increases h(h) and decreases K(h). The net
effect on ﬂows is not unidirectional; both a decrease by 25% or an
increase by 50% lead to an earlier onset of runoff and drainage. An
increase of n by 25% (not shown) would show a very similar onsetof runoff and drainage as the standard scenario.
When Manning’s roughness parameter, nM, is decreased so as torepresent a smoother surface, this leaves drainage onset unaffected(Fig. 7q), but causes a more rapid rise of runoff to its maximum value (Fig. 7r). The recessions in both runoff and drainage are accelerated by a decrease in nM (Fig. 7r and q).
Finally, in this particular scenario, the tortuosity parameter, s,
had practically no effect on both K(h) and simulated ﬂows. Nevertheless, s was included in the global optimization, since it cannot
be asserted that s will not show signiﬁcant effects when the valuesof the other VGM parameters differ. In the inverse analysis, a wide
range was chosen for all parameters in order to sample all possible
values. For s this range was set from 3 to +3.

1

0.5
0

0

0.5

1

1.5

2 0

0.5

F1

1

F1

2

1.5

F3

150

1

0.5
0

0

0.5

1

F2

1.5

2 0

0.5

1

F2

Fig. 10. Uniform-3 scenario (three objective functions): Pareto trade-off planes
between three combinations of objectives for BoxL (36,000 evaluations) and BoxR
(50,000). Pareto optimal solutions – red line of narrowly spaced dots; a cross –compromise solution. Objectives F1 and F2 refer to cumulative drainage at the midand down-stream drain, respectively; F3 refers to surface runoff. (For interpretationof the references to colour in this ﬁgure legend, the reader is referred to the web
version of this article.)

151

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

observed maximum values for pressure heads below zero revealthat there was no buildup of water pressure with depth duringconditions of ponding. By contrast, the simulation shows positive
pressures during ponding. Inserting a thin, highly conductive
drainage layer to mimic the geotextile below the soil would remove positive pressures from the simulation (not shown).
Optimized parameter values are shown in Table 6. For comparison, the manually calibrated values (BoxL) and those obtained in
the experimental measurements are displayed also in Table 6.
Most VGM parameters do not deviate much between the F0 compromise solutions and independent measurements. The optimized
Ks values ranged between those measured on furrow and ridge
samples (Table 6). The Manning’s roughness parameter was verysmall, nM < 0.001, suggesting an even surface. The surface of thisartiﬁcially prepared soil was indeed smooth compared to arable
soils in their natural settings. However, the effects of funneled runoff pattern due to the interrill micro-topography may have been
lumped into the Manning’s nM parameter. In the experiment, surface runoff was somewhat accelerated by being conﬁned toslope-parallel mini-furrows, whereas the model assumed sheet
ﬂow covering the entire soil surface. The conﬁnement of surface
ﬂow to a smaller fraction of the soil surface may lead to a smaller
effective value of the Manning’s nM parameter. We tried to mimic

this effect by assuming only half the width of the box. The optimization then produced a nM value increased by some 10%, while hs
doubled, and other VGM parameters changed considerably (notfurther shown).
3.2.2. Uniform-3 scenario
The Uniform-3 scenario was deﬁned by a search of three Pareto
extremes (see Eq. (10)). The bi-criterion F1–F2, F1–F3, and F2–F3 surfaces of the three-dimensional Pareto trade-off space (objective
space) are shown in Fig. 10 for BoxL and BoxR, respectively. Thereis no common minimum for all three objectives for both replicates.
For BoxR, the bi-criterion fronts exhibit a curved shape, with thecompromise solution getting equal shares of normalized RMSE
from all objective functions (Fig. 10).
The observed and simulated cumulative runoff and drainage atthree drains are compared in Fig. 11. Cumulative data were used in
the optimization in an attempt to force it to match the mass balance of the observed volumetric ﬂows. The compromise solution
F^0 for BoxL (Fig. 11, left) fails to match midstream drainage, but ﬁts
downstream drainage and surface runoff. The F^1 Pareto extreme
matched only up- and midstream drainage well. The results demonstrate that for this model scenario there is no F^0 for BoxL that ﬁts
all three objectives simultaneously.

BoxL

BoxR

upstream (no fit)

20

20

10

10

0

0

Drainage (dm³)

0

50

100

150

midstream

20

50

100

150

0

50

100

150

0

50

100

150

0

50

100

150

20

10

10

0

0
0

50

100

150

downstream

20

20

10

10

0

0
0

Runoff (dm³)

0

50

100

150

200

200

100

100

0

0
0

50

100

150

Time (min)

data

F^0

Time (min)

F^1

F^2

F^3

Fig. 11. Uniform-3 scenario (three objective functions): observed cumulative surface runoff and drainage at up-, mid-, and down-stream drains, and model simulation resultsfor different Pareto extreme solutions for drainage at midstream (F^1 ) and downstream (F^2 ) drains, surface runoff (F^3 ), and the compromise solution (F^0 ), for BoxL (left) andBoxR (right).

152

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

BoxL

Table 7
Optimized parameters for the Uniform-3 and MIM + Crust Scenarios with threeobjectives optimization. Fm, Fl, Fr – Pareto extreme solutions for the simulation ofdrainage from the midstream and downstream drains, and for runoff, respectively.

BoxR

2

F2

1.5

Pareto result

1

0.5

0

0

0.5

1

1.5

2 0

0.5

1

1.5

2

F1

F1

a

2

F3

1.5

0.5

0

0.5

1

1.5

2 0

0.5

F1

1

1.5

2

F1

2

F3

1.5

1

0.5
0

0

0.5

1

1.5

2

0

0.5

F2

1

1.5

Uniform-3 (BoxR)

F1

F2

F3

FR

F1

F2

F3

FR

Fm
Fl
Fr
F0

0.21
1.65
0.27
0.27

2.71
0.33
1.10
1.10

3.65
0.66
0.18
0.18

6.57
2.64
1.55
1.55

0.36
0.70
0.88
0.37

0.73
0.35
1.12
0.64

1.07
1.29
0.36
0.63

2.16
2.34
2.36
1.63

Fm
Fl
Fr
F0

MIM + Crust (BoxL)a
0.27
0.99
3.65
2.01
0.46
1.31
1.84
0.66
0.37
0.77
0.72
0.90

4.90
3.78
2.87
2.40

MIM + Crust (BoxR)a
0.37
0.65
1.01
0.76
0.35
1.52
0.90
1.07
0.39
0.46
0.71
0.59

2.02
2.62
2.36
1.76

With constraint Ks(crust) < Ks(soil).

In contrast, F^0 better represents the different BoxR data (Fig. 11,
right). This is because the increase of drainage towards the downstream end is not as extreme as in the BoxL replicate. The Pareto
extreme for midstream drainage, F^1 , agrees with the observations
nearly as well. The minimal upstream drainage is difﬁcult to explain and points to heterogeneity in the material properties atthe surface or subsurface. Upstream drainage was not included in
the ﬁtting procedure, but is shown for completeness.
The inverse simulation reveals that even a relatively good
match of surface runoff and lumped drainage (Uniform-2 scenario)
does not guarantee a correct model description of the spatial pattern of subsurface ﬂow.

10

Uniform-3 (BoxL)

2

F2

Fig. 12. MIM + Crust Scenario: Pareto trade-off planes between three combinationsof objectives for BoxL (21,200 evaluations) and BoxR (50,040). Front of red dots –Pareto optimal solutions; a cross – compromise solutions with the constraint
Ks(crust) < Ks(soil). Objectives F1 and F2 refer to cumulative drainage at the mid- anddownstream drains; F3 refers to surface runoff.

3.2.3. MIM + Crust scenario
The ‘MIM + Crust’ scenario considers mobile and immobile regions in a soil with a surface crust. This scenario, motivated by
the aggregated structure of the soil and observed crust formation,
again includes the three objectives for overland ﬂow and mid- anddownstream drainage. The bi-criterion fronts of the 3D Pareto
trade-off space (Fig. 12) look very similar to those obtained for
the Uniform-3 scenario (Fig. 11). We note that a constraint was imposed to keep the saturated hydraulic conductivity of the crust below that of the soil (Ks(crust) < Ks(soil)), which was obtained for
only 1053 runs out of 20,000 for BoxL and 14,030 runs out of
50,040 for BoxR. This constraint ﬁltered out some slightly lower
Pareto extremes for the MIM model, effectively without the surfaced crust. Overall, the Pareto extremes and the compromise solution did not improve the results compared to the Uniform-3

2

2

1.8

1.8

1.6

1.6

1.4

1.4

1.2

1.2
1.0 F3
0.8

F2

1.0

0.8
0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.5

1

F1

1.5

20

0.5

1

1.5

2

0

F1

Fig. 13. MIM + Crust scenario: the objective function (F1–F3) space; the color code indicates F3. Left: BoxL, Right: BoxR. The black solid line shows the F1–F2 front of the Pareto
optimal solutions, the dark solid line depicts the F1–F3 Pareto front, the hatched purple line indicates the F2–F3 Pareto front, and the stars represent the Pareto compromise
solution, F^0 .

153

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

BoxR

BoxL
20

upstream (no fit)

20

10

10

0

0

Drainage (dm³/min)

0

50

100

150

50

100

150

0

50

100

150

0

50

100

150

0

50

100

150

midstream

20

20

10

10

0

0
0

50

100

150

downstream

20

20

10

10

0

0
0

50

100

150

200

Runoff (dm³)

0

200

100

100

0

0
0

50

100

150

Time (min)

data

F^0

Time (min)

F^1

F^2

F^3

Fig. 14. MIM + Crust scenario (three objective functions): observed cumulative surface runoff and drainage at up-, mid-, and down-stream drains, and model simulation
results for different Pareto extremes for drainage at midstream (F^1 ) and downstream drain (F^2 ), surface runoff (F^3 ), and the compromise solution (F^0 ), for BoxL (left) and BoxR
(right).

example (Table 7). A more compact visualization of the objective
space is shown in Fig. 13 for the F1–F2 Pareto trade-off plane withthe color coded objective F3. Visually, cumulative midstream drainage and runoff were ﬁtted slightly better for BoxR (Fig. 14) than in
the Uniform-3 scenario (Fig. 12), but the differences were rather
small.
Since no compromise solution could be found that would matchall data equally well, the optimized model parameters may not bemeaningful, and were thus omitted. Hence, the added model complexity did not result in an improved match. If a local optimization
approach had been chosen, such as Levenberg Marquardt, it would
not have been possible to rule out a failure in the inverse procedurebehind this lack of agreement. Using multi-objective optimization
proved that there was a discrepancy between model structure andexperimental data.
Fig. 15 compares simulated surface runoff depths, ho, for the different scenarios at nine observation points along the surface. Maximum values for ho are less than 2.5 mm. For a given overlanddischarge rate, Qo, the depth ho increases with nM (Eqs. (2) and(3)). Hence, the differences between simulated ho originate mostly
from differences in the ﬁtted nM values, since the optimized values
for Qo did not differ that much. For each scenario, ho increases

downhill. The onset of runoff, in response to rainfall, occurs almost
simultaneously at different positions along the slope. Some scenarios differ with regard to the recessing limbs after the stop of rainfall. Only those scenarios with relatively high ho show fasterdissipation in the upslope part.
4. Discussion
In our case study, the task was to match drainage ﬂows at midand downstream hillslope positions and runoff. A ‘validation’ test
was added to predict (without ﬁtting) the pressure heads, therecession part of drainage ﬂows, and the drainage hydrographs at
the upstream position. However, an increase in drainage in thedownslope direction could not be reproduced by any of the considered model scenarios. This could have been caused by soil layerswith different soil bulk densities that may have generated saturated lateral ﬂow on layer interfaces. Such a situation, which could
be captured by explicit consideration of layers in the model, was
not simulated because of the considerable increase in the numberof calibration parameters.
The use of a global optimization technique rules out the ﬁtting
procedure itself as a possible pitfall. Since model assumptions were154

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

1

0

2

3

4

1

5

6

2
Box-L

0.2

8

3
Box-R

Uniform-2

9

7

9

4

5
m

0.2

5
0.1

0.1

Overland flow depth (cm)

1
0

0
0

50

0

50

100

0

50

100

0

50
Time (h)100100

Uniform-3
0.2

0.2

0.1

0.1

0

0
0

50

100

MIM+Crust
0.2

0.2

0.1

0.1

0

0
0

50
Time (h)100Fig. 15. Overland ﬂow depths simulated at different observation nodes (numbers in the top ﬁgure) along the upper boundary of HYDRUS-2D for different compromise
scenarios for BoxL (left) and BoxR (right).

explicitly known, the difﬁculty of simultaneously matching all different data suggested the presence of some sort of heterogeneity ofsoil hydraulic or surface properties that developed despite of thehighly controlled environment. The aggregated and layered soil
showed erosion and crust formation at the surface. The consideration of a crust, although observed, was not highly relevant in
the ﬂow simulation. Had the development of the crust been more
pronounced upstream than downstream, this could have explained
the increase in drainage ﬂux towards the lower half of the box. The
spatial and temporal development of a surface crust and its effecton inﬁltration was indeed measured and simulated recently(Augeard et al., 2007), although on a much smaller soil system. Thiskind of analysis was beyond the scope of this study. However, if
available, such information could be fed into the model.
The AMALGAM multi-objective global search method was useful for proving that there was no common optimal parameter set
that described all observations at different spatial locations equally
well. The shape of the Pareto trade-off space illustrates clearly the
incommensurateness of the objective functions. The analysis thus
revealed the divergence between model structure or hypothesis,
and the experimental observations. In our case, the assumptionof having uniform and isotropic soil, or aggregated soil with crust,does not accurately represent the spatial nature of the ﬂow process. Further information with regard to parameter uncertainty
could be gained by studying the parameter ranges for Pareto solutions (e.g., Wöhling et al., 2008).
In vadose zone modeling, there is a large body of literature on
including different types of observations into the model calibrationand its effects on parameter identiﬁcation (e.g., Friedel, 2005;
Vrugt et al., 2008; Wöhling and Vrugt, in press; Zhang et al.,2003). However, to the best of our knowledge, similar studies
involving coupled surface–subsurface ﬂow are so far missing.

Similar results, which we obtained for different scenarios andparameter sets (Uniform-3 and MIM + Crust), point to the wellknown equiﬁnality problem, where several underlying model
structures may yield similar simulation results (Beven, 2006). Reasons are, among others, the incomplete knowledge of the natural
system and its inaccurate representation in the model. To accountfor this type of uncertainty, (Klaus and Zehe, 2010) suggested to include multiple structural set-ups of a physically based model, e.g.,
by employing the Generalized Likelihood Uncertainty Estimation
method (Klaus and Zehe, 2010).

5. Conclusions
A coupled surface–subsurface ﬂow model based on HYDRUS-2D(Šimůnek et al., 1999) was presented and linked to the multiobjective parameter search method AMALGAM (Vrugt andRobinson, 2007). Further, results from a highly controlled benchscale experiment of surface and drainage ﬂows in response to
rainfall was presented. The data were used to evaluate the model
performance and to test various hypotheses regarding different
ﬂow processes. Multi-objective parameter estimation with different normalized RMSE objective functions, Pareto trade-off space
analysis, and visual comparison between model results and experimental observations were jointly used in this procedure.
The following conclusions can be drawn. Coupled surface–subsurface ﬂow is a complex problem which requires model
calibration to match several types of observations related to runoff
and soil water dynamics. While there is a large body of literature
related to inverse estimation for vadose zone models, more researchseems necessary to extend these techniques also to coupled models. To this end, the AMALGAM technique, along with Pareto plot

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

analysis, was a useful tool for global optimization and for studying
trade-off patterns between multiple objectives. In our case study, it
showed that a common optimal parameter set to describe all typesof observations equally well did not exist.
The result could be used to rethink the design of experimental
set-ups and model structures. Moreover, the analysis revealed thata good match between simulated and observed surface runoff andtotal drainage does not guarantee an accurate representation of theﬂow process.
Overall, the suggested model appears useful for studying thehydraulics of overland ﬂow and subsurface ﬂow in sloping soils,
potentially up to a hillslope scale. The multi-objective model calibration adds a powerful method for understanding the role of different hydraulic processes in such systems.
Acknowledgements
The helpful comments and suggestions of two anonymous
reviewers are highly appreciated. Funding of this work was provided by the German Research Foundation (DFG) under contract
DFG KO2899-1. The authors wish to thank the technical assistanceof B. Renaud and L. Prud’homme for rainfall simulations. They also
thank F. Darboux and C. Legout for their scientiﬁc help in the rainfall experiments.
References
Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, M.J., Ma, L., 2000. Root Zone Water
Quality Model – Modeling Management Effects on Water Quality and Crop
Production. Water Resources Publications, LLC, Highlands Ranch, CO, p. 372.
Augeard, B., Assouline, S., Fonty, A., Kao, C., Vauclin, M., 2007. Estimating hydraulicproperties of rainfall-induced soil surface seals from inﬁltration experiments
and X-ray bulk density measurements. J. Hydrol. 341 (1–2), 12–26.
Beven, K., 2006. A manifesto for the equiﬁnality thesis. J. Hydrol. 320 (1–2), 18–36.
Corradini, C., Morbidelli, R., Melone, F., 1998. On the interaction between inﬁltrationand Hortonian runoff. J. Hydrol. 204 (1–4), 52–67.
Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective
genetic algorithm: NSGA-II. Trans Evol Comput 6, 182–197.
Deng, Z.-Q., de Lima, J.L.M.P., Singh, V.P., 2005. Transport rate-based model foroverland ﬂow and solute transport: parameter estimation and processsimulation. J. Hydrol. 315 (1–4), 220–235.
Friedel, M.J., 2005. Coupled inverse modeling of vadose zone water, heat, and solutetransport: calibration constraints, parameter nonuniqueness, and predictive
uncertainty. J. Hydrol. 312 (1–4), 148–175.
Furman, A., 2008. Modeling coupled surface–subsurface ﬂow processes: a review.
Vadose Zone Journal 7 (2), 741–756.
Gandolﬁ, C., Savi, F., 2000. A mathematical model for the coupled simulation of
surface runoff and inﬁltration. J. Agric. Eng. Res. 75 (1), 49–55.
Govindaraju, R.S., Kavvas, M.L., 1991. Dynamics of moving boundary overland ﬂows
over inﬁltrating surfaces at hillslopes. Water Resour. Res. 27 (8), 1885–1898.
Govindaraju, R.S., Kavvas, M.L., Jones, S.E., 1990. Approximate analytical solutions
for overland ﬂows. Water Resour. Res. 26 (12), 2903–2912.
Gupta, H.V., Sorooshian, S., Yapo, P.O., 1998. Toward improved calibration ofhydrologic models: Multiple and noncommensurable measures of information.
Water Resour. Res. 34 (4), 751–763.
Herbst, M., Diekkrüger, B., Vanderborght, J., 2006. Numerical experiments on thesensitivity of runoff generation to the spatial variation of soil hydraulicproperties. J. Hydrol. 326 (1–4), 43–58.
Ippisch, O., Vogel, H.J., Bastian, P., 2006. Validity limits for the van Genuchten–
Mualem model and implications for parameter estimation and numericalsimulation. Adv. Water Resour. 29 (12), 1780–1789.
Klaus, J., Zehe, E., 2010. Modelling rapid ﬂow response of a tile-drained ﬁeld siteusing a 2D physically based model: assessment of ‘equiﬁnal’ model setups.
Hydrol. Process. 24 (12), 1595–1609.
Kodesová, R., Vignozzi, N., Rohosková, M., Hájková, T., Kocárek, M., Pagliai, M.,
Kozák, J., Simunek, J., 2009. Impact of varying soil structure on transport
processes in different diagnostic horizons of three soil types. J. Contam. Hydrol.
104 (1–4), 107–125. doi:10.1016/j.jconhyd.2008.10.008 (Special Issue ‘‘Flow
Domains’’).
Köhne, J.M., Köhne, S., Mohanty, B.P., Šimůnek, J., 2004. Inverse mobile–immobile
modeling of transport during transient ﬂow. Vadose Zone J. 3 (4), 1309–1321.
Kulli, B., Stamm, C., Papritz, A., Fluhler, H., 2003. Discrimination of ﬂow regions on
the basis of stained inﬁltration patterns in soil proﬁles. Vadose Zone J. 2 (3),338–348.
Le Bissonnais, Y., Benkhadra, H., Chaplot, V., Fox, D., King, D., Daroussin, J., 1998.
Crusting, runoff and sheet erosion on silty loamy soils at various scales andupscaling from m2 to small catchments. Soil Tillage Res. 46 (1–2), 69–80.

155

Leguédois, S., Bissonnais, Y.L., 2004. Size fractions resulting from an aggregatestability test, interrill detachment and transport. Earth Surf. Proc. Land. 29 (9),1117–1129.
Leguédois, S., Benoit, P., Pot, V., Le Bissonnais, Y., 2006. A 10-m2 Laboratory
Lysimeter to Assess Lateral and Vertical Pesticide Transfers in a Plowed Horizon,
Geophysical Research Abstracts. European Geosciences Union, Vienna, p. 1.
Léonard, J., Perrier, E., de Marsily, G., 2001. A model for simulating the inﬂuence of a
spatial distribution of large circular macropores on surface runoff. Water
Resour. Res. 37 (12), 3217–3225.
Li, Z., Zhang, J., 2001. Calculation of ﬁeld manning’s roughness coefﬁcient. Agric.
Water Manage. 49 (2), 153–161.
Madsen, H., 2000. Automatic calibration of a conceptual rainfall-runoff model using
multiple objectives. J. Hydrol. 235 (3–4), 276–288.
Madsen, H., 2003. Parameter estimation in distributed hydrological catchment
modelling using automatic calibration with multiple objectives. Adv. Water
Resour. 26 (2), 205–216.
Mallants, D., Mohanty, B.P., Vervoort, A., Feyen, J., 1997. Spatial analysis of saturated
hydraulic conductivity in a soil with macropores. Soil Technol. 10 (2), 115–131.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math. 11 (2), 431–441.
Mertens, J., Madsen, H., Feyen, L., Jacques, D., Feyen, J., 2004. Including priorinformation in the estimation of effective soil parameters in unsaturated zone
modelling. J. Hydrol. 294 (4), 251–269.
Merz, B., Plate, E.J., 1997. An analysis of the effects of spatial variability of soil andsoil moisture on runoff. Water Resour. Res. 33 (12), 2909–2922.
Mohamoud, Y.M., 1992. Evaluating Manning’s roughness coefﬁcients for tilled soils.
J. Hydrol. 135 (1–4), 143–156.
Morris, E.M., Woolhiser, D.A., 1980. Unsteady one-dimensional ﬂow over a plane:
partial equilibrium and recession hydrographs. Water Resour. Res. 16 (2), 355–
360.
Motha, J.A., Wigham, J.M., 1995. Modelling overland ﬂow with seepage. J. Hydrol.
169 (1–4), 265–280.
Nahar, N., Govindaraju, R.S., Corradini, C., Morbidelli, R., 2004. Role of run-on for
describing ﬁeld-scale inﬁltration and overland ﬂow over spatially variable soils.
J. Hydrol. 286 (1–4), 36–51.
Panday, S., Huyakorn, P.S., 2004. A fully coupled physically-based spatiallydistributed model for evaluating surface/subsurface ﬂow. Adv. Water Resour.
27 (4), 361–382.
Schwefel, H.P., 1993. Evolution and Optimum Seeking: The Sixth Generation. John
Wiley & Sons, Inc., p. 456.
Sepaskhah, A.R., Bondar, H., 2002. SW – soil and water: estimation of manningroughness coefﬁcient for bare and vegetated furrow irrigation. Biosyst. Eng. 82(3), 351–357.
Šimůnek, J., 2003. HYDRUS-2D Code Modiﬁcation: Modeling Overland Flow.
<http://www.pc-progress.com/en/Default.aspx?h3d-overland>.
Šimůnek, J., Sejna, M., van Genuchten, M.T., 1999. The HYDRUS-2D Software
Package for Simulating Two-dimensional Movement of Water, Heat, andMultiple Solutes in Variably-saturated Media, Version 2.0. International
Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado.
Simunek, J., Jarvis, N.J., van Genuchten, M.T., Gärdenäs, A., 2003. Review andcomparison of models for describing non-equilibrium and preferential ﬂow andtransport in the vadose zone. J. Hydrol. 272 (1–4), 14–35.
Simunek, J., van Genuchten, M.T., Sejna, M., 2008. Development and applications of
the HYDRUS and STANMOD software packages and related codes. Vadose Zone
J. 7 (2), 587–600.
Singh, V., Bhallamudi, S.M., 1998. Conjunctive surface–subsurface modeling ofoverland ﬂow. Adv. Water Resour. 21 (7), 567–579.
Singh, V.P., Jain, S.K., Sherif, M.M., 2005. Errors of kinematic wave and diffusion
wave approximations for time-independent ﬂows with inﬁltration andmomentum exchange included. Hydrol. Process. 19 (9), 1771–1790.
Tamari, S., Bruckler, L., Halbertsma, J., Chadoeuf, J., 1993. A simple method for
determining soil hydraulic properties in the laboratory. Soil Sci. Soc. Am. J. 57(3), 642–651.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898.
Vogel, H.J., Cousin, I., Ippisch, O., Bastian, P., 2005. The dominant role of structure forsolute transport in soil: experimental evidence and modelling of structure andtransport in a ﬁeld experiment. Hydrol. Earth Syst. Sci. Discuss. 2 (5), 2153–
2181.
Vrugt, J.A., Robinson, B.A., 2007. Improved evolutionary optimization fromgenetically adaptive multimethod search. Proc. Natl. Acad. Sci. 104 (3), 708–
711.
Vrugt, J.A., Gupta, H.V., Bastidas, L.A., Bouten, W., Sorooshian, S., 2003. Effective and
efﬁcient algorithm for multiobjective optimization of hydrologic models. Water
Resour. Res. 39 (8), 1214.
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V.V., 2008. Inverse
modeling of subsurface ﬂow and transport properties: a review with newdevelopments. Vadose Zone J. 7 (2), 843–864.
Wallach, R., Grigorin, G., Rivlin, J., 1997. The errors in surface runoff prediction by
neglecting the relationship between inﬁltration rate and overland ﬂow depth. J.
Hydrol. 200 (1–4), 243–259.
Wallach, R., Grigorin, G., Rivlin, J., 2001. A comprehensive mathematical model fortransport of soil-dissolved chemicals by overland ﬂow. J. Hydrol. 247 (1–2), 85–99.
Wind, G.P., 1966. Capillary conductivity data estimated by a simple method. In:
Rijtema, P.E., Wassink, H. (Eds.), Water in the Unsaturated Zone. IASH and
UNESCO, Gentbrugge and Paris, pp. 181–191.

156

J.M. Köhne et al. / Journal of Hydrology 403 (2011) 141–156

Wöhling, T., Mailhol, J.C., 2007. Physically based coupled model for simulating 1D
surface-2D subsurface ﬂow and plant water uptake in irrigation furrows. II:model test and evaluation. Journal of Irrigation and Drainage Engineering 133(6), 548–558.
Wöhling, T., Schmitz, G.H., 2007. Physically based coupled model for simulating
1D surface-2D subsurface ﬂow and plant water uptake in irrigation
furrows. I: model development. J. Irrigation Drainage Eng. 133 (6), 538–
547.
Wöhling, T., Vrugt, J.A., 2008. Combining multiobjective optimization and Bayesianmodel averaging to calibrate forecast ensembles of soil hydraulic models. Water
Resour. Res. 44 (12), W12432.

Wöhling, T., Vrugt, J., in press. Multi-Response Multi-layer Vadose Zone Model
Calibration using Markov Chain Monte Carlo Simulation and Field Water
Retention Data. Water Resources Research, 47, WR009265.
Wöhling, T., Vrugt, J.A., Barkle, G.F., 2008. Comparison of three multiobjective
optimization algorithms for inverse modeling of vadose zone hydraulicproperties. Soil Sci. Soc. Am. J. 72 (2), 305–319.
Woolhiser, D.A., Smith, R.E., Goodrich, D.C., 1990. KINEROS, A Kinematic Runoff and
Erosion Model: Documentation and User Manual.
Zhang, Z.F., Ward, A.L., Gee, G.W., 2003. Estimating soil hydraulic parameters of aﬁeld drainage experiment using inverse techniques. Vadose Zone J. 2 (2), 201–211.