Journal of Hydrology 254 (2001) 67±81

www.elsevier.com/locate/jhydrol

Simulating unsaturated ¯ow and transport in a macroporous soil totile drains subject to an entrance head: model development andpreliminary evaluation
A. Kohler a, K.C. Abbaspour b,*, M. Fritsch c, M. Th. van Genuchten d, R. Schulin eaSwiss Federal Institute of Technology, Department for Land and Water Management, ETH HoÈnggerberg, 8093 Zurich, Switzerlandb
Swiss Federal Institute for Environmental Science and Technology (EAWAG), 8600 DuÈbendorf, Switzerlandc
Bob Partners GmbH, Blaufahnenstrasse 14, CH-8093 Zurich, SwitzerlanddU.S. Salinity Laboratory, USDA, ARS, 450 West Big Spring Road, Riverside, CA, 92507 USAe
Swiss Federal Institute of Technology, Department of Soil Protection, Grabenstrasse 3, 8952 Schlieren, Switzerland
Received 6 October 2000; revised 8 May 2001; accepted 3 July 2001

Abstract
Accurate prediction of water ¯ow and chemical transport in agricultural soil pro®les requires the use of a simulation modelthat considers the most important physical, hydrological and chemical processes. Two important ¯ow-related processes in tiledrained ®eld systems are macropore ¯ow and water discharge from the tile drains. To better account for these two processes, we
extended an existing two-dimensional model (SWMS_2D) by adding a macropore ¯ow component as well as a Hooghoudt typeboundary condition that considers the presence of an entrance head at the tile drain. The macropore component is necessary to
account for water and solutes short-circuiting the soil matrix, while the drainage entrance head is needed to account for thecontraction of streamlines around the drains, a feature that causes delayed discharge. The applicability of the new model to aland®ll problem was examined. The simulation results, which included water ¯ow and solute transport, compared well withother models. q 2001 Elsevier Science B.V. All rights reserved.
Keywords: Unsaturated soil; Preferential ¯ow and transport; Macroporous soil; Tile drain; Entrance head; Modeling

1. Introduction
The gradual contamination of surface and subsurface water resources has important social, environmental, economic, and political repercussions.
Although public awareness about environmental
issues has increased sharply over the past two
decades, many situations and practices keep contributing to contamination of our soil and water* Corresponding author. Tel.: 141-1-823-5359; fax: 141-1-8235375.
E-mail address: abbaspour@eawag.ch (K.C. Abbaspour).

resources. These include agricultural management
practices such as fertilization and pesticide applications, the use of septic ®elds, and the construction ofland®lls, radioactive waste disposal sites, wastewaterand sewage lagoons, and mine tailing embankments.
Reliable analyses of these contamination problems
require models that account for the most important
physical, hydrological, and chemical processes. The
focus of this study is on two important ¯ow processes
in tile-drained agricultural ®elds: preferential ¯ow anddischarge through tile drains.
The Swiss Federal Of®ce of Environment, Forest
and Landscape (FOEFL, 1997) estimated that the0022-1694/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved.
PII: S 0022-169 4(01)00499-1

68

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

Nomenclature
a
A
B
c
cs
C
Cd
d
Dij
dhi
dzi
E
Ea
Epot
h
hT
he
K(h)
KijA
Ks
Ks,macKdrainLcmnn*q
Qmac,i
Qe
s
S
Se
t
xj
z
a
b
e
gw
gs
ms
mw
u
ur
us
u mac
r

Distance between the centre of the soil aggregates to the centre of the pores (l)Coef®cient of Hooghoudt equation (t 21)
Coef®cient of Hooghoudt equation (l 21 t 21)
Solution concentration (ml 23)
Concentration of the sink term (ml 23)
Coef®cient of Hooghoudt equation
Drain coef®cient
Number of days after a 5 mm rainfall event
Dispersion coef®cient tensor (l 2 t 21)
Difference between the water potentials of micro- and macro-node i (l)Vertical distance between macro-nodes (i 2 1) and i (l)Effective porosity of the macro-domain
Actual evapotranspiration (l t 21)
Potential evapotranspiration (l t 21)
Pressure head (l)Elevation of the groundwater table above the tile depth, midway between the drains (l)Entrance head (l)Unsaturated hydraulic conductivity (l t 21)
Components of a dimensionless anisotropy tensor K A
Saturated hydraulic conductivity (l t 21)
Saturated hydraulic conductivity of the macropore (l t 21)
Adjusted conductivity of the elements surrounding a drain (l t 21)
Characteristic length representing average spacing pf macropores (l)1 2 1/nvan Genuchten shape parameter
Empirical exponent accounting for the pore size
Water ¯ux (l t 21)
Water ¯ow from macro-node i down to macro-node i 1 1 (t 21)
Exchange ¯ow between a macro-node and a corresponding micro-node (t 21)
Absorbed concentration
Sink term (t 21)
Effective water saturation
Time (t)Spatial coordinates (l)Vertical spatial coordinate (l)van Genuchten shape parameter (l 21)
Geometric factor
Empirical factor
Zero-order rate constants for the liquid phase (ml 23 t 21)
Zero-order rate constants for the solid phase (ml 23 t 21)
First-order rate constants for solutes in solid phases (t 21)
First-order rate constant for solutes in the liquid phases (t 21)
Volumetric water content (l 3 l 23)
Residual water content (l 3 l 23)
Saturated water contents (l 3 l 23)
Water content of macropores (l 3 l 23)
Soil bulk density (ml 23)

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

69Fig. 1. Illustration of the two ¯ow domains used in the extended M-2D model.

