Simple Method for Predicting Drainage from Field Plots1J. B. SISSON, A. H. FERGUSON, AND M. TH. VAN GENUCHTEN
9 — 0(z,t) is the volumetric moisture content,t = time,z = depth (positive downward),
K = K(0) is the hydraulic conductivity, andH - H(6,z) = h(6)-z; i.e., H(hydraulic head) =h (pressure head) — ^(gravitational head).

ABSTRACT
When the one-dimensional moisture flow equation is simplified by applying the unit gradient approximation, a first-orderpartial differential equation results. The first-order equation
is hyperbolic and easily solved by the method of P. 0. Lax.
Three published K(6) relationships were used to generate three
analytical solutions for the drainage phase following infiltration. All three solutions produced straight lines or nearlystraight lines when log of total water above a depth was plottedversus log of time. Several suggestions for obtaining the required
parameters are presented and two example problems are included to demonstrate the accuracy and applicability of themethod.

When a unit gradient in the total head H is assumed,
9H/3z = -1, Eq. [1] becomes[2]dz
When only the drainage phase is considered, Eq. [2]may be solved subject to the conditions

Additional Index Words: Cauchy problem, redistribution,
characteristic value problem, hydraulic conductivity, infiltration.

16C if z < 0"<*•"> =Sisson, J. B., A. H. Ferguson, and M. Th. van Genuchten. 1980.
Simple method for predicting drainage from field plots. Soil
Sci. Soc. Am. J. 44:1147-1152.

is responsible for the annual destruction of thousandsof hectares of cropland by contributing to the formation and growth of saline seeps (Brown and Ferguson,
1973; Ferguson et al., 1972). In some irrigated areas,
subsurface drainage water contributes to river pollution (Wierenga and Patterson, 1972).
Estimation of the rate and quantity of drainage
water contributing to these problems is essential tofinding feasible solutions. But estimating these variables requires predicting the hydrologic behavior oflarge areas and frequently the characterization ofmany soils under field conditions. This necessitates
the use of simple, yet accurate, models which containparameters that can be obtained on site as quickly
as possible.
This paper considers a special class of models basedon the assumption of a unit gradient of the total potential head. Several studies (cf. Black et al., 1969,and Davidson et al., 1969, among others) have shownthat a unit gradient often exists during the redistribution and drainage phases when a uniform profile
is draining freely in the absence of a shallow watertable. Three solutions are presented, each based upona different conductivity equation. Two example problems are further included to demonstrate the use of
the present approach.

dz
dK/d0

<«

_
Q'

[4]

The right hand expression implies that d0 must be
zero or 9 remains constant for certain values of z andt. Since 0 is constant, dK/dO will be constant and thefirst two terms in Eq. [4] can be integrated to give

z =

dK

•t = At

[5]

where
A =

dK
de

These results imply that a set of curves (i.e. characteristics) propogate from the initial condition Oi- If

THEORETICAL CONSIDERATIONSThe equation for predicting the one-dimensional
flow of water in porous materials is (Taylor and Ashcroft, 1972):
937= U if* >owhere the subscripts c and m denote minimum andmaximum obtainable values, respectively.
When the profile is saturated initially, Om equals
the moisture content at saturation, 0S. In its generalcase, however, 6m may be less than 0S; for example,
following an irrigation with a flux less than the saturated hydraulic conductivity.
The initial value problem given by Eq. [2] and[3], also known as a Cauchy or a characteristic value
problem, has been the subject of many studies inmathematical and engineering literature (cf. Lax, 1972;
Aris and Amundson, 1973, where Aris and Amundson
present examples and are an excellent introduction
to this subject). The characteristics of Eq. [2] areobtained by solving the following system of ordinary
differential equations written in standard form.

ATER DRAINING from soil profiles is an importantW
factor in many contemporary and environmental problems. In the Northern Great Plains this water

80
3T

W

1
Contribution from Rockwell International, Rockwell Hanford Operations, Richland, Wash. Document RHO-SA-138. Partof a thesis submitted by senior author in partial fulfillment of
the M.S. Degree at Montana State University. Received 13 Nov.
1979.
Approved 6 Aug. 1980.
2
Senior Soil Physicist, Rockwell International, P.O. Box 800,
Richland, WA 99352; Professor, Agronomy Department, Montana State Univ.; and Soil Physicist, U.S. Salinity Lab., Riverside, Calif., respectively.

[1]

where
1147

