WATER RESOURCES RESEARCH, VOL. 34, NO. 5, PAGES 1061–1076, MAY 1998

Preferential transport of nitrate to a tile drain
in an intermittent-flood-irrigated field:
Model development and experimental evaluation
B. P. Mohanty
U.S. Salinity Laboratory, Riverside, California

R. S. Bowman and J. M. H. Hendrickx
Department of Earth and Environmental Sciences, New Mexico Institute of Mining and Technology, Socorro

J. Simunek and M. T. van Genuchten
U.S. Salinity Laboratory, Riverside, California

Abstract. A comprehensive field experiment was conducted near Las Nutrias, New
Mexico, to study field-scale flow and transport in the vadose zone. The field data wereanalyzed in terms of a two-dimensional numerical model based on the Richards equationfor variably saturated water flow, convection-dispersion equations with first-order chemical
decay chains for solute transport, and bimodal piecewise-continuous unsaturated hydraulicfunctions to account for preferential flow of water and nitrate-nitrogen (NO3-N; loosely
used as NO2
3 ) following flood irrigation events at the experimental site. The model wastested against measured NO2
3 flux concentrations in a subsurface tile drain, several
monitoring wells and nested piezometers, and against resident NO2
3 concentrations in thesoil profile (obtained at 52 spatial locations and four depths along a transect). NO2
3
transport at the field site could be described better with the bimodal hydraulic functions
than using the conventional approach with unimodal van Genuchten-Mualem typehydraulic functions. Average resident nitrate concentrations measured across the soil
profile were predicted reasonably well. However, NO2
3 flux concentrations in thesubsurface tile drain and piezometers at the field site were occasionally underestimated oroverestimated depending upon the irrigation sequence in three field benches, probably
reflecting unrepresented three-dimensional regional flow/transport processes. Limiting thecapture zone to a region closer to the tile drain did lead to a better match with observed
sharp increases and decreases in predicted NO2
3 flux concentrations during the irrigation
events. On the basis of this result we inferred that the preferential flow intercepted by thetile drain was generated in close proximity of the drain and essentially oriented vertically.
In summary, our study suggests that irrigation scheduling in adjacent field plots, drainage
design (e.g., spacing between tiles, drain depth, drain diameter) and effectiveness (e.g.,drain blockage), preferential flow in (horizontal) surface-opened shallow cracks and(vertical) macropores, and transient regional groundwater flow can add significant
uncertainty to the predictions of (local-scale) flow and transport to a tile drain.

1.

Introduction

Field-scale transport of agricultural chemicals from the soil
surface to and into groundwater through the vadose zone isincreasingly recognized as an extremely complicated process.
Available computer models generally lack one or several of themost important physical, chemical, or biological processes operative at the field scale. One such important field-scale process is the preferential flow of water and dissolved constituents
through a complex network of soil macropores, worm holes,
root channels, and structural cracks, especially during and after
the application of large amounts of water. Furthermore, processes occurring at larger (e.g., regional) scales may add significant uncertainty to the predictions based on local-scale
Copyright 1998 by the American Geophysical Union.
Paper number 98WR00294.
0148-0227/98/98WR-00294$09.00

physical, chemical, and biological process information. Deterministic identification and representation of the most relevant
flow and transport processes under site-specific soil, crop, water, and chemical management conditions and for the prevailing atmospheric, topographic, and hydrologic situations should
optimally be done through well-designed field experiments andmodeling in an interactive manner [Tod and Grismer, 1996].
Site-specific findings should be assimilated later to deduce
more general rules.
Tile drains are popularly used in many agricultural fields fordraining excess water from the vadose zone to sustain optimum
soil water contents for crop production and/or maintaining low
levels of soil salinity, sodicity, or other contaminants. Tile
drains have been used recently to advantage for studying fieldscale transport behavior [Hallberg et al., 1986]. In different
studies [Shalit and Steenhuis, 1996; Singh et al., 1996; Munster etal., 1996, 1994; Booltink, 1995; Kladivko et al., 1991; Czapar and

1061

1062
MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Plate 1. Areal view of the Las Nutrias field site. Numbers along the bench boundaries indicate monitoring wells andpiezometers (asterisks). Collector drain along the west boundary of the west bench, piezometer nests (32 and 34), andmonitoring well (33) are outside the snap shot boundary.

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Kanwar, 1991; Richard and Steenhuis, 1988; Bowman and Rice,
1986; White et al., 1983], researchers have interpreted sharp
increases in tile drain outflow and its chemical concentrations,
during and following heavy irrigation or rainfall, as being indicative of field-scale preferential flow. Moreover, many studies claim that tile drains are good spatial integrators of fieldscale flow and chemical transport processes (including
preferential flow) and hence that they are excellent tools for
estimating effective flow and transport properties of relatively
large fields. A limited number of two-dimensional numericalmodeling studies [e.g., Khan and Rushton, 1996a; Munster et al.,
1996; Singh et al., 1996; Singh and Kanwar, 1995a, b; Kamra etal., 1991a, b; Sanoja et al., 1990] have tested the use of tile
drains as spatial integrators of matrix-dominated flow in thepresumed absence of preferential flow. Also, Jury [1975a, b]
developed an analytical method to calculate travel times forflow originating at different distances from a tile drain using
the assumption of Darcian flow. Attempts to test models withfield data [e.g., Ragab et al., 1996; Singh and Kanwar, 1995a, b;
Kamra et al., 1991b] often were not successful because of a
failure to account for preferential flow and/or transient regional groundwater flow. These studies all assumed that thecapture zone of the tile drain remained constant with time,
irrespective of the presence or absence of matrix, preferential,
or regional flow regimes at the experimental site. Other studies
by Luxmoore [1991] and Thomas and Barfield [1974] concluded
that tile drains can be unreliable for quantifying chemicaltransport from a field since they intercept only a relatively
small proportion of the total flow in the field. The abovefindings are partly due to our still limited capability to accurately measure and model multidimensional field-scale waterflow and solute transport under variably saturated conditions,especially during conditions of preferential flow.
In a recent paper, Mohanty et al. [1997] identified and modeled preferential flow to a tile drain in a macroporous fieldunder intermittent flood irrigation practice near Las Nutrias,
New Mexico. Preferential flow was modeled by using bimodal
type piecewise-continuous unsaturated hydraulic functions.
Results of the study showed that the bimodal functions didlead to better predictions of the drainage flux in the subsurface
tile drain following any irrigation event as compared to theunimodal van Genuchten-Mualem hydraulic functions. In thispaper we will further evaluate the bimodal hydraulic functions
with regard to predicting field-scale preferential solute transport under intermittent flood irrigation. Specific objectives are(1) to monitor and interpret transient resident NO2
3 concentrations at different spatial locations and depths across thevadose zone, (2) to monitor and analyze NO2
3 flux concentration in a subsurface tile drain, in several monitoring wells, andin nested piezometers at the experimental field, (3) to develop
a two-dimensional numerical model for preferential flow andtransport following flood irrigation using the piecewisecontinuous hydraulic functions, (4) to test model performanceusing the observed NO2
3 flux and resident concentrations, and(5) to investigate the appropriateness of tile drains as field
integrators of flow and transport for the specific preferential
and regional groundwater flow conditions characteristic of ourexperimental site.

2.

Field Site

The experimental field site is located on a 24-ha commercial
agricultural farm, situated on an alluvial flood plain of the Rio

1063