annual total nitrogen loss in catchments of the Rhine
to be as high as 50 kg ha 21. Environmentalists often
identify drained lands as a major source of surface
water quality problems (Evans et al., 1989). These
concerns are supported by a recent study by Stamm(1997) who found that solute leached rapidly througha tile-drained ®eld soil and ended up quickly in acreek that captured the drainage water. Even suchsorbing constituents as phosphorus and ammonium
appeared quickly in the drain water after application
to the ®eld. This rapid transport to the drains was
attributed to the presence of preferential ¯ow paths.
Similar observations of preferential ¯ow and transport
in tile-drained ®elds were also made by Richard andSteenhuis (1988); Klavidko et al. (1991); Mohanty etal. (1997, 1998); and Lennartz et al. (1999), among
many others.
In order to more effectively use simulation models
for improving farming operations and fertilization
practices, minimizing tile-drain losses of nutrientsand pesticides, and optimising drainage design andwater table management, it is essential that thesemodels account for the most important ¯ow and transport processes. Van Genuchten and Wierenga (1976)
proposed a two-domain model to describe solute
displacement in a structured soil under steady-state¯ow condition. They partitioned the water phase
into a mobile (macropore) region subject to convec-

tive-dispersive transport, and an immobile, quasistagnant liquid (micropore) region that was coupled
to the mobile region by diffusive solute exchange
following ®rst-order kinetics. Gerke and van Genuchten (1993) later extended this approach to transientwater ¯ow by applying Richards equations to bothregions and coupling them in analogy to solute transport by using a ®rst-order rate equation. Jarvis et al.
(1991) developed a similar model, MACRO, todescribe water ¯ow and transport in macroporous
soils. Contrary to Gerke and van Genuchten (1993),
they considered water movement in the macropore
region to be a non-capillary laminar ¯ow process
driven by gravity only, and assumed solute transport
in the macropores to be a purely convective process.
Van Genuchten and Sudicky (1999) reviewed theseand other conceptualisations of macropore ¯ow insoils and unsaturated fractured rocks.
Jarvis (1991) included tile drainage in their onedimensional MACRO model by means of a special
boundary condition. Water ¯ow to the drains was
described using the seepage potential theory ofYoungs (1980) in which for each layer of the onedimensional column a lateral loss to the drain was
estimated. Loss of solute to the drains was calculated
assuming complete mixing within a ¯ow domain in
the horizontal dimension.
The objectives of this study were to add a macropore

70

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

Fig. 2. Steps in numerical procedure to solve the water ¯ow and solute transport equations. In step 1, the macropore ¯ow and the interaction
between the ¯ow regions are calculated. In step 2, the calculated ¯ow is used as a boundary condition for the next time step of the microdomain.

¯ow component and a tile drain boundary condition
to the SWMS_2D two-dimensional variably-saturated ¯ow/transport model of Simunek et al.
(1994), in order to describe water ¯ow and solutetransport in a tile-drained soil exhibiting preferential¯ow. In this paper we introduce the new model,
henceforth referred to as M-2D, discuss an example
comparing M-2D simulation results with thoseobtained with SWMS_2D and MACRO, and present
calibration of the model for the simulation of hourly
discharge through a land®ll in Switzerland. The calibrated model was then used to simulate electricalconductivity of the discharge through the land®ll.
In a second paper (Abbaspour et al., 2001) we used
the new model to analyse an experimental data set
collected at a tile-drained agricultural ®eld north ofZuÈrich, Switzerland.

2. Description of M-2D
2.1. General description
Fig. 1 depicts the two domains of the M-2D model.
The micropore domain (further referred to as the
micro-domain) represents the soil matrix in whichwater ¯ow is governed by the Richards equation andsolute transport by Fickian-based convection-dispersion equation, respectively. The macropore domain(macro-domain) in Fig. 1 represents surfaceconnected cracks and macropores in which water¯ow within the macropores occurs one-dimensionallyas non-capillary, laminar ¯ow and in which solutetransport is purely convective. Mohanty et al. (1998)
from their comprehensive study of the preferential

transport of nitrate to a tile drain in an intermittent¯ood-irrigated ®eld in New Mexico concluded that
preferential ¯ow intercepted by a tile drain was generated in close proximity of the drain and was essentially oriented vertically. A dye in®ltration test at ourstudy site also showed that in a ¯at drained ®eld and ina ¯at land®ll the mechanism of preferential ¯ow andtransport was directed vertically downward. Hence,
consideration of a two-dimensional ¯ow within
macropores, while adding considerable computational
effort in the program, did not seem to be warranted.
As depicted in Fig. 1, the macro-domain has thesame geometry as the micro-domain. A macropore isassumed to consist of a sequence of macropore nodes(macro-nodes) that extend vertically from a surface
micropore node (micro-node) down to a desired
depth. Each macro-node coincides with a certain
micro-node, but not vice-versa (i.e., not all microporesare close to or connected with a macropore). The widthof a macro-node is assumed to be equal to the width ofthe corresponding micro-node at the soil surface. Afterthe ®nite element grid of the system is implemented,
the user determines if a certain location of the micropore domain should correspond to a macropore
domain, for example near a drain.
Following Gerke and van Genuchten (1993),
exchange of water between a micro-node and corresponding macro-node is calculated based on the pressure head difference between the two nodes. We assume
that the macro-node has a pressure head equal to theheight of the water column in the macropore abovethat node. The surface nodes of macropores can also
store surface water. If the micro-domain is not able toabsorb all of the applied water at the surface, as indicated by a speci®ed parameter hCritS, excess water is

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

71

