Agricultural Water Management 146 (2014) 105–114

Contents lists available at ScienceDirect

Agricultural Water Management
journal homepage: www.elsevier.com/locate/agwat

Numerical simulation of water ﬂow in tile and mole drainage systems
Vilim Filipović a,b,c,∗ , Fábio Joel Kochem Mallmann d , Yves Coquet c , Jirka Šimůnek eaDepartment of Soil Amelioration, Faculty of Agriculture, University of Zagreb, Svetošimunska 25, 10000 Zagreb, Croatia
INRA, AgroParisTech, UMR 1091 EGC, F-78850 Thiverval-Grignon, France
Univ d’Orléans, ISTO, UMR 7327, 45071, Orléans, France; CNRS/INSU, ISTO, UMR 7327, 45071 Orléans, France; BRGM, ISTO, UMR 7327, BP 36009,
45060 Orléans, FrancedFaculdade Ingá/UNINGÁ, Rodovia PR 317, n. 6114, CEP 87035-510, Maringá, PR, BrazileDepartment of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USAbca r t i c l e

i n f oArticle history:
Received 17 December 2013
Accepted 16 July 2014
Keywords:
Tile drainage
Mole drainage
Numerical simulation
Heavy soil
Three-dimensional modelinga b s t r a c t
Tile drainage systems are sometimes not sufﬁcient to provide favorable unsaturated conditions in the
rootzone. These drainage systems then need to be supplemented with an additional high conductivity
material in the trenches above the tiles or by implementing mole drainage. The HYDRUS (2D/3D) model
was used to evaluate the impact of such additional measures for heavy clay soil. Three types of drainage
systems were simulated: (i) tile drains, (ii) tile drains with gravel trenches, and (iii) tile drains with gravel
trenches and mole drains, using either two-dimensional (the former two systems) or three-dimensional(the latter one) transport domains. Three scenarios were considered to test the efﬁciency of each system:(i) time to drain an initially saturated system, (ii) high intensity rainfall, and (iii) a real case scenario.
Different horizontal spacings between tile drains with or without gravel trenches were also comparedwith the system which included mole drainage. The results showed that the drainage system that included
mole drains and gravel trenches was the most efﬁcient. This system provided the largest drainage rate,
was the ﬁrst to reach steady-state in the time to drain scenario, and also efﬁciently reduced surface
runoff. Adding mole drains to a system with tile drains and gravel trenches resulted in a large reductionof surface runoff (75%). Simulations showed that the spacing of tile drains with or without gravel trenches
would have to be 40% or 55% smaller, respectively, in order to reproduce the same water table levels as
those observed for the drainage system with mole drains. Therefore, introducing mole drains in drainage
systems is an efﬁcient practice for reducing waterlogging and runoff.
© 2014 Elsevier B.V. All rights reserved.

1. Introduction
Soil drainage systems aim at limiting saturated conditions in thesoil proﬁle that can arise due to hydrological processes from above(downward water percolation) and from below (elevated groundwater). The main goal of drainage systems is to remove excesswater and maintain favorable unsaturated conditions in the rootzone. The most common agricultural drainage system consists ofperforated PVC tile drains installed in soil at various spacings anddepths, depending on soil hydraulic properties, climatic conditions,and cultivated crops. When designing drainage systems, it is crucial
to increase the hydraulic functioning of the entire system. These∗ Corresponding author at: University of Zagreb, Department of Soil Amelioration,
Faculty of Agriculture, Svetošimunska 25, 10000 Zagreb, Croatia.
Tel.: +385 1 2393711; fax: +385 1 2315300.
E-mail address: vﬁlipovic@agr.hr (V. Filipović).
http://dx.doi.org/10.1016/j.agwat.2014.07.020
0378-3774/© 2014 Elsevier B.V. All rights reserved.

systems are mostly installed in soils with high clay content, i.e.,heavy soils with very low hydraulic permeability (Tuli et al., 2005).
In such soils, subsurface drainage systems can substantially reducesurface runoff, shorten periods of surface ponding, and lower watertable (Konyha et al., 1992; Skaggs et al., 1994). Under certain
conditions, including high groundwater table, large and intensive precipitation, and heavy-textured soils, the installation of tile
drains may not be sufﬁcient to provide favorable conditions for
growing crops. Additional measures may then be needed, such as
using a backﬁll material (gravel) with a high hydraulic conductivity above tile drains or performing mole drainage. The presence of agravel layer above tile drains up to the tilled layer may increase the
efﬁciency of the entire system by promoting by-pass ﬂow of waterfrom the tilled layer directly into the drains. Mole drains are closelyspaced unlined channels of limited duration that are formed in clay
subsoil using a ripper blade with a cylindrical foot, often with an
expander, which helps to compact and stabilize the channel walls.
Mole drainage has been recommended for heavy soils with low

106

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

