Axisymmetric infiltration in soils, I. Numerical techniques for solution

Journal of Hydrology, 23 (19741,111--130 © North-Holland Publishing Company, Amsterdam -- Printe¢ ~ in The Netherlands

AXISYMMETRIC INFILTRATION NIQUES FOR SOLUTION

IN SOILS~ L NUMERICAL

TECH-

R.W. JEPPSON

Utah Water Research Laboratory, Utah State Trniversity, Logan, Utah (U.S.A.) (Accepted for publication January 16, 1974)

ABSTRACT Jeppson, R.W., 1974. Axisymmetric infiltration in soils, I. Numerical techniques for solution. J. Hydrol., 23: 111--130. Numerical solutions are obtained to the transient, three-dimensional axisymmetric, unsaturated moisture movement through bomogeneous soil due to infiltration over a horizontal circular surface area. The paper is divided into two parts. The present part deals with techniques for numerical solutions, the difficulties in achieving such solutions due to the strong nonlinearities of the equation of flow, and procedures which minimize such difficulties. In this part three adaptations of the Crank-Nicolson method, and three variations of the alternating direction implicit (ADI} method are studied and their solations compared. Only those methods which evaluate the nonlinear coefficients on the advanced time planes fully implicitly yield solution in close agreement. When dealing with the CrankNicolson adaptations, fully implicitly implies solving a system of nonlinear equations by the Newto,,1-Raphson iteration, and in the ADI methods the coefficients causing the nonlinearities must be iteratively corrected a number of times. A single predictor is not adequate. One valid adaptation of t~e Crank-Nicolson method produces a system of nonlinear equations, which is much easier to solve than the other, and the third adaptation is valid only for constant coefficients, and it produces solutions grossly in error, and thus illustrates the influence of the nonlinearity on the equation of flow. Part II (Jeppson, 1974} compiles solutions from a number of problems in graphs relating some items of interest in infiltration to such problem specifications as size of circle of application and parameters describing the hydraulic properties of unsaturated soils.

iNTRODUCTION

For solving unsaturated moisture movement problems the use of numerical methods, prominent among which are finite differences, has experienced rapid development during the past decade as the high speed digital computer ha~ become generally accessible. The publication of the recent book by Remson et al. (1971) suggests that, if not already, the future schooling of groundwater hydrologists, soil physicists, engineers with specialities in water resources, in~gation and drainage, or soils, and hydrogeologists at the graduate level will Include substantial training in numerical methods for solving initial-boundary-value problems associated with partial differential equations. Yet as Brutsaert (1971) points out~ all such solutions to unsaturated two-dimensional

112 problems, with only some unpublished exceptions prior to his paper, are restricted to situations in which sharp wetting fronts are not present and, consequently, functions vary smoothly in both space and time. The equation of flow becomes strongly nonlinear when applied to problems with the hydraulic properties of real soils under unsaturated conditions. A number of authors have eluded to numerical difficulties in solving this equation. In studying the transient saturated--unsaturated flow in a rectangalar region, Verma and Brutsaert (1970) found the performance of common implicit methods at best marginal and unacceptable as more of the region becomes unsaturated. Yet because of the lack of accuracy due to the first order approximation of the time derivative, and stability considerations, explicit methods have never been considered seriously. In his model of a groundwater basin, Freeze (1971) notes numerical difficulties due to oscillations in the predicted values of capillary pressure used to evaluate coefficients which he temporarily assumes are not functions of the unknowns in order to solve a linear system of finite difference equations. In solving the one-dimensional infiltration problem, Smith and Woolhiser (1971) note similar oscillations. In addition to these difficulties which are associated with convergence of an iterative solution method, the nonlinearities can cause scatter in such solution results as infiltration curves, may require that ad hoc procedures be adopted such as reducing the time step (or space increments to give convergence of the iterative method) or to compute the time increment from moisture accumulations, as in the pioneering work of Hanks and Bowers (1962) instead of establishing its size as is commonly done for independent variables, and worst of all can result in solutions grossly in error if inappropriate finearizations are used. Undoubtedly, many workers have encountered difficulties that in some cases prevented solutions, which are no~ mentioned in published or unpublished literature° FORMULATION OF PROBLEM The equation of flow describing moisture movement through soils is obtained by substituting Darcy's law into the differential form of the conservation of mass equation. In simplified form this equation is: aS v. (KVh) = n

(1)

at in which h is the hydraulic head and is the sum of elevation and pressure heads (L)~ K is the hydraulic conductivity (L/T), ~ is the soft porosity, S is the soil saturation and t is time {T). The following assumptions are required to a~vive at eq.lo ( i ) Darcy~s law is valid for saturated and unsaturated flOWo (2) The gas flow occurs under such small gioadients compazed to the water flow that it c a l be ignored.

113 (3) The fluid (water) is incompressible and, consequently, of constant; density. (4) The solid particles do not move or consolidate, and consequently ~ is constant. {5} On a macro-scale, the functions which describe the flow and their derivatives,are continuous so that the differentialform of the continuity equation is valid. The space coordinates and the hydraulic head in eq.l, as well as other length variables,will be nondimensionalized by dividing by a characteristiclength ~. The required functional relationship of K and S to Ph {the pressure head) in eq.l will be defined by the Brooks-Corey equations {1966), in this paper, even though other versions of the computer program have been developed which utilizetabular values of saturation versus capillarypressure and the Burdine theory (1953} for evaluating K. Reasons for selectingthe BrooksCorey equations are that good fitsof data for some soilsunder imbibition, as well as drainage, are achieved with these equations and the hydraulic properties of the soilsare completely defined by values of only three parameters, an important consideration in part I][{Jeppson, 1974) in establishingrelationships of these parameter values to pertinent solution features. From considerations of similitude,the bubbling pressure head Pb/Pg (one of the parameters in the Brooks-Corey equations) has been suggested as the appropriate length variable for nondimensionalization (see Corey and Corey, 1967; Brooks et al., 1971 }. For infiltrationproblems the physical interpretation of bubbling pressure as the maximum negative pressure at which the soil first begins to desaturate, may not be meaningful. Rather Pb is a parameter to give the best fit of the data to the equation. Nonetheless, in this paper ~ has been taken equal to Pb/Pg. The mathematics would be only slightly different if ~ is a different value than Pb/Pg. In terms of the dimensionless pressure Pt = - P / P b , the Brooks-Corey equa~ons are: S- Sr

S~ -

1 -

I- Sr

(2)

p~

and: I ~:~

=

2+3~

(3)

Pt in which Sr is the residual saturation, Se is the effective saturation, A is the pore size distributionexponent, and Kr is the relativehydrau!ic conductivity, which when multiplied by the saturated conductivity K0 equals the hydraulic conductivity K. This paper will be restrictedto considerations of homogeneous soilsand, consequently(K0 will be assumed constant. A/so %he Kirchhoff %raneformation (Ames, 1965) has been used %0 eliminate

114 the square of the first derivatives which occur in the expanded form of eq.1. This transformation is:

= fPt Kr (P~) Pt d

I

(4)

,

,

in which Pt' is a d u m m y variable of integration. Rants and Gardner (1971) and Others refer to ~ as the matric flux potential. The writer's experience indicates that Without a transformation of this nature an adequate solution in the vicinity of the sharp wetting front cannot be achieved without either going to very small space and time increments or resorting to some "round-about" method such as that Used by Hanks and Bowers (1962) of first computing the diffusivity and then estimating the coefficients which depend upon the relative hydraulic conductivity from this compute d diffusivity. The second-order polynomial approximations used in the finite differences are inadequate in defining the derivatives of the pressure head (or Kr which changes more abruptly) across the wetting front. The change in ~ is more gradual across the wetting front and, consequently, the polynomial approximations duplicate the variation of the continuous function more adequately. If Kr is defined by eq.3, then eq.4 gives the following relationship between and Pt:

1 1 p--C= I 1 - (1+3~,)~] '+3~

(5)

Substituting eqs.2 through 5 into eq.1 and ex~,anding in cylindrical coordinates for axisymmetric flows gives:

~r2

+ ..... + {2+3X) [1- (1+3~)~] ~÷3~ .... + - " ~ ~}z2 ~z r c~r

(6) I+2X

[1- (1+3X)~ ] 1+3 a in which r is the dimensionless radial coordinate in the horizontal plane, z is the dimensionless vertical coordinate, and the dimensionless time parameter r is given by: Ko t r =

(7)

~,(Pb /Pg)rt(1- Sr ) The initial condition as:~ociated with eq.6 for all prob!e~s considered here° in will assume that. static equilibrium exists at the initiation of the numerical solution. The hydraulic head is constant under Static equilibrium and~ there-

115 fore, from eq.5 the initial condition is:

~(r,z, 0)= [1-1/(z-ho) 1÷3~ ] / (i+Sx)

(S)

in which h0 is the value of this dimensionless hydraulic head, i.e., the actual head divided by 9. = Pb ]Pg. Boundary conditions for the axisymmetric infiltration problem consisting of moisture applied over a circle of dimensionless radius ra and moving through a soft of dimensionless depth D = depth/~ which is underlain by an impervious layer are: (A) Top surface of moisture application (1) Flux rate, v3 (r) specified: ~}~(r,D, r)

