Role of Groundwater Flow in Tile Drain DischargePeter J. Vaughan,* Donald L. Suarez, Jirka Simunek, Dennis L. Corwin, and James D. Rhoades
ABSTRACT
Tile systems drain water applied to agricultural fields as irrigationand precipitation but also may intercept regional groundwater flow.
Identification and characterization of the potential sources of tile
water is essential for informed management of salinity and contaminants. Factors influencing tile discharge including depth of water applied, evapotranspiration, water storage, drain blockage, and interception of regional groundwater flow were evaluated to determine which
may be related to a fivefold variation in cumulative tile discharge
among six sumps located 100 km west of Fresno, CA. Cumulative
depths drained were calculated for 5 yr of weekly irrigation, precipitation, and discharge data. Evapotranspiration and water storage were
estimated using the UnsatchemGeo variably-saturated water flow
model. Well water levels measured on 19 dates were spatially-averaged
providing spatial variation of depth-to- water among the drained areas.
Variability in depth of water drained (0.18-0.95 m) was large and was
not correlated with either water applied (3.26-4.58 m, r2 = 0.03) orwith computed water flux from the bottom of the soil column (0.050.31 m, r2 = 0.00). Groundwater interception by tile drains was afactor because depth-to-water was negatively correlated with discharge(r2 = 0.42) and drawdown of groundwater levels by drains was relatively larger for those drained areas encountered first during regional
groundwater flow. For all six sumps, drained water is likely derived
from locally applied water and interception of regional groundwaterflow implying that standard two-dimensional models of flow to drains,representing only water applied locally, would not be applicable tomodeling of drain flows or drain-water solute concentrations.

U.S. Salinity Laboratory, 450 W. Big Springs Road, Riverside, CA
92507. Received 18 Nov. 1997. *Corresponding author (pvaughan@
ussl.ars.usda.gov).
Published in J. Environ. Qual. 28:403-410 (1999).

T

ILE DRAINAGE in agricultural areas provides amethod of regulating the depth of shallow watertables. The water drained from tile systems is frequently
not reusable as irrigation water because of high salinity
and contamination. Two contaminants often present in
the San Joaquin Valley of California are B and Se.
Disposal of tile water in this area is a necessary component of a management system that includes tile drainage.
But tile water disposal raises environmental concerns
such as the well-known experience of Se contamination
at Kesterson reservoir in the San Joaquin Valley. The
source of tile water is shallow groundwater which may
be derived from water applied during irrigation of afield, from local precipitation, or it may have moved
laterally to the tile drain by regional groundwater flow.
The determination of the relative importance of thesecomponents is significant because drainage of irrigationwater represents an economic loss. Drainage of groundwater which has moved laterally from other areas is,however, not locally controllable and cannot be factored
into the local economics of production. Regulation ofdrain water quality is therefore a complex issue because
there are, potentially both local and regional components controlling water quality.
Tile drain tubes are normally installed in parallel setswith a fixed spacing. This arrangement can be modeledby analytical techniques which lead to expressions forflow rates within tubes based on drain diameter, spacing
Abbreviations: CIMIS, California Irrigation Management Information System; CIS, geographic information system.

404

J.

ENVIRON. QUAL., VOL 28,

and various boundary conditions for water flux and hydraulic head (Kirkham, 1949, 1958; Jury, 1975a; Fipps
and Skaggs, 1991). The lower boundary condition ofthese models is an impervious layer which permits asolution using the method of images. This assumed
boundaryimplies that all water arriving at a drain originated at the soil surface between the drains. In theanalytical solutions, certain flow paths pass throughlower elevations than the elevation of the drain resulting
in upwards flow to the drains. Such paths may have
travel times on the order of 10 yr or more (Jury, 1975b).
The utility of the impervious layer assumption is questionable under certain field conditions as, for example,
in the western San Joaquin Valley. In this area, a highly
permeable sand layer underlies clay-rich soils that havemuch lower permeability (Deverel and Gallanthine,
1989; Belitz and Phillips, 1995). Groundwater flow isfrom southwest to northeast and the historic record ofartesian wells at the margin of the alluvial fan near the
San Joaquin River (Mendenhall et al., 1916) suggested
that this groundwater flow was occurring in conditions
where the groundwater level was close to the surface.
Groundwater levels in this area were lowered by pumping during the first half of the 20th century, but the
introduction of irrigation water from northern California, through the Delta-Mendota canal, and the installation of tile drainage have again stabilized the water table
at a depth of 1.5 to 3 m(Deverel and Gallanthine, 1989).
Water flowing to tile drains in this area is a mixture of
local irrigation water and water that has movedupwards
from the sand layer (Deverel and Fujii, 1988). Evidence
for this partitioning includes isotopic enrichmentof g180
with depth, and decreases in tritium concentrations
when compared with irrigationwater sampled at thesame site (Deverel and Fio, 1991). Based on water chemistry data, the proportion of drained water that arrived
from greater depths was 30%for a drain at 1.8-m depth
and 60%for a 2.7-m drain. Furthermore, Se concentrations indicated that the deep groundwater component
had arrived from the west (Deverel and Fio, 1991).
4-2