directed to the surface macro-nodes from which water
may enter a macropore.
Within the soil pro®le we have the same ¯owexchange mechanism as for the surface nodes. If thepressure head at a certain micro-node reaches acritical value of close to saturation (speci®ed by theparameter hCritS), and a corresponding macro-node
exists that is not already saturated, water will ¯owfrom the micro-node to the macro-node according toan interaction term to be discussed later.
The numerical procedure for successive calculationof ¯ow and transport between the micro- and macrodomains is depicted in Fig. 2. Assuming that the
micro-domain was solved for time Tk, the model
calculates ¯ow and transport in the macro-domain as
follows. During step 1, pressure heads of the microdomain are compared to those of the macro-domain
and the interaction term for ¯ow between the two ¯owdomains is calculated. Based upon the actual watercontents in the macropores and the interaction term,water ¯ow and solute transport for time step Tk21 to Tk
are calculated for each macropore. During step 2, ¯ow
and transport in the micro-domain are solved for the
next time step from Tk to Tk11. During this step theinteraction term between the two ¯ow domains, which
was calculated in step 1, is treated as a boundary
condition by means of prescribed in- or out¯ow
rates from the corresponding micro-nodes.

water retention function:

2.2. Variably-saturated ¯owwhere Ks,mac is the saturated hydraulic conductivity ofthe macropore, u mac is the water content, E is themacroporosity of the soil (i.e., the total volume ofmacropores per bulk soil volume), and n * is an empirical exponent re¯ecting the macropore size distribution. More details about the formulation of Eq. (5),or how the parameters can be measured, can be foundin the manual of model MACRO (Jarvis, 1991), or in
Beven and Germann (1981).
Water ¯ow in the macro-domain is calculated bymeans of the following steps: (i) determination of thepressure head at each macro-node, (ii) calculation of
the interaction term between micro- and macro-nodesand (iii) redistribution of water within the macropores.
These steps are repeated several times within eachtime step of the micropore model. The time step of
the macropore sub-model is prescribed by the user andkept ®x during the model run.

2.2.1. Micropore domain
In SWMS_2D (Simunek et al., 1994), water ¯ow ina two-dimensional, isothermal, variably-saturated,
rigid porous medium is described with the following
form of the Richards equation:
"
!#
2u
2
A 2h
A

2S
k h Kij
1 Kij
2t
2x i
2x j

1

where u is the volumetric water content, h is thepressure head, S is a sink term, xj (j  1,2) are thespatial coordinates, t is time, KijA are components ofa dimensionless anisotropy tensor K A, and K(h) is theunsaturated hydraulic conductivity function.
Soil hydraulic properties are described using van
Genuchten's model (van Genuchten, 1980) for theSe h u h 2 ur1
h,0
us 2 ur
1 1 uahun mu h  uS

h$0

2
3

and the van Genuchten-Mualem model (Mualem,1976) for the hydraulic conductivity function:8
1=m m 22
< KS S0:5h,0e 1 2 1 2 Se  K h 
4
:
K
h$0
S

in which a and n are shape parameters, Se is effective
saturation, m  1 2 1/n, u r and u s are the residual andsaturated water contents, respectively, and Ks is thesaturated hydraulic conductivity.
2.2.2. Macropore domain
Water within a macropore is assumed to move bygravitational ¯ow from the surface node down to thebottom of the macropore. The hydraulic conductivity
at each node is calculated as a function of the relative
saturation of the macro-node according to:8>umac n* 0>>
0 , umac , us;mac
K>< s;mac E
Kmac umac  >
Ks;mac0mac  us;mac
>
>
>
:
0
umac  05

72

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

The exchange of water between the macro- andmicro-domains is calculated for each macro-node bycomparing the pressure head of the macro-node withthat of corresponding micro-node. The exchange rateof water, Qe,i, between the two ¯ow regions at macronode i is described as proposed by Beven andGermann (1981) in analogy to Darcy's law:
dhi
L2c

6

dhi  hmic;i 2 hmac;i

7

Qe;i  2Ki hmic in which:where Ki(hmic) is the hydraulic conductivity of micronode i, and Lc is the characteristic length representing
the average spacing of the macropores as well as thegeometry of the macroporous system. Gerke and van
Genuchten (1993) used a similar equation in theirmodel and derived the following expression for thecharacteristic length:
L2c 

a2 e
b

8

where b is a geometric factor (equal to 3 for rectangular slabs and 15 for spherical shapes), a is the
distance between the centre of the soil aggregates orsoil matrix domain to the centre of the macropores,
and e is an empirical factor set equal to 0.4.
For some cases the exchange of water is set to zero.
For example, when dhi , 0, i.e., pressure headgradient is directed towards micropore i, but at thesame time micropore i is saturated or macropore i isunsaturated. Also, Qe,i is set to zero when dhi . 0, i.e.,pressure head gradient is directed towards macropore
i, but micropore i is unsaturated.
After calculating the water exchange rate between
the two domains, the water content of each macronode is updated using the expression:

umac;new;i  umac;old;i 1 Qe;i dt

10

If u mac,new,i, calculated in this manner exceeds themacroporosity, Ei, of macro-node i, then its value isset equal to Ei and Qe,i is correspondingly adjusted.
Similarly, if u mac,new,i, is negative, its value is set equal
to zero and Qe,i is correspondingly reduced. The newwater contents, as discussed next, are starting values
for calculating water ¯ow within the macropores.

2.2.3. Redistribution of water within macropores
Redistribution of water in the macropores is calculated in two steps. The ®rst step involves solution of
the equations:

2umac;i2K umac;i 
2 Qe;i
2t
2z i

11