Grande, near Las Nutrias, New Mexico. The basin lies within
the Rio Grande rift, a series of north-south trending half grabens oriented parallel to the Rio Grande. The climate is arid tosemiarid and has a wide range in temperature and rainfall. The
parent material for the entire valley is alluvium. Plate 1 [Mohanty et al., 1997] shows the major soil classes across the 6-ha
center bench of the field, including our study section. In general, the field site consists of silty clay loam sediments withmoderate to poor drainage properties underlain by fine sands,
with no impeding strata to a depth of about 7 m [National
Resource Conservation Service (NRCS), 1992]. Surface horizons, however, contain visible root channels, worm holes, andcracks, thus potentially providing a complex network of preferential flow paths. The field (Plate 1) is relatively flat anddivided into three benches: east (240 m 3 270 m), center (240m 3 270 m), and west (420 m 3 270 m), separated by berms.
To sustain agricultural productivity, the field site was equippedwith surface irrigation and subsurface (tile) drainage systems.
The subsurface tile drainage system consisted of four lateraldrain lines installed at a depth of ;1.2 m. The tile drainsconnected to a main drain line leading to an off-site openchannel drain on the west. Irrigation return flows draining
from the field, as well as periodic buildup of shallow groundwater flowing in from off-site areas, were collected by thedrainage system. A southern section of the center bench drainage system was isolated by installing two manholes along asingle drain line at the east (up) and west (down) ends of thebench (Plate 1 and Figure 1). The manholes allowed regularsampling of water quantity and quality entering and leaving thestudy section of the drain line. The manholes were equippedwith Signett flow measurement systems including data loggers
for monitoring flow rates at 5-min intervals. Automated ISCOt
2700 samplers were used to collect drain water samples inthese manholes. We used the isolated section of the center
bench for collecting soil hydraulic parameters and for monitoring flow and transport in the vadose zone. Horizontal andvertical pressure gradients and chemical concentrations ofgroundwater at different (2, 3, 5, and 7 m) depths in theunderlying aquifer were determined regularly using piezometers and observation wells along the periphery of the field site.
For this purpose, eight nested piezometers at 2-, 3-, 5-, and 7-m
depths and 34 monitoring wells at 2-m depth (Plate 1) were
installed around the field. Three rain gauges were located on
the center bench to monitor the average precipitation input ofindividual rainfall events.
A concrete-lined canal on the north side of the field supplied
irrigation water to all three benches through a number ofcontrol gates using flood irrigation. Irrigation water flowedfrom the north end of the field to the south end in about 4 –7
hours and remained ponded in the field for several hours withsimultaneous infiltration. After the desired amount of water
was delivered to the field, irrigation supply was cutoff, thus
allowing continued variably saturated flow and redistribution
within the soil profile. A 25-cm-diameter circular weir of thetype developed by Samani [1993] was used in the irrigation
canal to measure water input to the center bench during eachirrigation event. The local farmer scheduled the water input
events to different benches depending on the consumptive
needs of his crops at the different growth stages. Table 1 shows
the irrigation dates and amounts, including the bench sequenceduring 1994 and 1995. All three benches normally received
water on the same day, although irrigation in the three benches
was spread occasionally over 2-day periods. Table 2 shows the1064

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Figure 1. Three-dimensional view of the field site including hydrology, soil and water quality monitoring
facilities, and locations of surface and subsurface site characterization studies.

crop rotation for 1994/1995. Details of the measurement andmathematical description of the surface and shallow subsurface hydraulic properties of the field site are given by Mohanty
et al. [1997].
2.1.

Solute Flux Concentration Sampling

Observation wells at 2-m depth and piezometers at 2-, 3-, 5-,
and 7-m depths (Plate 1 and Figure 1) were sampled intermittently for chemical concentrations to determine the spatiotemporal distributions and vertical concentration gradients of solutes in the shallow groundwater zone. Sampling intervals
ranged from a few weeks during summer to a few months

during winter. To avoid measuring stagnant well water, approximately five well volumes were bailed before samples wereactually collected. A total of 544 and 416 water samples werecollected during 1994 and 1995, respectively. ISCO 2700 automatic water samplers were used to collect tile effluents at theeast and west manholes (Plate 1 and Figure 1). During irrigations, water samples were taken every 2 hours for the first 48
hours, then once every 4 hours for the next 48 hours, andsubsequently once every 6 hours for 72 hours. This scheme wasfollowed by sampling at 24-hour intervals until the next irrigation. Water samples collected from the tile drains, monitoring
wells, and piezometers were later analyzed for NO2
3 and other

Table 1. Irrigation Application Dates and Volumes Applied to the Center Bench
Date

Day*

Total Elapsed
Time, hours

Total Volume,L

Total Depth,
mm

Bench
Sequence†

June 8, 1994
June 27, 1994
Aug. 8, 1994
Sept. 28, 1994
April 10, 1995
April 26, 1995

721
740
782
833
1026
1043

4.5
7.0
6.0
7.5
6.0
6.5

4.14E 1 06
4.56E 1 06
5.10E 1 06
5.88E 1 06
6.47E 1 06
5.37E 1 06

61.22
67.43
75.41
86.93
95.76
79.48

E-C-W
E-C-W
E-C-W
W-C-E
E-W-C
E-W-C

Bench sequence indicates the order of water application during a particular irrigation event. Irrigation
events prior to June 1994 were not recorded. Read 4.14E106 as 4.14 3 106.
*Sequential day from the start of the project.
†E, east bench; C, center bench; W, west bench.

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

solutes by using High Performance Liquid Chromatography(HPLC) in the New Mexico Institute of Mining and Technology soil water chemistry laboratory.
Table 3 presents mean NO2
3 flux concentrations for the 34
observation wells located at the experimental field site. No
general trend was found in the mean concentrations in horizontal space. The NO2
3 distributions were somewhat influenced by the irrigation sequence and timing of the three fieldbenches (to be discussed in section 4) and by spatial variability
in the flow/transport process. However, mean NO2
3 flux concentration showed a somewhat decreasing trend with time.
Except for a few hot points having high values (;7 mg/L),
NO2
3 concentrations were generally quite low across the fieldat 2-m depth on all sampling dates. No apparent reason could
be found for the localized hot points. Figure 2 shows thevertical profile in the NO2
3 flux concentration on different
dates measured in four nested piezometers (11, 14, 21, and 24,
Plate 1) located on the four corners of our study section.

1065

Table 2. Crop Rotation Schedule for the Center Bench at
the Las Nutrias Field Site
Crop

Planting Date

Harvest Dates

Winter wheat
Sorghum-Sudan crossAlfalfa

Oct. 10, 1993
May 13, 1994
Sept. 25, 1994

not harvested
July 19, 1994
May 29, 1995
June 2, 1995
July 5, 1995
Aug. 6, 1995
Sept. 6, 1995
Oct. 15, 1995

Except for one outlier, no significant vertical concentration
gradient (between 2- and 7-m depths) was found in the shallow
groundwater zone. Moreover, NO2
3 flux concentration at these
depths remained low and approximately unchanged through-Figure 2. Intermittent NO2
3 flux concentrations observed in the piezometers on the four corners of the study
section in the center bench.

1066

Table 3. NO2
3 Flux Concentration in the Monitoring Wells at the Las Nutrias Field Site in 1994
April
10

April
24

April
29

May
17

May
25

June
1

June
9

July
1

July
19

Aug.
19

Sept.
11

Sept.
24

Oct.
8

Oct.
23

Nov.
8

Nov.
15

354.8
221.6
96.21
377.2
119.9
396.8
279.3
150.7
397.5
316.1
278
245.8
242.1
208.8
173.7
151.6
420.9
183.8
436.6
365.9
328.9
296.3
292.2
261.6
223.5
207.7
438.5
328.1

0.926
0.000
2.438

2.725
0.149
1.084

1.370
0.200
0.130

0.923
0.056

3.801
0.050

0.176
0.065

0.086
0.097
0.939
0.047
0.514

2.660
0.170
0.080
0.060
1.280
0.030
0.030
0.030
0.030
0.000
0.030
0.020
0.030
0.030
0.030
0.020
0.030
0.000
0.030
0.040
0.090
0.040
0.000
0.000
0.050
0.050
0.040
0.000
0.070
0.040
0.040
0.070
0.160

1.560
0.080
0.080
0.000
7.000
0.000
0.000
0.020
0.000
0.000
0.000
0.000
0.000
0.000
0.020
7.000
0.000
7.000
0.690
0.000
0.000
0.040
0.040
0.000
0.030
0.040
0.000
0.000
0.000
0.000
0.030
0.000
0.738

1.650
0.010
0.070
0.490
0.310
0.020
0.410
0.120
0.200
0.040
0.030
0.070
0.050
0.050
0.060
0.050
0.070
0.050
0.050
0.050
0.060
0.080
0.080
0.130
0.090
0.090
0.070
0.080
0.400
0.090
0.100
0.120
0.164

0.46
0.13
0.03
0.03
0.30
0.01
0.04
0.03
0.24
0.02
0.01
0.02
0.01
0.00
0.01
0.00
0.01
0.15
0.51
0.04
0.01
0.01
0.01
0.01
0.03
0.01
0.05
0.01
2.76
0.04
0.01
0.01
0.16

