This article appeared in a journal published by Elsevier. The attachedcopy is furnished to the author for internal non-commercial researchand education use, including for instruction at the authors institution
and sharing with colleagues.
Other uses, including reproduction and distribution, or selling or
licensing copies, or posting to personal, institutional or third party
websites are prohibited.
In most cases authors are permitted to post their version of thearticle (e.g. in Word or Tex form) to their personal website or
institutional repository. Authors requiring further information
regarding Elsevier’s archiving and manuscript policies areencouraged to visit:
http://www.elsevier.com/copyrightAuthor's personal copy
water research 44 (2010) 1050–1061

Available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/watres

Modeling water flow and bacterial transport
in undisturbed lysimeters under irrigationsof dairy shed effluent and water using HYDRUS-1D
Shuang Jiang a,*, Liping Pang b, Graeme D. Buchan a, Jiřı́ Šimůnek c,Mike J. Noonan a, Murray E. Close baFaculty of Agriculture and Life Sciences and CSEQ, PO Box 84, Lincoln University, Lincoln 7647, Canterbury, New ZealandInstitute of Environmental Science & Research, PO Box 29-181, Christchurch, New Zealandc
Department of Environmental Sciences, University of California Riverside, CA, USAbarticle info

abstract

Article history:HYDRUS-1D was used to simulate water flow and leaching of fecal coliforms and bromide

Received 29 March 2009

(Br) through six undisturbed soil lysimeters (70 cm depth by 50 cm diameter) under fieldReceived in revised form

conditions. Dairy shed effluent (DSE) spiked with Br was applied to the lysimeters, which

23 July 2009

contained fine sandy loam layers. This application was followed by fortnightly spray or

Accepted 27 August 2009

flood water irrigation. Soil water contents were measured at four soil depths over 171 days,

Available online 2 September 2009and leachate was collected from the bottom. The post-DSE period simulations yieldeda generally decreased saturated water content compared to the pre-DSE period, and an

Keywords:

increased saturated hydraulic conductivity and air-entry index, suggesting that changes inModeling

soil hydraulic properties (e.g. via changes in structure) can be induced by irrigation andDual porosity

seasonal effects. The single-porosity flow model was successful in simulating water flow

Bacterial transport

under natural climatic conditions and spray irrigation. However, for lysimeters under flood

Undisturbed lysimeterirrigation, when the effect of preferential flow paths becomes more significant, the good

HYDRUS-1Dagreement between predicted and observed water contents could only be achieved byusing a dual-porosity flow model. Results derived from a mobile-immobile transport model
suggest that compared to Br, bacteria were transported through a narrower pore-networkwith less mass exchange between mobile and immobile water zones. Our study suggeststhat soils with higher topsoil clay content and soils under flood irrigation are at a high riskof bacteria leaching through preferential flow paths. Irrigation management strategies
must minimize the effect of preferential flow to reduce bacterial leaching from land
applications of effluent.
ª 2009 Elsevier Ltd. All rights reserved.

1.

Introduction

Applying dairy shed effluent (DSE) onto land to fertilize
pasture and dispose of effluent is a permitted activity in New
Zealand. However, the use of inappropriate irrigation

practices could result in an increased risk of contaminatingwaterways and groundwater with DSE. For example, highfecal coliforms were detected in the groundwater underlying
soils treated with DSE and pasturing dairy cows plus flood
irrigation at an area in South Canterbury (Close et al., 2008).

* Corresponding author. Tel.: þ64 21 2179922; fax: þ64 3 3423689.
E-mail address: shuang.jiang@msn.com (S. Jiang).
0043-1354/$ – see front matter ª 2009 Elsevier Ltd. All rights reserved.
doi:10.1016/j.watres.2009.08.039

Author's personal copy1051

water research 44 (2010) 1050–1061

Table 1 – Soil properties of the lysimeters.
Depth (cm)Clay (%)

Silt (%)

Sand (%)

Bulk density(g cm-3)

Porosity

Texture classa

Organic matter (%)

pH

A

0–20
20–35
35–50
50–70

20
17
10
6

50
50
58
47

30
33
32
48

1.09
1.19
1.42
1.46

0.58
0.55
0.46
0.45

Silt loam
Silt loam
Silt loam
Sandy loam

6.72

4.41

B

0–20
20–35
35–50
50–70

7
8
6
3

56
52
40
33

37
40
54
64

1.27
1.23
1.44
1.45

0.51
0.53
0.46
0.45

Silt loam
Silt loam
Sandy loam
Sandy loam

6.21

4.97

C

0–20
20–35
35–50
50–70

20
9
9
5

55
56
56
33

25
35
35
62

1.22
1.28
1.43
1.54

0.53
0.52
0.46
0.42

Silt loam
Silt loam
Silt loam
Sandy loam

7.24

5.75

D

0–20
20–35
35–50
50–70

9
14
10
10

57
66
55
49

33
20
35
41

1.19
1.20
1.59
1.53

0.54
0.55
0.40
0.42

Silt loam
Silt loam
Silt loam
Loam

7.18

4.72

E

0–20
20–35
35–50
50–70

10
9
6
6

57
56
42
48

33
35
52
46

1.17
1.16
1.40
1.32

0.55
0.56
0.47
0.50

Silt loam
Silt loam
Sandy loam
Sandy loam

6.30

4.66

F

0–20
20–35
35–50
50–70

9
8
10
7

55
54
55
46

36
38
36
47

1.03
1.25
1.54
1.48

0.60
0.53
0.42
0.44

Silt loam
Silt loam
Silt loam
Sandy loam

6.07

4.11

Lysimetera McLaren and Cameron (1996). The organic matter and pH are for the top 10 cm only.