v3(r) = Kr

,0
~z

(9)

K0

in which v3 (¢) is flux per unit area (with dimensions of velocity) and is positive when directed downward in the negative z direction. (2) Surface saturation S(D, ~') specified (Dirichlet condition), (S(D, ~) is denoted by S(1) in subsequent graphs):

[ 1 {S(D, ,) - S~) ~(r,D,r) =

-~- i - Sr

(,+3~)/~ ] / (1+3)~)

(lO)

(B) Top surface beyond radius of application: ~}~(r,D, r) = Kr

, r

> ra

(11)

~z (C) Impervious layer: = x~

(12)

(D) Axis of symmetry: =

o

(13)

~r

(E) Outer boundary beyond radius of influence, rf, (Divichlet condition}: ~(rf, z0 r) = ~{r,z, 0) FINITE DIFFERENCE SOLUTION METHODS

Three separate differencing schemes are used in adapting the Crank-

(14)

I16

Nicolson method for solution of the previously defined initial boundary value problem and three variations of the alternating direction implicit (ADI) method have been used. Each of these six schemes for solution is described below.

Crank-Nicolson methods Finite difference equations The Crafik-Nicolson method weights difference approximations of file derivatives with respect to the space coordinates at the current and advanced time steps equally as the derivative with respegt to time is approximated by a second order central difference centered midway between these two time steps. The discretization errors of the Crank-Nicolson method are second order for both space and time. The three different schemes which utilize the Crank-Nicolson method will be denoted by the finite difference Operators Fij, Gij, and Hi]. The Fij operator is obtained by first multiplying eq.6 through by:

