Consider a one-dimensional, vertical, homogeneous flow field with its upper boundary at the ground surface and its lower boundary below the water table. As outlined in Section 6.4, the equation of flow for such a system is
and the boundary conditions are
at the top, and
at the base. The solution is in terms of the pressure head, . The position of the water table at any time is given by the value of at which . The required input includes the initial conditions, , the rainfall rate , the groundwater recharge rate , and the unsaturated characteristic curves and . For , and , being the air entry pressure head.
The numerical finite-difference scheme used by Rubin and Steinhardt (1963), Liakopoulos (1965b), and Freeze (1969b) is the implicit scheme of Richtmyer (1957). With this method the plane is represented by a rectangular grid of points along the axis, and along the axis. The distance between the vertical nodes is, and between the time steps, . The solution at any given grid point is . Under this notation the finite-difference form for Eq. (A8.1) for an interior node , after suitable rearrangement, is
where the values of , , , which determine the and values to be applied for a particular node at a particular time step, are determined by extrapolation from previous time steps as described by Rubin and Steinhardt (1963).
Similar finite-difference equations, which incorporate the boundary conditions (A8.2) and (A8.3), can be written for the nodes at the upper and lower boundary .
For each time step, the finite-difference approximations constitute a system of linear algebraic equations in unknowns. The general form is
where the coefficients , , , and vary with the node , with the time step , and with the saturation . The variables on which , , , and depend are the values from the previous time step, the boundary values and , and the functional relationships and .
The values of are calculated from the recurrence relation
The and coefficients are calculated from to , and the ’s are back-calculated from to using Eq. (A8.6). For , and . Complications arise if the upper node is saturated.
If and are hysteretic, they are usually built into the computer program in the form of a table of values representing the main wetting, main drying, and principal scanning curves. The program must locate the correct curve on the basis of whether the node in question is wetting or drying, and by scanning the past history of the node to determine the value at which a change from wetting to drying or vice versa occurred.