2
Sumps
(numbered)¯o Wells[~ TileDrains~ Drainage
10-2 ~

Areas

StudyArea
Boundary/~ GroundwaterMovement

16,
15"~
i

I

1 km

I

N

Fig. 1. Map of six drainage areas that drain to the numbered sumpin the northeast corner of each area. Groundwater movement direction was a mean value computed for the entire water districtusingwater table elevation data for 1991. Well locations are 24 of the
57 wells used in water table surface calculations.

MARCH-APRIL 1999

These geochemical studies of groundwater and tiledrained water on the west side of the San Joaquin Valley
indicated that a significant portion of the water appearing in tile drains was not applied to the same field,
but was water that had flowed upwards from greater
depths and horizontally from the west.
In this study, we examined the water budget of anirrigated agricultural area located in the same generalarea of the western San Joaquin Valley as the aforementioned groundwater geochemical studies (Deverel andGallanthine, 1989; Deverel and Fio, 1991). ConsideringWet-]- I42ro ~- Wdp~- AWst-I- Wtd-- Wip-- Wig= 0 [1]a water balance in terms of water depth ILl for a control
volume consisting of a soil column and a tile drain.
Water enters the volumeas irrigation, precipitation ortile interception of groundwater and leaves as evapotranspiration, runoff, deep percolation, and tile drainage. The terms are: w~p depth of water applied throughirrigation and precipitation, wet evapotranspiration, Wrorunoff, Wdpdeep percolation, w,d tile drainage, Awst
change in soil water storage (expressed as a depth ofwater within a soil columnof specified height), and W~ggroundwaterinterception by tile drains. Spatial variability in the depth of water drained by tile systems could
be caused by spatial variability in any of the amounts
represented in Eq. [1] or by some combination theseterms. The question addressed here was the origin ofvariability in depth of water drained given these spatially-variable processes.
MATERIALS

AND METHODS

Field Area
The BroadviewWaterDistrict is located on the west sideof the San Joaquin Valley of California approximately 100
kmwest of Fresno. Mostof the water district is divided intoquarter-sections whichare square fields with area of 64.8 ha.
Onthe bordersof the waterdistrict there are also somesmaller
fields with irregular shapes. This study was concernedwith
only quarter-sections, and henceforthwe will refer to quartersections as fields. Thedistrict maintainsirrigation turnoutsin
the southwestcorner of each field. Watersupplied fromtheseturnouts is meteredwith a precision of 123m3 (0.1 acre-feet).
Sumpsare located in the northeast corner of certain fields forcollection and disposal of waterfromtile-drain systems. Water
draining into the sumpsis pumpedinto drainage ditches andmetered to the same precision as the irrigation water. The
sumpsmaydrain a single field, twoadjacent fields, or all four
fields that makeup an entire section (Fig. 1). Becausesumps
maydrain morethan one field, the term drainagearea is used
to refer to the entire area drained by a sump.
Approximately80%of the water district has tile drainage
at approximately 1.8-m depth. The most commonspacing ofdrain laterals is 91.4 m, but several other spacings are alsopresent (Fig. 1). Mapsof drains were compiledby the water
district staff using data obtainedfromindividual landowners.
Six sumpswereselected for analysis of irrigation and drainage,
each located in the northeast corner of an entire section.
Every2 d, cumulativeflow meters for irrigation were read
and weekly readings of meters on sumpswere taken during
the period of i Jan. 1991to 1 May1996 (the study period).
Daily precipitation data, collected at the California Irrigation
Management
Information System (CIMIS)station 10 kmeastVAUGttANET AL.: ROLEOF GROUNDWATERFLOWIN TILE DRAINDISCHARGEof the study area, in Firebaugh, were combinedwith the irrigation data to generate a complete record of water applied toeach field. There were probably minor differences in precipitation at the two locations but these differences were not significant for the analysis of spatial variability within the water
district.
For the summergrowing seasons between 1991 and1995, the land-use distribution for crops in the drainage areas
was: 47% cotton (Gossypium hirsutum L.), 15% tomato (Lycopersicon lycopersicum L.), 11%alfalfa (Medicagosativa L.),
9% melon (Cucurnis melo L.), and 18% fallow. No wintercrops were grown during this period.
Simulationof Drainage from the Root Zone