0.100
0.234
0.185
1.484
0.043
0.102
0.125
0.017
0.033
0.017
0.085
0.124
0.134
0.121
0.117

0.181
3.102
1.853
1.088

2.740
0.390
0.450
0.390
1.720
0.150
0.560
0.300
0.040
0.530
0.220
0.250
0.520
0.270
2.180
0.080
0.300
0.550
0.120
0.120
0.120
0.300
1.450
1.840
0.780
0.010
0.140
7.000
1.490
0.040
0.550
0.170
0.805

0.547
0.188
0.146
0.243
0.109
0.129
0.110
0.074
0.138
0.129
0.325
0.116
0.000
0.170
0.000
0.014
0.030
0.148
0.572
0.097
0.007
0.023
0.000
0.038
0.414
0.057
0.105
0.044

526.6
287.8
298.3

1.160
0.190
0.000
7.000
0.350
0.000
0.100
0.070
0.080
0.000
0.000
0.000
0.080
0.000
0.070
0.000
0.140
0.000
0.080
0.000
0.000
0.000
0.080
0.000
0.090
0.070
0.430
0.080
0.180
0.000
0.140
0.070
0.327

0.904
0.048
0.051
0.654
0.023
0.164
0.032
0.042
0.030
0.026
0.065
0.012
0.052
0.113
0.088

0.059
1.454

0.900
0.210
0.070
7.000
0.330
0.090
0.150
0.100
0.090
0.090
0.100
0.100
0.110
0.080
0.080
0.090
0.160
0.090
0.090
0.090
0.070
0.080
0.080
0.070
0.760
0.090
0.090
0.080
0.610
0.100
0.150
0.110
0.385

1.860
0.030
0.050
0.030
3.130
0.030
0.070
1.950
0.070
0.040
0.040
0.040
0.050
0.060
0.000
0.140
0.060
0.060
0.040
0.030

0.140
0.824

6.280
0.210
2.560
0.560
0.180
0.160
0.190
0.220
0.190
0.210
0.210
0.110
0.270
0.210
0.190
0.310
0.420
0.170
0.280
0.290
0.190
0.210
1.280
0.190
1.940
0.260
0.900
0.050
0.643

0.950
0.200
0.080
7.000
0.320
0.060
1.370
7.540
0.070
0.080
0.060
0.060
0.110
0.100
0.120
0.100
0.120
0.090
0.090
0.110
0.090
0.110
0.120
0.120
0.120
0.110
0.580
0.120
0.000
0.120
0.120
0.110
0.636

2.770

0.000
0.000

1.230
0.260
0.160
7.000
2.460
0.070
2.810
0.210
0.180
0.090
0.130
0.110
0.100
0.160
0.130
0.130
0.050
0.090
0.050
0.190
0.080
0.090
0.050
0.140
0.050
0.080
1.720
0.120
0.440
0.090
0.050
0.230
0.586

Easting

Northing

A1
A2
A3
A4
A5
A6
A7
A8
A9
A10
A11
A12
A13
A14
A15
A16
A17
A18
A19
A20
A21
A22
A23
A24
A25
A26
A27
A28
A30
A32
A33
A34
Mean

1583.6
1558.31532.31460.71428.81329
1314.41298.71325.61314.91310.41306.21305.71301.71297.31294.61208.81153.51088.6
1077.81072.41067.51066.81062.31056.71054.21083.21067.9
583.8
631
660.2

0.233
4.621
0.072
0.221
0.027

0.603
0.072
0.117
0.079

Concentrations are in milligrams per liter. Blank cells indicate that the data are missing.

0.320
0.030
0.000
0.000
0.230
4.000
0.000
1.700
0.000
0.040
0.000
0.455

1.680
0.000
0.310
0.000
0.000
0.000
0.780
0.040
0.010
0.200
0.000
0.000
0.000
0.020
0.000
0.000
0.000
0.030
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.201

0.000
0.000
0.000
0.128

0.047
0.072
0.062
0.007
0.041
0.048
0.062
0.158
0.060
0.102
1.156
0.050
0.027
0.022
0.145

0.345
0.111
0.385
0.128
0.117
0.108
0.036
0.109
0.153
0.033
0.065
6.619
0.033
0.031
0.036
0.374

Mean
1.44
0.16
0.40
2.24
1.46
0.07
0.56
0.88
0.13
0.09
0.09
0.10
0.10
0.10
0.21
0.93
0.11
0.53
0.18
0.09
0.08
0.10
0.15
0.19
0.19
0.09
0.64
0.53
1.16
0.07
0.37
0.17

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

April
1

Well

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

1067

out the sampling period (August 1994 to approximately May
1995).
NO2
3 flux concentration in the tile drain water collectedfrom the east and west manholes are shown in Figure 3. Evidently, in both manholes, sharp increases in the drainage waterNO2
3 concentrations occurred within a few hours following all
irrigation events. The concentrations remained high for a few
days, followed by sharp drops (with little tailing) to the relatively low background preirrigation concentration level. In
general, NO2
3 concentration in the tile water became as high as
10 ; 15 mg/L before returning to background levels of about
0.3 ; 0.5 mg/L. The sharp increases and decreases in NO2
3
concentration qualitatively suggest that transport to the tile
drain was mostly convection-dominated following the irrigation events.
2.2.

Solute Resident Concentration Sampling

From late February thru early March 1994 and 1995, soil
samples were collected from the center bench study section
along a north-south transect perpendicular to the tile drain, as
shown in Figure 1. Fifty sampling sites along this transect werearranged 1.5 m apart on both sides of the tile drain. In addition, several soil samples were collected very close to the tile
drain. An ;2-m separation distance was maintained between

Figure 4. Profile soil moisture content across the N-Stransect in February 1994 and February 1995 at the field site.

Figure 3. Tile flux NO2
3 concentrations at the east and west
manholes. Breaks in the line indicate missing data because ofhuman or instrumental factors.

the 1994 and 1995 sampling transects. This sampling scheme
was used to investigate possible spatial patterns (trends) in soil
attributes (e.g., resident NO2
3 concentration) created as a result of flow to the tile drain. Soil samples were collected in
30-cm increments from the soil surface to a depth of 120 cm
using Veihmeyert soil sampling tubes. The soil samples werelater analyzed for gravimetric soil water content, soil waterNO2
3 concentration (by using HPLC), soil bulk density, andsoil texture. Note that the same sites (within 1-m-diameterarea) were previously also used for determining water retention and hydraulic conductivity functions using in situ andlaboratory methods [Mohanty et al., 1997].
Figures 4 and 5 show spatial distributions of the soil watercontent and the resident NO32 concentration at different
depths across the north-south transect for February 1994 and1995, respectively. Except for a narrow sharp dip in concentration near the tile drain in 1995, no plot-scale tile flow induced trend was found in the resident NO2
3 concentration at
all four depths. This feature may indicate in part that theradial/vertical capture zone for the tile drain was quite narrow.
Moreover, geostatistical analyses [Journel and Huijbregts, 1978]of the spatial data resulted in white noise dominated semivariograms (Figure 6) indicating no spatial structure. Table 4
gives descriptive statistics including mean, range, standard deviation, and W statistics [Shapiro and Wilk, 1965] of the 1994
and 1995 resident nitrate concentration and soil water content
data. The concentration and soil water content data were bothfound normally distributed at all depths and times. The high1068

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Figure 5. Profile resident NO32 concentrations across the
N-S transect in February 1994 and February 1995 at the fieldsite.
resident nitrate concentrations near the soil surface (0 –30 cm)
in 1994 were partly due to failure of winter wheat in thepreceding season, which left large amounts of unused urea
fertilizer in the soil. Furthermore, freezing of the surface soil
during the winter months of 1993–1994 may have helped topreserve the NO2
3 until spring 1994. A statistical comparison
showed significant differences between resident NO2
3 concen-

Figure 6. Semivariograms of profile resident NO2
3 concentrations in 1994 and 1995.
trations in 1994 and 1995, indicating its depletion from the soil
profile over several crop growing seasons during the 1-year
period. Depletion of NO2
3 was likely due to several simultaneous processes including convective-dispersive transport togroundwater, immobilization to biomass (i.e., organic nitrogen), denitrification to nitrogen gas, and uptake by plants.
Total nitrate depletion was maximum near the soil surface (at0 –30 cm) and decreased gradually with depth. Low and statis-

