|
<br />halfway between grid points i and (i + I), respec-
<br />tively. The hypothetical lower boundary condition at
<br />z = ""l Eq. 5, win not be utilized in the numerical
<br />computation, therefore it is not included in the
<br />numerical infiltration model that consists of Eqs, 20
<br />through 24, At a glance, one might suspect that the
<br />proposed numerical model is a combination of
<br />explicit and implicit schemes in a sense that causes
<br />computational oscillation, but actuaDy it is not that
<br />case, It is explicit all the way, namely unknowns are
<br />being solved for, one at a time. The numerical
<br />computation will begin with the finite-difference
<br />initial condition (Eq, 21), then proceed to determine
<br />unknowns for i = 2, 3, 4, ... at next time level by
<br />using Eq, 20, and finally end up for each time ievel
<br />with the solution of Y,or 0 for i = I by means of the
<br />finite-difference boundary condition, Eq, 22 (before
<br />ponding) or 23 (after ponding), with Eq, 24 to be
<br />used as a criterion for determining the time of
<br />ponding. The fact that Rubin and Steinhardt (1963)
<br />also formulated a boundary condition similar to Eq,
<br />22 in their linearized implicit.scheme model for rain
<br />infIltration does not necessarily mean that Eq, 22
<br />cannot be used in the explicit,scheme model. Because
<br />Eq, 20 is valid for all i's except for i = I at time level
<br />(j + I), the boundary condition such as Eq, 22 or 23
<br />must be formulated at time level (j + I) in order to
<br />solve for the remaining unknown at i = 1.
<br />
<br />Finite-difference approximation of Eq, 3, as
<br />expressed by Eq, 22, has two approximate expres-
<br />sions of K, In Eq, 3, the first term in which K is
<br />multiplied by 3y,/3z should also be designated on the
<br />same soil surface (i = I) as the second term, but
<br />could not be done so because an approximation of
<br />3 # az could only be accomplished between grid
<br />points i = I and i = 2, This in effect forces K in the
<br />first term to be specified halfway between the two
<br />grid points, Equation 22 so formulated is somewhat
<br />different from Rubin and Steinhardt's (1963)
<br />formulation in which both K's were evaluated at i =
<br />I, but without giving enough account of how
<br />31/13z was approximated. Other investigators such as
<br />Smith and Woolhiser (1971) did more or less the
<br />same as Rubin and Steinhardt.
<br />
<br />Ponding occurs at time level j = n, the value of
<br />which can be determined by use of Eq. 24, If l'I. is
<br />invariable, ponding mayor may not occur exactly at
<br />time level j = D. For. convenience, in the present
<br />computation, only 6.z is specified and varying lit's are
<br />computed from the following stability criterion
<br />(Richtmeyer, 1957; Gupta and Staple, 1964):
<br />
<br />D At
<br />max
<br />(Az)2
<br />
<br />5 A
<br />
<br />, , . . , , , . , ' .(25)
<br />
<br />where Dm ax is the maximum effective diffusivity for
<br />a given moisture profIle and ^ is a constant. With a
<br />
<br />linearized form of Eq. I without a 3K/3z term,
<br />Richtmeyer (1957) obtained ^ = 0.5, Since K(O)
<br />varies with the soil depth, z, in the unsaturated zone,
<br />Eq, I is not linear in the present form, Despite
<br />Richtmeyer's K = 504, it is questionable to use Eq,
<br />25 as a criterion for stability of the explicit solution
<br />in the case that K( 0) varies considerably with the
<br />specifiedtJ.z. In other words, one can expect appreci-
<br />able oscillation in the solution if the change in the
<br />soil moisture profIle is rapid. As will be shown later,
<br />this situation actually occurred when the specified
<br />values of t::.z, r, and h became too large and that of 00
<br />was too small. Despite this possible computational
<br />oscillation due to the preceding factors, Eq, 25 was
<br />applied to Eq, 20, for there seems no other available
<br />criterion which can be applied herein.
<br />
<br />In general, use of Eq. 25 for all cells of the
<br />finite-difference mesh (Figure 2) at each time level
<br />permits the computation of the respective Lit values
<br />for each cell and from them the one, whichever is
<br />smaller, will be selected as a I'i for next calculation,
<br />Therefore, given a ^ value which is arbitrarily chosen
<br />to be close to 0.5, the 6.t can be computed from Eq,
<br />25 as
<br />
<br />ll.t :S
<br />
<br />A(AZ)2
<br />Dmax
<br />
<br />. . ' , , ' , . , . . . (26)
<br />
<br />After the,6.t is determined, the computation of
<br />the unknown et 1 for various i and j values by use of
<br />Eqs, 20 through 24 can proceed in the following
<br />orderly way,
<br />
<br />For grid points other than i = I
<br />
<br />The values of y, and K at grid points i . I, i and i
<br />+ I at time level j, as shown in Figure 2, can be used
<br />to compute the 8 value at grid point i and time level j
<br />+ I from Eq, 20, Let (RHS)1 represent the right-hand
<br />side of Eq, 20, Then, from Eq, 20
<br />
<br />6J+1 K 6J + (RHs)j At
<br />i i i
<br />
<br />, , , .(27)
<br />
<br />Because the y, and K values in (R.~S~ are actually
<br />used in the computation of the oj+ value, any slight
<br />error in the y, computation at time j, especially
<br />around the saturatio,! front, would reflect in the
<br />computation of the oj+! value that may sometimes
<br />exceed 8s' If this situation happens, the 0 value at a
<br />grid point under consideration will then be set at 8 s
<br />and it will be assumed that the saturation front
<br />already arrived at or passed that grid point, The exact
<br />location of the saturation front can be determined by
<br />means of Eq, 17, If the saturation front so computed
<br />stays in between two grid points snch as i . I and i, it
<br />would be more accnrate to weigh the K value at
<br />respective grid points according to the exact location
<br />of the saturation front. A linear weighing techniqne
<br />
<br />20
<br />
|