1+2k [1-(1+3X)}] 1+ 3~ {which will be denoted by a~ ) and then differencing the space derivatives by second order central differences. For a grid network with Ar = Az = As, ~his operator is:

Fi]

=

;:k+l

.-k+l

alk+l(b~i+lj+b2~i_li)+ - a2

) ~i.i+ i -

(ak+l

+

+ c) ~i]

ak+l

~k+l ) ~ij-1 i+ i j

~- lj )

(15) 2,3

1=

.Nz

1

in which the subscripts i and ] denote the space grid points such that i = 1 + r/As and ] = 1 + (D-z)/As, the superscript k denotes the time step such that i+ = 1 + r/Ar~

2As 2 c = ~

(16)

hr 0.5 ~ +~ i -1

As

b~ = ~ + ~ = 2r As b~ = !

0.5 -

2r

(17)

1-

(18) ~ -I

117

1+2~

(~9)

a~ = [1-(l+3X)t~]] ~+3~ and: 2+2k

(2 +3h)As a2 = - - - - - - - - - [1- (1 +3X)~ii] ~+a~" 2

(20)

(In eqs.19 and 20, the superscripts k of ~ correspond to the superscript of a in eq.15). T h e Gii operator is obtained by differencing eq.6 without prior multiplication by al, and the coefficient of the time derivative is evaluated by assuming that ~ + ~ / 2 which is midway between the two time steps, equals the average of ~h and ~ + t . This operator is: =

j+b2

i-l/+(l+aak+l)~i/-i +

)~/]+l k

+ (1 +a~) ~kij-l + (1 a~) ~i.i+l k

~k = o

- 4 ij

-

(21)

in which the b's are the same as above,

a3 = A s ( l + l . 5 h ) [ 1 - ( l + 3 h ) ~ i j ]

l+3~k

(22)

and: 4+ 1/2 a4 =

2&s2

(23)