where 2u mac,i, represents the change in water content
of macro-node i. Eq. (11) is calculated from thesurface node down to the bottom node of the macropore. After solving Eq. (11), the macropores arechecked for over-saturation (u mac,i . Ei) and thesurplus of water is redistributed from bottom to thetop. Eq. (11) is solved using an implicit, iterativeinterval halving method (Gerald and Wheatley, 1989).
2.3. Solute transport
2.3.1. Solute transport in the micro-domain
Solute transport in the micro-domain is describedby the convection-dispersion equation including
source/sink terms for transformation processes that
follow ®rst- and zero-order kinetics (Simunek et al.,1994):


2uc
2
2c
2q c

uDij
2 i 1 mw uc 1 ms sr
2xi
2t
2x i
2x i
12
1 gw u 1 gs r 2 Scs
where c is the solution concentration, s is the sorbed
concentration, qi is the ith component of the volumetric ¯ux, mw and ms are ®rst-order rate constantsfor solutes in the liquid and solid phases, respectively,g w and g s are zero-order rate constants for the liquidand solid phases, respectively, r is the soil bulk
density, S is the sink term, cs is the concentration of
the sink term, and Dij is the dispersion coef®cienttensor.
2.3.2. Solute concentration of macro-nodes and soluteexchange between domains
Solute transport within the macro-domain andbetween the domains is calculated by assuming apurely convective process. Water with a certain soluteconcentration entering a macro-node is immediately
mixed with water associated with that macro-node.
After mixing, water is moved down to the next
macro-node as follows. First, water ¯ow, Qmac,i,A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

73Fig. 3. Drain discharge versus groundwater table above the tile drain, midway between drain as a function of varying parameters in Hooghoudt
drain boundary condition.

from macro-node i vertically downward to node i 1 1
is calculated according to:
Qmac;i 1u
2 umac;new;i  1 Qmac; i21 1 Qe;idt mac;old;i13

The procedure starts at the surface node and proceeds
downward. Using the water ¯uxes calculated in thismanner for each macro-node, a solute mass balance iscarried out next. To perform the mass balance, theconcentration of the exchange ¯ow is set equal to theconcentration of the macro-node if the exchange ¯ow is
directed toward the micropore, and set equal to theconcentration of the micro-node if the ¯ow is directed
toward the macropore. The system of solute massbalance equations for all nodes of a macropore is subsequently solved using Gaussian elimination.
2.4. Hooghoudt boundary condition
The traditional way to model an ideal drain in a®eld-scale ®nite element model is to represent thedrain as a node with a system-dependent boundary
condition in which the pressure head is set to zero
when the drain node is saturated, and the ¯ux is set
to zero when the drain node is unsaturated. Thisapproach does not account for the contraction of thestreamlines near the drain, the presence of a nonuniform pressure along the drain boundary, and thepresence of an entrance head around the drain.