1148

SOIL SCI. SOC. AM. J., VOL. 44,

1980

9, cmVcm 3

REGION WHERE ~~~——x CHARACTERISTICS ARE 'K GIVEN BY LAX METHOD -o
NTFig. 2—Typical behavior of Lax solutions with a monotonicallyincreasing K(e). Points denoted by ta and ta are also denoted in Fig. 1.
Fig. 1—Characteristics defined by Eq. [5]. Solid lines fillregion where water content is unchanging and dashed lines

fill region where desorbtion is taking place.

one could travel along one of these curves with an
instrument that sensed moisture content, the instrument readings would not change.
The characteristics of Eq. [2] subject to [3] areshown in Fig. 1 as parallel solid lines propagatingfrom t—0 and lie in the region given by z > At. The
wedge-shaped region defined by 0 < z < At and shownin Fig. 1 as being filled with dashed lines is devoid ofcharacteristics when the usual mathematical methods
are applied. Lax (1972) studied Eq. [2] and showedthat unless 0i(z) is a monotonically increasing function of z then no continuous solution exists in thisregion (0 < z < At). If discontinuous solutions wereallowed Lax demonstrated that several distinct solu-

tions exist to Eq. [2] with the same initial condition[3]. Lax noted that only one solution exists in theregion (0 < z < At) such that the characteristics do
not intersect (Note: Intersecting characteristics indicate the development of a shock wave or a mathematical discontinuity). The Lax solution is given byd(z,t) =f(z/t)
nfor 0 < z < At

for z > At

[6]

where the function / is defined such that..dAT,It should be noted that the values of Oc (i-e., the[7]provided thatd*K/de2 ?& 0 for all 0.

increasing functions K(9) is schematically shown in
Fig. 2 where the desorption of a uniform profile during the early drainage phase is considered (the desorption zone is still confined to the upper part of the
profile, i.e., z < At). The lower boundary of the desorption zone is given by Eq. (5), with Oi = Om- Thisboundary propagates downward at a constant rate A= dK/dO\<im- The method for obtaining solutions canbe summarized in the following three-step algorithm
based on the Lax solution scheme:
1) Differentiate a chosen K(9) with respect to 0 obtaining dK/de.
2) Replace the resulting dK/dd term with z/t.
3) Solve the resulting expression explicitly for 0 =f(z/t) = 6(z,t).
Three frequently used K(Q) relationships found incurrent literature are presented in Table 1 with all
steps in the Lax algorithm. Similar solutions havebeen published previously for the K(0) relationships
referenced as Davidson et al. (1969) and Watson(1967). [See Davidson et al. (1969) and Gardner etal. (1970), respectively.] But unlike the previous solutions, found by methods apparently analogous to
the separation of variables, the Lax solutions have no
arbitrary constants arising from integration. More important, however, the Lax solutions are valid for alltimes t ^ 0, i.e., also for the early drainage phase, aphase which has to date been largely ignored.

[8]

Note here that 6 defined by Eq. [6] is continuous
throughout z but that the derivative of 9 with respect
to z is two valued at z = 0 and z = At, depending
whether the limit is taken from the "left" or the"right." In general, if Eq. [8] is violated Eq. [6]becomes multivalued, or imaginary, imaginary solutions are easily obtained by choosing K(0) as a cubic
polynomial, solving Eq. [6] and determining valuesof the polynomial coefficients where Eq. [8] is violated.
The behavior of the Lax solution for monotonicallyminimum obtainable value for the moisture content)
in models of Watson (1967) and Davison et al. (1969)(Table 1) are assumed to be zero. Although this does
not present any problems for the Watson model, themodel based on the exponential conductivity curve
will exhibit a mathematical peculiarity in this respect.
Because both K and its slope dK/dO are non-zero
when 6 = 0, another characteristic will emanate fromthe region (z ^ 0, t — Q) describing the downward
propagation of the zero water content. The positionof zero moisture content is again given by Eq. [5],whereby it is understood that 6\. is now zero. Although
this particular case is physically unobtainable (since
no one has observed a dry layer of measurable depth
immediately following an irrigation), the complete
and correct solution is

1149

SISSON ET AL.: SIMPLE METHOD FOR PREDICTING DRAINAGE FROM FIELD PLOTS

Table I—Algorithm steps for obtaining Laz's solution for the given hydraulic conductivity (K) water content (9) relations.