Table 4. Descriptive Statistics of Resident NO2
3 Concentration and Soil Moisture Content at the Field Site in 1994 and
1995
1994
Statistics

0 –30 cm

30 – 60 cm

1995

60 –90 cm

90 –120 cm

0 –30 cm2
3

30 – 60 cm

60 –90 cm

90 –120 cm

NumberMeanStandard deviation
Minimum
Maximum
W statistics*

54
163.7572.23
73.03
470.82
0.775

54
45.28
23.10
0.00
100.87
0.97

Resident NO Concentration, mg/L
54
54
54
17.74
23.14
29.29
12.54
35.25
15.88
2.79
0.17
0.04
53.27
153.29
78.51
0.874
0.664
0.964

54
12.42
10.06
0.31
35.74
0.889

53
3.11
2.98
0.23
13.20
0.795

54
2.01
2.13
0.00
11.03
0.365

NumberMeanStandard deviation
Minimum
Maximum
W statistics*

54
16.61
1.77
12.36
23.15
0.938

54
23.09
3.10
15.41
38.73
0.829

Gravimetric Soil Moisture Content, %
54
54
54
26.89
20.11
15.56
3.71
2.56
1.74
17.31
15.09
12.68
32.78
27.68
22.40
0.936
0.954
0.943

54
22.42
2.74
15.81
27.03
0.957

54
27.59
2.09
20.64
31.77
0.950

54
22.38
2.97
16.04
28.45
0.964

*Normally distributed at a confidence level .0.9.

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

tically similar mean soil water contents during the 1994 and1995 soil sampling events (Table 4) facilitated an unbiased
year-to-year comparison of resident NO2
3 concentrations inthe soil profile.

3.

Model Description

The field experimental data were analyzed in terms of a
two-dimensional numerical model accounting for variably saturated flow, convective-dispersive transport of several nitrogenspecies, and preferential flow/transport using bimodal type hydraulic functions. The model is a modification and extension of
the CHAIN_2D code documented by Simunek and van Genuchten [1994]. This section gives a brief overview of the model.
3.1. Governing Equations for Water Flow and Solute
Transport
Considering two-dimensional isothermal Darcian flow ofwater and assuming that the different pore regions (micropores and macropores) are at hydraulic equilibrium in a variably saturated rigid porous medium, the governing flow equation
is given by the following modified form of the Richards equation:
­u
­
5
­t
­ xi

FS

K K Aij­h
1 K Aiz
­ xj

DG

2S

(1)

where u is the volumetric water content [L 3 L 23 ], h is thepressure head [L], S is a sink term [T 21 ], x i (i 5 1, 2) are
the spatial coordinates [L], t is the time [T], K A
ij are components of a dimensionless anisotropy tensor K A , and K is theunsaturated hydraulic conductivity function [LT 21 ] given by
K~h, x, z! 5 K s~ x, z! K r~h, x, z!

(2)

where K r is the relative hydraulic conductivity and K s is thesaturated hydraulic conductivity [LT 21 ]. The anisotropy tensor K A
ij in (1) is used to account for an anisotropic medium.
The diagonal entries of K A
ij equal 1 and off-diagonal entries
zero for an isotropic medium. If (1) is applied to planar flow ina vertical cross section, x 1 5 x is the horizontal coordinate andx 2 5 z is the vertical coordinate, the latter taken to be positive
upward. Einstein’s summation convention is used in (1). Hencewhen an index appears twice in an algebraic term, this particular term must be summed over all possible values of the index.
Piecewise-continuous hydraulic functions [Mohanty et al.,1997] were found to perform better than the normal unimodal
hydraulic functions when simulating water flow at our experimental field site. For a multimodal pore-size distribution containing p capillary-dominated flow domains, the piecewisecontinuous soil water retention and hydraulic conductivity
functions can be written as

u ~h! 5O

w pu p~h! 5pO F

w p u r,p 1pu s,p 2 u r,p@1 1 ~ a ph! np# mp

G(3)h # h *u

K~h! 5Ou 5 u s~ p51!

h . h *u

(4)

O
p

np21

kp

h # h *K(5)

K np~h! 5 K* 1 K*@exp~h2h*!d 2 1#

h *K , h # 0

(6)

h.0

(7)

K np~h! 5 K* 1 K*@exp~2h*!d 2 1#

where K p is the hydraulic conductivity of the pth capillarydominated flow domain [LT 21 ], K np is the hydraulic conductivity for noncapillary-dominated flow domain [LT 21 ], u s,p isthe saturated water content of the pth capillary-dominatedflow domain [L 3 L 23 ], u r,p is the residual water content of thepth capillary-dominated flow domain [L 3 L 23 ], h is the equilibrium soil water pressure head of the bulk soil across all flow
domains [L], h *K ' h *u 5 h* is the critical or break point soil
water pressure head where flow changes from capillarydominated to noncapillary-dominated flow or vice versa [L],K* is the hydraulic conductivity corresponding to h*[LT 21 ],d is a fitting parameter representing effective macroporosity orother structural features contributing to noncapillary dominated flow [L 21 ], a p , and n p are fitting parameters [vanGenuchten, 1980] for the pth capillary-dominated flow domain[L 21 , dimensionless], w p is the weighting factor of the pthcapillary-dominated flow domain such that ¥ w p 5 1, and 0 #w p # 1 [dimensionless], and k p is the saturated hydraulicconductivity for capillary-dominated flow domain p subjectedto ¥ k p 5 K*[LT 21 ].
Conceptually, the above functions indicate that the van Genuchten-Mualem models (equations (3) and (5)) should be
used only for tensions greater than h* when capillary flow
dominates, in conjunction with other more appropriate models(e.g., equations (4) and (6)) for soil water retention and especially the hydraulic conductivity at tensions less than h* whengravity flow dominates. Note that when p 5 1, the sum typemultimodal van Genuchten-Mualem hydraulic functions(equations (3) and (5)) reduce to the original unimodal van
Genuchten-Mualem functions.
The partial differential equations governing two-dimensional equilibrium transport of solutes involved in a sequential
first-order decay chain during transient water flow in variably
saturated rigid porous medium are taken as [Simunek and van
Genuchten, 1994]­ u c 1 ­ r s 1 ­a v g 1
­
1
1
5
­t
­t
­t
­ xi
2

S

v
u D ij,1

D

­c 1
­
1
­ xj
­ xi

S

g
a v D ij,1

­g 1
­ xj

D

­q ic 1
2 Sc r,1 2 ~ m v,1 1 m 9v,1! u c 1 2 ~ m s,1 1 m 9s,1! r s 1­ xi

2 ~ m g,1 1 m 9g,1!a v g 1 1 g v,1u 1 g s,1r 1 g g,1a v­ u c k ­ r s k ­a v g k
­
1
1
5
­t
­t
­t
­ xi
2

S

v
u D ij,k

D

­c k
­
1
­ xj
­ xi

k51

S

g
a v D ij,k

(8)
­g k
­ xj

D

­q ic k2 Sc r,k 2 ~ m v,k 1 m 9v,k! u c k 2 ~ m s,k 1 m 9s,k! r s k­ xi