Someterms in Eq. [1] can be evaluated easily from available
data but others must be estimated. Evapotranspiration was
an important componentof the water balance but direct measurements of soil evaporation and transpiration were not available, so model calculation of evapotranspiration was necessary. In the San Joaquin Valley, crops such as cotton andalfalfa are often grown under conditions of periodic water
stress, as a managementtool, to increase crop production. A
significant effect of water stress is reduction of transpirationfrom that of a fully-watered crop implying that daily transpiration can be considered as coupled to the soil water content.
One model that includes coupling between these processes is
the UnsatchemGeomodel which computes variably-saturated
water flow, plant root water uptake, multicomponent solutetransport, and other processes occurring in the vadose zone
(Vaughan et al., 1997). UnsatchemGeois a generalization
the original Unsatchemmodel and operates in the context ofa geographic information system (GIS), which stores both
the data required by the model and the results. The original
Unsatchemmodel is a finite-element solution of the Richards
equation coupled to various numerical models of heat, multicomponent chemical, and carbon dioxide transport (Suarez
and Simunek, 1994; Simunek and Suarez, 1994; Simunek etal., 1996; Suarez and Simunek, 1997). UnsatchemGeomodels
one-dimensional water flow and other processes at point locations within a geographic area. Data required by this model,such as initial water content as a function of depth, are specified for point locations. Other data, such as material proper*
ties, are specified for areas rather than points. The hierarchicalnature .of data organization in the model allows storage ofdata for points or areas and avoids unnecessary duplicationof data. For the calculations described in this paper, the UnsatchemGeomodel was set up to compute only water flow, plantroot water uptake, heat flow, and root growth. The computed
water flux at the lower boundary of the modeled soil column
was of primary importance in this paper because it can becompared with the measured drainage.
Other terms in Eq. [1] that can be studied using UnsatchemGeo are the change in water storage (Aw~t) and runoff (Wro).
The surface boundary condition for water flow for the modeling described here was time-dependent, prescribed flux. At
certain times, this prescribed flux exceeded the computedinfiltration rate, resulting in ponding at the soil surface and an
automatic reset to a constant head condition. Runoff wasassumed for ponding depths >0.05 m. For a one-dimensionalmodel this means ponding depths were limited to 0.05 m, andany further water application was considered runoff and was
ignored. No measurements comparing runoff with surface water depth were available for the study area, so the ponding
depth limitation of 0.05 m was simply a reasonable assumption.
At the lower boundary of the modeled soil column thespecified flux boundary condition was the "deep drainage"
condition (Simuneket al., 1996; Hopmansand Stricker, 1989),

405

q = -A exp(-Blh- h0[)[2]
where q is the flux [L3 L-2 T-~], A [L3 L-~ -Z]
T-~], and B [L
are parameters, and h0 is the pressure head at the bottom
node when the entire soil profile is fully-saturated [L]. For
the calculations discussed here, A = 0.169 cm d-~ and B =0.027 cm-L
Hydraulic properties required by the model were measured
or estimated for the two main soil types that occur in thedrainage areas. The soil types are the Lillis series (very fine,montmorillonitic, thermic Entic Chromoxererts) and the Cerini series (fine-loamy, mixed (calcareous), thermic Typic Torrifluvents)
mapped by the Natural Resources Conservation
Service (a preliminary, unpublished mapbased on soil surveys
taken during 1990-1992). The measured hydraulic property
was saturated hydraulic conductivity which was a mean value
for two samples collected for each soil texture (Table 1). The
two cylindrical soil samples (5.40 cm diam.) were collectedusing a bulk density sampler from a depth of 0.3 mat a spacingof 0.1 m. In 1991, a soil sampling operation collected soil
samples from 315 locations within the Broadview Water District. From this sampling, an average saturation percentage
for each soil type was determined by taking the mean for all
samples collected from locations within the mapped area ofthat soil type. Locations for the soil sampling reported here,
were selected from the sets of locations in the 1991 survey,
as those which most closely matched the mean saturation percentage.
UnsatchemGeo uses the van Genuchten-Mualem model ofwater retention with m = 1 - 1/n,0(h) = 0r +(1 + Ic~hln) 0-~j’)h < hsh >- hs[3]