permeability, which would otherwise require a small drain spacing(Hudson et al., 1962). During the construction of the mole channels
the foot and expander create a ‘leg slot’ directly above the moleas well as ﬁssures in the soil upwards from the mole toward theplough layer. These cracks and ﬁssures can promote preferential
ﬂow toward mole drains (Leeds-Harrison et al., 1982). Mole drains
are intended to improve lateral ﬂow to tile drains and are usuallyused only in soils with a high clay content (>35%). If establishedproperly and under the right conditions, they may still be operating
after ﬁve years (Harris, 1984). Tile and mole drains are nonconductive under unsaturated conditions because positive pressure mustoccur before water can start ﬂowing into them (Stormont and Zhou,2005).
In recent years, numerical models have been developed tosimulate water ﬂow in tile-drained soils. One-dimensional (1D)models (e.g. MACRO, Larsbo and Jarvis, 2003) use an approximate
analytical solution of Darcy’s equation to account for water loss
through the drainage system (e.g., Hooghout, Ernst, or Boussinesq’s
equations), while two- or three-dimensional (2D or 3D) models (e.g. HYDRUS (2D/3D), Šimůnek et al., 2006) use the Richards
equation to explicitly model the dynamics of the water table.
In 1D models, the ﬂow to tile drains is implemented as a sinkterm in the mass balance equation. The Houghoudt’s equation isbased on the Dupuit–Forchheimer assumptions with corrections
for convergence of radial ﬂow near the tiles. When combining theHoughoudt’s equation (for the saturated zone) with the Richards
equation (for the unsaturated zone), the Houghoudt’s equation
gives instant lateral ﬂuxes to tile drains, while the Richards equation gives transient ﬂuxes in the unsaturated zone, leading to a riseor drop of the water table. An example of such model is DrainMod, which has been used in many applications. For example,
Skaggs et al. (2012) used DrainMod to simulate water ﬂow in a
subsurface-drained agricultural ﬁeld in eastern North Carolina. The
performance statistics indicated that the model with calibrated
input data accurately predicted daily water table depths, daily
drainage rates, and monthly drainage volumes. Singh et al. (2006)
calibrated and tested DrainMod on two types of soils in Iowa andused it to simulate the impacts of different designs of subsurface
drainage systems. Simulation results suggested that a drainage system designed for a drainage intensity of 4.6 mm d−1 with a drain
depth of 1.05 m and a drain spacing of 25 m was sufﬁcient tomaximize crop production under the prevailing local agricultural
conditions.
Water ﬂow and solute transport to tile drains has been also evaluated using the HYDRUS family of codes (Šimůnek et al., 2008). For
example, the HYDRUS-2D software package was used by De Vos
et al. (2000, 2002) to simulate nitrate transport in a tile-drained layered silt loam soil in a reclaimed Dutch polder, or by Castanheiraand Serralheiro (2010) to evaluate the impact of mole drains on
salinity of a vertisol under irrigation.
Numerical modeling of drainage systems involving both tileand mole drains is much more limited. Snow et al. (2007) used
the APSIM-SWIM model to predict drainage rates and runoff ina mole-tile drained silty loam soil of New Zealand. There was an
excellent agreement between simulated and measured drainage,
as well as a reasonable agreement between measured and simulated cumulative surface runoff. Armstrong et al. (2000) performeda modeling study on a macroporous clay soil with mole and tile
drains, in which they compared four “preferential ﬂow” models(MACRO, CRACK-NP, SIMULAT and PLM) in their ability to simulate
isoproturon leaching. MACRO model gave the best results, although
globally the simulations showed the difﬁculty of deriving adequate
parameters, even where relatively complete soil physical data wereavailable. A similar study was performed by Besien et al. (1997) ona structured heavy clay soil with tile (0.75 m depth and 50 m spacing) and mole drains (0.5 m depth and 3 m spacing) in which the

