An Hermitian finite element solution of the two-dimensional
saturated-unsaturated flow equation
M. Th. V A N G E N U C H T E N

US Salinity Laboratory, 4500 Glenwood Drive, RiversMe, California 92501, USA
This paper describes a Gaierkin-type finite element solution of the two-dimensional saturatedunsaturated flow equation. The numerical solution uses an incomplete (reduced) set of Hermitian
cubic basis functions and is formulated in terms of normal and tangential coordinates. The formulation leads to continuous pressure gradients across interelement boundaries for a number of welldefined element configurations, such as for rectangular and circular elements. Other elements generally
lead to discontinuous gradients; however, the gradients remain uniquely defined at the nodes. The
method avoids calculation of second-order derivatives, yet retains many of the advantages associatedwith Hermitian elements. A nine-point Lobatto-type integration scheme is used to evaluate all local
element integrals. This alternative scheme produces about the same accuracy as the usual 9- or 16point Gaussian quadrature schemes, but is computationally more efficient.

INTRODUCTION
The frequent use of land for the disposal of a wide variety
of industrial, municipal and agricultural wastes emphasizes
the need for an accurate quantitative description of flow
and transport in the unsaturated zone. Initially, mostly
finite difference methods were used to solve the governing
flow equations, x-4 Finite element techniques became
available in the late 1960s and have since been applied to
two-dimensional saturated-unsaturated flow problems by
Neuman, s Reeves and Duguid, 6 Segol, 7 and many others.
The finite element solution presented here uses first-ordercontinuous Hermitian cubic basis functions. If properly
applied, this method generates continuous pressure headgradients across interelement boundaries, and hence produces a continuous flow velocity field over the entire
domain. Earlier work dealing with one-dimensional flow
and transport showed that the Hermitian approach generated better results than was possible with either standard
finite differences or linear finite elements, particularly withrespect to the material balance for the transport equation. 8
This study is limited to solutions of the flow equation only.
GOVERNING FLOW EQUATION

the time and x~ is the a-th spatial coordinate. The soil
water capacity C is0

c= ~s+~

Presenied at the Fourth International Conference on I:inlte Elements in Water Resources, Inst. f. Stromungsmechanik, Universityof Hannover, Fed. Rep. of Germany, 21-25 June 1982.
Received July 1982. Written discussion closes August 1983.

106

Advances in Water Resources, 1983, Volume 6, June

(2)