Steps to the Lax solution
Step IfdKdeztDavidson et al(1969)

dKz

Brooks & Corey(1964)

dK

Source of K(6)Watson (1967)

0

Step 2deStep 3z v«i-fl

t
z

de

t

t Subscripts m and c denote maximum and minimum obtainable values, respectively.
% Parameter A defined with Eq. [5].
0for 2 < dK/d6\e=o-t6(z,t) = 0m + -ln(z/At)0m

o-t < z <At.

[9]otherwise

PARAMETER ESTIMATION
Any solution of Eq. [2] involves parameters thatmust be estimated from experimental data. While 0m
and Km can be estimated in principle from direct observation during the latter stages of infiltration (inthe case of ponding), or from measurements on soilcores, the exponential parameters are probably mosteasily obtained from nonlinear regression of field datataken sometime after the entire profile is undergoingdesorption. Although this procedure holds strictly only for uniform profiles, it may be extended tosituations wherein the profile is reasonably uniform
or exhibits some weak layering. This may be done byconsidering the total water, W, above a depth, z:

W =

[10]

This function is generally smoother with depth andhence more amenable to curve fitting than the moisture contents.

Results of forming W from the solutions in Table
1 (see Table 2) are rather surprising since they all
yield nearly straight lines when W vs. t is plotted on
log-log paper (provided 0C/0S is small). This straightline relationship was apparently first noted by Ogata
and Richards (1957). Wilcox (1959) fitted parameters
to a wide range of soils and also found that the log
W-log t relationship provided a reasonable approximation of the drainage process. However, Wilcox also
reported that the exponential parameter generally
decreased with depth. It is easily verified that the
total water equations in Table 2 predict the Wilcox
anomaly exactly by ignoring the early drainage phaseand using times 1 and 10 to obtain the exponents inthe Watson model. The resulting exponents will appear to decrease with depth.
Wilcox (1959) modified the original Ogata-Richards model by writing the flux at a depth explicitly
as a function of the total water above that depth, andused the resulting expression to estimate deep drainage losses. For the conductivity equations presented
here, the Wilcox flux model results in expressions of
the form K = K(W) (see Table 2), where K(W) was
obtained by eliminating the term z/At from W(z,t)
and K(z,t).
To estimate the exponents required by the models
in Table 1, z/At was eliminated from W(z,t) and 0(z,t)
and the results given in Table 2. When total water
above a depth is treated as a function of 0, it is ap-

Table 2—Total water (W) above a depth as W(z, t) and W(z, 9) and hydraulic conductivity (K} as K(z, t) and K( W).

WatsonSource of K(0)Davidson et al.

Brooks & Corey

W=W(z,t) =Exponent parameter in termsof discernible quan titles t
W=W(z,6) =- vn/U-n)( At I
wm

AW

AW

wm-we

z(6-lla)nz6c + (l-n)z0

K = K(z,t) =K = K(W) =f See text for definition of A W.

Kmexp(l-a(Wm-W)lz)

W-WC

1150

SOIL SCI. SOC. AM. J., VOL. 44,

1980

100,0001—•ff

100

1000

10

100
o

jf
1.0

10

0.1

0.2

0.3

0.4

0.6

0.5

9, cm3/cm3
Fig. 4—Moisture release curve for Glendale clay loam.

.01

absence of evaporation. The experimental setup of
the column and its soil is described by Dane and Wierenga (1975). The initial and boundary conditions
for the numerical model are:

.001
0.2

0.3

0.4

0.5

0.6

9, cm3/cm3
Fig. 3—Hydraulic conductivity K(6) for Glendale clay loam.

parent that the exponents can be estimated from simple regression of average total water above a depth
versus water content at that depth. The true implication is that parameters can be obtained while conducting the usual evapotranspiration studies without
doing an actual infiltration drainage study. Once theexponent is estimated and Km estimated from infiltration rates, drainage from cropped profiles can be
estimated. The added advantage of this type of analysis is that variation among soils from plot to plot canbe compensated for by using a drainage estimate foreach individual plot.
To obtain a check on the value of exponents computed, the total water equations in Table 2 were
evaluated at z = At and exponents solved for explicitly in terms of physical measurements. The quantityof water lost between t = 0 and the time when desorption starts at a depth, i.e., Wm — W(z = At, t), is denoted as AFF and the results given in Table 2.

0(z,0) = 0m 0 < z < 150 t = 0f^ = 08zz = 0^1± = -1z

3z

Example 1
In this example a large (150 cm long) soil column,
uniformly filled with Glendale clay loam and initiallysaturated (0S = 0.52), is allowed to drain freely in the= 150 t0

The hydraulic functions K(0) and h(0) of the soil are
given in Fig. 3 and 4, respectively. The hydraulic conductivity function, K(Q), is well described by the
Brooks and Corey model (Table 1), i.e.:Co

A 94fi\ 4.25

^27T^)[12](0m = 0.52; 0C = 0.246; n = 0.2353). Figure 5 presents
results obtained with the numerical (finite differenceCSMP) solution of Eq. [1] and [11] (Dane andWierenga, 1975) with results obtained with the unit
gradient method. The moisture distributions based
upon the latter method are obtained directly by following the three-step algorithm as outlined previously(see also Table 1):
1. Obtain dS:/d0. Hence from [12]:0 - 0.246\ 3.25
0.274 jAPPLICATIONS
Two example problems are now considered to demonstrate the use and accuracy of the unit gradientmethod. The first example compares the unit gradientmethod with a numerical solution of the complete
flow equation [ 1 ], while the second example demonstrates the application of the method to an actual field
experiment.