MACRO model was calibrated and used to investigate the leaching of isoproturon. Madvar et al. (2007) used the SEEP/W model
to simulate the hydraulic performance of mole and tile drains inheavy-textured soil. Based on their numerical study the combination of mole and tile drains resulted in economical and hydraulicimprovement of the agricultural system.
Most of modeling studies considered water ﬂow in one- or
two-dimensional soil proﬁles, either taking into account only tile
drainage without mole drainage or mole drainage without tile
drainage (e.g., Castanheira and Serralheiro, 2010), but not both in a
fully three-dimensional system. This is mainly because mole drains,
being installed in a perpendicular direction to tile drains, makethe entire system three-dimensional. Since we could not identify
any study that simulated water ﬂow in systems with tile and moledrains as a three-dimensional problem, we conducted 3D numerical experiments to evaluate the performance of these systems. In
our numerical simulations we considered environmental conditions that are typical for the eastern part of Croatia. Hydromorphic
soils in Croatia cover an area of 1618,500 ha, which is approximately one third (29%) of the total area occupied by agriculturalsoils (Husnjak, 2007). Subsurface tile drainage systems are installedon 161,530 ha (Petošić et al., 2004), with mole drainage often used
to improve drainage efﬁciency.
The main objective of this paper was to numerically evaluate theperformance of different subsurface drainage systems of increasing complexity in hydromorphic soils using a three-dimensionalmodel that can explicitly account for any drainage system. (i) The
simplest system to be considered consists of a layered soil proﬁle with parallel tile drains. Such system can be simulated using
a two-dimensional simulation domain that is perpendicular to tile
drains. (ii) In the second set of simulations it is assumed that a
backﬁll material (gravel) of higher hydraulic conductivity than theoriginal soil is placed above the drains. Such system can also be
evaluated assuming a two-dimensional simulation domain perpendicular to drains. (iii) Finally, the most complex system consists oftile drains with backﬁlled gravel above drains, combined with perpendicular mole drains. Such system has to be analyzed using afully three-dimensional ﬂow model.
2. Materials and methods
2.1. Theory
Subsurface water ﬂow in the unsaturated and saturated
zones is governed by the mass balance equation and theDarcy–Buckingham’s law:∂∂J
= − wi − S(h)∂t∂xiJwi = −K(h)KijA

(1)



∂H
∂h
= −K(h) KijA
+ KijA
∂xj
∂xj


(2)

respectively, where  is the volumetric water content [L3 L−3 ],Jwi is Darcy’s ﬂux [LT−1 ] in the i-th direction, S is a sink/source term
[T−1 ], K(h) is the saturated/unsaturated soil hydraulic conductivity function [LT−1 ], Kij A are the components of the dimensionlessanisotropy tensor for hydraulic conductivity KA [−], H is the total
hydraulic head [L] deﬁned as the sum of the pressure head and thegravitational head, H = h + z, xi are the spatial coordinates [L], withx and y being horizontal coordinates and z the vertical coordinate
directed upwards, and t is time [T]. The governing ﬂow equation
for these conditions, resulting from combining (1) and (2), is givenby the following modiﬁed form of the Richards’ equation:

 

∂
∂
K
=
∂t
∂xi

∂h
KijA
+ KijA
∂xj



−S

(3)

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

107

Table 1
van Genuchten soil hydraulic parameters for all materials of the soil proﬁle used in the numerical simulations.
Material

Depth (cm) r (cm3 cm−3 ) s (cm3 cm−3 )˛ (cm−1 )nKs (cm day−1 )

l

1
2
3
4
5

0–30
30–70
70–100100–200
Gravel (trench
and mole)
material

0.095
0.095
0.095
0.102
0.005

0.42
0.41
0.41
0.50
0.42

0.001360.002120.001360.0121
0.1

1.20
1.18
1.20
1.41
2.1

12
14
12
12.8
3000

0.5
0.5
0.5
0.5
0.5

The unsaturated hydraulic conductivity function K(h) is givenby:K(h, x, y, z) = Ks (x, y, z)Kr (h, x, y, z)(4)where Kr is the relative hydraulic conductivity [−] and Ks is thesaturated hydraulic conductivity [LT−1 ].
Soil hydraulic functions were described using the van
Genuchten–Mualem model (van Genuchten, 1980), which isdeﬁned as follows:s − r(h) = r +(h) = s(1 + |˛h|n )

m





1/m

 m 2 − rs − rm=1−(5)

for h ≥ 0

K(h) = Ks Sel 1 − 1 − Se
Se =for h < 0

1;n > 1n(6)(7)

and Parrish (1988) for sand, with an increased Ks value to imitate
the hydraulic properties of gravel (Table 1). As only one seepage
face boundary condition could be considered in HYDRUS (2D/3D),
mole drains were simulated as a gravel material with hydraulicproperties identical to those of the material used for tile trenches.
2.2.2. Drainage systems
Three different drainage systems were evaluated in our numerical simulations. The simplest system (System 1) consisted ofparallel tile drains and no other additional measures (denoted
below in graphs and tables as “tile”). In simulations for System
2 (denoted as “tile gravel”), it was assumed that a backﬁll material (gravel) of higher hydraulic conductivity was placed above
drains. Finally, the most complex system (System 3) consisted oftile drains with a backﬁlled gravel above drains combined withperpendicular mole drains (denoted as “tile gravel mole”). Fig. 1
shows a schematic of the simulated drainage systems. System 1 and2 were evaluated assuming a two-dimensional transport domain
perpendicular to tile drains (Fig. 1a), while System 3 needed a fully
three-dimensional water ﬂow model (Fig. 1b).

(8)