&r [1- (.5+1.5~) --ii{~+1 + ~kii) ] cx , 2~.)1(x + 8~,~ The Hij operator is similar to Gij with the exception that the right side of ecl.6 is replaced by the invalid approximation"

(ask*~ ~*'.. _ ~ j ~ ~)/2as~

c a5 = - - = al

(24)

2A s'~

(25) i +2 h.

A ~ [ 1 - (1 +3~) ~ij~ 1 +3

118 which replaces the term associated with a4 in eq.21. If the coefficient of the time derivative were not a function of g, then the operators Gij and Hij would be identical. The finite dil.~erence operators for the top surface (if eq.9 is used as the condition over the circle of application) and the impervious layers are obtained by approximating the derivatives in the boundary conditions with second order central differences centered on the boundary, and subsequently eliminating the value of ~ at the nonexistent grid p o i n ~ u t s i d e the boundary by combining with the appropriate finite difference operator for interior grid points. These difference equations are not given herein for the sake of brevity. The boundary condition at the axis of symmetry where r = 0 has been handled by setting ~~d = ~~,J in all the above finite difference equations. This approach was used instead of developing an operator by combining the central difference approximation of eq.13 with the regular operator, because b~ and b~ are undefined along the line singularity at r = 0.

Solving finite difference equations Writing the finite difference operators for all grid poin~ (interior and boundary) produces a system of nonlinear algebraic equations for the unknowns ~.+i (which w i 11be denoted by the vector ">k+l g below). A solution of this system advances the solution of the infiltration problem through one time step A r. The general Newton-Raphson iterative method has been utilized in solving this sytem. This method provides a better approximation to the un° known vector ~k+~ after each iteration by means of the formula (see Saaty and Bran, 1964, for example): 1

m+

1

m

in which the sup.erscript m outside the parentheses denotes the iteration number, the vector ~ consists of the elements composed of the finite difference operator Fij, Gij, or Hij depending upon the finite difference approximation used, each of which wit1 equal zero when the solution has been obtained, and D is the Jacobian which consists ~f a banded matr~. The actual implementation of eq.26 in a computer progr .am will replace the quantity behind the minus sign by the solution vector ~ to the linear system: (D)m

~= (~)m

(27)

The solution vector ~ was initially obtained by using the ~gcrithm described by Thumau {1963}. (The actual listing of the ALGOL program was obtained from Burrough Corporation.} This algorithm is more efficient in computation and storage requiremenSs than standard Hnear algebra algorithms, because it opeloates only on the band portion of the matrix Do Subsequently0

119 upon noting that the diagonal elements of D at the zero of the functions are considerably larger in magnitude than the off-diagonal nonzero elements for values of ~ which cause the elements of F to be zero, partkularly during the first time steps of the solution, it became apparent that considerable reduction in computer storage requirements could be achieved and, at the same time, possibly decrease the amount of computer execution time required for a solution, by utilizing an inner iterafive scheme. This scheme in essence combines the line-successive-relaxation iterafive method with the Newton-Raphson method; that is, an iteration is created within an iteration. The innermost iteration solves for the ~/*+~ 's along consecutive horizontal lines from the system of equations resulting under the assumption that the ~k+~ 's on the previous and the next line are known. In other words, the two outer bands of D are assumed zero and, consequently, a tridiagonal system of equations, equal in number to the grid points along this line, is solved during each inner iteration. This inner iteration is continued until the sum c.f absolute changes in g's along the line beeomes less in magnitude than a specified error (i.e., approximately 10 -7 ), before repeating the same inner ite~ative process at the next line. A pas.s through all lines constitutes one outer iteration, and provides values of ~ throughout the flow region which are close to those that would be obtained from one iteration by the Newton-Raphson method, (eq. 28). When the accumulated absolute changes in ~ from all inner iterations become less than a second error parameter, the outer iteration is terminated. The above solution processes will be referred to as the Newton-line-relaxation method. Over-relaxation, or over-adjustment, of individual lines has n¢~t been studied, but would be expected to decrease the number of iterations required. ADI methods

ADI methods (Douglas, 1961, 1962) complete the advance of the solution through eacb. time increment by a two step operation. The first step operates along consecutive lines in one direction and the second along lines in the other direction. Three variations of the ADI method are studied herein. These consist of: (1} evaluating all coefficients in the finite difference equations from the dependent variable at the current time step ]z; this me;hod is classified by Blair and Weinang {1969) as a mixed finite difference equation and will be referred to herein as the ADI method; {2) use of a predictor to estimate ~'s at advanced time steps from which coefficients are evaluated; this method was used by ~ubin (1968) and is referred to herein as the ADIP method without iterative corrections; (3) use of a predictor which is iteratively corrected; this method is referred to herein as the ADIPC method with iterative cor~'ections. ~he finite difference eqaations for interior grid points for all three ADI methods can be written as:

120

~k+ ~t~ k+~/~ ~k+~P ~k+~l~ -b;o~_~j +(2+a~ )~ij - b~ - i + q =

(l+a~)

k

~j-~

. k+~/4 +~a~

-

2) k +(1 ~ij

-

a~) k ~ij+~

(28)

for the first portion of the time step and:

-(1 +a3 )~ij-I +(2+a4 . k+l/2

= b2 ~ i - l j

. k+3/4

+ ~a4

- 2)

)

.. - ( 1 ~3

a3

)-~j-1

k+ll2 + bl ~k+l/2

~i1

~i+~j

(29)

for the second portion of the time step. (The difference equations for the boundary grid points can be obtained from eqs. 28 and 29 by the procedure described for the Crank-Nicolson operators.) The superscripts of the a's as given in eqs.28 and 29 apply only for the ADIP and ADIPC methods. In the ADI method, a k+~/4 ~ set equal to a k and ak+ 3/4 is set equal to a k+~12. (The actual fini~ difference equations used in method i were obtained by an alternate procedure described by Douglas (1962) which results in slightly simpler equations than eqs.28 and 29.) The predictor used in the ADIP and ADIPC methods extrapolates from values in the current time plane by the formula: ~ k Ar (~h+ ~/2 )0 = ~k + _ _ (30) ~}r 2 A sii~:ilar fermula~ but with superscripts advanced by 1/2 is used to predict at the second half time step. The a4 's are computed using the average of the predicted (~'s) ° and those at the current time plane or half time plane respectively depending upon whether eq.28 or 29 is used. In eq.30 (~}~/~r) is evaluated, after solving for it from eq.6, by replacing the derivatives on the other side of the equal sign by second order central differences. In the ADIPC method, the predicted values for the second iteration (i.e., those with a 1 superscript outside the parentheses) are obtained from seMng the resulting system of equations given by evaluating the a's with the preo dictor with a 0 superscript° This process is repeated iteratively at both the first and second portions of each time step until values of (~)m+~ and (~)m are equal to within a specified error. CHARACTERISTICS

OF FINITE

DIFFERNNCE

EQUATIONS

In attempting to obtain a solution using the F operaters~ eqol 5, etco, the ~iter discovered he could not program an adequate initi~ guess so that the NewtonoRaphson method would converge by means of a few sta~ments~ deo spite the fact that the nature of the problem permits re~onab!e estLraates for

121

at advanced time steps to be generated easily. Even after achieving success in getting a proper solution for a few time steps by meticulously adjusting the initial guesses, the Newton-R~phson iteration failed in solving the system at later time steps. Determining why the New~on-Raphson ~tera$ion failed reveals the strong effect of the nonlinear coefficients on the functional beharbor of the d~fference equations. To simplify the examina~on o~ this beharbor, finite d~fference operators for the one-~mensional vertical moisture movement problem will be used. These latter operators can be obtained either by differencing the one-dh~nens~onal equation of flow or setting ~i* q and ~i-g equal ~o ~q in eq.15. These operators will be denoted by Fi, j = 1~2o.. Nz.

|

=

iNDEPENDENT VARIABLE 3=~

F~g.1. Functional relationship of F= for several application rates, Q.

The variation o~" F1 with ~1 for six different application rates Q = v~/K0 and for ~, = 1.5 is shown in Fig.1. All of the curves on this figure have been obtained using a value of ~+~ which satisfies the static equilibrium condition to 16 significant digit. The zero for the function FI occurs when the curve intersects the horizont~A dashed line approximately midway through the graph° In plotting the carves on this figure each curve ends at the right of the graph at a point beyond which the function F: becomes undefined because the quantity [ 1 - (1 + 3 ~ ) ~ +~ ] which is raised to an exponent in the equa° tion for FI becomes negative° For the smaller appIicati0n rates it appears ~ a t ~he point where F, becomes undefined and where its zero exis~ are co° incident. This is not the case~ however~ as shown on Figs.2 and 3 in which the

plotting in the immediate ~cinity of this point has been greatly magnified, for ~he curves of Q = 0.001 and Q = 0.01~ respectively. Rather the following are true.

122

_o ~,

O

"

-

" /--=

I

!

0

P,

g

F

i

t" .18181738

+2"qO

'~10

7

+6xlO 7 .181818iB

INDEPENDENT

VARIABLE

91

O~ I" ,1BI61818

4xlO-~

8~tO-.~

INDEPENDENT

12~I0-~ VARIABLE

IGxlO-,~

.I 8 1 8 1 8 i 8

~

Fig.2. F u n c t i o n a l relationship of F~ to ~ for dimensionless application rate, Q = O.0Ol. Fig.3. F u n c t i o n a l relationship of F~ to ~ with dimensionless application rate, Q = 0.01.

(1) the zero of F~ occurs at a value of ~ which is very close to the w~lue whe~'e F~ becomes undefined, but which is always less. (2) The function F1 is positive for ~1 less than its value at which FI = 0 and negative for ~ greater than this value. (3) Between the zero of F~ and its undefined region, F~ reaches a rainimum. For ~ less than the value causing the minimum ;~F1/~}t~ is negative and for ~ grea~er than this value ~}F~/~}~ is positive. (4) As the application rate decreases the phenomena described in (1) through (3) occur essentially at the same point, and the separate phenomena can be detected only by examining the values to a relatively large number of significant digits. This erratic variation of F~ is due to the strong+Ieffect that small changes in k+l ~ have on the magnitude of the coefficients a{ and a k+~ 2 in F~. This erratic functional behavior is not restricted to the finite difference operator F~ for the surface ~ i d point, but is characteristic of the F's at all grid points. The relationship of F~ at the next grid point with ~2 varying over a relatively small range is given in Fig.4. The separate curves on Fig.4 were obtained by setting ~+~ equal to a value which is less than its static equilibrium value by the amount indicated on each curve. It is not difficult to appreciate wh~+the NewZon-Raphson method by itself is !naaequate. Should the guess for ~7 be larger than that which gives Fj a mmn~num, subsequent iterations would cause Fj to very rapidly~ if not immediately~ become undefined. On the other hand, if ~his g~ess is veiny much sm~Jler than the zero of Fj, the first iteration would project into regions where F/is undefined. Furthermore, if any iteration should cause Fj to take on its minimum value or close to it, the next vector of ~ j " from the New°

123 0

(~

~ -

tu 0

o,o

~

below the surface

/

I [

l r.~62

.~7o

INDEPEHDEHT

.~so .~z

V~RIABLE

~z

Fig.4. Functional relationship of F~ at the first grid point below the surface with the value of t~~. The value of ~, at the surface grid point has been decreased by the indicated increment to define the separate c u r v e s .

ton-Raphson iteration would be very unrealistically large or small. Clearly, to solve the problem using the F operators requires some approach in addition to the Newton-Raphson method, at least until the values of ~ are smaller than those associated with no moisture movement, but such values always exist just in head of the wetting front. The approach which has been used to successfully solve problems using the F finite difference operators takes advantage cf the fact that ~}Fj/b~j in the vicinity of the root is larger than either ~Fj./~}~]-, or aFj/~j+~ and, therefore, if the correct root for each separate ~ can be deterrained, then an iteration which determines in sequence the root of each F defined by using the most current value of ~j_, and ~j+~ will converge. The root to each Fj is obtained by first squaring and then utilizing a Fibonacci search {Wilde and Beighter, 1967} to obtain the minimum of the squared function, A small interval for the search is relatively easy to determine from the characteristics of Fj described earlier. After using the Fibonacci search-iterative scheme for a few iterations,the solution process is turned over to the Newton-Raphson itera~on. A much better approach to a solution is avaflable~,however, by using the G finitedifference operators° Fig°5 shows that the (~j'sdo not: behave erraticallyin the vicinity of the root as do the F/s. Cc.nsequent~.y, no particular difficultiesoccur in using the Ne~on°Raphson, or Newton°line-relaxation method in so!ring the system given by %he G finite difference operators. Solution of the equations from the F operators are considerably more difo

124

\.\.

\, 5:

J

ra

O

INDEPEIVD?.N T VA RIABL E

Fig.5. Functional relationship of G~ a~ ~he first grid point below the surface with ~ , ~he 'vah~e of ~ at ~he surface has bee~ decreased from ~t.s equilibrium value by the indicated araoun~ to define the separate curves.

fieu!t to achieve than those from Lhe (2 operators. Contrasting F and G gives the message to fl~ose contemplating development of a finite difference solution to an unsaturated flow problem tbat the functional behavior of the finite difference equations should be exa~fined, and if this behavior is undesirable heom the view point of obtaining a soiution, ~ha~ ~J~erna~ive developments should be studied. Wha~ m~ght appear to be inconsequential changes in the development of the equagion can cause major changes in the ease with which the resulting system can be solved,

"FABLE

Variat.ion o~' ( ~ +~t2 )m and {g'~+l)m at ~hree ffeid points* Iteration, m

0

1

2

8

4

5

6

7

8

.02~0

.0942

.1253

.1408

.1491

.1541

.1572

.1593

.1607

~a,~

.1538

o1752

.1793

.1805

1810

.18_a

oi8i3

o1814

.1814

~v,a

.1627

.1778

.1804

.t8!1

o1814

oi815

.1815

k+,/a }m

':~The p~°ob~enaproducing ~hese variations is given ~n Fig°7°

125 The predicted values of ~ used in the AD][P method are n o t very close to those resulting f r o m a solution of the G operators, particularly in the vicinity of the wetting front where the functions vary rapidly. Table I illustrates how the predicted values (~k+~p)m and (~k+~)m vary over a few iterations for the first time step in obtaining a solution by the ADIPC m e t h o d . COMPARISON OF SOLUTIONS

Solutions based on operators F and G agree closely. The detectable differences are restricted to the region surrounding the wetting front. For example, solutions with ;~ = 1.5~ ho = - 3.0, Sr = 0.15, 77 = 0.4, D = 3.0, v3/Ko = 0 . 6 , As = 0.1, Ar = 0.025 and an error parameter of 5 . 1 0 -8 , show a difference of 1.71% at the surface grid point when r = 0.05, with the solution based on G giving the larger saturation. After advancing through 18 time steps (r = 0.45), the saturations at the surface agree to within 0.07%, but at the wetting front, which at this time has advance:l to a depth of approximately 0.t~ units~ ~he m a x i m u m difference is 4.51%~with the solution based on G giving the highest saturations. Consequently, for practical purposes these two operators produce identical results. Clearly, since the G functions are better behaved tl-.an the F's, their use is preferable. By comparing solutions to the same problem obtained using different time and space increments, an indication is given whether these increments need to be decreased in size or not. The problem specified by h = 0.9, Sr = 0°22° = 0.4, b0 = - 3.0, D = 2.0, and ,°(D~r) = 0.85 was solved thrice using first AT = 0.05 and h s = 0.2, second h r "= 0.025 and As = 0.2, ~nd third A • = 0.025 and As = 0.05. Very small differences exist between these three solutiol~s. These differences are confined to the vicinity of the wetting front and soon disappear. For instance, the largest difference in saturations between the first and second solutions occurs at the time step r -- 0.15 at the fir~ grid point beneath the surface and near the centerline where ~he difference is 0.41%. When r = 0.6 the difference at this point diminishes to 0.04%, b u t the m a x i m u m difference occurs at a depth of 0.7 units neur the wetting front and equals 1.5%. A comparison of the first two solutions with ~he third shows similar very small differences, with first and third in closer agreement than the second and third. These small differences lead to the conclusion that reducing the grid spacing or time increment does n o t significantly improve the solution. In contrast to this excellent agreement, solugons based on the H finite difference operators show no quantitative agreement with those from the other methods. This lack of agreement is illustrated on Fig.6, which shows both the volume of water infiltrated and t.be changes in saturation at a p o i n t ~ given by ~olutions of the identical problem by the G and H operators and the ADIP and ADIPC methods. The source of error in the solution from the H operator is an inappropl~ate differencing of the right side of eq.6 as given

126

0

25

50

75 1.0 T~me Parameter, "r

1.25

! 5

Fig.6. I n t a k e volume of water and saturation at the c e n t e r h n e at a d e p t h of 0.1 units as given by solutions based on t h e Gi] and Hij finite difference e q u a t i o n s and obtained bg the ADIP and ADIPC methods.

by eq.24. Actually, eq.24 is a valid approximation of (Ar/2As ~ }~(es ~)/Or, bu~ no~ of (Ar/2As =) as O~/0r. However: ~ = ~ 2

3r

5 ........... + ~

~)r

~)z

The vast differences in the solution based on the H operators are caused by the term ~{aas/a~'}° This term takes on relatively large values, particularly for moisture conditions near static equilibrium (see eq°25). These values are opposite in sign to as {~}~/ar) and cause the terms in the H operators from ~he fight side o5 eq.6 to exert too small an influence in changing values of ~ between consecutive time steps. Since the H operators would be identical to the G operators if ~he coefficient,as, were constant, the vast difference in the solution from the H operators dramatizes the influence of one of the nonlinear terms in the unsaturated equation of flowo Any method o5 ]ineafi~ zation~ it appea~.~,m a y result in a solution which qualitatively appears valid. bu~ which, in fact, is quantitatively grossly in error. Fur~her evidence o~ ~h~s ~ac~ comes from noting on Fig.6~ ~ha~ the solution by Zhe AD][P method does no~ agree closely wi~h ~hose obtained by ~he C~°ankoN~co~son me~,hod (G-operaZor) and ~he ADIPC rae~hod. ©bviousl:/, even eraS.unSUng ~he nonolinear coe_~c~en~s from predicted values o~ ~ is no$ adeo q~,~a~,eo Figo7 provides addi~iona~ evidence of ~his ina~equacyo A~ ~ees~ ~,o~

127

o

~

:

3

5

6

[ , m t P.~ramele¢ "~

Fig.7. Saturations at two points obtained from solutions using: (1) the Crank-N~colson Method and Gi] finite difference equations; (2) the ADI Method; (8) the ADIP Method; and (4) the ADIPC Method.

problems with abrl:pt we~ting fronts, a fully implicit method, requiring the solution of a system of nonlinear difference equations, such as the CrankNicolson me~hod, or the ADIP(A (~¢ith a number of iterative corrections) is needed for valid solutions. Fig.7 shows saturation changes at two poin~ as given from solutions by the G operators and ~he three ADI methods described earlier, to identical problems specifying a constant infiltration rate of v 3 / K o = 0.6. Only the solutions from ADIPC and the Crank-Nicolson methods are in agreement. The reason why the ADI method produces saturation values which are too small and tbc ADIP method gives saturations which are too large, particularly near ~he wetting fron%js that the first predicted values of are too smail (over adjusted) as shown ir:~~able I. CONCLUSIONS REGARDING SOLUTION METHODS

While the exact results of the different solution methods described herein apply only to the particular formulagon used, similar results would be expected in sobAng other fGrmula~ions of unsaturated flow in porous media. These c h a r a c ~ s t i c s are a consequence of the s~rongly nonlinear nature of the physical problem, in ~he absence of genera! theoi-y of methods for solving nonlinear i n i t ~ - b o u n d a ~ °value problems~ it appears %hat those methods of finite differences which have been developed for linear problems, can be adapted to provide adequate so!ugons to nonlinear equatk,ns as evidenced by

!28

the close agreement between solutions from the Crank-Nicolson A D I P C methods for different time and space increments. The solutions will only be adequate, however, if the method takes into account that the coefficients which cause the nonlinearity vary rapidly and, consequen~dy, must be evaluated appropriately at advanced time steps. Appropriate evaluation of these coefficients cannot be accomplished by using only data from current time steps. Either the evaluation of these coefficients must take place by solving the nonlinear systern of finite difference equations in Which ~he coefficients are unknown, as well as the dependent variable (as was done by use of the Newton-Raphson method), or by iteratively correcting the coefficient (as was done in the ADIPC method). The latter approach will not be convergent for all iterative correction methods, however. Solutions produced by adaptations of explicit methods used in solving Hnear parabolic partial differential equations should, therefore, be viewed with suspicion regardless of how small a time increment is used. The Crank-Nicolson method requires slightly more computer execution time than the ADIPC for comparable problems solved. This comparison is based on solutions that required the sum of changes between consecutive Newton-Raphson iterations be less than 3 • 10 "s before terminating the solution and terminating the corrective iteration in the ADIPC method when the largest change between (~)m and (~)m-~ was less than 0.0001 (to decrease this error in the ADIPC method causes considerable increases in execution time)o No particular numerical difficulties occurred in obtaining the solutions to the Gij finitedifference operators, but in using the A D I P C method values of ~ at grid point at radial distances beyond the wetting front, particularly near the surface, tended to become slightly larger than the initialstatic equio librium values. Logic was therefore added into the computer solution to set any value of ~ which was computed larger than the static equilibrium value equal to the static equilibrium value. Since no such constraints were programmed into solving the Crank-Nicolson method using the Gij operator, the ~ i t e r favors this method slightly, but also because there is le~s question regard~ng the order of approximation from use of the Crank-Nicolson method. ACKNOWLEDGEMENT

The majority c.f the support for this work was by the Agricultural lte~earch Service, Northwest Watershed Research Center, Boise, Idaho, through a cooperative agreement with the Utah Water Research Laboratory~ Utah State UnbJersityo APPENDIX

X- NOTATION

a s , a ~ ~'~3, a4, as = v a r i a b l e c o e f f i c i e n t s b~ = variable coefficient c~ = constant = 2As 2/4

D

= dimensionless depth~ or depth/~

179

g ho K K~ Kr h P Pb Pt r

ra rf S

Se Sr s(1) As l) Z

k

P T

= = = = =

designation for a nonlinear system of finite difference equations designation for a nonlinear system of finite difference equations acceleration due to gravity designation for a nonlinear system of finite difference equations dimensionless hydraulic head at which the soils exist initially under static equilibrium, i.e. this head/~ = hydraulic conductivity = saturated hydraulic conductivity = relative hydraulic conductivity = hydraulic head or energy per unit weight of fluid being the sum of the elevation and pressure heads = capillary pressure = bubbling pressure and parameter in Brooks-Corey equation = dimensionless pressure head, or pressure head/¢ ffi dimensionless radial (horizontal) coordinate, or this coordinete/~ = dimensionless radius of circle of application, or this radius/~ = dimensionless radius of outer boundary of region of interest, or this radius/~ = saturation, or volume of water divided by volume of voids = effective saturation = residual saturation, and parameter in Brooks-Corey equatior, s = saturation on surface of application = dimensionless space increment, or this increment/~ = application flux = dimensionless axial (vertical) coordinate, or this coordinate/~ = pore size distribution exponent, and parameter in Brooks-Corey equations = porosity, or volume of voids divided by total volume = dependent variable obtained from Kirchhoff transformation = dei~zity of fluid = dimensionless time parameter = K o t P g / t ~ , P b ~(1-Sr)]

#

T

Ar

= time increment

REFERENCES Ames, W.F., 1965. Nonlinear Partial Differential Equations in Engineering. Academic Press, New York, N.Y., 511 pp. Blair, P.M. and Weinang, C.F., 1969. Solution of two-phase flow problems using implicit difference equations. Soc. Petrol. Eng. AIME, Trans., 246: 417--424. Brooks R.H. and Corey, A.T., 1966. Properties of porous media affecting fluid flow. J. Irrig. Drain. Div., ASCE, 92 (IR2): 61--87. Brooks, R.H., Ng, B., Corey, G.L. and Corey, A.T., 1971. Drainage of Soil Profiles. J. Irrig. Drain. Div., ASCE, 97(IR3): 455--467. Brutsaert, W.F., 1971. A Functional Iteration Technique for Solving the Richards Equations Applied to Two-Dimensional Infiltration Problems. Water Resour. Res., 7(6): 1583~-1596. Burdine, H.T., 1953. Relative permeability calculations from pork-size distribution data. Petrol. Trans. Am. Inst. Mining, Metali. Eng., 198: 71--77. Corey, G.L. and Corey, A.T., 1967. Similitude of drainage of soils. J. Irrig. Drain. Div., ASCE, 93(IR3): 3--23. Douglas Jr., J., 1961. A Survey of Numerical Methods for Paraboiic Differential Equations. Advan. Computers, 2: 1--15. Acad. Press, New York~ N.Y. Freeze, R.Ao~ 1971. ~ree-dimensiona~transient, saturated--unsaturated flow in a groundwater basin. Water Resour. Re~., 7(2): 347--366.

130 Hanks, R.J. and Bowers, S.A°, ~962. Numerical solution of the moisture flow equation for infiltration into layered soils. Soil Sci. Soc. Am.. Proc., 26(6): 5 3 0 - 5 3 4 . Jeppson, R.W., 1972. Relationships of Infiltration Characteristics to Parameters Describing the Hydraulic Properties of Soils. Utah Water Research Laboratory, Utah S~.ate University, Logan, Utah, PRWG-59c-7, 57 pp. Jeppson, R.W., 1974. Axisymmetrie infiltration in soils, II. Summary of infiltration characteristics related to problem specifications. J. Hydrol., in press. Raats, P.A.C. and Gardner, W.R., 1971. Comparison of empirical relationship between pressure head and hydraulic conductivity and some observations on radially symmetric flow. Water Resour. Res., 7(4): 921--928. Remson, I.R., Hornberger, G.M. and Molz, F.J., 1971. Numerical Methods in Subsurface Hydrology. Wiley-Interscience, New York, N.Y., 389 pp. Rubin, J., 1968. Theoretical anaiysis of two-dimensional, transient flow of water in unsaturated and partly unsaturated soils. Soil Sci. Soc. Am.,Proe,32(5): 607--615. Saaty, T.L. and Bran, J., 1964. Non-Linear Mathematics. McGraw-Hill, New York, N.Y., 381 pp. Smith, R.E. and Woolhiser, D.A., 1971. Mathematical Simulation of Infiltrating Watersheds. Hydrol. Pap., 47, Colo, State Univ., Fort Collins, Colo. (See also: Overland flow on an infiltration surface. Water Resour. Res., 7(4): 899--913. Thurnau, D.H., 1963. Algorithm 195 BANDSOLVE. Association for Computing Machinery, 6(8). Verma, R.D. and Bru~aert, W.F., 1970. Unconfined aquifer seepage by capillary flow theory. J. Hydraul. Div., ASCE, 96(HY6): 1331--1334. Wilde, D.J. and Beigbter, C.S., 1967. Foundations of Optimization. Prentice-Hall, London, 480 pp.