where 0~ and 0r are the saturated and residual water contents(cm~ cm-~), ~ (cm-~) and n are empirical parameters, h (cm)is the pressure head, and h~ (cm) is the pressure headsaturation (van Genuchten,1980; Simuneket al., 1996). 0~, 0r,~, and n were selected from a list of properties for a variety ofsoils (Carsel and Parrish, 1988). Error estimates for hydraulicproperties given in Table 1 were not available because thevalues were either assumed,or, in the case of hydraulic conductivity, only two measurements were available.
Soil evaporation and plant transpiration
were modeled
through a combinedpotential evapotranspiration. Daily reference evapotranspiration,
ET0 (ram -~) was o btained f rom
stored records representing conditions measured at the CIMIS
station in Firebaugh. ET0in these records was calculated using
the Penmanequation (Penman, 1963). A potential crop evapotranspiration, Tp (cm -~) was o btained f rom
T~ = 0.~ E~0~
where K~is the crop coefficient (Doorenbos and Pruitt, 1977).
The evapotranspiration was split into soil evaporation and
Table1. Soil hydraulicproperties.
Saturated
Saturated Residual
hydraulic
Soil
water
water Parameters: conductivity’~
texturecontent
content[Eq.3lnK, (m/d)-3) 0, (cmS cm-~)
0~fern~ cm~ {m~)
Sandyclay~
0.45
0.10
2.7
1.23 0.016
Clayloam§
0.41
0.09
1.9
1.31 0.lll
Meanof measured
values for twosamples.
Lillisseries.
Ceriniseries.

406

J.

ENVIRON. QUAL., VOL. 28,

MARCH-APRIL1999

plant transpiration under the assumptionthat the fractionalrooting depthalso represents the fraction of the evapotranspiration that is transpiration (Simuneket al., 1996). Whenthe
actual rooting depth reaches the maximum
specified value, thecanopycoverageof the soil is complete,the soil evaporationis
zero, and the potential transpiration is equal to the potential
evapotranspiration. Root water uptake at depth z was calculated using

wasa digital line graphfile of topographiccontour locations
obtained from the U.S. Geological Survey. Elevation of thewater table wasdetermined,for each well, by subtracting thethe measureddepth-to-water from the computedsurface elevation. Theresulting watertable elevations werefitted using athird-order trend surface. Residualsfor water table elevations
were calculated at each well location. These residuals werekriged to provide an improvedrepresentation of each surface~c(Z) -- ,~(h)~(Z)7;[5]where Tc(Z) is the rate of water removal(cm -1 d- 1) with
13(z) the normalizedroot distribution function (cm ~).
water stress function is

Zest = Zt ~- Zk