Knowledge of the fate and transport of bacteria in structured soils is crucial for assessing the risks for groundwater
quality of land-applied animal effluent and wastes. Contaminant transport models can be very useful tools for thequantitative assessment of microbial transport in soils and to
elucidate the important factors and processes that control
microbial transport. With some assumptions, complex
processes are reduced to a set of mathematical equations withparameters characterizing the system.
Although bacterial transport in intact soils has been
investigated (Smith et al., 1985; Hekman et al., 1995; McMurry

5
Aug-0

5
Sep-0

Oct-0

5

5
Nov-0

5
Dec-0

Jan-0

6

Feb-0

Rainfall
Evapotranspiration
DSE application
Spray irrigation
Flood irrigation

160
140

6
120

100

120

80

100
60
80
60

40

I r r ig a t i o n ( m m )M o n t h l y r a in f a ll o r e v a p o t r a n s p i r a t i o n ( m m )

180

40
20
20
0

0
-60

-40

-20

0

Days before DSE application

20

40

60

80

100

120

Days after DSE application

Fig. 1 – Comparison of rainfall, evapotranspiration (ET) andirrigation between pre- and post-DSE application periods.

et al., 1998; Aislabie et al., 2001; McLeod et al., 2001, 2003, 2004;
Guber et al., 2005), most of these studies do not include
modeling work. Most reported studies on modeling bacterialtransport are with repacked homogenous porous media(Hijnen et al., 2005; Bradford et al., 2006; Gargiulo et al., 2008;
Torkzaban et al., 2008). Comparatively fewer studies (Shelton
et al., 2003; McGechan and Vinten, 2003; Pang et al., 2008) havecarried out modeling of bacterial transport through undisturbed soils. However, data used for modeling in these studies
were largely obtained from experiments carried out indoors
under simulated rainfall (Shelton et al., 2003; Pang et al., 2008).
This may not be fully representative of behavior under natural
climatic conditions and irrigation.
Classically, water flow through variably-saturated soils is
described by the Richards’ equation with a uniform flow domain(Šimůnek and van Genuchten, 2008). Recently, efforts have been
made at simulating contaminant transport under the influenceof preferential flow using dual-porosity models (Gerke and van
Genuchten, 1993, 1996; Šimůnek et al., 2001, 2003; Jarvis, 2007;
van Genuchten and Šimůnek, 2004; Šimůnek and van
Genuchten, 2008). These dual-porosity models assume thecoexistence of two separate pore domains: fractures (or interaggregate pores, cracks and macropores); and matrix pores,
with water exchange between the two domains. A few types ofdual-porosity model have been incorporated into HYDRUS-1D(Šimůnek et al., 1998, 2005), the most extensively used model for
simulating water and contaminant transport in variably-saturated porous media.
Our objective in this paper is, using HYDRUS-1D as a tool, todevelop our understanding of (a) the dynamic changes of soil

Author's personal copy1052

water research 44 (2010) 1050–1061

0.50

D

A

0.45
0.40
0.35
0.30

10 cm depth
25 cm depth

0.25

45 cm depth
60 cm depth

0.20
0.50

-1
W a te r c o n t e n t (v v )

0.45

B

E

C

F

0.40
0.35
0.30
0.25
0.20
0.50
0.45
0.40
0.35
0.30
0.25
0.20
-50

-40

-30

-20

-10

0 -50

-40

-30

-20

-10

0

Days before DSE application
Fig. 2 – Comparison of water contents observed during the pre-DSE application period and predicted using the singleporosity flow model. Predicted values are in solid lines and observed values in empty circles. Black-depth 10 cm, red-depth
25 cm, green-depth 45 cm, and yellow-depth 60 cm. Same for Figs. 3–5. Note that a large portion of data was missing forlysimeter B due to instrument failure.

structure (as reflected in parameters of soil hydraulic properties)
under different irrigation regimes and climatic conditions, and(b) the effect of preferential flow on water flow and bacterialtransport. Based on our findings, recommendations were also
provided for management strategies to minimize the risk ofbacterial leaching through soils into groundwater.

2.

Materials and methods

2.1.

Leaching experiments and experimental data

The six lysimeters (50 cm diameter by 70 cm depth) used inthis study were originally collected in 1996 from a dairy farmnear Lincoln. The soil is an alluvial Templeton fine sandy loamwith lithologic units of silt loam and sandy loam, with texture
mostly coarser towards the lysimeter base (Table 1). Soilproperties determined from destructive sampling at the endof this study are also listed in Table 1. These lysimeters werepreviously used for a study of nitrogen leaching from DSE until
2002. The lysimeters were located outdoors in the groundalongside an open trench at Lincoln University’s Centre forSoil and Environmental Quality site. The lysimeters hada typical New Zealand pasture mix of ryegrass and whiteclover. Pasture was cut periodically to simulate typical grazing
practice.
Rainfall was recorded automatically at the site using
a tipping bucket gauge connected to a datalogger (Campbell

Author's personal copy1053

water research 44 (2010) 1050–1061

Table 2 – Hydraulic parameters optimized from observed water contents using the single-porosity flow model.
Lysimeter

Pre-application
1

1

1

qr (v v ) qs (v v ) a (m )nPost-application
1

Ks (m d ) SSRs

R

2

A

0.147
0.056
0.007

0.361a
0.440
0.457

0.080
0.210
0.274

1.942
1.176
1.123

7.119
7.000
0.252

0.023 0.978

B

0.003
0.088
0.158

0.465a
0.466
0.439

0.783
2.672
4.203

2.831
1.352
1.241

0.014
0.199
0.507

0.124 0.938

C

0.076
0.051
0.001

0.452a
0.415
0.375

2.479
5.623
5.512

1.697
1.338
2.423

0.121
0.223
0.000

D

0.010
0.004
0.085

0.424a
0.415
0.395

0.682
0.118
0.189

1.596
1.641
1.909

E

0.133
0.150
0.000

0.415a
0.377
0.402

0.114
0.112
0.011

F

0.001
0.039
0.025

0.380a
0.418
0.485

0.771
0.691
0.928

1

1

qr (v v ) qs (v v ) a (m1)n

Ks (m d1) SSRs

R2

0.042
0.103
0.000

0.308
0.348
0.349

1.340
0.800
0.286

1.300
1.305
1.415

6.339
7.128
1.993

0.435 0.565

0.165 0.836

0.058
0.065
0.128

0.382
0.366
0.350

0.634
0.600
0.267

1.788
1.829
1.606

1.168
7.000
3.370

0.377 0.626

0.008
2.672
1.414

0.084 0.916

0.027
0.173
0.202

0.418
0.430
0.404

0.601
0.343
0.617

1.788
1.959
1.654

7.000
0.268
0.313

0.288 0.728

1.243
1.221
2.120

0.160
0.565
0.000

0.075 0.929

3.125
2.280
1.872

0.006
1.316
0.006

0.037 0.965

0.049
0.007
0.014

0.358
0.368
0.393

1.428
1.045
1.255

2.638
2.298
3.714

1.405
1.989
0.484

0.202 0.803

SSRs and R2 values are for the whole lysimeter. Calibration of the single-porosity flow model was failed for post-application of lysimeters B and E.
a is not optimized but use maximum observed value.

Scientific, USA). Potential evapotranspiration data wereobtained from the Broadfield climate station, located about
2 km north of the lysimeter site. In this study each lysimeter
had a single soil moisture (TDR) sensor and temperature
sensor at each of four depths (10, 25, 45 and 60 cm). Sensor
data were collected with a datalogger and monitoring undernatural climatic conditions commenced on 14 August 2005.
More detailed information is contained in Jiang et al. (2008).
On 30 September 2005, 25 mm of DSE (spiked with NaBr)
collected on the same day was applied to all lysimeters in
5 min (Jiang et al., 2008). The chemical characteristics of the
DSE were Br 312 mg L1, fecal coliform 1.65  104 cfu mL1,
total N 220 mg L1, total C 1612 mg L1, electric conductivity
1.56 mS cm1 and pH 8.5. Fecal coliforms were not detected in
the leachate collected prior to the leaching experiments. The
lysimeters were irrigated with water on the 6th day after DSE
application and every two weeks afterwards. The experiments
were carried out for a period of 124 days. Flood irrigation(100 mm applied in 0.5 h on a 14-day cycle, extreme case) wasapplied to lysimeters B and E, which had very similar soil
properties (Table 1) with sandy loam appearing at muchshallower depth (35 cm) than in other lysimeters (50 cm
depth). Spray irrigation (50 mm applied in 1.5 h on a 14-day
cycle, a typical practice used in the region) was applied to theother four lysimeters. Irrigation and climatic data are given in
Fig. 1.
Leachate samples were collected intensively in the first
12 h after DSE application. After 12 h, samples were taken
once a day if available. During heavy rain, leachate was
sampled as it became available. Leachate was collected understerilized conditions. The bacterial samples were analyzed
within 12 h of sampling using the Membrane Filtration Technique (APHA, 1998) with mFC agar (Difco). The Br samples