t > 0

[13]

2. Equate the resulting expression to z/t,t

= 155110 - 0.246\ 3.25
0.274 }

[14]

3. Solve for 0,
0 = 0.246 + 0.02858 (z/«)°-3077 (0< z < 1551t).
[15]
Figure 5 demonstrates that Eq. [15] generates results which compare well with those obtained with
the numerical solution of Eq. [1], except near the

SISSON ET AL.: SIMPLE METHOD FOR PREDICTING DRAINAGE FROM FIELD PLOTS
9, cmVcm3
.30

1151

40

.40

25
25 _

9c

15

10

50

7
u

Kf

5

75

.01

.10

1.0

10

100

t (days)100

125 NUMERICAL
150

Fig. 5—Comparison of unit gradient solution with numericalsolution for predicting water contents.

surface where the approximate solution approaches
6C — 0.246 at z = 0. The largest differences between
the two methods occur at intermediate drainage times(between 0.3 and 2 days), the maximum deviation(for z > 25 cm) being about 0.01 cm3/cm3. It is doubtful whether such differences are important when drainage from deep profiles is being estimated.
Figure 6 further compares the amount of drainable
water (W — 150- Oc) as obtained with the numericaland the unit gradient solutions. The solid line wasobtained by integrating Eq. [15] over the entire profile, leading to the expression (see also Table 2):
W(150,t) = 36.90 + 15.318 t~OMfn (t > 150/1551) [16]Because of the restriction z < At on Eq. [15], Eq.
[16] does not hold for the early drainage phase whenpart of the profile is still saturated. For this period
W(l5Q,t) is given by
fzi
f !50
W(150,i),t) = I6dz + I
9mdz
«/0
*/

ri50
6mdz - Kmt = 78 - 100 t (t < 150/1551).
Jo[17]
The use of Eq. [17] constitutes an important improvement over earlier drainage models in that it accurately describes drainage at small times, as demonstrated in Fig. 6. Overall, the unit gradient method
provides results which approach the numerical predictions; although, the method seems to overpredictthe drainage rate somewhat, especially at intermediatetimes (0.10 < t < 0.30).
Thus far in this example, the unit gradient method
was used to predict moisture distributions and amountsof drainable water in the profile. A demonstrationof how the method can also be used to analyze drainage equations like those shown in Table 2 will now
be carried out. Suppose the total water curve in
Fig. 6 given by the CSMP solution is given and one
knows Oc- For illustrative purposes, 0TO and 6C are
assumed to be the same as before (0.52 and 0.246, respectively). Hence the drainage data can be plottedin the same way as before (log (W — 36.9) vs. log t).
Ignoring the early drainage phase, the dashed line in
Fig. 6 is given byW(l5Q,t) = 36.9 + 16.60 r"-303 (t > 0.3 days).
[18]

Fig. 6—Comaprison of numerically predicted total drainable
water above 150 cm with unit gradient predictions. During
early drainage phase (times less than value shown with
open circle) only part of the profile is undergoing desorption.