where  r and  s denote residual and saturated volumetric watercontents [L3 L−3 ], respectively; Se is the effective saturation [−], ˛[L−1 ], and n [−] are retention curve shape factors, and l is a pore connectivity parameter [−]. A numerical model for simulating waterﬂow and solute transport in two- and three-dimensional transportdomains, HYDRUS (2D/3D) (Šimůnek et al., 2006), was used it thisstudy.
2.2. Numerical experiments and input data
2.2.1. Soil hydraulic properties
Numerical experiments were set up for a Eutric Haplic Calcaric
Siltic Gleysol (horizons: Ap; Bg; Cr; Cg), which is located in eastern
Croatia (45◦ 09 N and 18◦ 42 E). Disturbed soil samples were takenfrom each horizon and particle size distribution was determinedusing the pipette method (Gee and Or, 2002). Undisturbed soil samples (100 cm3 ) were used to measure bulk density and soil hydraulicproperties. The saturated hydraulic conductivity, Ks , was measured
using the constant head method (Klute and Dirksen, 1986). The saturated water content,  s , was measured using the ISO 11274:1998
sandbox method (Clement, 1966). Data points of the soil water
retention curve were measured using the pressure plate apparatus (Dane and Hopmans, 2002) for applied pressures of 3, 10, 33,
100, 625, and 1500 kPa. While the saturated water content,  s , wasmeasured, the remaining parameters of the soil water retention
curve (5) ( r , ˛, and n) were optimized using the RETC software(van Genuchten et al., 1991) by ﬁtting measured data. R2 values
for soil water retention curve ﬁtting were in range from 0.91 to
0.99. The pore connectivity parameter, l, was assumed to be equal
to an average value for many soils (l = 0.5) (Mualem, 1976). The
hydraulic properties of the material that was assumed to be backﬁlled to the tile trenches were taken from the soil catalog of Carsel

2.2.3. Transport domain
Because of the symmetry of water ﬂow between two parallel
tile drains, it was possible to consider only a half domain between
two tile drains. A two-dimensional transport domain for Systems
1 and 2 (Fig. 1a) thus had a width of 750 cm (tile drain spacing was
assumed to be 15 m) and a height of 200 cm. The 8 cm diameter tile
drain, which was located on the left side of the transport domain,
was assumed to be centered at a depth of 96 cm from the bottomof the soil proﬁle. The soil proﬁle was divided into 4 soil horizons(Table 1) based on a soil survey of the area. In simulations of System
2, an additional vertical layer of gravel, representing a backﬁll material (material 5; Table 1), was considered. The three-dimensionaltransport domain for System 3 (Fig. 1b) with tile and mole drains
and a backﬁll material above the tile drains had the same ﬁrst two
dimensions as the two-dimensional transport domains, i.e., a widthof 750 cm and a depth of 200 cm. The third dimension, perpendicular to tile drains, was equal to 100 cm (mole drain spacing was
assumed to be 2 m). The mole drain, which was perpendicular totile drains, was located at the front of the three-dimensional transport domain (Fig. 1b) and was assumed to have a diameter of 6.5 cm,
an inclination of 0.5%, and was 40 cm deep (from the soil surface)
above the tile drain.
2.2.4. Simulated scenarios
Numerical simulations were carried out for each drainage
system for three different scenarios with different initial andboundary conditions. In Scenario 1, a “time to drain” scenario wasused to compare the drainage response of the three simulated
systems, i.e., tile drains, tile drains with gravel trenches, and tile
drains with gravel trenches and mole drains, starting from wet
soil conditions. The initial pressure head was assumed to be atequilibrium conditions with a positive pressure head of 175 cm
at the bottom of the transport domain which left the ﬁrst 25 cm

108

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

Fig. 1. Transport domains with 4 or 5 soil materials considered in the numerical simulations (a) a two-dimensional transport domain was considered for the tile drains andthe tile drains plus gravel trenches systems, (b) a three-dimensional transport domain was considered for the tile drains with gravel trenches plus mole drains system.

(corresponding to the tilled layer) unsaturated. No ﬂow boundary
condition was considered for all boundaries except for the small
semicircle on the left of the transport domain which represented
the tile drain (Fig. 1a). A seepage face boundary condition was used
to represent the tile drain. Simulation time was equal to 2400 h(100 days). However, as there was still some limited outﬂow after
2400 h, the time when drainage was assumed to stop was deﬁned
as the time when the drainage rate decreased below 0.001 mm h−1 .
Simulations for Systems 1 and 2 were also performed in 3D tocheck if the results were identical to those provided by 2D
modeling.
In the “rainfall” scenario (scenario 2), the initial conditions were
set up as the pressure head distribution from the end of the ﬁrst
scenario (steady state conditions), and a rainfall event (with anintensity of 2 mm h−1 and a duration set up to get a total of either
50 mm, 100 mm, or 150 mm) was applied starting at the ﬁrst day ofsimulation. In addition to a seepage face boundary condition representing a tile drain, an atmospheric boundary condition was appliedat the top of the transport domain to include rainfall. Evaporation
was neglected in this scenario.
In the “real case” scenario (scenario 3), 4-year simulations(2009–2012) were performed with meteorological data collected
at the Gradište station (45◦ 09 N and 18◦ 42 E) using the same
initial conditions as in scenario (2). An atmospheric boundary
condition was applied at the top of the transport domain using
evapotranspiration values calculated by the Penman–Monteith
approach (Monteith, 1981). Transpiration was calculated for maize(Zea Mays L., 2009, 2010, 2012) and barley (Hordeum vulgare L.,2011) which were grown at the actual location. Again, a seepage face boundary condition was selected to represent the tile
drain.