were analyzed using a Bromide Ion-Selective Electrode (ISE,Orion 250A 96-35 BN, Thermo Scientific, Inc.).

2.2.

HYDRUS-1D model simulations

2.2.1.

Boundary conditions and root water uptake

For Hydrus-1D model input, measured potential evapotranspiration (PET) was converted to potential evaporation rate (PE)and potential transpiration rate (PT) based on the following
equation (Ritchie, 1972):
PE ¼ PET  exp ðrExtinct  LAIÞ(1)

PT ¼ PET  PE

(2)

Here rExtinct is a canopy extinction coefficient. A rExtinct
value of 0.6 is adopted from Ludlow (1985) for ryegrass, which
was the predominant vegetation for our lysimeters. Similarly,
Hay and Porter (2006) reported a rExtinct value of 0.59–0.60 for
wheat. LAI is leaf area index, generally LAI ¼ 3 for grassland,
adopting from Hay and Porter (2006) and Ludlow (1985).
However the LAI value was adjusted according to the healthstatus of the grass. For example, we used LAI ¼ 2 for lysimeterF because there was some vegetation loss due to grass grub.
A time-variable upper boundary condition was applied at theplant-atmosphere interface, with irrigation incorporated into
rainfall data.
The bottom of each lysimeter had an added 20–30 mm thicklayer of gravel. The base was sealed except for a central drainage
hole, open to the atmosphere, for collecting drainage samples.
No suction was applied to the base when samples were taken.
Water drained from the bottom only when it was saturated, thus
mathematically it was a seepage face at the lower boundaryAuthor's personal copy1054

water research 44 (2010) 1050–1061(Flury et al., 1999). The saturated depth at the bottom, a few mm,
was effectively the depth of the capillary fringe. The root zone
depth observed, during destructive sampling at the end of thestudy, was about 40 cm, with most roots between 5 and 10 cm
depths. During destructive sampling at the end of the study, weobserved that more roots were distributed near the top and root
density gradually reduced with depth. Therefore we assumed
a linear distribution of root water uptake (from 1 to 0 from top tobottom of the root zone), with uptake decreasing linearly withdepth. This allows more loss at the surface and no root uptake atthe bottom of the root zone. The function of Feddes et al. (1978)
was employed for root water uptake in HYDRUS-1D (Šimůnek
et al., 1998, 2005).

10 cm depth
25 cm depth

0.50

45 cm depth

Water content (v v-1)

60 cm depth

0.40

0.30

0.20

2.2.2. Water content simulations using the single-porosity
flow model

0.10
0

20

40

60

80

100

120

Days after DSE application

Fig. 3 – An example (lysimeter A) of the departure ofpredicted soil water contents from those observed for the
post-DSE application period, due to direct application ofparameter values derived from the pre-DSE application
period. A single-porosity model was used for both periods.

Measured soil moisture contents were simulated using themodified form of the Richards’ equation (Šimůnek et al., 1998,
2005) for water movement in unsaturated soils under onedimensional uniform flow:



vq
v
vh
¼
KðhÞ
þ1 S
vt vz
vz

(3)

Here q is volumetric water content [], h is the soil waterpressure head [L], t is time [T], z is the soil depth [L], and S is thesink term for root water uptake [T1]. The unsaturated
hydraulic conductivity K [LT1], as a function of h, is given in
the van Genuchten Equation (van Genuchten, 1980)

0.50

D

A
0.40

0.30

10 cm depth

Water content (v v-1)

0.20

25 cm depth
45 cm depth
60 cm depth

0.10
0.50

C

F

0.40

0.30

0.20

0.10
0

20

40

60

80

100

120 0

20

40

60

80

100

120

Days after DSE application
Fig. 4 – Comparison of soil water contents predicted from the calibrated single-porosity flow model and those observed
during the post-DSE application period, for lysimeters A, C, D and F (spray irrigation).

Author's personal copy1055

water research 44 (2010) 1050–1061

Water content (v v-1)

0.50

B

E

0.40

0.30

10 cm depth

0.20

25 cm depth
45 cm depth
60 cm depth

0.10
0

20

40

60

80

100

120

0

20

40

60

80

100

120

Days after DSE application
Fig. 5 – Comparison of observed soil water contents with those predicted from the dual-porosity flow model during the postDSE application period, for lysimeters B and E (flood irrigation).

