Appendix VIII

Appendix VIII
Numerical Solution of the Boundary-Value Problem Representing One-Dimensional Infiltration Above a Recharging Groundwater Flow System

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

\frac{\delta}{\delta z}\left[K(\psi)(\frac{\delta\psi}{\delta z}+ 1)\right] = C(\psi)\frac{\delta\psi}{\delta t} (A8.1)

and the boundary conditions are

\frac{\delta\psi}{\delta t} = \frac{R}{K(\psi)} - 1 (A8.2)

at the top, and

\frac{\delta\psi}{\delta z} = \frac{Q}{K_0} - 1 (A8.3)

at the base. The solution is in terms of the pressure head, \psi(z, t). The position of the water table at any time is given by the value of z at which \psi = 0. The required input includes the initial conditions, \psi(z, 0), the rainfall rate R, the groundwater recharge rate Q, and the unsaturated characteristic curves K(\psi) and C(\psi). For \psi \geq \psi_a, K = K_0 and C = 0, \psi_a 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 (z, t) plane is represented by a rectangular grid of points j = 1, 2, …, L along the z axis, and n = 1, 2, … along the t axis. The distance between the vertical nodes is\Delta z, and between the time steps, \Delta t. The solution at any given grid point (j, n) is \psi_j^n. Under this notation the finite-difference form for Eq. (A8.1) for an interior node (j = 2 \text{ to } j = L - 1), after suitable rearrangement, is

C(\psi_{\text{III}})\left(\frac{\psi_j^n - \psi_j^{n-1}}{\Delta t}\right) = \frac{1}{\Delta z}\left\{\left[K(\psi_{\text{I}})(1 + \frac{1}{2\Delta z})[\psi_{j+1}^{n-1} + \psi_{j+1}^{n} - \psi_{j}^{n-1} - \psi_j^n]\right]\right\} (A8.4)

where the values of \psi_{\text{I}}, \psi_{\text{II}}, \psi_{\text{III}}, which determine the K and C 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 (j = 1, j = L).

For each time step, the finite-difference approximations constitute a system of L linear algebraic equations in L unknowns. The general form is

-A_j\psi_{j+1}^n + B_j\psi_j^n - C_j\psi_{j-1}^n = D_j (A8.5)

where the coefficients A, B, C, and D vary with the node (j = 1, j = 2 \text{ to } L - 1, j = L), with the time step (n = 1, n = 2, n > 2), and with the saturation (\psi  \geq \psi_a, \psi < \psi_a). The variables on which A, B, C, and D depend are the \psi values from the previous time step, the boundary values R and Q, and the functional relationships K(\psi) and C(\psi).

The values of \psi_j^n are calculated from the recurrence relation

\psi_j^n = E_j\psi_{j+1}^n + F_j (A8.6)


E_j = \frac{A_j}{B_j - C_j E_{j-1}}

E_j = \frac{D_j + C_j F_{j-1}}{B_j - C_j E_{j-1}}

The E and F coefficients are calculated from j = 1 to j = L, and the \psi’s are back-calculated from j = L to j = 1 using Eq. (A8.6). For j = 1, E_1 = A_1/B_1 and F_1 = D_1/B_1. Complications arise if the upper node (j = L) is saturated.

If K(\psi) and C(\psi) 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 \psi value at which a change from wetting to drying or vice versa occurred.