2 ~ m g,k 1 m 9g,k!a v g k 1 m 9v,k21u c k21 1 m 9s,k21r s k21
1 m 9g,k21a vg k21 1 g v,ku 1 g s,kr 1 g g,ka v(9)k [ ~2, · · ·, K!

k pK p~h!

p

5~m p 5 1 2 1/n p!

1069

np 2mp 2~1 2 ~ a ph! @1 1 ~ a ph! #@1 1 ~ a ph! np# mp/ 2!

where c, s, and g are solute concentrations in the liquid[ML 23 ], solid [MM 21 ], and gaseous [ML 23 ] phases, respectively; q i is the ith component of the volumetric flux density[LT 21 ]; m v , m s , and m g are first-order rate constants for1070

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Table 5. Fertilization Applications at the Las Nutrias Field Site
Date

Day

Fertilizer

West Bench,kg/haCenter Bench,kg/ha

East Bench,kg/haJune 25, 1993
Oct. 15, 1993
Oct. 25, 1993
March 7, 1994
May 13, 1994
April 3, 1995

374
486
496
628
695
1019

ammonium sulphate
urea
urea
ammonium sulphateurea
N(8)-P(36)-K(4)-S(4)

448
224
zzz
448
224
280

448
zzz
224
zzz
zzz
280

zzz
zzz
zzz
zzz
zzz
280

solutes in the liquid, solid, and gas phases [T 21 ], respectively;
m 9v , m 9s , and m 9g are similar first-order reaction rate constants
providing connections between the individual chain species;g v , g s , and g g are zero-order rate constants for the liquid[ML 23 T 21 ], solid [T 21 ], and gas [ML 23 T 21 ] phases, respectively; r is the soil bulk density [ML 23 ], a v is the aircontent [L 3 L 23 ], S is the sink term in the water flow equation, c r is the concentration of the sink term [ML 23 ], D ijv isthe dispersion coefficient tensor [L 2 T 21 ] for the liquid phase;
D gij is the diffusion coefficient tensor [L 2 T 21 ] for the gas
phase; k represents the kth chain number; and K is the numberof solutes involved in the chain reaction. The indicial notation
used in this paper assumes summations over indices i and j(i,j) but not over index k.
3.2.

Nitrogen Dynamics and Modeling

Novotny and Chesters [1981] summarized the complex behavior and transformation processes of nitrogen in soils. Fixation of atmospheric N, accumulation of ammonium-N, ammonification of organic-N, hydrolysis of urea, nitrification ofammonium-N, denitrification of nitrate-N, fixation of ammonium-N, plant uptake of nitrate-N, and ammonium volatilization constitute the most important transformation process ofnitrogen in the soil-air-water-plant environment. Bacterial
populations, soil temperature, soil water content, soil aeration,and soil pH are some of the factors that govern the transformation rates of nitrogen fertilizers to soluble nitrate-nitrogenand nitrogen gases. While in reality most of these nitrogentransformation processes will occur in a simultaneous and nonlinear fashion, we consider here only the most important forward processes, thereby keeping our model as simple as possible. Consistent with (7) and (8), nitrogen transformations
were simulated by means of an apparent first-order decay
chain in different phases as follows:
g2
N2
1
m
2
2
urea/(NH4)2SO4 3 NH1
4 3 NO2 3 NO3
c1, s1c2, s2 c3
c4 n
❘3
❘3
N2O
where c, s, and g represent concentrations in the liquid, solid,
and gaseous phases, respectively; straight arrows represent thedifferent zero-order (g) and first-order (m, m9) rate reactions,and curved arrows indicate equilibrium distribution coefficients between phases. The approach lumps the processes ofhydrolysis, nitrification, denitrification, volatilization, and plant
uptake into a set of first- and zero-order reactions, while immobilization, mineralization, and biological fixation of nitrogen are considered unimportant during the study period.
Table 5 lists the fertilization schedule at our experimental
field for 1993–1995. Ammonium sulphate and urea were thetwo most commonly used nitrogen fertilizers. Ammonium sulphate contains ;21% nitrogen by weight and is quite readilyoxidized to NO2
3 . Urea, on the other hand, contains 45%
nitrogen in organic form that first must be converted to ammonium before oxidation to nitrate can occur. This process ismuch slower, and hence urea-derived nitrogen will persist in
the soil profile for much longer periods of time. Moreover, due
to damage to the winter wheat crop in the fall of 1993, we
anticipated that most of the urea remained unused at the fieldsite until the spring of 1994.
The above set of flow and transport equations ((1), (8), and(9)) were solved using the Galerkin finite element method[Simunek and van Genuchten, 1994]. We used upstream
weighting [Yeh and Tripathi, 1990] and an appropriate combination of space and time discretizations to adjust the local
Peclet numbers and Courant numbers [Perrochet and Berod,1993] for reducing undesired oscillations or numerical dispersion during convection-dominated transport. Convection was
assumed to dominate transport processes during conditions ofpreferential flow following flood irrigation at the field site. For
the numerical simulations we used a du convergence criterion[Huang et al., 1996] in conjunction with the “mass conservative” modified Picard iteration method [Celia et al., 1990] toenable relatively fast and robust convergence of the numericalsolution.
3.3.

Model Application to the Las Nutrias Field Site

Two-dimensional numerical simulations were carried out fora N-S vertical cross section across the center bench study
section. The cross section was divided into 9350 nodes and9126 quadrilateral elements as shown in Figure 7. On the basisof available surface infiltration and hydraulic property data,
the soil profile was divided into three soil horizons characterized by three sets of field-averaged hydraulic functions. The
three horizons were discretized into a total of 169 layers fornumerical simulation purposes. Finer discretizations were usednear the soil surface, across the horizon interfaces, and aroundthe subsurface tile drain to handle relatively large local fluxes
and pressure gradients in these areas. In addition, we kept our
maximum vertical discretization relatively small (Dz 5 1 cm)
in the upper two horizons (to a depth of 100 cm). This finer
discretization improved the numerical solutions, in particular,when the bimodal hydraulic functions were used. For the simulations we used a time-variant evapotranspiration rate assuming a normalized water uptake distribution across the top 40
cm of the soil profile in accordance with the approaches ofFeddes et al. [1978] and Simunek and van Genuchten [1994].
Empirical functions [Nielsen and MacDonald, 1978] for NO2
3
uptake by different crops were used for the simulation. Measured resident NO2
3 concentrations before the 1994 crop growing season (i.e., February 1994) at different nodal points andMOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

1071

Figure 7. Typical/adapted finite element mesh including system-independent (time-dependent) and systemdependent flux/head-type atmospheric boundary, no-flux bottom and hydrologic boundaries, tile drain boundaries, and root water uptake sink for the two-dimensional flow/transport simulation in the vertical cross sectionof the field (study section) at Las Nutrias, New Mexico. The cross section of unit length (in Y direction) is
assumed to be representative over the 240-m-long study section. Also shown is the step-zoom-in procedure
adapted to study the capture zone of the tile drain.

interpolated water table elevation and NO32 concentrationfrom the groundwater monitoring wells along the field boundary were used to define the initial conditions in the twodimensional flow-transport domain V. In a few instances weused prior or posterior information to interpolate or extrapolate the data to our simulation starting date. Time-dependent
boundary conditions, such as intermittent irrigation/precipitation rates and chemical application rates at atmosphericboundary nodes, groundwater head values along hydrologic
boundary nodes, and no-flow conditions along impermeablebottom boundary nodes, were invoked. We simplified the twohydrologic boundaries (GAD and GBC, Figure 7) by assuming
them to be impermeable (i.e., no regional flow contributions)
during short-term irrigation-event-based simulations. The tile
drain at the center of the study section was treated as a boundary node surrounded by four regular square elements withadjusted hydraulic conductivities using the electric analog approach of Vimoke et al. [1963] and Fipps [1986]. Mathematicaldetails of the initial and boundary conditions for water flow at
the field site are given by Mohanty et al. [1997].
For the long-term year-long simulation, we invoked an initialtime step of one day with a starting date of March 1, 1994,

which coincided with our first soil sampling campaign at theLas Nutrias field site. As initial conditions we assumed theobserved resident concentration for NO2
3 (C1994), and zeroconcentrations for the other nitrogen species for all finite element nodes in the simulation domain V. For the short-term
irrigation-event based simulations, we adopted 0.1-hour initialtime steps and used the predicted soil water pressure heads
and resident solute concentration (c k ( x, z)) levels at prespecified times from the long-term simulation as initial conditions.
Given the specific field-hydrologic set-up and results of ourexperiment, this approach was found to be the best available
option considering our objective of short-term high-frequency
preferential flow predictions in the tile drain following flood
irrigation.
Nitrogen fertilizer concentration of the irrigation/precipitation water was input as a Cauchy boundary condition acrossthe top (soil-atmosphere) boundary, GAB, of the flow-transportdomain:2u D ij­c kDg~k gc k 2 g atm!
n 1 q in ic k 5 q in ic 0,k 1­ xj id~ x, z! [ G AB

k51 ··· K(10)t [ t flux

1072

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Table 6a. Field-Averaged Parameters of Bimodal Piecewise-Continuous and Unimodal van Genuchten-Mualem Hydraulic
Functions for Different Soil Horizons at the Field Site near Las Nutrias, New Mexico
Bimodal Piecewise-Continuous Hydraulic Functions
Depth,cm^ u s,i51 &,cm3 cm23^ u r,i51 &,cm3 cm23^ a &,
cm21

^n&

^K*&,
cm/s

^h*&,
cm

^ d &,
cm21

0–40
40–100100–700

0.475
0.459
0.430

0.110
0.050
0.090

0.015
0.021
0.083

1.60
1.40
2.12

5.55E 2 04
4.16E 2 04
2.73E 2 03

3
3
0

0.92
0.70
n/a

Unimodal van Genuchten-Mualem Hydraulic Functions (High/Low)
Depth,cm^ u s &,
cm3 cm23^ u r &,
cm3 cm23

^a&,
cm21

^n&

^K s &,
cm/s

0–40
40–100100–700

0.484/0.475
0.464/0.459
0.430/0.430

0.080/0.120
0.045/0.050
0.090/0.090

0.004/0.015
0.010/0.021
0.083/0.083

1.44/1.55
1.25/1.38
2.12/2.12

4.08E 2 03/9.50E 2 04
3.50E 2 03/6.80E 2 04
2.73E 2 03/2.73E 2 03

From Mohanty et al. [1997]. Read 5.55E204 as 5.55 3 1024.

where q i n i represents the outward fluid flux, n i is the outward
unit normal vector, c 0,k is the concentration of the kth nitrogen species in the incoming fluid [ML 23 ], D g is the molecular
diffusion coefficient of nitrogen gas [L 2 T 21 ], d is the thickness of a stagnant boundary layer on the soil surface accounting for gas diffusion [L] [Jury et al., 1983], g atm is the nitrogengas concentration in the atmosphere above the stagnant
boundary layer [ML 23 ], and other variables are as definedearlier. For the impermeable bottom boundary, GCD, and thehydrologic boundaries between adjacent tile lines (GBC, GAD),
when (q i n i 5 0), we used a Neumann boundary condition:­c k
u D ijn 50­ xj i~ x, z! [ G CD ~ G BC ~ G AD

(11)
k51 ··· K;t

Since a major focus of this study was on preferential flowthat would occur during or relatively soon after flood irrigation, we limited the short-term simulation time to ;100 hours
Table 6b. Transport Parameters of Different Nitrogen
Chain Products for Different Soil Horizons at the Field Site
Parameter

Chain
Product

Depth,
cm

Values of theParameters

NH1
4
NH1
4
NH1
4
NH1
4
NO2
2
2
NO3
NH1
4
NO2
2
NO2
3
NH1
4
NO2
2
2
NO3
NO2
3

0–40
40–100100–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700
0–700

1300
1200
1400
0.05
0.01
7.1 3 1026
0.0
0.00050.0005
0.0
0.005
0.1
0.0
0.005
0.0
0.0
0.0
0.0
0.0
0.002

r, kg/m3
D L, m
D T, m
D v , m2 h21
D g , m2 h21
m v , h21
m s , h21
m g , h21
m 9v , h21
m 9v , h21
m 9v , h21m 9s , h21m 9s , h21m 9s , h21
m 9g , h21
m 9g , h21
m 9g , h21
g v , h21

from the start of an irrigation event. Also we adopted 0.1-hour(56 min) initial time steps for our numerical simulations toclosely follow the 5-min intervals of flow measurements in thefield. The hydraulic parameters given in Table 6 were measured in the field [Mohanty et al., 1997]. The transport parameters for different solutes and soil types (Table 6) used in the
numerical simulation were estimated from available literature[van Genuchten, 1985].

4.

Simulation Results and Discussion

We present simulation results for four state variables. Figures 8 and 9 show predicted and observed water flow rates andNO2
3 flux concentrations, respectively, in the tile drain for atypical irrigation event at the field site. The flow simulations
were previously described by Mohanty et al. [1997]. In general,
when the bimodal hydraulic functions were used, model predictions and observed data for tile flow matched reasonably
well in magnitude. However, for different irrigation events, anearly or delayed (e.g., Figure 8) predicted flow peak relative tothe measured data was observed [Mohanty et al., 1997]. Furthermore, zero (westerly) flow was observed for several hours(e.g., hours 23 ; 59, Figure 8) following every irrigation event.
Also, observed water table levels in the adjacent monitoring

Figure 8. Comparison of predicted and observed tile flows
following an irrigation event [after Mohanty et al., 1997].

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

wells and piezometers (not shown here) during and following aflood irrigation event revealed certain variations in the groundwater flow fields near the manholes warranting further investigation. We will discuss possible reasons for this unusual flowtransport behavior at our field site in the following sections.
Many earlier tile drainage studies neglected, for simplicity,
the effects of a regional groundwater table gradient on the tile
flow calculations. Tod and Grismer [1996] and Singh and Kanwar [1995a, b] indicated that the topographic and regionalhydrologic conditions could well be a significant factor determining tile flow behavior. To investigate the significance of aregional ground water table gradient (at the field or benchscale), we analyzed data from two different dates with contrasting irrigation sequencing: the east-center-west bench (asdefined in Table 1), August 24, 1995 (Figure 10a), and thewest-center-east bench, September 7, 1995 (Figure 10b). Several factors (as shown in Figure 11) may contribute to the flowrate in the tile drain at a given time including factor a, thedegree of submergence of the underdesigned tile drain below
the water table (similar to observations by Khan and Rushton[1996b]); factor b, a regional (field-scale) groundwater table
gradient across the three benches; factor c, a local groundwatertable gradient immediately adjacent to the drain; and factor d,a bench-scale groundwater table gradient across the neighboring benches.
In a normal day between consecutive irrigation events,
groundwater table showed a general northeast to approximately southwest gradient across the entire field and a local
east to approximately west gradient near our study section(e.g., Figure 11b). Figure 10a shows that as water infiltrated
during and following flooding of the east bench, the east tile
drain became more submerged (factor a) and the east-tocenter bench groundwater table gradient increased (factor d),thus increasing flow through both the east and west manholes.
Apparently, this also resulted in a slight net loss of flow fromthe center bench due to the decreased local water table gradient there (factor c). Net center bench flow rates were determined simply by subtracting the flow measured in the east
manhole from the flow measured in the west manhole. Afterward, subsequent to flooding and infiltration in the center
bench, net flow rates increased for the same reasons as during

Figure 9. Comparison of predicted and observed tile flux
NO2
3 concentrations following the irrigation event on June 8,
1994, at the field site. Time 0 hour matches the midnight prior
to the irrigation.

1073

Figure 10. Typical tile flow behavior at the field site following (a) an east-center-west and (b) a west-center-east irrigation
sequence. Vertical dashed lines represent irrigation startingtimes for the respective benches. Positive values indicate westerly flow. Occasional easterly flow are shown as zero flow.

flooding of east bench. However, flow rates from the east
bench decreased as the east-to-center bench groundwater table
gradient decreased (factor d). Eventually, flow rates from theeast bench decreased to unmeasurable rates, and temporarily
stagnated or even reversed direction (i.e., flow occurred from
center to east bench). One-way valve of our flow meter could
not register the amount of reverse flow thus simply indicating
zero westerly flow. Subsequent to flooding and infiltration on
the west bench, the center-to-west bench water table gradient
was reduced (factor d), and net flow rates from the center
bench declined rapidly to unmeasurable rates, which again
stagnated or reversed direction (i.e., flow occurred from west
to center bench). Unmeasurable flow rates in both the manholes persisted for a period of time until the east-to-west gradient (factor b) in the drain recovered sufficiently to establish
measurable westerly flow rates.
For the west-to-east irrigation sequence (Figure 10b), following flooding and infiltration in the west bench, flow rates
from the center and later the east bench rapidly stagnated or
temporarily reversed because of west-east groundwater table
gradient (factor b). Westerly flow was reestablished following
flooding of the center and east benches. The observed peaks inthe net center bench flow rate following the irrigation event are
due to the recovery of measurable flow rates in the west manhole prior to their recovery in the east manhole. These findings
suggest that tile flow and conservative (convection-dominated)
solute transport in this field perhaps should have been betterinterpreted and modeled by means of a three-dimensionalmodel that considers regional water and chemical management
schemes (e.g., irrigation sequencing and fertilizer scheduling inneighboring benches) rather than the simplified two-dimensional approach used here for flow and transport to a tile drain.
Observed and predicted NO2
3 flux concentrations in the tile1074

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Figure 11. Factors contributing to transient tile flowtransport at our field site including (a) degree of submergenceof the tile drain, (b) regional groundwater table gradient (e.g.,groundwater table contour determined from 34 monitoring
wells along the experimental field boundary on May 30, 1995),(c) local groundwater table gradient adjacent to the tile drain,and (d) bench-scale groundwater table gradient based on irrigation sequence of different field benches (e.g., west bench
irrigation prior to center or east bench irrigation).
drain showed far better agreement when the bimodal hydraulicfunctions were used as compared to the unimodal van Genuchten-Mualem hydraulic functions (Figure 9). The most striking feature of Figure 9 is the presence of three peaks in theobserved nitrate data. We believe that the peak concentration(;12 mg/L at hour 18) for this east-to-west irrigation event
depicts field-scale preferential transport, while the peak (;9mg/L at hour 60) is likely an artifact of the irrigation sequence
and reverse flow as discussed in section 3. A smaller peak (;5mg/L at hour 12) in NO2
3 concentration within ;3– 4 hours
from the start of the irrigation event likely occurred through anarrow zone above the tile drain, possibly because of soil
disturbance during tile installation. Our piecewise-continuous
hydraulic functions, however, could not fully distinguish between the three fast-flow regimes because of their lumped(continuum) nature. Moreover, following an irrigation event,simulated and observed NO2
3 flux concentration in the tile flow
showed discrepancies of different degrees with time because of
the same reasons of irrigation sequence and timing as for water
flow. Piezometer NO2
3 concentrations (Figure 2) were somewhat over predicted by the model. In general, using the uni-

modal or bimodal hydraulic functions, the model predicts 2–3
times higher concentrations than the observed values. This
result may be due to the assumption of impermeable hydrologic boundaries which neglect any mixing, dilution, and transport toward or from adjacent sites. Besides, uncertainties associated with the transport parameters of nitrogen species andother unrepresented nitrogen transformation processes may
also contribute to these discrepancies.
Figure 12 shows predicted and measured field-averaged resident NO2
3 concentrations in the soil profile for the long-term
simulation. Predicted field-averaged resident NO2
3 concentrations in the soil profile compared well to the corresponding
field measured values. The model predictions were better forthe shallow depths than for the deeper depths. We suggest thatchanges in the groundwater flow fields, as described earlier,
may have been partly responsible for this depth-dependent
behavior of resident NO2
3 concentration.
To verify the contribution of groundwater flow fields to theobserved discrepancies in flux and resident NO2
3 concentrations, as well as in the tile flow volume, we conducted a stepzoom-in analyses for determining the tile drain capture zone
during (preferential) peak and other peak flow periods. Also ina somewhat similar study in Indiana, Kladivko et al. [1991]
found evidence of preferential flow to the tile drain originating
in its close proximity. Keeping these findings in perspective, thestep-zoom-in procedure was adopted by reducing the flow domain width in discrete steps, as shown in Figure 7. In each step,
the flow domain was shrunk by approximately one half on bothsides of the tile drain and the hydrologic boundaries between
consecutive tile drains were assumed to match the modified
flow transport domain boundaries. Hence hydrologic boundaries between consecutive tile drains spanned a finite length
recharging to shallow groundwater that did not converge to
any local tile drain (i.e., followed the regional groundwaterflow path). For the step-zoom-in procedure we kept the depth

Figure 12. Comparison of predicted and observed (1995)
field-averaged profile resident NO2
3 concentrations during thelong-term simulation. Measured 1994 data were used as initialcondition. Predictions (results) use bimodal hydraulic functions only.

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Figure 13. Comparison of tile flux NO2
3 concentrations fordifferent effective capture lengths of the tile drain. The test
case is for the irrigation event on June 8, 1994, at the field site.
of the flow domain constant (i.e., 7 m, up to the impermeable
layer). This procedure resulted in five simulation runs with theflow domain width decreasing from 76 m up to 3.5 m.
By zooming in, flow fields gradually changed from threedimensional (composite effect of local vertical seepage andregional groundwater flow) to two-dimensional local flow inclose proximity of the tile drain (i.e., flow domain width of
3.5 m). Predicted tile flow NO2
3 flux concentration using thefive selected domain widths (for June 8, 1994, east-center-west
irrigation event) are presented in Figure 13. Apparently, adecrease in domain width to 3.5 m leads to somewhat better
matches of the short-term model predictions to the observed
nitrate concentration data in the tile drain at initial times (i.e.,
between hours 15 and 25). This finding suggests that during
and immediately after an irrigation or heavy rainfall event, thetile drain located at 1.2-m depth at the field site is able to catchthe flow and soluble NO2
3 from a narrow region adjacent to thetile drain. Assuming a somewhat similar mixing layer concept
as used by Shalit and Steenhuis [1996], water during irrigation
short-circuited horizontally from the irrigation ports via a network of surface-opened (shallow) cracks toward the tile drain
and then infiltrated through the (disturbed) soil zone above thetile drains. Any chemical transported beyond this effective
capture zone reached the groundwater table through someother pathways and eventually may or may not have reached
the tile drains since nitrate reaching groundwater will undergo
some denitrification by converting to N2, N2O, NO, and NO2
under conditions of reduced oxygen supply.
Measured NO2
3 flux concentrations following the irrigationwater application (between hours 25 and 40, Figure 13)
matched reasonably well with the model prediction when using
the full flow domain width of the tile drain (i.e., 76 m). In otherwords, during this period the total field width started contributing to the tile flow. Between hours 40 and 60 a dip in the tile
flow NO32 concentration data reflects the typical reduced/
reverse flow process during the east-center-west irrigation sequence as described earlier. In Figure 13, underpredictions
and overpredictions in nitrate flux concentration beyond hour
60 are suggested to be associated with unrepresented nitrogentransport/transformation processes. Similar trends for tile flow
NO2
3 concentrations and their match with model predictions
were found for other irrigation events (not further shown

1075

here), although peaks gradually reduced over time reflecting
depletion (Figure 3).
The above findings relating to reduced/reverse flow also
qualitatively agree with observations near the manholes andadjacent piezometric levels that the tile drain remained submerged during wet seasons. One of the reasons for this behavior is likely an underdesigned tile drain. In addition, the observed three-dimensional flow at the site as caused by thesequence and timing of irrigation of the three field benches(west to east or east to west) governed the field-scale hydrology
and solute transport to the tile drain, resulting in some temporal shifts in the model predictions. Regional groundwater
gradients in our two-dimensional flow-transport model wererepresented by monthly head values (for long-term simulations) at hydrologic divides across the study section. Boundary
conditions of this type were found to be too crude to catch theeffect of field-scale hydrology on tile flow following any shortterm flood irrigation event. Hence, on the basis of our fieldmonitoring program, our mathematical modeling, and our data
analyses, we infer that tile drains do not necessarily provide agood integrated representation of field-scale flow/transport behavior, in particular when the tiles are underdesigned and/or
the hydrology at the field site is three-dimensional. Regionalhydrology may also be a significant factor, in addition to the
local hydrology. However, tile drains still should be useful forcatching preferential flow phenomena from areas close to thedrains, following any heavy water input event. We conclude
that the proposed model based on CHAIN_2D and using
piecewise-continuous hydraulic functions appears to be an appropriate tool for making flow-transport prediction in tiledrained fields.
Acknowledgments. This project was partly funded by the USDACSRS through special water quality grant 9301322. Thanks to graduate
students Tracy Roth, Bob Reedy, Gunter Hilmes, and Jeff Chaves at
New Mexico Institute of Mining and Technology, Socorro, New Mexico, for field and laboratory data collection.

References
Booltink, H. W. G., Field monitoring of nitrate leaching and water flow
in a structured clay soil, Agric. Ecosyst. Environ., 52(2–3), 251–261,
1995.
Bowman, R. S., and R. C. Rice, Transport of conservative tracers in
the field under intermittent flood irrigation, Water Resour. Res., 22,1531–1536, 1986.
Celia, M. A., E. T. Bouloutas, and R. L. Zabra, A general massconservative numerical solution for the unsaturated flow equation,
Water Resour. Res., 26, 1483–1496, 1990.
Czapar, G. F., and R. S. Kanwar, Field measurement of preferential
flow using subsurface drainage tiles, paper presented at Preferential
Flow Symposium, Am. Soc. of Agric. Eng., St. Joseph, Mich., 1991.
Feddes, R. A., E. Bresler, and S. P. Neuman, Field test of a modified
numerical model for water uptake by root systems, Water Resour.
Res., 10, 1199 –1206, 1974.
Feddes, R. A., P. J. Kowalik, and H. Zaradny, Simulation of Field Water
Use and Crop Yield, John Wiley, New York, 1978.
Fipps, G., R. W. Skaggs, and J. L. Nieber, Drains as a boundary
condition in finite elements, Water Resour. Res., 22, 1613–1621, 1986.
Hallberg, G. R., J. L. Baker, and G. W. Randall, Utility of tile-line
effluent studies to evaluate the impact of agricultural practices ongroundwater, paper presented at Agricultural Impacts on Ground
Water Conference, Natl. Water Well Assoc., Dublin, Ohio, 1986.
Huang, K., B. P. Mohanty, and M. T. van Genuchten, A new convergence criterion for the modified picard iteration method to solve thevariably saturated flow equation, J. Hydrol., 178, 69 –91, 1996.
Journel, A. G., and C. J. Huijbregts, Mining Geostatistics, 600 pp.,
Academic, San Diego, Calif., 1978.

1076

MOHANTY ET AL.: PREFERENTIAL TRANSPORT OF NITRATE TO A TILE DRAIN

Jury, W. A., Solute travel-time estimates for tile-drained fields, I,
Theory, Soil Sci. Soc. Am. J., 39, 1020 –1024, 1975a.
Jury, W. A., Solute travel-time estimates for tile-drained fields, II,
Application to experimental studies, Soil Sci. Soc. Am. J., 39, 1024 –
1028, 1975b.
Jury, W. A., W. F. Spencer, and W. J. Farmer, Behavior assessment
model for trace organics in soil, I, Model description, J. Environ.
Qual., 12, 558 –564, 1983.
Kamra, S. K., S. R. Singh, K. V. G. K. Rao, and M. T. van Genuchten,
A semidiscrete model for water and solute movement in tile-drainedsoils, 1, Governing equations and solution, Water Resour. Res., 27,2439 –2447, 1991a.
Kamra, S. K., S. R. Singh, K. V. G. K. Rao, and M. T. van Genuchten,
A semidiscrete model for water and solute movement in tile-drainedsoils, 2, Field validation and applications, Water Resour. Res., 27,2449 –2456, 1991b.
Khan, S., and K. R. Rushton, Reappraisal of flow to tile drains, II,
Time-variant response, J. Hydrol., 183, 367–382, 1996a.
Khan, S., and K. R. Rushton, Reappraisal of flow to tile drains, III,
Drains with limited flow capacity, J. Hydrol., 183, 383–395, 1996b.
Kladivko, E. J., G. E. Van Scoyoc, E. J. Monke, K. M. Oates, and W.
Pask, Pesticides and nutrient movement into tile drains on a silt
loam soil in Indianan, J. Environ. Qual., 20, 264 –270, 1991.
Luxmoore, R. J., On preferential flow and its measurement, paper
presented at Preferential Flow Symposium, Am. Soc. of Agric. Eng.,
St. Joseph, Mich., 1991.
Mohanty, B. P., R. S. Bowman, J. M. H. Hendrickx, and M. T. van
Genuchten, New piecewise-continuous hydraulic functions for modeling preferential flow in an intermittent-flood-irrigated field, Water
Resour. Res., 33, 2049 –2063, 1997.
Munster, C. L., R. W. Skaggs, J. E. Parsons, R. O. Evans, J. W. Gilliam,and M. A. Breve, Simulating Aldicrab transport in a drained field,
Trans. ASAE, 37, 1817–1824, 1994.
Munster, C. L., R. W. Skaggs, and V. R. Pemmireddy, Effect of watertable management on the fate of the pesticide Aldicrab, Trans.
ASAE, 39, 55– 66, 1996.
National Resource Conservation Service (NRCS), Soil survey map of
Las Nutrias, New Mexico, Lincoln, Neb., 1992.
Nielsen, D. R., and J. G. MacDonald (Ed.), Nitrogen in the Environment, vol. 2, 528 pp., Academic, San Diego, Calif., 1978.
Novotny, V., and G. Chesters, Handbook of Nonpoint Pollution, 555
pp., van Nostrand Reinhold, New York, 1981.
Perrochet, P., and D. Berod, Stability of the standard Crank-NicolsonGalerkin scheme applied to the diffusion-convection equation:
Some new insights, Water Resour. Res., 29, 3291–3297, 1993.
Ragab, R., D. M. Cooper, G. L. Harris, and J. A. Catt, Simulating
nitrate under winter wheat grown on a structured clay soil considering bypass flow, J. Hydrol., 182, 157–173, 1996.
Richard, T. L., and T. S. Steenhuis, Tile drain sampling of preferential
flow on a field scale, J. Contam. Hydrol., 3, 307–325, 1988.
Samani, Z., Measuring water in trapezoidal canals, J. Irrigat. Drain.
Eng., 119, 181–186, 1993.
Sanoja, J., R. S. Kanwar, and S. W. Melvin, Comparison of simulated(DRAINMOD) and measured tile outflow and water table elevations from two field sites in Iowa, Trans. ASAE, 33, 827– 833, 1990.

Shalit, G., and T. Steenhuis, A simple mixing layer model predicting
solute flow to drainage lines under preferential flow, J. Hydrol., 183,
139 –149, 1996.
Shapiro, S. S., and W. B. Wilk, An analysis of variance test for normality (complete examples), Biometrika, 52, 591– 611, 1965.
Simunek, J., and M. T. van Genuchten, The CHAIN_2D code forsimulating the two-dimensional movement of water, heat, and multiple solutes in variably-saturated porous media, Res. Rep. 136, U.S.
Salinity Lab., Riverside, Calif., 1994.
Singh, P., and R. S. Kanwar, Modification of RZWQM for simulating
subsurface drainage by adding a tile flow component, Trans. ASAE,
38, 489 – 498, 1995a.
Singh, P., and R. S. Kanwar, Simulating NO3-N transport to subsurface
drain flows as affected by tillage under continuous corn using modified RZWQM, Trans. ASAE, 38, 499 –506, 1995b.
Singh, P., R. S. Kanwar, K. E. Johnsen, and L. R. Ahuja, Calibrationand evaluation of subsurface drainage component of RZWQM v.
2.5, J. Environ. Qual., 25, 56 – 63, 1996.
Thomas, G. W., and B. J. Barfield, The unreliability of tile effluent formonitoring subsurface nitrate-nitrogen losses from soils, J. Environ.
Qual., 3, 183–185, 1974.
Tod, I. C., and M. E. Grismer, Issues confronting salinity and drainage
management in irrigated agriculture, USCID Newsl., pp. 4 –5, U.S.
Comm. on Irrigat. and Drain., Denver, Colo., Jan.–July 1996.
van Genuchten, M. T., A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892– 898, 1980.
van Genuchten, M. T., Convective-dispersive transport of solutes involved in sequential first-order decay reactions, Comput. Geosci., 11,
129 –147, 1985.
Vimoke, B. S., and G. S. Taylor, Simulating water flow in soil with anelectric resistance network, Res. Rep. 41-65, Soil and Water Conserv.
Res. Div., U.S. Dep. of Agric., Agric. Res. Serv., Columbus, Ohio,
1962.
Vimoke, B. S., T. D. Tura, T. J. Thiel, and G. S. Taylor, Improvements
in construction and of resistance networks for studying drainage
problems, Soil Sci. Soc. Am. Proc., 26, 203–207, 1963.
White, R. E., S. R. Wellings, and J. P. Bell, Seasonal variations in
nitrate leaching in structured clay soils under mixed land use, Agric.
Water Manage., 7, 391– 410, 1983.
Yeh, G. T., and V. S. Tripathi, HYDROGEOCHEM: A coupledmodel of HYDROlogic transport and GEOCHEMical equilibria inreactive multicomponent systems, Publ. 3170, Environs Sci. Div.,Oak Ridge Natl. Lab., Oak Ridge, Tenn., 1990.
R. S. Bowman and J. M. H. Hendrickx, Department of Earth andEnvironmental Sciences, New Mexico Institute of Mining and Technology, Socorro, NM 87801.
B. P. Mohanty, J. Simunek, and M. T. van Genuchten, U.S. Salinity
Laboratory, 450 W. Big Springs Road, Riverside, CA 92507. (e-mail:bmohanty@ussl.ars.usda.gov)(Received August 15, 1997; revised January 13, 1998;
accepted January 22, 1998.)