(qðhÞ ¼

qs  qr
qr þ 
n m
1 þ jahj
qs

h<0

(4)

h0

h

m i2
KðhÞ ¼ Ks Sle 1  1  Se1=m

(5)

and
m ¼ 1  1=n ðn > 1Þ(6)

qr and qs denote the residual and saturated volumetric watercontents, respectively; Ks is the saturated hydraulic conductivity [LT1]; Se is effective saturation; l is the pore connectivity
parameter and is about 0.5 on average for many soils (Mualem,1976); n is an empirical parameter related to the pore-size
distribution, that is reflected in the slope of water retention
curve; and a is an empirical parameter, often assumed to berelated to the inverse of the air-entry suction (L1). Water
retention curves of coarse-textured soils are generally sharpwith steep slopes, thus high n values. A high a value (e.g. forsands and gravels) implies that soils de-saturate quickly at
small suctions, while a low a value (e.g. for clays) indicatesa slow de-saturation with increasing suctions.
Values of qr, qs, a, n and Ks were optimized by inverse
modeling of water contents measured at four depths using thevan Genuchten-Mualem single-porosity model from HYDRUS1D. As HYDRUS allows a maximum of 15 parameters to be
optimized, we reduced four soil layers (Table 1) to three layers
by combining two soil layers that had similar texturalcompositions and water content patterns. The combined layer
varied for different lysimeters. The initial inputs of qr, qs, a, n
and Ks were estimated from measured water retention data
(i.e. h vs. q) using RETC (van Genuchten et al., 1991), or from thepedotransfer functions (PTFs) based on soil texture (% of sand,silt and clay) and bulk density using the Rosetta package(Schaap et al., 2001), which is embedded in HYDRUS-1D. Initialwater contents measured were used as the inputs for initialwater contents in the model. Values of qr, qs, a, n and Ks optimized from the pre-DSE application model were used as initial
inputs for the post-DSE application model. To reduce theuncertainty of fitted parameters, we fixed saturated watercontents using measured maximum water contents for the
topsoil for the pre-DSE period simulations, reducing to a totalof 14 parameters to be optimized.

2.2.3. Water content simulations using the dual-porosity
flow model
For lysimeters B and E, the above single-porosity flow model
failed to simulate the water contents measured during theTable 3 – Hydraulic parameters optimized from water contents observed during the post-DSE application period using thedual-porosity flow model for lysimeters B and E (flood irrigaiton).
Layerqs (v v1)

Ks (m d1)a (m1)nu (d1 m1)

SSRs

R2

B

L1
L2
L3

0.255
0.160
0.039

2.659
2.963
2.028

0.367
0.085
0.002

1.168
1.173
1.526

0.014
0.003
0.000

0.403

0.627

E

L1
L2
L3

0.069
0.093
0.048

1.768
1.428
5.462

0.988
1.007
0.367

1.591
1.294
1.797

0.008
0.003
0.001

0.293

0.730

Lysimeter

SSRs and R2 values are for the whole lysimeter.

Author's personal copy1056

water research 44 (2010) 1050–1061

0.05

-1

B r (C C 0 )

A

Obs.
Pred.

0.04
0.03
0.02
0.01
0.00

0

10

20

30

40

50

-1

B r (C C 0 )

1.2

C

1.0
0.4

Here subscripts m and im refer to the mobile and immobile
water regions, respectively; q ¼ qm þ qim is the volumetricmoisture content [], Sim and Sm are sink terms (root water
uptake) for both regions [T1], and Gw is the transfer rate forwater exchange between macropores and matrix [T1]. We
assume that root water uptake is preferentially from macropores, thus Sim ¼ 0.
In the dual porosity flow model based on mass transfer
driven by differences in soil water pressure head, theexchange rate of water between the macropores and matrix
regions, Gw, is assumed to be proportional to the difference inpressure heads between the two pore regions (Gerke and van
Genuchten, 1993; Šimůnek et al., 2003):

0.2

Gw ¼ uðhm  him Þ
0.0

0

20

40

60

80

100

120

0.14

D

0.12

-1

B r (C C 0 )

0.10
0.08
0.06
0.04
0.02
0.00
0

20

40

60

80

100

120

0.12

F

0.10
0.08
-1

B r (C C 0 )

(9)
1 1

0.06
0.04
0.02
0.00
0

10

20

30

40

50

Days after DSE application

Fig. 6 – Br concentrations observed and predicted from the
MIM transport model and the single-porosity flow model,
in drainage from lysimeters A, C, D and F (spray irrigation).

post-DSE application period. During destructive sampling at
the end of the field study, we observed that these two lysimeters contained much larger macropores than the otherlysimeters. To consider the significant effect of preferential
flow on these two lysimeters, a dual-porosity flow model(based on mass transfer driven by differences in soil waterpressure head) was selected from HYDRUS-1D.
The dual-porosity formulation for water flow is based ona mixed formulation of the Richards equation to describewater flow in the macropores (mobile water region) anda mass balance equation to describe moisture dynamics in thematrix (immobile water region), as follows (Šimůnek
et al., 2003):

Here u is a first-order mass transfer coefficient (L T ). Sincepressure heads are now needed for both regions, thisapproach requires estimating retention curves for both poreregions. That means each region has its own values of qr, qs, a,and n. As a result, soil hydraulic properties are now describedby six parameters for macropores (qr, qs, a, n, Ks, l ), fourparameters for the matrix (qr-im, qs-im, aim, nim), and a parameter (u) for mass transfer between the two zones (Šimůnek
et al., 2003). Note that the soil total saturated water content(e.g. qs þ qsim) can be obtained from measurement. To reduce
the number of parameters optimized (maximum 15), we
assumed that there is no residual water in macropores andmatrix (i.e. qr ¼ 0, qr-im ¼ 0) and that the matrix is made of lesspermeable media with known values of qs-im, aim, and nim,
which were determined from soil texture and bulk density
using Rosetta package. Therefore five parameters (qs, Ks, a, n,
and u) were optimized for each soil layer (total 15 parameters).
Sometimes, qs or Ks were set at their measured values, or as
reference values, as a constraint to reduce the simulation
uncertainty.

2.2.4. Simulations of bacteria and Br transport using themobile-immobile contaminant transport model (MIM)
Regardless of whether a single- or dual-porosity model ischosen for the water flow, HYDRUS-1D allows simulation ofcontaminant transport using a two-region mobile–immobile
model, MIM (van Genuchten & Wagenet, 1989), The MIM
assumes that contaminant transport is limited to the mobile
water region and that water in the immobile water region isstagnant, with a first-order diffusive exchange process
between the two regions (Šimůnek et al., 2003, 2005).
As there was no instantaneous sorption of bacteria in ourdata and Br is non-reactive, the following simplified MIM(Šimůnek and van Genuchten, 2006, 2008) was used in thisstudy for simulating bacteria and Br transport:
vqm Cm
v
vCm
vqCm
¼

 fm  Gs
qm Dm
vt
vz
vz
vz

qim




vqm
v
vh
¼
KðhÞ
þ 1  Sm  Gw
vt
vz
vz

(7)

vqim
¼ Sim þ Gw
vt

(8)

vCim
¼ fim þ Gs
vt

(10)

(11)

Here C is the concentration in the liquid phase [ML3]; D is thedispersion coefficient [L2T1], which is the product of dispersivity x[L] and pore-water velocity V (LT1), i.e. D ¼ x v; q is
the water flux [LT1]; 4m and 4im are reaction terms in theAuthor's personal copy
water research 44 (2010) 1050–1061

Fecal coliform (C C0 )1.0e-3
-1

A

Obs.
Pred.

8.0e-46.0e-4

when a single-porosity flow model is used with MIM since inthis case there is no air-phase present in the immobile zone.
When a single-porosity flow model is used with MIM, masstransfer between the two regions is through diffusion

4.0e-4

Gs ¼ as ðCm  Cim Þ

2.0e-40.0-1Fecal coliform (C C0 )-2.0e-4-1

Fecal coliform (C C0 )

0

10

20

30

40

50

C

1.0
0.4

0.2

0

20

40

60

80

100

120

D

3e-4
2e-4
1e-4
1e-4
5e-5
0
0

20

40

60

80

100

120

Fecal coliform (C C0 )1.6e-4F

1.4e-41.2e-41.0e-48.0e-56.0e-54.0e-52.0e-5
0.0
0

10

20

30

40

50

Days after DSE application

Fig. 7 – Fecal coliform concentrations observed andpredicted from the MIM transport model and the singleporosity flow model, in drainage from lysimeters A, C, D
and F (spray irrigation).

mobile and immobile domains [ML3T1], respectively(4m ¼ 4im ¼ 0 for Br), and Gs is the term for mass transfer
between the two regions [ML3T1].
The reaction terms (4m and 4im), which represent any
mechanisms responsible for removal of bacteria from theliquid phase (i.e. inactivation, or sorption to the solid phase or
air-water interfaces), are defined as follows:
fm ¼ ka qm Cm þ kaa qm Cm þ mqm Cm ¼ ðka þ kaa þ mÞqm Cm ¼ lqm Cm(12)
fim ¼ ka qim Cim þ kaa qim Cim þ mqim Cim ¼ ðka þ kaa þ mÞqim Cim¼ lqim Cim

ð13Þ

Here ka is the attachment coefficient to the solid phase (T1),
kaa is the attachment coefficient to the air-water interfaces(T1), m is the inactivation rate (T1), and l ¼ ka þ kaa þ m is the
lumped total removal rate (T1). Note that kaa is equal to zero

(14)

where as is the first-order solute mass transfer coefficient(T1). When a dual-porosity flow model is used with MIM,
contaminant is assumed to be able to move in and out of thematrix domain by diffusive and convective fluxes, and theterm for mass transfer between the two regions, Gs, is written
as (Šimůnek and van Genuchten, 2008):
Gs ¼ as ðCm  Cim Þ þ Gw C

0.0

3e-4

-1

1057

(15)

where C* is equal to Cm for Gw > 0 and Cim for Gw < 0.
The observed patterns of Br and bacterial transport weredistinctively different (Figs. 6 and 7), due to size-exclusion inbacterial transport. Therefore transport parameters forbacteria were estimated independently from those for Br.
Values of x, qim and as were optimized for Br data, while x, qim, land as were optimized for bacterial data. The values ofhydraulic property parameters optimized from water flow
were fixed during simulations of Br and bacteria. As theirconcentrations for model calibration were obtained from only
one sampling location (i.e. the column base), we assumed thatparameter values were the same for all soil layers. We also
assumed when using the MIM model that while the immobile
water content, qim, could differ for different soil layers, theporosity of the mobile domain was the same in all layers.
HYDRUS-1D code was modified in this study to include theseassumptions.
For lysimeters B and E, the dual-porosity flow model wasused and properties of mobile and immobile water domains
have been already defined from water content simulations.
Thus we no longer needed to optimize qim for Br and bacterialtransport. Therefore values of x and as were optimized for Br
data, while x, l and as were optimized for bacterial data.

3.

Results and discussion

3.1.

Water contents

3.1.1.

Pre-DSE application period

Fig. 2 shows that soil water contents predicted from the singleporosity flow model fit well with those observed during the preDSE application period for all lysimeters. This is reflected in thevery low values of SSRs (the Sum of Squared Residuals) (0.02–
0.17) and very high values of R2 (0.84–0.99) as shown in Table 2.
Although 15 parameters were optimized simultaneously, there
was little correlation between parameters, indicating they are
relatively independent. The dual-porosity model was also
applied for the pre-DSE application period but the results (notshown) were much worse. During this period, no irrigation wasapplied and lysimeters received only natural precipitation.
Minor preferential flow is expected without irrigation (Chenet al., 2002). These results suggest that the effect of preferential
flow was not significant for water movement during the preDSE application period and a single-porosity flow model isAuthor's personal copy1058

water research 44 (2010) 1050–1061

flood irrigation, much of the water could go through preferential flow, bypassing the matrix and draining deeper into the
profile. Hotter and drier weather would promote the formation of macropores (Brown et al., 2000; Jarvis, 2007). This isespecially true for soils with higher topsoil clay content (Jiang
et al., 2008; Lin and Zhou, 2008) as clay may shrink and crack.
All these facts would reduce the effective soil water holding
capacity.
The single-porosity flow model was also successful for
simulating water contents for the post-DSE application forlysimeters A, C, D, and F. Although it gave greater SSRs values(0.20–0.44) and lower R2 values (0.57–0.80) than those from thepre-DSE application period (SSRs ¼ 0.02–0.17, R2 ¼ 0.84–0.99),
the post-DSE application period (124 days) was much longer
than the pre-DSE application period (47 days). Statistically the
goodness-of-fit for both data series has the same significance(P ¼ 0.000).
The calibrated model for the post-DSE application period
yielded generally decreased values of qs but increased valuesof a and Ks (Table 2) compared to those calibrated from thepre-DSE application period. This further supports the above
evaluation on the change of soil hydraulic properties. As more
macropores are developed, soil conductivity increases (higherKs), more air-entry (greater a) is allowed, and less water is held
(smaller qs).

0.12
0.10

B

Obs.
Pred.

-1

Br (C C0 )

0.08
0.06
0.04
0.02
0.00
0.16

E

0.14
0.12

-1

Br (C C0 )

0.10
0.08
0.06
0.04
0.02
0.00
0

20

40

60

80

100

120

Days after DSE application

3.1.2.

Post-DSE application period

Using parameter values obtained from the pre-DSE application
period for the post-DSE application period, the predicted watercontents (without any calibrations) from the single-porosity
model are all significantly higher than those observed.
An example is given in Fig. 3 for lysimeter A. The discrepancies
indicate that soil hydraulic properties had significantly
changed and soils did not hold as much water as they did in the
pre-DSE application period. If soil hydraulic properties had
remained the same, irrigation should have resulted in a muchwetter soil profile, which is opposite to what was observed.
Compared to the pre-DSE application period, the majorchanges in the post-DSE application period (refer to Fig. 1)
were fortnightly irrigation events (100 mm for flood irrigationand 50 mm for spray irrigation) and drier and hotter climate
conditions (spring and summer) than the pre-DSE application
period (winter). Irrigation is more likely to generate preferential flow compared to rainfall (Chen et al., 2002). Whenapplied in dry conditions (especially on the surface) andbefore the soil can wet up and swell, irrigation makes soil
prone to have more preferential flow due to prior shrinkage
and water repellence (Edwards et al., 1993; Jarvis, 2007; Lin andZhou, 2008). Thus when irrigation water is applied, especially

Obs.
Pred.

2.0e-3

Fecal coliform (C C0-1)

adequate for describing water movement in the absence of
irrigation for the soils investigated.

2.5e-3B

1.5e-31.0e-35.0e-40.0
1.6e-3E

1.4e-3

Fecal coliform (C C0-1)

Fig. 8 – Br concentrations observed and predicted from the
MIM transport model and the dual-porosity flow model, indrainage from lysimeters B and E (floor irrigation).

1.2e-31.0e-38.0e-46.0e-44.0e-42.0e-4
0.0
0

20

40

60

80

100

120

Days after DSE application

Fig. 9 – Fecal coliform concentrations observed andpredicted from the MIM transport model and the dualporosity flow model, in drainage from lysimeters B and E(floor irrigation).

Author's personal copy1059

water research 44 (2010) 1050–1061

Table 4 – Transport parameters optimized from observed Br and bacterial concentrations using the mobile-immobiletransport model and single-porosity flow model for lysimeters A, C, D and F.
Lysimeter

A
C
D
F

Br

Fecal coliformx(m)qim (v v )as (d )

SSRs

R

x(m)

qim (v v )l (day1)as (d1)

SSRs

R2

0.016
0.547
0.179
0.031

0.245
0.335
0.316
0.254

0.026
0.037
0.014
0.403

0.437
0.501
0.537
0.324

0.581
0.559
0.468
0.892

0.006
0.254
0.005
0.042

0.181
0.330
0.353
0.100

0.753
0.359
0.300
0.715

0.026
0.013
0.001
0.027

0.470
0.518
0.508
0.567

0.530
0.498
0.493
0.436

1

1

2

The changes in soil hydraulic properties appear to haveoccurred over a period of time during the post-DSE application
period. This is reflected in the fact that predicted soil watercontents do not fit observations equally well for the post-DSE
application data and are commonly under-predicted in theinitial stage (Figs. 4 and 5). If the soil hydraulic properties wereconstant, the model should fit all data equally well as it does
for the pre-DSE application data (Fig. 2).
While the single-porosity flow model failed to simulate
water contents for the post-DSE application period for lysimeters B and E, application of the dual-porosity model to thesedata was successful. Unlike the other four lysimeters underspray irrigation, these two were under flood irrigation, thus
preferential flow is expected to be more significant. A dye
tracer experiment and destructive sampling at the end of the
field experiments (Jiang et al., 2008) revealed that there weremore macropores in these two lysimeters than in otherlysimeters. The dual-porosity model also yielded decreased
values of qs and increased values of Ks for the post-DSE
application period (Table 3) than those predicted from the preDSE application period (Table 2). This is in consistent with thefindings for other lysimeters, as discussed earlier. For Ks it isappropriate to compare results from both flow models as Ks
applies to the macropore domain only. However, it is difficult
to compare other parameters derived from two different
models. Unfortunately, the dual-porosity model failed tosimulate water contents obtained from pre-DSE application
for lysimeters B and E. Thus apart from Ks, we are unable todirectly compare model results for assessing the changes insoil hydraulic properties of lysimeters B and E during the postDSE application period. Nevertheless the success vs. failure inusing different flow models for the pre- and post-DSE application periods suggests a significant change in soil hydraulicproperties.

3.2.

Br and fecal coliform transport

MIM-predicted Br and bacterial concentrations match wellwith most of those observed (Figs. 6–9). However the predicted1

peak bacterial concentrations for all lysimeters are significantly lower than those observed. This is because an equal
weight was assigned in the model for all observed data during
optimization. As there were many more data with low
concentrations, modeling results are very much affected by
the data with lower concentrations.
Bacterial breakthrough was either earlier than or at thesame time as Br concentrations. The time for the peak Br andbacterial concentrations is 20 min (same time as DSE application) for lysimeter C and 6 days (same time as 1st irrigation)
for lysimeters A and D. For lysimeters B, E and F, the time forthe peak bacterial concentration is 6–7 days while for Br it is
18–20 days, about three times slower. This suggests thatbacteria were not retarded in the soils and our assumption of
negligible instantaneous sorption made in MIM is reasonable.
Instead, bacterial transport could be faster than Br as a resultof pore-size exclusion.
We believe that the rapid leaching of Br and bacteria inlysimeter C immediately (20 min) after DSE application andbefore any water irrigation is through a preferential flow paththat penetrated most or all of the lysimeter depth. At thattime, bacterial concentration in the leachate was 85% of the
DSE applied (see C/Co in Fig. 7). This indicates that applied DSE
leached through preferential flow paths with little or no
dilution. The topsoil of this lysimeter had a higher clay
content, which readily forms preferential flow paths byshrinking and cracking. The greater risk of microbial leaching
through structured clayey soils is also highlighted in Panget al. (2008) and Pang (2009) who analyzed studies carried out
both in New Zealand and overseas on microbial transport
through intact soils.
For all lysimeters, bacterial leaching occurred at the start of
the first irrigation event and their concentrations dropped
rapidly as irrigation continued. In contrast, Br leaching had
a relatively slower response to irrigation events. The MIM
simulated well these different leaching patterns. The explanation for these distinct patterns is that the initial part of theflow event is dominated by drainage from large pores(Hallberg et al. 1986), which were followed by both Br andTable 5 – Transport parameters optimized from observed Br and bacterial concentrations using the mobile-immobiletransport model and dual-porosity flow model for lysimeters B and E.
Lysimeter

B
E

Br

Fecal coliform

x(m)

as

SSRs

R2

x(m)

l (day1)

as

SSRs

R2

0.049
0.017

0.051
0.077

0.403
0.421

0.561
0.680

0.151
4.88E-04

0.589
0.706

0.040
0.729

0.288
0.908

0.713
0.187

Author's personal copy
1060

water research 44 (2010) 1050–1061

bacteria. As drainage slows down, much of the water being
delivered to the lysimeter base is from smaller pores, anddominated by Br as the bacteria would be excluded or
removed by the small pores. Compared to other lysimeters, Br
leaching from lysimeters B and E responded well to irrigation
events as more macropores were developed, and MIM simulated well these individual responses (Fig. 8).
Unlike in other lysimeters in which bacteria concentrations declined over time, there were abnormal rises of bacterial concentration in lysimeter E with irrigation events (Fig. 9).
The reason is unclear but the possible explanations are (a) an
external fecal source (e.g. from bird faeces), (b) release ofentrapped bacteria due to soil structure change, and (c) undersuitable nutrient state, bacteria growth during summer. The
MIM is unable to simulate these abnormal increases ofbacterial concentrations. As a result, the goodness of fit
parameters (SSRs and R2) for modeling of bacterial transport inlysimeter E is poor (Table 5). In trying to simulate theseabnormal rises in bacterial concentrations, MIM produced
a very high and probably unrealistic value of mass exchange
coefficient (as).
Similar to the situation for water contents mentioned in
Section 3.1.2, there are some under-predictions of Br andbacterial concentrations (Figs. 6–9) for the initial-stage simulations. Again we believe this discrepancy is a result ofchanges in soil structure, while the MIM assumes constantmobile and immobile regions for the whole simulated process(Mitchell and van Genuchten, 1991).
Compared to Br, values of dipersivity (x), immobile watercontent (qim), and mass exchange coefficient (as) determinedfrom bacterial data are generally smaller (Table 4), indicatinga significantly narrow pore-network and less mass exchange
between mobile and immobile water zones in bacterialtransport. This suggests that bacteria took short-circuit paths
by preferential flow, consistent with pore-size exclusion
theory (Ginn, 2002; Sinton et al., 2005).

4.

Conclusions and recommendations

Our modeling results suggest that the single-porosity flow
model is sufficient and adequate for simulating water movement through structured sandy loam soils under natural
climatic conditions and spray irrigation (50 mm fortnightly inthis study). However, a dual-porosity flow model is more
appropriate for simulating water movement under flood irrigation (100 mm fortnightly in this study) as the effect ofpreferential flow paths becomes more significant under flood
irrigation.
This study has demonstrated that soil structure and thus
soil hydraulic properties can change with irrigation (especially
under high hydraulic loading like flood irrigation), and withclimate conditions. Thus parameter values that best describewater flow and contaminant transport under one set of
conditions may not be applicable for another set of conditions.
This dynamic change of soil structure needs to be considered.
Compared to soils with less topsoil clay and soils underspray irrigation, soils with higher topsoil clay content andsoils under flood irrigation have a greater risk of bacterial
leaching through preferential flow paths, potentially

impairing groundwater quality. Management strategies toreduce bacterial leaching from land application of effluent
must aim to minimize the effect of preferential flow byadopting appropriate irrigation schemes, e.g. using spray
instead of flood irrigation and decreasing DSE application
depth but increasing irrigation frequency. Special attention isneeded to use of a low irrigation rate at the beginning of eachirrigation event (particularly if the soils are very dry) to allow
sufficient time for soil to wet up and swell for more effective
water retention so that bacterial leaching can be minimized.

Acknowledgements
We thank the New Zealand Foundation for Research, Scienceand Technology and the Institute of Environmental Science &Research for funding this research with Scholarships, and theCentre for Soil and Environmental Quality for providing
lysimeters and technical support.

references

Aislabie, J., Smith, J., Fraser, R., McLeod, M., 2001. Leaching ofbacterial indicators of faecal contamination through four New
Zealand soils. Australian Journal of Soil Research 39 (6),1397–1406.
American Public Health Association, American Water Works
Association and Water Environment Federation, 1998.
Standard Methods for the Examination of Water and Waste
Water. American Public Health Association, Washington, DC.
Bradford, S.A., Šimůnek, J., Walker, S.L., 2006. Transport andstraining of E. coli O157:H7 in saturated porous media. Water
Resources Research 42, W12S12,. doi:10.1029/2005WR004805.
Brown, C.D., Hollis, J.M., Bettinson, R.J., Walker, A., 2000. Leachingof pesticides and a bromide tracer through lysimeters fromfive contrasting soils. Pest Management Science 56, 83–93.
Chen, C., Roseberg, R.J., Selker, J.S., 2002. Using microsprinkler
irrigation to reduce leaching in a shrink/swell clay soil.
Agricultural Water Management 54, 159–171.
Close, M., Dann, R., Ball, A., Pirie, R., Savill, M., Smith, Z., 2008.
Microbial groundwater quality and its health implications fora border-strip irrigated dairy farm catchment, South Island,
New Zealand. Journal of Water and Health 6 (1), 83–98.
Edwards, W.M., Shipitalo, M.J., Owens, L.B., Dick, W.A., 1993.
Factors affecting preferential flow of water and atrazine
through earth-worm burrows under continuous no-till corn.
Journal of Environmental Quality 22, 453–457.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field
Water Use and Crop Yield. Wageningen, The Netherlands.
Flury, M., Yates, M.V., Jury, W.A., 1999. Numerical analysis of theeffect of the lower boundary condition on solute transport inlysimeters. Soil Science Society of America Journal 63, 1493–1499.
Gargiulo, G., Bradford, S.A., Simunek, J., Ustohal, P.,
Vereecken, H., Klumpp, E., 2008. Bacteria transport anddeposition on under unsaturated flow conditions: the role ofwater content and bacteria surface hydrophobicity. Vadose
Zone Journal 7, 406–419.
Gerke, H.H., van Genuchten, M.T., 1993. A dual-porosity model for
simulating the preferential movement of water and solutes instructured porous media. Water Resources Research 29 (2),
305–320.
Gerke, H.H., van Genuchten, M.T., 1996. Macroscopic
representation of structural geometry for simulating water

Author's personal copy
water research 44 (2010) 1050–1061

and solute movement in dual-porosity media. Advances in
Water Resources 19 (6), 343–351.
Ginn, T.R., 2002. A travel time approach to exclusion on transportin porous media. Water Resources Research 38 (4), 12.01–12.12.
Guber, A.K., Shelton, D.R., Pachepsky, Y.A., 2005. Transport andretention of manure-borne coliforms in soil. Vadose Zone
Journal 4 (3), 828–837.
Hallberg, G.G., Baker, J.L., Randall, G.W., 1986. Utility of Tile-line
Effluent Studies to Evaluate the Impact of Agricultural Practiceon Ground Water. National Water Well Assoc., Dublin, OH/
Omaha, NE, pp. 298–326.
Hay, R.K.M., Porter, J.R., 2006. The Physiology of Crop Yield Ames.
Iowa: Blackwell Pub., Oxford, UK.
Hekman, W.E., Heijnen, C.E., Burgers, S.L.G.E., van Veen, J.A., vanElsas, J.D., 1995. Transport of bacterial inoculants throughintact cores of two different soils as affected by waterpercolation and the. FEMS Microbiology Ecology 16 (2),
143–157.
Hijnen, W., Brouwer-Hanzens, A., Charles, K., Medema, G., 2005.
Transport of MS2 phage, E. coli, Clostridium perfringens,
Cryptosporidium parvum and Giardia intestinalis in a gravel anda sandy soil. Environmental Science and Technology 39,
7860–7868.
Jarvis, N.J., 2007. A review of non-equilibrium water flow andsolute transport in soil macropores: principles, controlling
factors and consequences for water quality. European Journalof Soil Science 58, 523–546.
Jiang, S., Buchan, G.D., Noonan, M.J., Smith, N., Pang, L., Close, M.,
2008. Bacterial leaching from dairy shed effluent applied toa fine sandy loam under irrigated pasture. Australian Journalof Soil Research 46, 552–564.
Lin, H., Zhou, X., 2008. Evidence of subsurface preferential flowusing soil hydrologic monitoring in the Shale Hills catchment.
European Journal of Soil Science 59, 34–49.
Ludlow, M.M., 1985. Photosynthesis and dry matter production in
C3 and C4 pasture plants, with special emphasis on tropical C3
legumes and C4 grasses. Australian Journal of Plant Physiology
12, 557–572.
McGechan, M.B., Vinten, A.J.A., 2003. Simulation of transport
through soil of E. coli derived from livestock slurry using the
MACRO model. Soil Use and Management 19 (4), 321–330.
McLaren, R.G., Cameron, K.C., 1996. Soil Science: Sustainable
Production and Environmental Protection. Oxford University
Press, Auckland, NZ.
McLeod, M., Aislabie, J., Smith, J., Fraser, R.H., Robert, A.,
Taylor, M.D., 2001. Viral and chemical tracer movement
through contrasting soil. Journal of Environmental Quality 30(6), 2134–2140.
McLeod, M., Aislabie, J., Ryburn, J., McGill, A., Taylor, M.D., 2003.
Microbial and chemical tracer movement through two
Southland soils, New Zealand. Australian Journal of Soil
Research 41 (6), 1163–1169.
McLeod, M., Aislabie, J., Ryburn, J., McGill, A., 2004. Microbial andchemical tracer movement through Granular, Ultic and Recent
Soils. New Zealand Journal of Agricultural Research 47 (4),557–563.
McMurry, S.W., Coyne, M.S., Perfect, E., 1998. Fecal coliform
transport through intact soil blocks amended with poultry
manure. Journal of Environmental Quality 27 (1), 86.
Mitchell, A.R. and van Genuchten, M.T., 1991. Preferential Flow.
Gish, T.J. and Shirmohammadi, A. (eds), pp. 278–287, Am. Soc.
Agric. Eng., St. Joseph, MI.
Mualem, Y., 1976. A new model for prediction of the hydraulicconductivity of unsaturated porous media. Water Resources
Research 12 (3), 513–522.
Pang, L., Mcleod, M., Aislabie, J., Simunek, J., Close, M., Hector, R.,
2008. Modeling transport of microbes in ten undisturbed soils
under effluent irrigation. Vadose Zone Journal 7 (1), 97–111.

1061

Pang, L., 2009. Microbial removal rates in subsurface media
estimated from published studies of field experiments and large
intact soil cores. Journal of Environmental Quality 38, 1531–1559.
Ritchie, J.T., 1972. Model for predicting evaporation from a row
crop with incomplete cover. Water Resources Research 8 (5),1204–1213.
Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. Rosetta:
a computer program for estimating soil hydraulic parameterswith hierarchical pedotransfer functions. Journal of Hydrology
251, 163–176.
Shelton, D.R., Pachepsky, Y.A., Sadeghi, A.M., Stout, W.L.,
Karns, J.S., Gburek, W.J., 2003. Release rates of manure-borne
coliform bacteria from data on leaching through stony soil.
Vadose Zone Journal 2 (1), 34–39.
Šimůnek, J., Jarvis, N.J., van Genuchten, M.T., Gärdenäs, A., 2003.
Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadosezone. Journal of Hydrology 272, 14–35.
Šimůnek, J., Sejna, M., Van Genuchten, M.T., 1998. The HYDRUS1D Software Package for Simulating the One-Dimensional
Movement of Water, Heat, and Multiple Solutes in Variablysaturated Media: Version 2.0. IGWMC-TPS-70. International
Ground Water Modeling Center. Colorado School of Mines,
Golden.
Šimůnek, J., van Genuchten, M.T., 2006. The Handbook of
Groundwater Engineering. In: Delleur, J.W. (Ed.). CRC Press
Boca Raton, FL, pp. 22.21–22.46.
Šimůnek, J., van Genuchten, M.T., 2008. Modeling nonequilibrium
flow and transport processes using HYDRUS. Vadose Zone
Journal 7 (2), 782–797.
Šimůnek, J., Van Genuchten, M.T., Sejna, M., 2005. The HYDRUS1D Software Package for Simulating the One-dimensional
Movement of Water, Heat and Multiple Solutes in Variablysaturated Media. Department of Environmental Sciences,
University of California Riverside, Riverside, California.
Šimůnek, J., Wendroth, O., Wypler, N., van Genuchten, M.T., 2001.
Non-equilibrium water flow characterized by means of
upward infiltration experiments. European Journal of Soil
Science 52, 13–24.
Sinton, L.W., Braithwaite, R.R., Hall, C.H., Pang, L., Close, M.E.,
Noonan, M.J., 2005. Tracing the movement of irrigated effluent
into an alluvial gravel aquifer. Water Air and Soil Pollution
166, 287–301.
Smith, M.S., Thomas, G.W., White, R.E., Ritonga, D., 1985.
Transport of Escherichia coli through intact and disturbed
soil columns. Journal of Environmental Quality 14 (1),87–91.
Torkzaban, S., Tazehkand, S.S., Walker, S.L., Bradford, S.A.,
2008. Transport and fate of bacteria in porous media:
Coupled effects of chemical conditions and pore spacegeometry. Water Resources Research 44. doi:10.1029/
2007WR006541.
van Genuchten, M.T., 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Science
Society of America Journal 44, 892–898.
van Genuchten, M.T., Leij, F.J., Yates, S.R., 1991. Version 1.0. The
RETC Code for Quantifying the Hydraulic Functions of
Unsaturated Soils. U.S. Salinity Laboratory USDA, ARS,
Riverside, CA.
van Genuchten, M.T., Šimůnek, J., 2004. Unsaturated-Zone
modeling; Progress, Challenges and applications, Wageningen
UR frontis series. In: Feddes, R.A., de Rooij, G.H., van Dam, J.C.
(Eds.). Kluwer Academic Publishers, Dordrecht, The
Netherlands, pp. 37–69. x–xi.
van Genuchten, M.T., Wagenet, R.J., 1989. Two-site/two-region
models for pesticide transport and degradation: theoreticaldevelopment and analytical solutions. Soil Science Society ofAmerica Journal 53 (5), 1303–1310.