By comparing Eq. [18] with the total water equation
based on the Brooks and Corey model (Table 2), it
is easily shown that n equals 0.2353 and A equals
1249.2. Hence from the definition of A = Ks/n(9m —Oc), one obtains Ks = 79.6. This value for the saturated hydraulic conductivity is certainly a reasonableapproximation of the "correct" value, i.e., of Ks = 100
as originally employed in the numerical solution of
Eq. [ 1 ]. Using the value of m, the complete conductivity equation is (see Table 1)- 0.246 \ 4 - 30
[19]

*» =

which is very close to the original curve shown in Fig.
3. Equation [19] in turn can now be used to predict
moisture distributions within the 150-cm profile, as
discussed before.
Example 2This example considers drainage from a field site
containing Gerber silty clay loam, one of the dominant soils of a saline seeped area on the Northern
Great Plains. Some physical characteristics of this soil
are given in Table 3. The values given are representative of the glacial soils on which many of the saline
seeps occur. The experimental conditions of this example are extensively discussed by Sisson (J. B. Sisson,
1972. Hydraulic properties of the Gerber soil. M.S.
Thesis, Montana State Univ., Bozeman) and are briefly summarized in Table 3.
In the fall of 1970, two field plots on the Gerbersoil were established, and neutron access tubes were
installed to allow for an over-winter stabilization
period. The plots were irrigated in the spring of 1971,
one for 4 hours and the other for 24 hours, and covered with plastic to prevent evaporation. Moisture
contents were determined at 15-cm intervals, at first
Table 3—Physical properties of Gerber silty clay loam.
Depth

Texture

cm
0-20

20-36
36-53
53-74
74-94
94-122
122-198
198-244

sicl
sic
sicl
sicl
sicl
cl
cl
cl

Bulk density
ofpedsWater content
0.3 bar

15 bars

g/cm-

g/g

g/g

1.70
1.86
1.59
1.59
1.59
1.62
1.86
1.86

31.5
33.7
28.2
32.1
28.0
26.7
24.8
26.0

14.0
20.2
15.0
14.7
13.2
12.3
11.4
10.3

1152

SOIL SCI. SOC. AM. J., VOL. 44, 1980

on a daily basis, later every 2 or 3 days. The observed
values were assumed to be representative for each
15-cm depth increment. Depth to the water table was> 13.8 m.
The drainage equations for this particular experiment were based on the Watson model (see Tables
1 and 2), Table 4 presents observed and calculated
total water contents, W(z,t), for two times. The drainage equations used are:
1 026

W(z,t) = 0.320 z -

f

[20]

and

W(z,t) = 0.312 z1-033 t -0.033[21]

Both equations were obtained by fitting data to the
Watson model, using a multiple regression program.
The coefficients in Eq. [20] and [21] were based on
soil water content data for z < 105 and z < 180 cm,
respectively, and utilizing data for days 2 through 10
after irrigation only. The greatest difference between
observed and calculated values is about 1 cm. Equation [20] appears to give somewhat better results
than Eq. [21], especially in the upper part of the profile. This is probably a result of the textural changes
which occur at about 100 cm (Table 3). It is emphasized here that the predictions assume a uniform profile, so that some inaccuracies may ge expected whensoil horizons with different hydraulic properties are
lumped into one profile, which is then assumed to beuniform. Notwithstanding the apparent nonuniformity of the profile, the correspondence between measured
and predicted values after 54 days in Table 4 is certainly acceptable for most purposes where drainage
has to be estimated for large areas using a simple
mathematical model.
SUMMARY
Although only two examples were used for assessing
the accuracy and applicability of the prediction equations, it does appear that the unit gradient method
holds considerable promise as a means for estimating
long-term drainage from soils. The unknown parameters are easily determined, whether by direct observation or by regression analysis on observed data,
and the determinations based on shallower depths canbe used to predict drainage from deeper depths, provided the profile is reasonably uniform.

Table 4—Measured and predicted total water above a given
depth in Gerber silty clay loam.
Total water above a depth
Measured

Predicted day 54

Initial
Depth

day 2

Day 54

Based on180 cmt

75
90
105
120
135
150
165
180

26.8
32.0
37.4
42.6
47.9
63.4
59.0
64.8

23.4
28.1
33.0
38.0
43.1
48.2
53.8
59.3

24.2
29.2
34.2
39.2
44.2
49.3
54.3
59.4

t Prediction equation: W(z,t) = 0.320 z'°M r°".
t Prediction equation: W(z,t) = 0.312 z'-T-"3.

Based on105 cmt

23.6
28.6
33.5
38.5
43.4
48.4
53.4
58.5