3. Results and discussion
3.1. “Time to drain” scenario
To evaluate the speciﬁc drainage response of each of the threedrainage systems, a “time to drain” scenario was performed. The
ending time of drainage was reached after 681 h (28.4 days) for
the tile system, 594 h (24.7 days) for the tile gravel system, and506 h (21.1 days) for the tile gravel mole system. The value corresponding to 90% of the ﬁnal drainage volume from each system
was selected as the time when a quasi steady-state condition
was reached. According to this criterion, the systems reachedquasi steady-state after cumulated outﬂow amounts of 4.86 mm,
6.57 mm, 9.18 mm for the tile, tile gravel, and tile gravel mole systems, which corresponded to days 11, 7.1, and 6, respectively(Fig. 2). The difference in time response of the three systems may
have a large inﬂuence on the crop root system, i.e., the faster thesoil is drained the less consequences on the crop root system andyield can be expected (Oosterbaan, 1991). The soil horizons had arelatively low permeability and, in the absence of a highly conductive layer, the drainage is very slow. When a gravel trench
above the tile is added, the time response is faster (Fig. 2) since
the trench is ﬁlled up with a material that has a high saturated
hydraulic conductivity. The response of the system with moles isfastest since it includes both the trench and the mole drain constructed from the same highly conductive material as the trenches.
Since the domain was initially saturated well above the moledrain depth, the response of the system is increasing in the ordertile < tile gravel < tile gravel mole.
In addition to the time response of each system, a large
difference in the cumulative amount of drainage was found.

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

109

Fig. 2. Cumulative outﬂow rates for the “time to drain” scenario.

This difference is due to the substitution of soil materials
having different water retention properties within the simulation
domain. In System 2, parts of the volumes of materials 1, 2, and3 are replaced by material 5 to represent the gravel trench. In
System 3, in addition to the trench, a part of material 2 is replacedby material 5 to represent the mole drain. The difference in water
storage between System 1 and System 2 is 1.86 mm and 4.86 mm
between System 1 and 3. Material 5 compared to other materials
has very different water retention properties, which explain the
different cumulative outﬂows (Table 2). Inserting a gravel trench
above the tile drain increased outﬂow by 1.86 mm, 90% of it comingfrom the water initially stored in the trench itself. Adding a moledrain (System 3) further increases the outﬂow by 3 mm, withmore than half of it coming from the gravel material representing
the mole drain. From these results one can already expect asigniﬁcant effect of the mole drain on the recession of the watertable.
Fig. 3 shows the simulated cross sections of all three systems
during the ﬁrst 24 h. The results are displayed after 2 and 24 h. Pressure head values are lower in the tile gravel system than in the tile
system and the lowest in the tile gravel mole system. After 2 h ofdrainage, the inﬂuence of the trench in System 2 is limited to theright side of the domain while in System 3 the effect is expanding
more to the left side because of the presence of the mole drain. After
24 h, there is a large drop in the water table depth for all three systems in the order (high to low) tile > tile gravel > tile gravel mole.
The mean pressure head after 24 h in the domain above the tile
depth (from the surface down to 100 cm) was −75 cm, −78 cm,
and −84 cm for the tile, tile gravel, and tile gravel mole system,
respectively. These results indicate the higher efﬁciency of thetile gravel mole system in reducing the water table height compared to the tile gravel system.

3.2. High intensity rainfall scenario
High intensity rainfall events with an intensity of 2 mm per hourand various durations, to get total amounts of 50 mm, 100 mm, or
150 mm, were applied in the second scenario. The initial condition
was set equal to the ﬁnal pressure head distribution at the endof the ﬁrst scenario, i.e., the surface pressure head was −108 cm
and the bottom pressure head was 92 cm on the right side of thedomain. The cumulative outﬂow increased with rainfall in all threesystems (Fig. 4). In this scenario, the trench and mole drain were
not initially saturated as in the previous scenario, which led to theslower response time at the beginning of the experiment. The gravelmaterial (representing the trench and mole) was not initially saturated and the time needed for its saturation induced a delay in thesystem response, which can be best seen for the tile gravel molesystem during the ﬁrst 4 h (Fig. 4). There is no outﬂow during theﬁrst 4 h in System 3 since the mole drain needs to become saturated
to trigger drainage. The tile system without a gravel layer becamealmost instantaneously saturated since the rest of the domain had
a very low permeability (5 mm h−1 ).
The low permeability of the soil promotes runoff during highintensity rainfalls. In Fig. 5, the total amount of surface runoff ispresented for each system. One can see that the cumulative amountof runoff increases according to rainfall, as for cumulative outﬂow(Fig. 4). The largest runoff occurs in the tile system, then in thetile gravel system, and the smallest in the tile gravel mole system.
The mole drain was very efﬁcient in reducing surface runoff. These
results are in agreement with several studies, in which authors
found that artiﬁcial drainage increases the amount of inﬁltratingwater and reduces runoff (e.g., Bengtson et al., 1995).
Mass balance information for all three system and three rainfall
events is shown in Table 3. Table 3 shows that the effectiveness