GALERKIN APPROXIMATION
Equation (1) will be solved with the Galerkin-type finite
element method. The dependent variable h is assumed to beof the formnh(x,~, t) ~- y Cj(x,O ai(t)

(3)

/'=1
where the ¢i are the selected basis functions and where theaj are the associated time-dependent coefficients. The basis
functions are used also to minimize the residual L(h) of
equation ( 1) by requiring thatf L [h(xa, t)] ¢pi(xa) dI2 = 0
I2L(h)=-Tx * K ~ - ~ + K.~ - C - - + Q = O a t (~,¢=1,2)

where L is a differential operator, h is the pressure head,K~¢ is the hydraulic conductivity, Q is a fluid source, t is

ahl

where 0 is the volumetric soil water content, S is the specific storage coefficient, S w is the degree of fluid saturationand e is the porosity.

The governing equation for transient saturated-unsaturated
flow in a two-dimensional vertical cross-section is taken as

(1)

OSw~

(4)

where ~2 represents the two-dimensional domain. Combining equations (1), (3) and (4), and applying Green's first
theorem to the spatial derivatives leads to[P] {X} + [R]

d{X}
dt

+ {F} = 0

(5)

where
{Xi} = {ai}(6a)0309-1708/83/020107-06 $2.00© 1983 CML PubLications

Hermitian finite element solution o f two-dimensional saturated-unsaturated flow equation." M. Th. van Genuchten[Pu] = J K~ -ax~- -axa d~2(6b)[Ri/]= ( C(9/(9i d~

(6c)

z

(~'z.3)

I2

12

1
~"~

•

in which I" represents the boundary of ~2. The parameter qn
in equation (6d) holds for applied flux-type boundary
conditions:qn(Xa, t) = --K~~ + Ka2 na

aILI(x 1 ,Z 1 )(-I ,-I)( I ,-I )Xob

Figure 1. Deformed isoparametric Hermitian elements inlocal {a) and global (b) coordinates(7)
where the basis functions ¢ are given by

where nc~ is the outward normal to the boundary F.

¢oo/= -~(1 + ~o) (1 + 777/0)(2 + ~ o + r/% -- ~: -- r~2) (1 2a)Basis functions
When equation (3) is applied to a one-dimensional element (Xa, Xb) and first-order continuous Hermitian basis
functions are used, its form reduces toh(~,t)= Y~ ~,o/(0~/(t)+~1/(~)---2-~ j
/=1

(8)

The local coordinate ~ in equation (8) is related to theglobal x-coordinate through the expression
2= -- 1 q - - (X --Xa)AX(9)where Ax = x b - x a is the nodal distance of the element.
The one-dimensional Hermitian cubic basis functions are¢oj = --~(~ + ~o)2 ( ~ o - - 2)

(Go = + 1)

(lOa)

~bl/= ~ ~o(~ + ~o)2 (~Jo -- 1)(Go = -+ 1)(1 0b)

As is the case with zero-order continuous basis functions,the two-dimensional Hermitian basis functions can begenerated from the one-dimensional expressions in at least
two ways. In one approach, the two-dimensional functions
are obtained by simply multiplying the one-dimensional
expressions with each other. This approach leads to four
unknowns per node:the function H, its spatial derivatives
aH/bx and aH/az, and the cross-derivative
where x = xl and z = x2 are the global coordinates. Unfortunately, this approach also leads to the derivatives
a2H/ax 2 and a2H/az ~ in the global matrix equation whenapplied to deformed non-rectangular elements. These higher
order derivatives must be retained in the formulation if
first-order continuity in the pressure head across interelement boundaries is required. When the higher orderderivatives are neglected, the gradients still will be continuous at the nodes, but generally not across element boundaries
between nodes.
An. alternative approach would be to eliminate the crossderivative from the formulation altogether, thus avoiding
the calculation of unnecessary quantities and reducing thenumber of unknowns per node from four to three. Thisapproach leads to the following expansion for the unknown
function over the local element of Fig. I : 9

a2H/Oxbz,h(~, r/, t) = ~4 (/=1

OHJo~mlaHii¢oo/Hi+¢lO/-'~-~'.-.+~OIj'-~'--/ (ll)~xo/-- -- ~ o ( 1 + / ~ o ) (1 + rrOo) (I -- ~2)(12b)¢Olj = ---~%(1 + ~o) (1 + ~?¢/o)(1 --772)(12c)

in which Go= +- 1 and % + 1.

Normal and tangential coordinates
The finite etement scheme presented here is formulated
in terms of normal and tangential coordinates (0, z; see
Fig. lb). The advantage of this method is that flux-typeboundary conditions are much easier to incorporate into
the numerical scheme. In addition, the study is restricted todeformed elements that have 90 ° internal angles. This
restriction could limit in some cases the applicability of thesolution, particularly when highly irregular boundaries arepresent, n However, because most field problems can be
described easily with elements having 90 ° angles, andbecause the formulation could be simplified, it was felt
advantageous to maintain the square angles. A different
formulation is needed for those few cases where thisrestriction cannot be met.
When normal and tangential coordinates are used, theunknown nodal coefficients become H and its derivatives
aH/ao and aH/az. Hence, approximation ( 1 1 ) m u s t be
rewritten in terms of these unknowns. From the chain ruleof differentiation we haveattj

all~ ao/aN~ aria~ao a~az a t(13a)0M/= 0M/ao/+ 0.; a

(13b)
at7
aa a~7a.r at7
Because the distorted elements have 90 ° internal angles, itfollows immediately that all derivatives ar/a~ and aa/ar/are
zero. The other nodal derivatives in equation (13) can be
expressed in terms of the geometric properties of theelement. For example, for node l of the deformed elementin Fig. 1b we have
Oal Ll-- = -a~
2

0~'l

L4

(14)
an
2
where L1 and L4 are the lengths of element sides 1 and 4,
respectively. Substituting equation (13) into equation (11)and using the information above yields,h(~, r/, t) = ~
j=l

-

a.,

a.q

ooiHi+~oi--+q~riaoar /Advances in Water Resources, 1983, Volume 6, June

(15)

107

Hermitian finite element solution of two-dimensional saturated-unsaturated flow equation." M. Th van Genuchten
4+5QS

with the new basis functionsao/f~a] = (~lO/ ~

o~

3X3GP

ari
f~r] = (~O1] - -(16)
0

The following isoparametric transformation is used for
distorted elements: ~l

c

0

o

C

.

i
= Z

e oojX +

+ em

!

0

(17a)

0

0

0 I

t;k

nk

Wk

j=l
t~k

z (~, 7?) -- Y~ ¢oojZj + ¢ o / - - + ¢,~/
/=l

ao

-~r]

rtk

wk

(17b) ±1

where the coefficients represent the nodal values of the x
and z coordinates and their normal and tangential derivatives:

-+ 1

1/9

+-.x/'O~

±~

25/81

0

10/9

0

±%/~

40181

0

+-~,/E4

10/9

-+x/~

0

40/81

o

0

--8/9

0

0

64/81

aZi= sin(e/)
ao

(lSa)

- cos(s/)(18b)

Figure 2. Location of integration points and values ofweighting coefficients for the nine-point Lobatto-type(4 + 5 QS) and nine-point Gaussian quadrature (3 x 3 GP)
schemes

Finally, the side lengths L i in equation (14) can beobtained by numerically integrating equations of the form?

scheme is slightly less than the error for 3 x 3 Gaussian
quadrature.a°aXi = cos(e/)
ao
ax; _

sin(a/)

arar

1(°z/'tL=f{~a~/+\~/~d~( ~,~o=+- 1) (19)--1

Numerical integration
Numerical techniques are used to evaluate the integrals
in equation (6). The integrations are carried out directly on
the local (~, 7?) element between the limits - 1 and + 1. For
this purpose, the area elements dFZ in equation (6) are
replaced by
dI2 = det(J) d~ d77(20)

where J is the Jacobian of the transformation. Several
methods are available to carry out the numerical integrations over each local element. Typically, an integral isapproximated by a sum consisting of the values • of thefunction at selected points in the local element and multiplied by appropriate weighting coefficients. For example,
for the function f(~, r/) we have11

wkf(~k, rTk) + E(21)dX] +at/2{x}t+At--/X}t(22a){X} t+At/2 --- w{x}t+At + (1 -- w ) l X } t(22b)where At is the time step and w a temporal weighting
coefficient (0 ~ w ~< 1). Defining matrix equation (5) at thehalf-time level and substituting equations (22a, b) into
equation (5) leads to the following set of equations:[A] t + A t l 2 { x } t + A t = [ B l t + A t / 2 { x } t + {F} t+At/'2 (23)where1

[AI = w[Pl + Att [R](24a)[B] = ( w - - l ) [P] + ~ t

[R]

(248)

k=l

where (~k, l/k) are the coordinates of the integration point,
wk is the associated weighting coefficient, and E is the errorof the approximation. Gaussian quadrature schemes havebeen used nearly without exception in related finite element
studies. Although Gaussian quadrature probably is the mostaccurate scheme for a given number of integration points,
this approach does not take full advantage of the properties
of Hermitian basis functions. These basis functions, including their derivatives, are zero on several or all corner nodes.
Hence, the computational effort can be reduced by locatingsome of the integration points at the nodes. Consequently,
an alternative nine-point Lobatto-type quadrature scheme
was used for all numerical integrations. Figure 2 shows thelocation of the nine integration points: four at the nodes
and five inside the element.'The integration error E for this108

Equation (5) defines a set of n ordinary differential
equations in n unknowns. A finite difference scheme wasintroduced to approximate the time derivative in the matrix
equation. Define for that purpose:1f(~, ~) d~ dl? = ~--1 --1

Integration in time

Advances in WaterResources, 1983, Volume 6, June

The iterative scheme used to solve equation (23)was exactly
the same as the scheme employed earlier for similar onedimensional flow situations, sBoundary conditions
Constant pressure head (Dirichlet-type) boundary conditions were included in two ways. First. all algebraic equations that are associated with constant boundary nodes
were eliminated from matrix equation (23). Second,
Dirichlet boundary conditions also affect those line integralsof equation (6d) that are associated with the normal andtangential gradients of the pressure head. To calculate theseterms, one must first evaluate the flux qn normal to theboundary. This procedure leads to equivalent nodal values
for the line integrals; these nodal values must be added to
the vector {F}.11

Hermitian finite element solution of two-dimensional saturated-unsaturated flow equation: M. Th. van Genuchten
0

A similar approach applies also to flux-type boundary
conditions. When the pressure gradient or the flux normal
to a boundary is specified, either aH]ao or all~at alongthat boundary will be known. This also leads to a reductionin unknowns. In addition, the line integrals for both thepressure head and one of its spatial derivatives will be nonzero for the nodes along the boundary.

20
~"

40

RESULTS
Some initial results obtained with the Hermitian scheme are
given below. The examples are presented as a verification of
the accuracy of the numerical scheme, and also to illustrate
the type of problems that have been solved thus far.

A

One-dimensional vertical infiltration"1 ~In this example, 12 water is allowed to inf'dtrate into a
deep homogeneous soil profile that has the following
hydraulic properties:
0.6829 -- 0.09524 ln(I h [) h <~ -- 29.484

O(h) =

0.4531 - 0.02732 ln(Ih I)

(25a)

2 hr.

b
~ ~/.I nitiol
Distribution

E

6O

.

8.
LO

80

I00

2 x 2 GP3 x 3 GP--29.484 ~<h ~<-- 14.495~

41h[ -3"409s h~<--29.484

t19.34x10(25b)--29.484 ~<h ~<-- 14.495

K(h) = ~516.8 Ih j-0.97814

120

0 (x, z, 0) =t--60~z~O
(26)

10.20

-- 125 < z < - - 6 0

The equivalent initial pressure profile follows immediatelyfrom equations (25a) and (26). The medium was divided
into 25 rectangular elements (Ax = 10 cm, Az = 5 cm)
with a total of 52 nodes. No-flow boundary conditions(qn = 0) were applied at the sides, and a constant pressure
head was imposed at the soil surface (h = -- 14.495 cm). A
free draining soil prof'tle was assumed at z = -- 125 cm:ah-- (x,-az

125, t) = 0.

(27)

Figure 3 presents calculated moisture distributions
versus depth for four different quadrature schemes. The
solid line represents the assumed correct solution and was
obtained with several one-dimensional finite difference andfinite element solutions that used extremely small time andspatial increments, s Relatively poor results were generated
with the four-point Gaussian quadrature (2 x 2 GP) scheme.
The calculated distributions not only lag behind the correct
solution, but also exhibit serious oscillations. Much better
results were obtained with the nine-point Gaussian quadrature (3 x 3 GP) scheme, although the computed distributions
also in this case lag somewhat behind the correct solution.
The most accurate results were obtained with the 16-point
Gaussian (4 x 4 GP) and the nine-point alternative (4 + 5 QS)
integration schemes. This last scheme actually generated
results that were slightly better than those obtained withthe 4 x 4 GP scheme. The 4 + 5 QS scheme nearly duplicatedthe correct solution, except for some minor oscillations
near the toe of the moisture front.
Table 1 gives a summary of the various computer runs
for this example, including a comparison of the execution
times needed to complete the 0.4-day simulation on an IBM
360191 computer. The table also shows results that werel.................4 x 4 G P

0000

4 + 5 QS (nodal values)
correct solution

where K is given in cm/day and h in cm. The initial moisture
content distribution is0.15 --z/1200

J-""

9hr

0.20

030

MOISTURE

040

CONTENT (cm3cn'?)

Figure 3. Effect of different quadrature schemes on computed moisture distributions

obtained with equivalent one-dimensional Hermitian finite
element solutions, a For example, the four-point Gaussian
quadrature scheme in two dimensions reduces to a two-point
Gaussian scheme (2 GP) in one dimension. Similarly, the4 + 5 quadrature scheme reduces to (and produces the same
results as) a five-point Lobatto integration scheme (5 LP).
The infdtration example shows that the 4 + 5 quadratqre
scheme generates results that are at least as accurate as
those obtained with 4 x 4 Gaussian quadrature (Fig. 2), butwith a saving in computer time of about 45% (Table 1). The
equivalent savings for the 5ne-dimensional programs are
considerably less.

Two.dimensional horizontal infiltration
This example considers the two-dimensional horizontal
infiltration of water into an initially air-dry, semi-infinite
soil slab. A detailed description of the physical problem isgiven by Rubin. 13 In summary, water infffltrates from theleft side of the rectangular soil slab under the influence

Table 1. Summary o f various one-dimensional (l-D) and twodimensional (2-D) simulations o f example 1. Total simulation timeis 0.4 days
Integration method
1-D

2-D

Total numberof time steps
1-D and 2-D

2GP
3 GP
4GP
5 LP

2X2GP
3 X 3 GP
4X4GP
4 + 5 QS

132
170
149
156

CPU time (s)
1-D

2-D

6.6
9.9
10.2
9.4

63.6
130.8
160.0
89.8

Advances in WaterResources, 1983, Volume 6, June

109

Hermitian finite element solution o f two-dimensional saturated-unsaturated flow equation: M. Th. van Genuchten
SOIL SURFACE7REHOVOT SAND

0.200
'0[

0150 OlO00.050

~

~,.

.....

\',

'°
F

\:,
"%-,.

\

',

,~ ', \ ~z/

R'UBIN(,968)....... z,.oLos,,, vr ,,.

\s.,-,.,o:"'"--I z° •

X,>. \ \ \ ',X,,~

Xo

Q=Qo8(X'Xo)8(z- Zo)
,.
o

0

2

4

6

8

I0

12

14

16

HORIZONTAL DISTANCE (cm)Figure 4. Calculated moisture distributions durin& infiltration into a horizontal soil slab

x

(cm)

of a constant total head (h + z) o f - - 13 cm. No-flow boundary conditions are present along the top (z = 0) and bottom( z = - - 1 0 c m ) boundaries of the system. Soil-hydraulicproperties of the Rehovot sand used in the experiment are
given elsewhere} a'14 Figure 4 compares the Hermitian
simulation results after 10 min of infiltration with results
obtained earlier by Rubin ~a using finite differences and by
Zyvoloski et al. ~s using linear finite elements. Note that theHermitian-based moisture distributions are very close to the
distributions calculated by Rubin} a
Subsurface irrigation
This example describes the infiltration of water from ahorizontal line source in a two-dimensional cross-section.
Water is applied through buried laterals (point sources) thatare located at a depth of 15 cm below the soil surface(Fig. 5). The experiment was used earlier by Thomas et al. 16
to compare calculated and measured steady-state pressure
head distributions in a soil box designed to model a subsurface irrigation system. The calculated distributions werebased on an analytical Solution developed by Zachmann
and Thomas. 17 This solution holds for a semi-inf'mite soil
medium with no-flow boundary conditions at x =-+ Xo.
Here we will compare the analytical solution with Hermitian
finite element results using exactly the same experimental
conditions as described by Thomas et al. ]6
Figure 5 shows schematically the finite element meshfor the line source problem. The soil slab was extended
arbitrarily to a depth of 350 cm, thereby assuming that thesemi-infinite system was closely approximated by this deep
profile. The source strength Q equals 1.05 x 10-3 cma/cm/s,
which is equivalent to a constant application rate of 0.7436
cm/day when distributed evenly over the entire soil surface.
The hydraulic conductivity K (cm/day) of the clay loam
soil isg = 96.768 exp(--0.1258h)(h >t 0)(28)

Two methods can be used to evaluate the steady-statepressure distribution. One approach would be to calculate
the distribution directly by forcing the time-derivative inequation ( I ) to be zero. This method often requires numerous iterations and could at times even lead to convergenceproblems, especially for highly non-linear problems. The
direct method was followed first. However, several numerical adjustments were required in the program, such as an

110

Advances in Water Resources, 1983, Volume 6, Juneb

Figure5. Schematic representation of the line source
problem and the associated finite element mesheffective under-relaxation, to speed up convergence andhence to limit the number of iterations. Alternatively, theproblem can be solved indirectly by starting out with anarbitrary initial condition and solving the problem as if itwere a transient one. A fictitious soil water retention curve
was used for this purpose. This last approach proved to bemuch easier and less costly.
Calculated steady-state pressure head distributions areshown in Fig. 6. The curves duplicated exactly the analytical results obtained with the solution of Zachmann andThomas. 17 The pressure head deep in the profile was foundHermitian finite element solution o f two-dimensional saturated-unsaturated f l o w equation." M. Th. van Genuchten(

at z = -- 100 cm. Note that the pressure distributions close
to the soil surface are nearly identical to those shown in
Fig. 6.

2(
CONCLUSION
4(

The three examples shown in this study demonstrate thatthe Hermitian finite element solution can be used to calculate pressure head (and moisture content) distributions fora number of two-dimensional problems. Application of thealternative 4 + 5 Lobatto-type quadrature scheme resultsin considerable computer savings as compared to the use ofGaussian quadrature. Because the Hermitian solution generates a continuous flow velocity field for most finite element
configurations, it is expected that the scheme will be especi:
ally attractive when combined with numerical solutions of
the solute transport equation.

T
n

w¸

8(

t:)

..J

5
to

I0(

12(

REFERENCES
14(

-60

-40

-20

0

20

40

60

HORIZONTAL DISTANCE (cm)Figure 6. Calculated steady-state pressure head distribution
resulting f r o m a buried line source at z = - 15 cm (semiinfinite profile)

0

~

f

i

I

i

i

i

i

i

40
-

-J 60

~
-3()

80

-zo

-I0
I00.
'

"

.

. - 2 0.

.

.

0

2'o

'

4O
'

°'

6O

HORIZONTAL DISTANCE (cm)Figure 7. Calculated steady-state pressure head distribution
resulting from a line source at z = - 15 cm (fixed watertable at --100 cm)

to be --38.7 cm, and was independent of the horizontal
coordinate. The value o f - - 3 8 . 7 cm is consistent with thepressure head at which the hydraulic conductivity reaches
a value of 0.7436 cm/day, the latter being the application
rate as averaged over the soil surface. As a modification of
the last example, Fig. 7 shows calculated steady-statepressure head distributions resulting from the same line
source as before, but now with a fixed water table present

1 Hanks, R. J. and Bowers, S. A. Numerical solution of the moisture flow equation for infiltration into layered soils, Soil Sci.
Soc. Am. Proc. 1962, 26, 530
2 Freeze, R. A. The mechanism of natural groundwater recharge
and discharge. 1. One-dimensional, vertical, unsteady, unsaturated flow above a recharging or discharging groundwaterflow system, Water Resources Res. 1969, 5, 153
3 Brandt, A., Bresler, E., Diner, N., Ben-Asher, I., Heller, J. andGoldberg, D. Infiltration from a trickle source. I. Mathematicalmodels, Soil ScL Soc. Am. Proc. 1971, 35, 675
4 Pikul, M. F., Street, R. L. and Remson, I. A numerical model
based on coupled one-dimensional Richards and Boussinesq
equations, WaterResources Res. 1974, 4, 1069
5 Neuman, S. P. Saturated-unsaturated seepage by finite elements, J. Hydraulics Div., ASCE 1973, 99, HYI2:2233-2251
6 Reeves, M. and Duguid, J. O. Water movement throughsaturated-unsaturated porous media: a finite-element Galerkinmodel. ORNL-4927, Oak Ridge National Lab., Oak Ridge,
Tenn., 1975
7 Segol, G. A three-dimensional Galerkin-f'mite element modelfor the analysis of contaminant transport in saturated-unsaturated porous media. In: Finite Elements in WaterResources,
Gray, W. G., Pinder, G. F. and Brebbia, C. A., eds., Pentech
Press, London, pp. 2.123-144, 1977
8 van Genuchten, M. Tit. A comparison of numerical solutions of
the one-dimensional unsaturated-saturated flow and masstransport equations, Adv. Water Resources 1982, 5, 47
9 van Genuchten, M. Th., Pinder, G. F. and Frind, E. O. Simulationof two-dimensional contaminant transport with isoparametric
Emite elements, Water Resources Res. 1977, 13, 451
10 Gray, W. G. and van Genuchten, M. Th. Economic alternatives
to Gaussian quadrature over isoparametric quadrilaterals, lnt.
J. Num. Methods Eng. 1978, 12, 1478
11 Frind, E. O. An isoparametric Hermitian element for the solution of field problems, lnt. J. Num. Methods Eng. 1977, 11,945
12 Warrick, A. W., Biggar, J. W. and Nielsen, D. R. Simultaneous
solute and water transfer for an unsaturated soil, Water
ResourcesRes. 1971, 7, 1216
13 Rubin, J. Theoretical analysis of two-dimensional transient
flow in unsaturated and partly saturated soils, Soil Sci. Soc.
Am. Proc. 1968, 32, 607
14 Rubin, J. Numerical analysis of ponded rainfall infiltration. In:
Water in the Unsaturated Zone, Rijtema, P. E. and Wassink, H.,
eds., Wageningen, Netherlands, Vol. I, pp. 440-451, 1969
15 Zyvoloski, G., Bruch, J. C., Jr. and Sloss, J. M. Solution ofequation for two-dimensional infiltration problems, Soil Sct
1976, 122, 65
16 Thomas, A. W., Duke, H. R., Zachmann, D. W. and Kruse, E. G.
Comparisons of calculated and measured capillary potentials
from line sources, Soil Sci. Soc. Am. J. 1976, 40, 10
17 Zachmann, D. W. and Thomas, A. W. A mathematical investigation of steady infiltration from line sources, Soil ScL Soc.
Am. Proc. 1973, 37, 495

Advances in Water Resources, 1983, Volume 6, June111