|
<br />can be used for this purpose, The computation of
<br />81'1 by means ofEq, 27 is thus straightforward.
<br />
<br />For grid point i = 1 with j<n
<br />
<br />If the .p~+ I value was, already computed by
<br />using Eq, 20 at i = 2, the .pJtl value on the upper
<br />boundary can be determined by means of Eq, 22, in
<br />which the rainfall inten~ty at ~~me level j + Ii rl+ I, is
<br />also given, J:jecause Kl3{'h, KJI l, and .pJ+ are all
<br />related to ell 1 through he known phYsic~ pro~erty
<br />relationships of soil, Eq, 22 can be solved for .p\ I by
<br />trial and error. The iteration of the computation can
<br />'+1 .
<br />pr~crd by ,~rlst assuming ell ",ej. Then compute
<br />K)/2 and KJI from the known K( 8) relationship as
<br />follows: Substituting the expressions
<br />
<br />.j+l = (oj+1 + .j+l) /2
<br />3/2 1 2
<br />
<br />.(28)
<br />
<br />Kj+1 _ K(Oj+l)
<br />3/2 3/2
<br />
<br />,(29)
<br />
<br />Kj+1 = K(Oj+l)
<br />1 1
<br />
<br />, (30)
<br />
<br />into Eq, 22 lields .pit I for known .p~+1 The
<br />computed .pt value in turn >,yill be substituted into
<br />the .p( e) relationship fo~+t+ 1, If the difference
<br />between the computed ej from, the .p(e) relation-
<br />ship and the initially assumed eJ+I is found to be
<br />within the tolerable accuracy, Bq, 22 is solved,
<br />Otherwise assume the e1' I value just computed in the
<br />last step and follow the foregoing computation
<br />procedure until the accuracy is met. The iteration can
<br />be accomplished in a systematic way,
<br />
<br />For grid points in the saturated zone
<br />
<br />The Laplace equation, Eq. 14, in the fmite-
<br />difference form can be formulated to compute the
<br />.pt I values for i = 2, 3,4, ..', up to the last saturated
<br />nodal point as follows:
<br />
<br />,,)+1 .. 2,)+1 + ..,.1+1 =- 0 (31)
<br />'i-l ' i '1+1 ' . , . . . .
<br />
<br />Since .p(z)-proflle is linear in the saturated zone, Eq,
<br />15 can be used instead of Eq, 3 L Recalling the upper
<br />boundary condition after ponding, Eq, 23, one can
<br />formulate Eq. 15 in the finite-difference form as
<br />
<br />(fj+l - K )
<br />.t1. -I~s s (i - 1) Az + .{+1 ,(32)
<br />
<br />
<br />where the current infIltration rate, fj+l, is computed
<br />by using Eq, 10 in the finite-difference form:
<br />
<br />fj+l =- -; (1/2)(oji- +1 _ oj "j+l
<br />i + "'i+l
<br />1=1
<br />
<br />oj ) ^Z + K (33)
<br />1+1 ll.t 0
<br />
<br />Actually, f(t) computed by using Eq. 33 corresponds
<br />more closely to fJ+1/2 than to fJ+1 Accordingly, the
<br />.p value computed from Eq, 32 should be at time level
<br />G + 1(2) ralher than G + I) for given pondlng depth,
<br />I/J i + 1 ~ at the corresponding time level. There certain-
<br />ly will be a time lag, I:1t/2, for .p and f values so
<br />computed at each time level, but the computed values
<br />will in uo way be affected by the time lag,
<br />
<br />Before one computes the .p(z)-proflle in the
<br />saturated zone, the e (z)-profile in the unsaturated
<br />zone must be determined. As mentioned previously,
<br />the .p value in the saturated zone so CornjlUted could
<br />be in a significant error if the current fJ 1 computa~
<br />tion for the unsaturated zone by using Eq, 33 is not
<br />accurate enough. In other words, the inaccuracy in
<br />the computation of the current fi+ I value would
<br />result in the erroneous computatjon of the location
<br />of the current saturation frort, q+ 1 by using Eq. 17,
<br />and hence the current ",r values in the saturated
<br />zone by using Eq, 32, Thf apparent interaction in the
<br />computation of the .pr and f J+ I appearing in both
<br />Eqs, 32 and 33 may result in computational oscilla-
<br />tion which, howevcr, damps out quickly with the
<br />advancing saturation front.
<br />
<br />Use of Eq, 16 in the evaluation of the infiltra-
<br />tion rate, f(t), has a problem at the time of ponding,
<br />tp' when both h(1)) and Lf(tp) become zero, as
<br />pointed out previously, Therefore, Eq. 16 cannot be
<br />used in the computation of f(t), Instead, Eq, 10 or,
<br />more specifically a finite-difference form thereof, Eq,
<br />33 was used throughout the study, Use of Eq, 33
<br />does not require the exact location of the saturation
<br />front, Lf(t), for the total rate of change of the
<br />moisture content in the saturated zone is always
<br />assumed to be zero. Furthermore, because the present
<br />method requires that the moisture content in the
<br />unsaturated zone at the current time level G + I) be
<br />computed before proceeding to compute .p in the
<br />saturated zone, there should not be any technical
<br />difficulty in determining Lf (t) from Eq, 17 in which
<br />f(t) is approximated by Eq, 33, Of course, the
<br />accuracy of Lf(t) so determined depends on how well
<br />Eq, 33 can approximate f(t), which is in turn
<br />dependent upon the accuracy of 8 for all nodal points
<br />in the unsaturated zone.
<br />
<br />It is understood that in any numerical scheme.
<br />explicit or implicit, the roundoff error is generally
<br />caused by a combination of numerous factors includ-
<br />ing the step size in space (llz) and time (lit) and
<br />additional assumptions or conditions imposed in the
<br />computation, If the e value varies rapidly with the
<br />soil depth (te" the cases usually associated with very
<br />large 6.z, r, and h or very small 80, as mentioned
<br />before), determination of 8 by using Eqs, 20 and 22
<br />may not be sufficiently accurate, Since Eq, 33 was
<br />also used in the approximation of f(t) before pond.
<br />
<br />21
<br />
|