Table 2
Water balance (WS—change in storage) for the various material layers composing the three simulated systems at the end of the “time to drain” scenario (mm).
All Materials
WS in eachsystemDifference between
the systems

Tile
Tile gravel
Tile gravel mole
Tile vs Tile gravel
Tile vs Tile-gravel mole

−5.35
−7.21
−10.21
−1.86
−4.86

Mat 1

Mat 2

Mat 3

Mat 4

Mat 5

−1.12
−1.14
−1.14
−0.02
−0.02

−0.87
−0.96
−1.91
−0.09
−1.04

−0.68
−0.72
−0.72
−0.04
−0.04

−2.68
−2.68
−2.68
0.00
0.00

0.00
−1.71
−3.76
−1.71
−3.76

110

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

Fig. 3. Pressure head distributions after 2 and 24 h for the tile drains, tile drains with gravel trenches, and tile drains with gravel trenches and mole drains systems.

Fig. 4. Cumulative outﬂow in the second scenario with high intensity rainfall events of 50 mm, 100 mm, and 150 mm.

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

111

Fig. 5. Cumulative surface runoff in the second scenario with high intensity rainfall events of 50 mm, 100 mm, and 150 mm.

of mole drainage in reducing surface runoff is not only due to abetter efﬁciency in promoting drainage outﬂow, but also due toa better capacity for water storage than the other two drainage
systems.
3.3. Real case scenario
For the real case scenario, agroecological conditions that correspond to East Croatia were taken into account to be able to explain
the behavior of each system in real ﬁeld conditions. Fig. 6a shows
the cumulative outﬂow during the 2009–2012 period. At the beginning, one can see that the order of drainage relative to cumulative
outﬂow is a bit different from what was found in ﬁrst two scenarios
in which the tile gravel mole system had the largest outﬂow. Thiseffect is due to the small rainfall amounts at the beginning of thesimulated time period which were drained by the tile gravel system, but not by the tile gravel mole system, because mole drains
did not become saturated enough to trigger outﬂow. Long term simulations showed the opposite effect on the 1st of June 2010 whena high intensity rainfall event (72.8 mm) saturated mole drains
and promoted a large outﬂow from the tile gravel mole system,
while in the same time the other two systems generated low outﬂow and large surface runoff (Fig. 6b). Adding gravel in the trench
above a drain resulted in small decrease of surface runoff by (1%)
However, inserting a mole drain into the system decreased therunoff by an additional 75%. The tile drain system with gravel andmole drains effectively reduced surface runoff in real ﬁeld conditions. In such situation, the tile gravel system cannot conduct suchlarge amounts of water as the system with moles, which leads towaterlogging conditions at the surface or the presence of a highwater table.

3.4. Drain spacing
Finally, simulations with different tile drain spacings were performed in order to compare the water table variations for eachsystem and to try to obtain, either with only tile drains or tile
drains with backﬁlled trenches, the same water table level as for
System 3. The simulation set up i.e., initial and boundary conditions, was taken from the “time to drain” scenario. The tile withgravel trenches and mole system had the same spacing of 15 m asin above studied scenarios. The tile and tile gravel systems were set
up with various spacings from 6 to 12 m. Fig. 7a and b shows thewater table elevation during 96 h in the middle between two tile
drains for the three systems (numbers 6 to 12 corresponds to thedrain spacing in meters). At the beginning, there is a rapid decreasein water table depths in all three systems, clearly visible during
the ﬁrst 24 h. After that, the decrease is smaller and starts to stabilize as time approaches the end of simulation. The tile gravelsystem was more efﬁcient than the tile system at each spacing
in lowering the water table (by ∼2 cm). To reproduce the samewater table level after 96 h (−100.2 cm) as in the tile gravel molesystem, one would need to lower the spacing of the tile gravel system to 9 m, or to 7 m if the system without gravel is used. System
1 and 2 with a 6-m spacing were able to lower the water table
more efﬁciently than System 3 with a 15-m tile drain spacing anda 2-m mole drain spacing. However, a system with such small tile
drain spacing is not economically sustainable. A small tile drain
spacing is mostly used only in golf course construction or in civilengineering. Since the construction of mole drains requires only asimple equipment (i.e., a special kind of plough) (Hopkins, 2002),
this system could signiﬁcantly lower the overall cost of the drainage
system.