Entrance head is the minimum head that is required
to initiate ¯ow in a drain. A comparison of discharges
calculated with the analytical solution of Kirkham(1949) and numerically simulated drain ¯ows showedthat the discrepancy could be substantial. Fipps andSkaggs (1986) demonstrated that only a model representing the drainage tube at suf®ciently high resolution and in its full shape would lead to accurate
estimation of the groundwater table and the drainage
discharge rate.
Lacking suf®cient spatial resolution to explicitly
describe the local non-uniformity of the ¯ow ®eldaround and along slotted drainpipes, ®eld-scalemodels usually try to account for the head loss causedby streamline convergence by invoking an empiricalresistance. The dependence of this resistance on
perforation sizes and arrangements, as well as on thedimensions and conductivity of the surrounding material has been theoretically analyzed by Dierickx(1980), among others. Other factors such as cloggingof ¯ow paths by particle deposition and disturbance ofthe soil around drain pipes or the envelope material
may further add to the non-uniformity of the ¯ow ®eld
at a drain boundary. As it is usually neither possible
nor necessary to separate all these individual effects,
the concept has been extended to include all otherforms of resistance encountered by the ¯ow entering
a drain. This compound resistance as described in
Stuyt et al. (2000) is called the `approach ¯owresistance'.

74

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

Table 1
Soil hydraulic parameters used for the MACRO and M-2D model comparison testab

C b  boundary soil water tension in the MACRO which is similar to hCritS in M-2D.
l  Brooks and Corey water retention shape factor set equal to (n 2 1).

Simunek et al. (1994) represented drains in the
SWMS_2D code as single nodes with a systemdependent boundary condition, while accounting for
approach ¯ow resistance in the entrance conditions byadjusting the conductivity of the ®nite elements
surrounding the drainage node according to theequation:
Kdrain  KS Cd

14

where Kdrain is the adjusted conductivity of theelements surrounding the drain, Ks is the actual
hydraulic conductivity of the elements surrounding
the drain, and Cd (,1) is a coef®cient that depends
on the diameter of the drain and the dimensions of thesurrounding elements. This approach enables
SWMS_2D to describe the hydraulic behaviour ofideal drains accurately. An ideal drain is de®ned inthis paper as a drain with a homogenous ¯ow ®eldFig. 4. A comparison of the simulated drainage out¯ow and the solute leaching rate between M-2D, MACRO, and SWMS_2D.

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

75

equation proposed by Oosterbaan et al. (1989), i.e.:q  A hT 2 he  1 B h2T 2 h2e Fig. 5. Water retentions and ¯ow curves based on Brooks and Coeryand van Genuchten formulations used to compare M-2D with
MACRO.

around and along the drainpipe. Although the reduction factor Cd accounts for the contraction of the ¯ow®eld around the drainage pipe, it may not account for
the typically non-uniform entrance head of non-ideal
drains.
In a recent ®eld study (Kohler, 2000) we studied theperformance of a ®eld tile-drainage system bymeasuring the groundwater drawdown during periodsof high groundwater tables. We found that the ®elddrains were subjected to relatively large entrance
heads, he0 $ 0, that were a function of the watertable height above the tile depth, midway between®eld drains. To be able to simulate such non-ideal
drainage system with an entrance head, we implemented in the M-2D code a new type of boundary condition, further referred to as a Hooghoudt boundary
condition. This condition prescribes the discharge
rate q at the drain node as a function of the watertable height hT at some reference point, i.e.,q  f(hT). For the experimental ®eld studied we
obtained good results using the extended Hooghoudt15

in combination with the empirical ®nding that he 
ChT 1 he0, in which C and he0 are site-speci®c parameters, A and B are ®tting constants, and where hT was
taken midway between the drains. The q(hT)-function
is site speci®c and generally should be measured.
However, if no suitable q 2 hT data are available,
the relationship may also be determined numericallyby means of local high-resolution modelling of thedrain performance, after which the q(hT)-function
can be used in ®eld-scale models in which the drains
are represented as single nodes. Note that the above
empirical relationships can easily be replaced byanother locally calibrated model relating draindischarge to the water table height measured at aconveniently location.
The Hooghoudt boundary condition is hence a
Neumann type boundary condition in which tile
drainage ¯ow is prescribed according to the head at
a reference node. The new boundary condition has the¯exibility to describe a large number of drain situations since the parameters A and B in Eq. (15) can be
determined empirically. This is illustrated in Fig. 3
where drain discharge is plotted as a function of thegroundwater elevation midway between the drains
and changing values of A, B, and he0. Notice the risein the groundwater table needed in order to initiate the¯ow in the case where he0 is not equal to zero.
A more accurate way to account for the drains
would be to implicitly allow for ¯ow processes inthe drain by including the drain geometry in the¯ow domain. The only problem with this approach
would be that a high-resolution discretization would
be required around the drains, as already pointed out
by Fipps and Skaggs (1986), which would substantially increase the number of nodes and, hence theexecution time. Such a slow model would not thenlend itself easily to the type of inverse analysis carriedout in this study for parameter estimation.
3. Comparing M-2D and macro simulation results
We tested the macropore component of M-2D bycomparing M-2D simulation results with thoseobtained using MACRO (1991) for a one-dimensional

76

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

test problem. The problem considers the arrival of a
wetting front at the bottom of a 1.5-m long soil column.
An atmospheric boundary at the soil surface and aseepage face at the bottom of the column wereassumed. The set up hence resembles a free-draining
lysimeter. The simulation lasted 40 days, starting witha 10-day period of precipitation at a rate of 0.01 m d 21,
followed by 30 days of water redistribution without any
rainfall or evapo-transpiration. More detail about thistest site can be found in Abbaspour et al. (2001).
The initial soil water pressure head in the column
was assumed to be in hydraulic equilibrium with thebottom of the soil set at zero pressure head. The soilcolumn was initially free of solutes, while the rainfall
had a constant solute concentration of 1000 mg m 23.
Pertinent information regarding the soil column islisted in Table 1. The values in Table 1 were calculated by inverse modeling in a previous study by
Abbaspour et al. (2001). Other parameters include
dispersivity for solute transport, which was set equal
to 0.1 m, and the exchange length (macropore
spacing) between the ¯ow regions, which was set
equal to 0.2 m in both models.
Fig. 4 shows calculated drainage and solute leaching rates obtained with M-2D and MACRO. The same
system, but without macropores, was also simulatedusing SWMS_2D for comparison purposes. Drainage
as well as solute leaching started earlier for MACRO
as compared to M-2D, but no solute leached from thecolumn during the simulation period when
SWMS_2D was used. It should be noted that
MACRO uses the Brooks and Corey (1964) modelof soil hydraulic characteristics, while SWMS_2D
uses the van Genuchten (1980) model. Although thisproblem is not important for simulations under highly
unsaturated conditions, it may have some effect at andnear the air-entry pressure head where most of theexchange of water between the macro- and microdomains takes place. The two hydraulic characteristic
models are illustrated graphically for the ®rst soil
layer in Fig. 5. Differences between the two hydraulicfunctions may be important in terms of affecting thetiming (initiation) and extent of a bypass ¯ow throughthe macropores.
Another important difference between MACRO
and M-2D may be the method in which the waterexchange rate between the micro- and macropore
regions is calculated. In MACRO the exchange rateis calculated as a function of the difference between
the actual micropore water content and the boundarywater content, whereas in M-2D water exchange isbased on the pressure head differences between amicro-node and a macro-node. Geometric discretization in MACRO is also different from that in M-2D;
MACRO allows only a limited number of soil layers
to be considered in a simulation (set to 10 in thisexample), while M-2D allows essentially an unlimited
number of soil layers. Given these differences, thesimulation results for MACRO and M-2D in Fig. 5
are quite comparable, especially the peak tile-drainage and concentration results.

4. Testing M-2D in a ®eld application
4.1. Field site
In a real ®eld test of the M-2D model we simulated
the hourly drainage discharge of a municipal solid
waste incinerator (MSWI) bottom ash land®ll near
Buchs, Switzerland. Lostorf is situated in a disused
gravel pit. After calibrating the model parameters fordischarge, solute transport through the land®ll was
simulated using measured electrical conductivity ofthe discharge water. The land®ll is 6 m deep, and aliner consisting of 0.8 m opalinous clay supporting agravel drainage layer (0.2 m) serves to collect theleachate via high-density polyethylene (HDPE) tubing
into a shaft. Two geotextiles separate the MSWI
bottom ash from the drainage layer and the drainage
layer from the clay liner. The surface of the land®ll has
not yet been covered or cultivated so there is directcontact between the atmosphere and the disposed ash.
The average yearly precipitation at the local Swiss
Metrological Association (SMA) station (Buchs/Suhr,approximately 2 km west of the land®ll in a comparable geographical location) for the period of 1987 toand including 1996 is 1060 mm. The average rainfall
maximum occurs between May and June. The driest
period is January to April. The summer rains tend tooccur in storm events of high rain intensity.
We selected the above site for testing M-2D because
the site has already been modelled in a previous study(Johnson et al., 2001), and this would allow us acomparison of M-2D with other models. Johnson andco-workers compared the application of several models

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

77Fig. 6. Comparison of the M-2D simulation with measurements of: (a) the drainage discharge, (b) the cumulative discharge, and (c) the
conductivity of the seepage water for the 1996 calibration data set.

such as a linear storage model (Huwe et al., 1994), aneural network model, HYDRUS5 (Vogel et al., 1996),
and MACRO (Jarvis, 1991) for the simulation of ¯owthrough the Lostorf land®ll. The hourly nature of thesimulation plus the apparent existence of preferential¯ow makes modeling of ¯ow and transport through theland®ll especially challenging. The hourly simulationis necessary for an accurate accounting of the soluteload leaching from the land®ll.
4.2. Input data
Automatic registration of drainage discharge, electrical conductivity and temperature were made overthe periods November 1993 to February 1993,
November 1994 to November 1995, and May 1996
to December 1996. Flowtec (DI 652) instrumentation,
based on changes in the magnetic ®eld as a function of¯ow rate, was used to measure discharge. Conductivity and temperature were registered on-line with a
combined WTW LF 196 electrode. Average values
taken every second were saved to a data logger
every 15 min. Precipitation measurements weremade on site using a tilting-siphon rain gauge. The
on-site rainfall measurements were compared to datacollected by the Buchs/Suhr SMA station.
Potential evapotranspiration from the surface of theland®ll was estimated using the hourly form of the
FAO Penman-Monteith equation (Allen et al.,1994). Hourly values of global radiation, relative
humidity, pressure, wind velocity (at 2 m) andtemperature were obtained from the SMA at theBuchs/Suhr site and in addition from on-site measurements of temperature, relative humidity and windvelocity in 1996.
The actual evaporation Ea from the bare bottom ash

78

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

Table 2
The results of statistical comparisons between measured and simulated cumulative discharge for different models
Model

Calibration Objective function, gValidation objective function, g

Neural network model aLinear storage model a
HYDRUS5 a
MACRO aM-2D

7.42
40.4 (29.4.96±31.8.96), 7.2 (1.9.96±31.12.96)
27.6
10.1
13.1

15.1
8.2
32.1
41.1
19.4

a

Data published in Johnson et al. (2001).

surface was estimated from the potential evapotranspiration using an empirical approach adopted fromBlack et al. (1969) by applying the following correction:
Ea 

Epot
d0:2

16

where d is the number of days after a 5 mm rain event.
Using this empirical correction approach, the total

evaporation for the 1996 data sets was corrected from3960 m 3 to 2950 m 3 and agreed well with the difference
betweenrainfallanddischarge(5860±
2940  2950 m 3). The agreement for the 1995 data set
was not as good. The evaporation estimated from rainfall, minus discharge (3490 m 3) was higher than thecorrected Penman-Monteith value (3010 m 3).
Using the data set 29.4.96±31.12.96 (5928 h),Fig. 7. Comparison of the M-2D simulation with measurements of: (a) the drainage discharge, (b) the cumulative discharge, and (c) the
conductivity of the seepage water for the 1995 test data set.

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81

M-2D was calibrated and the parameters were optimized. The model was then used to predict thedischarge monitored in the period 1.12.94±7.2.95(2125 h). Missing discharge and conductivity datacaused by instrumentation failure were linearly interpolated. Such periods did not exceed 10 h and constituted less than 1% of the data sets.
4.3. Model parameterization
The initial land®ll condition was expressed by aknown pressure head distribution over the verticaldepth of the land®ll. The boundary condition at the
top of the land®ll was depicted as an atmosphericboundary condition with prescribed precipitationand evapotranspiration data, while at the bottom aseepage face boundary condition was imposed tosimulate the drained bottom of the land®ll. Initialand boundary conditions applied to the ¯ow throughthe land®ll were mathematically formulated as:² Initial condition:
h z; t  hi z t  017

² Atmospheric condition:

2k


2h
2 1  q 0 t z  0
2z18

² Seepage face condition:h L  0 if h L . 0;

19

and:

2h
 0 if h L , 02z
where hi is the initial pressure head (mm), L is thecoordinate of the bottom of the soil, and q 0 (t) is theprescribed ¯ux (mm h 21) at the surface.
Initial estimates of u r, u s, a and n were taken from
unpublished laboratory desorption curves of intact
MSWI bottom ash samples taken from the land®ll
Lostorf (Buchter, 1997). The values were u r  0.14
u s  0.42 mm 3 mm 23, a  0.0085
mm 3 mm 23,
21
mm and n  1.18. The SUFI program of Abbaspour

79

et al. (1997) was used to inversely ®t the unknown
parameters of the micropore region a , n, u r, u s, and Ks
as well as the parameters n *, E, Ks,mac, and LC describing the macropore region. The objective function, g,
was expressed as:v
uXu mu
Q 2 Q 0 2i
u
t i
g
m

20

where m  5928 is the number of data points, Q ismeasured cumulative discharge, and Q 0 is simulated
cumulative discharge.
4.4. Simulation results
It is important to emphasise that the simulation with
the M-2D model in this study was conducted with theminimum of required information. The initial pressure
head hi, the parameters a , n, u r, u s, and Ks as well as
the parameters describing the macropore region hCritS,n *, E, Ks,mac, LC were ®rst ®tted on the measured cumulative discharge assuming a uniform distribution of
the initial pressure head in the land®ll.
At the beginning of the simulation period in April
1996 land®ll was in a very dry state. In order toobtain values for the initial pressure head, estimates
were made using the parameters obtained by the
initial ®tting. The ¯ow 3000 h beyond the end of
the simulation period was modelled without rainfall
until the discharge rates were equal to the measured
discharge at the beginning of the simulation period
in April. Using the pressure head values at this stage
as the initial values, the inverse estimation of thehydraulic parameters was once again performed.
The results of simulations were greatly improvedby the new non-uniform-initial-pressure-head data.
Fig. 6a shows the measured and simulated hourly
discharge as a function of time, while Fig. 6b illustrates the measured and simulated cumulative
discharge. Although there are some differences in
the hourly discharge rates of the measured andsimulated data, the cumulative values are rather
closely matched (g  13.1). In Table 2, the resultsof simulations with M-2D are compared with othermodels reported by Johnson et al. (2001). The estimated values of the parameters of the micropore

80

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81
21

region were a  0.0015 mm , n  1.10, u r  0.14
mm 3 mm 23, u s  0.45 mm 3 mm 23, and Ks  0.58
mm h 21, and of the macropore region n *  5,
E  0.002 mm 3 mm 23, Ks,mac  6 mm h 21, hCritS 
2450 mm, and Lc  1 mm.
We tested the above parameters, obtained by ®tting
the 1996 discharge data, by simulating the hourly
discharge rates measured in 1995 (Fig. 7). The agreement between modelled and measured cumulative
discharge in Fig. 7b is excellent up to about 1600 h
and diverges beyond this time. The overall value ofg  19.4 is, however, comparable with other models(Table 2). The discrepancy in the hourly discharge
rates in Fig. 7 indicates an overestimation of the fast¯ow component conducted by preferential ¯ow paths,
while the simulated matrix ¯ow is not as responsive as
the measured values.
As shown in Table 2 the simulation results with
M-2D are comparable to those of the MACRO
model published in Johnson et al. (2001). The
peak ¯ow rates simulated with M-2D are toohigh and too sharp compared to the measured
hourly discharge, whereas the base ¯ow from thematrix is too damped. Some reasons for these
discrepancies could be the assumption of a homogeneous system for the land®ll, and lack of accurate knowledge of initial and boundary conditions.
For this study we assumed a homogeneous systembecause no information was available on system
heterogeneity, and assumption of a multiplayer
system would have drastically increase the computational effort.
In addition to the water ¯ow we also simulated thesolute transport by means of measured electricalconductivity of the seepage ¯ow at the bottom of
the land®ll. Assuming a conductivity of 4 ms mm 21
for the precipitation and a conductivity of
14 ms mm 21 for the water stored in the matrix ofthe land®ll, we obtained a simulated conductivity ofthe seepage water as shown in Fig. 6c for the calibration phase and Fig. 7c for the validation case. Since
the measured conductivity of the seepage water issimulated rather well by M-2D, it can be assumed
that the model describes correctly the dynamics of
the water ¯ow within the land®ll system. It should
be noted that the model was only ®tted to the
cumulative discharge and not to the solute transport
information.

5. Conclusions
We extended the SWMS_2D model of Simunek etal. (1994) by adding a macropore ¯ow component andan entrance head in the drainage boundary condition.
In this paper, simulation of ¯ow and transport througha land®ll system tested the macropore component of
the model.
An inverse program ®tted parameters describing
the ¯ow, during calibration phase based on cumulative discharge. The calibrated model was tested
using data from a previous year. Although themodel over estimated the peaky dynamics of themacropore ¯ow and under estimated the base ¯owthrough the matrix, but based on the relatively small
amount of available information, the model
performed satisfactorily and the discharge comparedwell with other models tested by Johnson et al.
(2001). Simulation of the solute transport componentof the model using electrical conductivity as a tracer
also provided relatively good results for both calibration and validation stages.
References
Abbaspour, K.C., Kohler, A., Schulin, R., Fritsch, M., van Genuchten, M.Th., 2001. Application of a two-dimensional unsaturated¯ow and transport model to a macroporous agricultural ®eldwith entrance resistance to tile drains. European J. of Soil
Sci., 52 (3).
Abbaspour, K.C., van Genuchten, M.T., Schulin, R., SchlaÈppi, E.,
1997, A., 1892. sequential uncertainty domain inverse procedure for estimating subsurface ¯ow and transport parameters.
Water Resour. Res. 33, 1879±1879.
Allen, R.G., Smith, M., Pereira, L.S., Perrier, A., 1994. An update
for the calculation of reference evapotranspiration. ICID Bulletin, 43, 35±92.
Beven, K., Germann, P., 1981. Water ¯ow in soil macropores. II: A
combined ¯ow model. J. Soil Sci. 32, 15±29.
Black, T.A., Gardner, W.R., Thurtell, G.W., 1969. The prediction of
evaporation, drainage and soil storage for a bare soil. Soil Sci.
Soc. Amer. Proc., 33, 655±660.
Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous
media. , Hydrology Paper no. 3. Colorado State University, Fort
Collins, Colorado.
Buchter, B., 1997. Personal communication.
Dierickx, W., 1980. Electrolytic analogue study of the effect of
openings and surrounds of various permeability on the performance of ®eld drainage pipes. Rijksstation voor landbouwtechniek, Publication 77, Merelbeke, Belgium.
Evans, R.O., Gilliam, J.W., Skaggs, R.W., 1989. Effects of agricultural water table management on drainage water quality. Report

A. Kohler et al. / Journal of Hydrology 254 (2001) 67±81
237 of the Water Resource Research Institute of the Universityof North Carolina.
Fipps, G., Skaggs, R.W., 1986. Drains as a boundary condition in®nite elements. Water Resour. Res. 22, 1613±1621.
FOEFL, 1997. Verminderung des NaÈhrstoffeintrags in die GewaÈsser
durch Massnahmen in der Landwirtschaft. Schriftenreihe
Umwelt Nr. 209 vom Bundesamt fuÈr Umwelt, Wald und Landschaft.
Gerald, C.F., Wheatley, P.O., 1989. Applied Numerical Analysis.
Addison-Wesley Publishing Company, New York.
Gerke, H.H., van Genuchten, M.Th., 1993. A dual porosity model
for simulating the preferential movement of water and solutes instructured porous media. Water Resour. Res. 29, 305±319.
Huwe, B., GoÈlz-Huwe, H., Eberhardt, J., 1994. ParameterschaÈtzung
und Madellrechnungen zum Gebietswasserhaushalt kleiner,
heterogener Einzugsgebiete mit einfachen Modellkonzepten.
Mitteilgn. Dtsch. Bodenkundl. Gesellsch., 74, 273±276.
Jarvis, N., 1991. MACRO A model of water movement and solutetransport in macroporous soils. Department of Soil Sciences,
Report and Dissertations 9, Swedish University of Agricultural
Sciences, Uppsala.
Johnson, C.A., Schaap, M.G., Abbaspour, K.C., 2001. Modelcomparison of ¯ow through a municipal solid waste incinerator
ash land®ll. J. of Hydrol. 243, 55±72.
Kirkham, D., 1949. Flow of ponded water into drain tubes overlaying an impervious layer. Am. Geophys. Union Trans. 30,
369±385.
Klavidko, E.J., van Scoyoc, G.E., Monke, E.J., Oates, K.M., Pask,
W., 1991. Pesticide and nutrient movement into subsurface tile
drains on a silt loam soil in Indiana. J. Environ. Qual. 20, 264±
270.
Kohler, A., 2000. Investigation of water and solute out¯ow frommacroporous agricultural ®elds with tile drains. PhD Dissertation, Swiss Federal Institute of Technology ZuÈrich. In preparation.
Lennartz, B., Michaelsen, J., Wichtmann, W., Widmoser, P., 1999.
Time variance analysis of preferential solute movement in a tiledrained ®eld site. Soil Sci. Soc. Am. J. 63, 39±47.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., van Genuchten,
M.Th., 1997. New piece-wise continuous hydraulic functions
for modeling preferential ¯ow in an intermittent ¯ood-irrigated®eld. Water Resour. Res. 33, 2049±2063.
Mohanty, B.P., Bowman, R.S., Hendrickx, J.M.H., Simunek, J., van

81

Genuchten, M.Th., 1998. Preferential transport of nitrate to atile drain in an intermittent-¯ood-irrigated ®eld: Model development and experimental evaluation. Water Resour. Res. 34,
1061±1076.
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour. Res.
12, 513±522.
Oosterbaan, R.J., Pissarra, A., van Alphen, J.G., 1989. Hydraulichead and discharge relations of pipe drainage systems withentrance resistance. p. 86-98. In Proceedings of 15th European
Regional Conference on Agricultural Water Management. Vol.
III ICID, Dubrovnik.
Richard, T.L., Steenhuis, T.S., 1988. Tile drain sampling of preferential ¯ow on a ®eld scale. J. Contam. Hydrol. 3, 307±325.
Simunek, J., Vogel, T., van Genuchten, M.Th., 1994. The
SWMS_2D Code for simulating water ¯ow and solute transport
in tow dimensional variably saturated media (Version 1.21).
Research Report No. 132. U.S. Salinity Laboratory Agricultural
Research Service, U.S. Department of Agriculture, Riverside,
California.
Stamm, C., 1997. Rapid transport of phosphorus in drained grassland soils. PhD Dissertation No. 12486, Swiss Federal Instituteof Technology ZuÈrich.
Stuyt, L.C.P.M., Dierickx, W., Martinez Beltran, J., 2000. Materials
for subsurface land drainage systems. FAO Irrigation andDrainage Paper 60, Rome, Italy.
van Genuchten, M.Th., 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc.
Am. J. 44, 892±898.
van Genuchten, M.Th., Wierenga, P.J., 1976. Mass transfer in sorbing porous media. I: Analytical solutions. Soil Sci. Soc. of Am.
J. 40, 473±480.
van Genuchten, M.Th., Sudicky, E.A., 1999. Recent advances invadose zone ¯ow and transport modeling. In: Parlange, M.B.,
Hopmans, J.W. (Eds.). Vadose Zone Hydrology: Cutting AcrossDisciplines. Oxford University Press, New York, pp. 155±193.
Vogel, T., Huang, K., Zhang, R., Van Genuchten, M. Th., 1996. The
HYDRUS code for simulating one-dimensional water ¯ow,solute transport, and heat movement in variable-saturated
media, version 5. U.S. Salinity Laboratory, ARS, Riverside,
California.
Youngs, E.G., 1980. The analysis of groundwater seepage in heterogeneous aquifers. Hydrol. Sci. Bull. 25, 155±165.