1[61
1 + (h/h5o)wherehs0 is the pressure head causing a 50%reduction whenthe exponent p = 1 (see Simuneket al., 1996, for further
details on root water uptake). The parameter hs0 was set to-20 mand p = 3 in the calculations discussed here.
The UnsatchemGeo
simulation of water flow in the various
drainageareas wasset to run over the study period. Theinitialconditionfor this simulationwasvolumetricsoil watercontent,
0(z). The gravimetric soil water content was determinedfor
samplescollected at 0.3-m intervals to a depth of 1.2 mfor
each location (Leschet al., 1992).Aseparateset of 75 measurementsof bulk density taken both within the study area (Fig.
f), and in fields nearby, were kriged to obtain estimates ofbulk density at sampling locations (Vaughanet al., 1996).
0(z) was computedusing these estimates and was linearlyinterpolated betweenthe samplingdepths. The finite-element
columnextendedto a maximum
depth of 1.84 m. The time step
was varied during the simulation depending on howrapidly
convergenceoccurredduring iteration of the water flow calculation. The minimumtime step was 8.6 s and the maximum
was1 d. All changesin prescribedwaterflux at the soil surface
were precededby a gradual reduction of the time step to theminimumvalue.
Calculations were performedfor two locations within eachfield. Cumulativedrainage from the bottom of the modeled
profile was comparedwith measureddepths of drainage whichwere recorded on a weekly basis. Computedresults werestored at time intervals varying from one to two days rather
than at the computationaltime step to minimizedata storage
requirements. A subset of computedresults was selected tomatch the dates whenmeasurementswere taken. Meancumulative drainage depths for each drainage area werecalculated
as a simple arithmetic meanof the coincident data from all
locations within the area.

-pwherezt is estimated from the trend surface and zk isestimated by kriging the residuals. Kriging permitted better
representationof the variability at shorter distances. All kriging calculations used the ordinary kriging program(OKB2D)
in the GSLIBpackage (Deutsch and Journel, 1992). For thewater table elevations and residuals, the surfaces were represented as grids with cell size of 20 m. Spatial covariancewas
modeledfor semivariogramdata for each of the 19 dates.
Thesemodelsnormally included a small nugget effect (meanof nugget values: 0.016 m2) and an omnidirectional sphericalmodelbut the residuals for 2 of the 19 dates required directional modeling.
(Zest)RESULTSAND DISCUSSION

During the 1948 d represented in the study period,
cumulative tile drainage for the six fields was a maximumdepth of 0.95 m in field 16-2 and a minimumof
0.18 min field 9-2, a factor of five variation (Fig. 2).
This variation in depth drained was also consistent overtime. The main question, addressed here, was the originof this variability in depth of water drained. Somepossible contributing factors were: different amountsof irrigation water applied to the fields, different degrees ofwater use caused by planting different crops, spatial
variability of Awst[Eq. 1], blockageof tile drains causing
reduced flow and correspondingly greater deep percolation to the groundwater, lateral groundwater flow to
the sumpdrainage area, and variation of hydraulic conductivity for different soil types causing variability inthe efficiency of the drains.
Variation in cumulative depth of applied water among
the drainage areas was 3.26 to 4.58 m depth, or a factorof 1.4. This was significantly less than the factor of five~-

SumpIDj~ "----’~’---- 4-2

Determination of Depth-to-Water

[7]

0.8

~

{ --+-~-2~ ~ ~0-2

/
~

/"
¢

Groundwaterlevels throughoutthe study area were calculated from depth-to-water measuredat 57 wells for 19 dates
during the study period. Theselevels providedan independent
data set for use in evaluating the role of groundwaterin tile
drainage. Becauseof equipmentfailures, the actual numberof measuredwells for any given date varied between34 and54 with a meanvalue of 48. Thesurface elevation of well sites
was determined from a digital elevation modelproduced by~ program(a modified version of an algorithm
the TOPOGRID
described by Hutchinson,1989). Input data for the calculation
0
tA program in the ARC/INFO softwarenames in this report is for identification
constituteendorsement by the USDA.

package. The use of brandpurposes only and does not

500

1500
1000
Oa~s Since danua~ ~, ~ 991~
2000

Fig. 2. Cumulalive depth of water drained from six drainage areas
by tile drain syslems belween 1 Jan. 1991 ~nd 1 Nay 1996.

VAUGHAN ET AL.:

ROLE OF GROUNDWATERFLOW IN TILE

variation that occurred for the cumulative drainage
depths (Fig. 3). Furthermore, the depth of water appliedand depth drained for specific fields were poorly correlated (coefficient of determination, 2 =0. 03). Th
drainage area with the greatest drainage depth had thelowest depth of applied water (16-2). Also, the greatest
depth of applied water (4.6 m) was almost identical for
two drainage areas (3-2 and 4-2), but the depth drained
for these areas differed by a factor of 1.8. The datafor drainage area 4-2 did suggest a positive correlation
because 4-2 had the greatest depth of applied water andthe second largest cumulative drainage. But, considering
the entire data set, there was only a small, poorly-correlated variation in irrigation amount, implying that variations in drainage depths are probably better explainedby other factors.
Spatial variability in water storage is likely to be significant in the analysis of cumulative drainage for shortmodeling periods. But variations in water storage are
restricted to a finite range of soil water content implyingthat such variations will be less significant as the modeled time period increases. Aw~,[Eq. 1] is expressed as
the change in depth of water over the 1.84 m modeledsoil column. For the 5-yr period, computed Aws~varied
from a mean value of -0.11 m in drainage area 4-2 to
0.10 m in 10-2. This variation was not negligible andsuggested that spatial variation of soil water storage wasa potential contributing factor to the observedvariations
in drainage. However, the variation was not large
enough to consider spatial variation in/~ws~ as a major
factor in the 0.8 m variation in cumulative drainage
depth.
Comparedwith irrigation and precipitation amounts,
the explanation of drainage variability based on spatial
variation of evapotranspiration rates was less certain
because there was no source of continuous evapotranspiration estimates from measured data within individual
fields. However, evapotranspiration was computed by
the UnsatchemGeoprogram using the plant root water
uptake module. The calculated cumulative drainage
from the base of the modeledprofile predicted the frac5
~
~
4

407

tion of applied water that was drained. The remaining
fraction was either removed by evapotranspiration or
was stored in the profile. Because this study attempted
to determine whether variability in evapotranspiration
could have resulted in the observed variability in drainage, a conservative assumption was that the entire computed flow entered the tile drains and appeared at thesump. This assumption was equivalent to Wdp= 0 [Eq.
1], and would tend to maximizethe effect of variability
in evapotranspiration. The difference between observed
and modeled drainage, termed the drainage residual inthis paper, was plotted to demonstrate the degree towhich variability in observed drainage can be explainedby the model (Fig. 4). The model made substantiallybetter predictions of observed drainage for sumps 3-2,
9-2, and 10-2 than for sumps 4-2, 15-2, and 16-2. Predicted and observed cumulative drainage for each setof sumps, taken separately, were positively correlated.
But, when the data were combined, the correlation became nonexistent (r 2 = 0.00). By subtracting an arbitrary 0.5-m depth from the measured depths for sumps
4-2, 15-2, and 16-2, the correlation for all six sumpswas
improved (r 2 = 0.45). This suggested that poor model
performance, indicated by the nonexistent correlation,
could have been caused by additional water, from someother source, that entered sumps4-2, 15-2, and 16-2. In
summary,while variations in evapotranspiration among
different fields might have been a secondary contributing factor to observedvariability in cumulative drainage,
such variations were not a prime factor.
Variations in cumulative drainage depths could also
be caused by drain blockage. Although such variations
were not directly quantifiable because no data wereavailable regarding performance of individual drain laterals, someobservations of relative drainage rates wererelevant. The least cumulative depths drained occurred
for sumps 9-2 and 10-2 whereas the greatest depths
drained occurred at sumps 4-2 and 16-2. If blockage
were an important factor, then the argument would bethat drains in 9-2 and 10-2 tended to have more
blockage, resulting in slower flow rates. Drainage areas
0.8

SumpID:3-2
4°2

DRAIN DISCHARGE

SumpIDf~ 10-2

~=

~ 0.2
1
0
0

500

1000
1500
DaysSinceJanuary1, 1991

Fig. 3. Cumulative depth of water (irrigation + precipitation)
to six drainage areas designated by the sump identification
between 1 Jan. 1991 and 1 May 1996.

2000
appliednumber

500

1000
1500
DaysSinceJanuary1,1991

2000

Fig. 4. Cumulativedrainagedepthresidual(observedcalculated).
Drainagedepthswerecalculatedassuming
all waterdrainedfrom
the soil column
enteredthe tile drains.

408

J.

ENVIRON. QUAL., VOL. 28,

9-2 and 10-2, however, have the greatest density of tile
drains and the slowest drainage rates (Fig. 5). It seemed
unlikely that tile drain blockage would have been moresignificant in fields where the density of drains was
greatest. Also, there was no specific evidence of substantial tile drain blockage in drainage areas 9-2 and 10-2
(David Cone, manager of the water district, 1998, personal communication). Finally, drainage area 16-2 only
had a main drain running east-west at the northern endof the western field (Fig. 1), there are no tile drainlaterals. This accounted for the low tile drain density inthis drainage area. The lack of tile drain laterals could
be thought of as a hypothetical blockage of all drainsin this field except the main drain. Drainage area 16-2
had the highest cumulative depth drained, however, so
the complete elimination of flow from the hypotheticalset of drain laterals in half of the drainage area resulted
in no reduction of cumulative drainage depths relative
to other areas. This was quite surprising, but it certainly
suggested that blockage of tile drain laterals was unlikely to have been an important factor in this drainage
area. A more reasonable conclusion is that, over longtime periods, flow within tile drains in the study area
was probably not the rate-controlling factor for sumpdrainage rates.
The fifth possibility, groundwater flow into drains,
was difficult to evaluate because the flux of groundwater
into the drains cannot be calculated accurately. Instead,the approach we chose was calculation of mean depthto-groundwater within drainage areas, and comparisonof these depths with mean daily drainage. If depth-togroundwater were a significant factor then there shouldhave been a negative correlation between depth-togroundwater and mean daily drainage at least forgroundwater depths that were shallower than the depthof the drains. Modeling of water table elevation was
done to obtain estimates of the mean depth-to-water
within each drainage area for each of the 19 sets ofmeasurements of depth-to-water. The addition of thekriged residual surface to the original third-order trend
surface substantially improvedthe fit [Eq. 7]. The tern-

MARCH-APRIL1999

porally-averaged sum of squares taken over the 19 dates
was 6.22 m2 for the trend surface, but only 0.80 m2 when
the kriged surface was added. The kriged surface was
not an exact interpolation because of the non-zero nugget effect.
Comparison of well levels and mean daily depth
drained was done for the period between 1 Jan. 1991
and 2 Jan. 1996 (Fig. 6). This final date was the date
the last well-level measurementduring the study period.
The estimated mean depth-to-water, Zest [Eq. 7], wasspatially averaged over each drainage area by taking
the arithmetic meanfor all grid cells located within thearea including an areal weighting to accommodatethosecells along the boundaries which are smaller than 20 by
20 m. This was done for each measurement date. These
results were temporally-averaged over the 19 dates toobtain a mean depth-to-water. Resulting mean depths
were likely biased from the true mean depth-to-water
because the dates of the measurements did not follow
a sampling strategy that would adequately characterize
the seasonal fluctuations in water levels. Nonetheless,if there was spatial variation in depth-to-water that was
negatively correlated with daily depth drained, such spatial variation should have been apparent even if thecomputed means did not represent true mean values
because the seasonal fluctuations were approximately
the same at all sites. The results were consistent withnegative correlation of mean daily depth drained andmean depth-to-water (coefficientof determination,r 2 = 0.42, Fig. 6). A regression model, such as a linear
model, was not plotted because the functional form of
the relationship is likely to exhibit a threshold depthto-water occurring whenever the depth-to-water is lessthan the depth of the bottom of the drain. A rapid
increase in drain flows would be expected as the meanwater table moves upwards above the drains. The data
were generally consistent with a threshold model ofdrain flow because the three greatest depths drained
occurred in drainage areas having the three smallest
values of depth-to-water.
However, there are notenough data to adequately characterize a threshold
0.5

0.7

~] 16-2

0.6
0.4~" 0.54-2

16-2 []
4-2[~._¢ 0.4’~- 0.3[]15-2ls.2[]

3-2 I~[]

3-21~~ 13.2[] 9-2~o-2[]

0.1

10-2[]

9-2 []

60

80-1)SpecificTile Length
(mha

100

120

Fig. 5. Meandaily depth ofwater drained vs. areal density ofdrains forseveral drainage areas. Non-numbered points represent drainage
areas that surrounded the study area.

2.2
2.1
MeanDepthto Water(m)Fig. 6. Mean daily depth drained was spatially-averaged
within eachdrainage area and temporally-averaged for the 19 dates when measurements were made.

VAUGHAN ET AL.:

ROLE OF GROUNDWATER FLOW IN TILE

model, and the mean depth-to-water was likely biasedby the sampling such that it could not be compared tothe actual depth of the drains.
Other evidence, based on groundwater movement,
also supported the idea that regional groundwater was
flowing into tile drains. The groundwater movement
direction (north-northeast) was calculated from the average slope of the water table surface in 1991 (Fig. 1).
Waterlevels in six wells located along an east-west road,
0.8 km south of the study area, were compared withdata from all other wells in the district (Fig. 7). The
meanwell depths for these six wells were slightly greaterthan mean well depths for the remaining wells in the
district during 1991but, starting in late 1992, these welldepths were consistently less than those in the remainder of the district. This trend continued through Juneof 1997, as indicated by the last point on the plot. Thus,
for most of the study period, depth-to-water was actuallyless at locations upslope from the drainage areas. Twoof the three drainage areas showingthe greatest cumulative depth drained were located near the southernboundary of the water district (1.5-2 and 16-2). Groundwater moving north-northeast from areas to the south
and west of the district wouldencounter these tile drain
systems first and the maintenance of groundwater levels
at tile drain depths wouldtend to require higher flows ofnonlocal groundwater in these sumps. The two drainage
areas with the lowest cumulative drainage amounts (9-2
and 10-2) were not in the direct path of groundwater
moving into the district from the south and would be
expected to drain less groundwater, as observed (Fig.
2). Thus, drainage in areas 9-2, 10-2, 15-2, and 16-2 canbe rationalized by this argument based on groundwater
movement. Mean depth-to-water in drainage area 4-2
was 2.0 m, essentially the same as the mean depth-towater in 15-2 and 16-2. Area 4-2, however, was located in
the interior of the water district so the high groundwaterlevel and large cumulative depth drained would not
be explained by groundwater movement. The depth ofwater applied in 4-2 was equal to that applied to 3-2

409

DRAIN DISCHARGE

and approximately 1.2 m more than the depth applied
to the other drainage areas during the 5-yr period.
Depth of water applied was, therefore, probably a factor
in the greater amountof drainage in this drainage area.
Cumulative depth drained in drainage area 3-2 wassomewhatgreater than for 9-2 and 10-2, but the simulation of depth drained for 3-2 brought the drainage residual close to those of 9-2 and 10-2 (Fig. 4). Thus, spatial
variation in evapotranspiration and/or hydraulic conductivity could account for the variation in depth
drained among3-2, 9-2, and 10-2.
Because regional groundwater flow to tile drains was
not easily quantifiable, another approach involved examining relationships between the factors that werelikely related to groundwater movement and the observed drainage. One such factor was saturated hydraulic conductivity which varied by a factor of 7 for thetwo soil textures present in the study area (Table 1).
Greater saturated hydraulic conductivities would be expected to result in more rapid drainage of irrigationwater or groundwater. For irrigation water, greater hydraulic conductivities would be expected to drain water
more rapidly out of the root zone to the groundwaterlevel, reducing evapotranspiration losses. More rapid
drainage of groundwater that has movedlaterally would
be expected for the clay-loam soil which occurs in northnortheast-trending stringers. Thus, variations in areal
fraction of clay-loam amongthe fields might have been
related to the effectiveness of removal of either type ofwater by tile drains. The areal fraction of clay-loam, ineach drainage area, was determined using a GIS overlay
procedure. The data set was expanded to include someother drainage areas for which modeling data were notavailable. The data for daily depth drained versus areal
fraction of clay-loam were not correlated based on theavailable data (coefficient of determination, 2 =0. 002,
Fig. 8). Thus, a straightforward conclusion that greaterdrainage rates occurred in areas that had a greater areal
fraction of clay-loam was not supported by these data.
0.7 -

0.6
2.5
0.5

16-2[~

2
0.4.

[~4-2

1.5
0.~

~15-2
3-2~]

[] []
0.5
All OtherWells

500

1000
1500
2000
DaysSince January1,1991
Fig. 7. Depth-to-water is an arithmetic mean for well depths that were
measured on each specific date. Triangles represent six wells ’along
the southern boundary of the water district.
Squares represent all
the remaining wells.

~o-2
[]

0.1
[] 9-20204060Areal % Clay-Loam

80

100

Fig. 8. Areal percent of each soil texture was determined for eachdrainage area from GIS overlays of a preliminary soil texture map
based on a 1990 to 1992 soil survey of the area by the Natural
Resources Conservation Service.

410

J. ENVIRON. QUAL., VOL. 28, MARCH-APRIL 1999

CONCLUSIONSIn this 5-yr study of Broadview Water District irrigation and drainage records we have shown that variability
in cumulative tile drain flows could not be uniquely
explained by specific local variations in any one of thefollowing factors: irrigation amount, evapotranspiration, variation in water storage, or tile drain blockage.
Drainage depths that were consistently higher for certain sumps were especially difficult to rationalize on
the basis of yearly irrigation amounts which were quiteconsistent among the various fields. All these factors,
however, may have possibly contributed to the observed
variability and there were probably interactions among
them. Another possible explanation was variation in thefraction of nonlocal groundwater appearing at different
sumps. Evidence for this was negative correlation ofgroundwater levels and cumulative depth drained. In
addition, greater cumulative depths of drainage occurred in the two drainage areas that would be encountered first by groundwater moving in a north-northeast
direction. The relatively large depth drained by one
sump near the northern boundary of the water district
was probably not related to lateral groundwater flow but
might be explained by the greater amount of irrigationwater applied to that drainage area. Because several of
the possible explanations that were presented might be
operable to different degrees and might exhibit interactions, it was not possible to complete a quantitativeevaluation of the fractional contribution by each process. In general, caution is warranted when modelingtile-drain flow or the chemical composition of tile drainwaters because mixing with groundwater that was transported laterally from other areas may be occurring. Furthermore, this study, and the earlier work cited (Deverel
and Gallanthine, 1989; Deverel and Fio, 1991), suggestthat any modeling of tile drain flows in this area which
requires the assumption of a shallow impermeable layer
and no mixing of drained irrigation water with non-local
groundwater will have poor predictive capability unless
it can be independently confirmed that such mixing is
not occurring.
ACKNOWLEDGMENTS
We thank David G. Cone, Manager of the Broadview Water
District, for providing the irrigation and drainage data utilized
in this study and for consultations concerning many aspectsof water management in the district.