Table 3Mass balance for all three simulated systems for scenario 2—subscripts 50, 100, and 150 correspond to the cumulative rainfall amount (mm).

RainfallOutﬂow
Runoff
WS

System 150

System 250

System 350

System 1100

System 2100

System 3100

System 1150

System 2150

System 3150

50.00
16.57
31.93
1.49

50.00
19.58
29.69
0.72

50.00
28.13
6.36
15.51

100.00
29.71
68.66
1.63

100.00
35.65
63.52
0.83

100.00
51.22
13.61
35.18

150.0042.81
105.391.79

150.00
51.69
97.35
0.96

150.00
74.26
20.85
54.89

112

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

Fig. 6. Cumulative outﬂow (a), cumulative runoff (b) and cumulative surface ﬂux (c) for the real case scenario (2009–2012) in East Croatia.

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

113

Fig. 7. Comparison of water table levels for different drain spacings between: (a) tile system and tile gravel mole system (15 m), (b) tile gravel system and tile gravel molesystem (15 m).

4. Conclusions
Numerical experiments have been performed for three different drainage systems: (i) tile drains, (ii) tile drains with gravel
trenches, and (iii) tile drains with gravel trenches and moledrains. The tile drain system showed in all scenarios the slowest
outﬂow response and the smallest cumulative outﬂow. When
compared with the other two systems, the drainage system withtile drains, gravel trenches, and mole drains provides the largest
efﬁciency. Under real case scenario, the system with tile drains
and gravel trenches had the fastest reaction time and provided
drainage outﬂow for small intensity rainfalls. The system withtile drains, gravel trenches and mole drains was able to provide
larger outﬂow during higher intensity rainfall (>15 mm). Adding
mole drains to a system with tile drains and gravel trenches
resulted in a very large reduction of surface runoff (75%), which
can have a very important effect during high intensity rainfalls.

Additional simulations showed that the spacing of tile drains withor without gravel trenches would have to be 40% or 55% smaller,
respectively, in order to reproduce the same water table levels as
those observed for the drainage system with mole drains. Since
the use of the drainage systems with mole drains is still limited to
special agro-ecological conditions, there is only limited extensiveﬁeld data available which allow to test fully three-dimensional
simulations of these particular drainage systems. Research needs tobe expanded through ﬁeld trials in order to have more information
on water ﬂow dynamics in such complex systems.

Acknowledgements
We acknowledge the ﬁnancial support by GENESIS project on
groundwater systems (FP7) and Agreen Skills project (INRA) for
V. Filipović postdoctoral fellowship. F.J.K. Mallmann acknowledges

114

V. Filipović et al. / Agricultural Water Management 146 (2014) 105–114

CNPq (Conselho Nacional de Desenvolvimento Cientíﬁco e Tecnológico – Brazil) for postdoctoral scholarship.
References
Armstrong, A., Aden, K., Amraoui, N., Diekkrüger, B., Jarvis, N., Mouvet, C., Nicholls, P.,
Wittwer, C., 2000. Comparison of the performance of pesticide-leaching models
on a cracking clay soil: results using the Brimstone Farm dataset. Agric. Water
Manage. 44, 85–104.
Bengtson, R.L., Carter, C.E., Fouss, J.L., Southwick, L.M., Willis, G.H., 1995. Agriculturaldrainage and water quality in Mississippi delta. J. Irrig. Drain. Eng.-ASCE 121,
292–295.
Besien, T.J., Jarvis, N.J., Williams, R.J., 1997. Simulation of water movement and isoproturon behaviour in a heavy clay soil using the MACRO model. Hydrol. Earth
Syst. Sci. 1, 835–844.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769.
Castanheira, P.J.N., Serralheiro, R.P., 2010. Impact of mole drains on salinity of avertisoil under irrigation. Biosystems Eng. 105 (1), 25–33.
Clement, C.R., 1966. A simple and reliable tension table. J. Soil Sci. 17, 132–133.
Dane, J.H., Hopmans, J.W., 2002. Water retention and storage: laboratory. In: Dane,
J.H., Topp, G.C. (Eds.), Methods of Soil Analysis: Part 4, Physical Methods. , third
ed. SSSA, Madison, WI, pp. 675–680.
De Vos, J.A., Hesterberg, D., Raats, P.A.C., 2000. Nitrate leaching in a tile-drained silt
loam soil. Soil Sci. Soc. Am. J. 64, 517–527.
De Vos, J.A., Raats, P.A.C., Feddes, R.A., 2002. Chloride transport in a recently
reclaimed Dutch polder. J. Hydrol. 257 (1–4), 59–77.
Gee, G.W., Or, D., 2002. Particle size analysis. In: Dane, J.H., Topp, G.C. (Eds.), Methodsof Soil Analysis: Part 4, Physical Methods. , third ed. SSSA, Madison, WI, pp.
255–293.
Harris, G.L., 1984. Effect of mole submergence on the life of mole channels. Agric.
Water Manage. 8, 361–374.
Hopkins, D., 2002. Managing wet soils: mole drainage. In: Agriculture Notes,
AG0949, Available on: http://www.cag.org.uk/docs/Mole%20drains.pdf
(accessed December 12th, 2013).
Hudson, A.W., Hopewell, H.G., Bowler, D.G., Cross, M.W., 1962. The Draining of Farm
Lands. Massey College, Palmerston North.
Husnjak, S., 2007. The Croatian agricultural soils and the requirements for theamelioration measures. In: Proceedings of the Scientiﬁc Symposium Ameliorative Measures Aimed at Rural Space Development With Focus on the National
Irrigation Project, Croatian Academy of Art and Science, Zagreb, pp. 21–37.
Klute, A., Dirksen, C., 1986. Hydraulic conductivity and diffusivity: laboratory methods. In: Klute, A. (Ed.), Methods of Soil Analysis, Part 1, Physical and Mineralogical
Methods. Agronomy Monograph 9. ASA-SSSA, Madison, WI, pp. 687–734.
Konyha, K.D., Skaggs, R.W., Gilliam, J.W., 1992. Effects of drainage and watermanagement practices on hydrology. J. Irrig. Drain. Eng.-ASCE 118, 807–819.

Larsbo, M., Jarvis, N., 2003. MACRO 5.0. A model of water ﬂow and solute transport in macroporous soil. Technical description. In: Emergo Report 2003:6.
Swedish University of Agricultural Sciences, Department of Soil Sciences, Uppsala, Sweden.
Leeds-Harrison, P., Spoor, G., Godwin, R.J., 1982. Water ﬂow to mole drains. J. Agric.
Eng. Res. 27, 81–91.
Madvar, H.R., Samani, J.M.V., Salami, M., 2007. Modeling hydraulic performanceof mole drains and combined mole-pipe drains. In: Salokhe, V.M., Jayasuriya,
H.P.W., Soni, P. (Eds.), Cutting edge technologies and innovations on sustainable
resources for world food sufﬁciency, Proceedings of the International Agricultural Engineering Conference, Bangkok, unpaginated.
Monteith, J.L., 1981. Evaporation and surface temperature. Q. J. Meteorol. Soc. 107,
1–27.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522.
Oosterbaan, R.J., 1991. Agricultural land drainage: a wider application through caution and restraint. In: ILRI Annual Report 21-36. International Institute for Land
Reclamation and Improvement, Wageningen.
Petošić, D., Tadić, L., Romić, D., Tomić, F., 2004. Drainage outﬂow in different pipedrainage variants on gleyic podzoluvisol in the Sava river valleyy. Irrig. Drain.
53, 17–27.
Šimůnek, J., van Genuchten, M., Šejna, M.Th., 2006. The HYDRUS software package
for simulating two- and three dimensional movement of water, heat, and multiple solute in variably-saturated media. In: Technical Manual, Version 1.0. PC
Progress, Prague.
Šimůnek, J., van Genuchten, M., Šejna, M.Th., 2008. Development and applications of
the HYDRUS and STANMOD software packages, and related codes. Vadose Zone
J. 7, 587–600.
Singh, R., Helmers, M.J., Zhiming, Q., 2006. Calibration and validation of DRAINMOD
to design subsurface drainage systems for Iowa’s tile landscapes. Agric. Water
Manage. 85, 221–232.
Skaggs, R.W., Breve, M.A., Gilliam, J.W., 1994. Hydrologic and water quality impactof agricultural drainage. Crit. Rev. Environ. Sci. Technol. 24, 1–32.
Skaggs, R.W., Youssef, M.A., Chescheir, G.M., 2012. DrainMod: model use, calibrationand validation. Trans. ASABE 55, 1509–1522.
Snow, V.O., Houlbrooke, D.J., Huth, N.I., 2007. Predicting soil water, tile drainage,
and runoff in a mole-tile drained soil. N.Z. J. Agric. Res. 50, 13–24.
Stormont, J., Zhou, S., 2005. Impact of unsaturated ﬂow on pavement edgedrain
performance. J. Transp. Eng.-ASCE 131, 46–53.
Tuli, A., Hopmans, J.W., Rolston, D.E., Moldrup, P., 2005. Comparison of air and water
permeability between disturbed and undisturbed soils. Soil Sci. Soc. Am. J. 69,
1361–1371.
van Genuchten, M.Th., 1980. A closed form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–1037.
van Genuchten, M., Leij, F.J.Th., Yates, S.R., 1991. The RETC Code for Quantifying theHydraulic Functions of Unsaturated Soils. In: Version 1.0. EPA Report 600/291/065. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA.