Appendix IX

Appendix IX
Development of Finite-Difference Equation for Transient Flow in a Heterogeneous, Anisotropic, Horizontal, Confined Aquifer

The partial differential equation that describes transient flow through a saturated anisotropic medium (Section 2.11) is

\frac{\delta}{\delta x} \left(K_x\frac{\delta h}{\delta x}\right) + \frac{\delta}{\delta y} \left(K_y\frac{\delta h}{\delta y}\right) + \frac{\delta}{\delta z} \left(K_z\frac{\delta h}{\delta z}\right) = S_s\frac{\delta h}{\delta t} (A9.1)

For a horizontal, confined aquifer of thickness, b, the two-dimensional form of Eq. (A9.1) reduces to

\frac{\delta}{\delta x} \left(T_x\frac{\delta h}{\delta x}\right) + \frac{\delta}{\delta y} \left(T_y\frac{\delta h}{\delta y}\right) = S_s\frac{\delta h}{\delta t} (A9.2)

where T_x and T_y are the principal components of the transmissivity tensor defined by T = Kb, and S is the storativity defined by S = S_sb. To find the finite-difference equation for an interior node in the nodal grid used to discretize the region of flow, we must replace the partial derivatives in Eq. (A9.2) by differences. Using the definitions developed in Appendix VI and the ijk notation of Section 8.8 and Figure 8.26(c), we can write the finite-difference expression:

\frac{\delta}{\delta x} \left(T_x\frac{\delta h}{\delta x}\right) \simeq \frac{1}{\Delta x}\left[\left(T_x\frac{\delta h}{\delta x}\right)_{i+1/2,j}^k - \left(T_x\frac{\delta h}{\delta x}\right)_{i-1/2,j}^k\right] (A9.3)

where the subscript (i + \frac{1}{2}, j) means that the bracketed quantity is evaluated at the midpoint between nodes (i, j) and (i + 1, j), and the superscript k means that the bracketed quantity is evaluated at time step k. We can further approximate terms on the right-hand side of Eq. (A9.3) by

\left(T_x\frac{\delta h}{\delta x}\right)_{i+1/2,j}^k = (T_x)_{i+1/2,j}(h_{i+1,j}^k - h_{i,j}^k)\frac{1}{\Delta x} (A9.4A)

\left(T_x\frac{\delta h}{\delta x}\right)_{i-1/2,j}^k = (T_x)_{i-1/2,j}(h_{i,j}^k - h_{i-1,j}^k)\frac{1}{\Delta x} (A9.4B)

If we evaluate (T_x)_{i+1/2,j} and (T_x)_{i-1/2,j} by simple averages of the form

(T_x)_{i+1/2,j} \simeq \frac{(T_x)_{i+1,j} + (T_x)_{i,j}}{2} (A9.5)

then these expressions can be substituted in Eqs. (A9.4), and Eqs. (A9.4) can in turn be substituted in Eq. (A9.3) to give

\frac{\delta}{\delta x} \left(T_x\frac{\delta h}{\delta x}\right) \simeq \frac{1}{2\Delta x^2} \left\{\left[(T_x)_{i+1/2,j} + (T_x)_{i,j}\right]h_{i+1,j}^k - \left[(T_x)_{i+1/2,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i,j}^k + \left[(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i-1,j}^k \right\} (A9.6)


\frac{\delta}{\delta y} \left(T_y\frac{\delta h}{\delta y}\right) \simeq \frac{1}{2\Delta y^2} \left\{\left[(T_y)_{i+1/2,j} + (T_y)_{i,j}\right]h_{i+1,j}^k - \left[(T_y)_{i+1/2,j} + 2(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i,j}^k + \left[(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i-1,j}^k \right\} (A9.7)

Finally, we can approximate the right-hand side of Eq. (A9.2) by

S\frac{\delta h}{\delta t} \simeq S_{i,j}\left(\frac{h_{i,j}^k - h_{i,j}^{k-1}}{\Delta t}\right) (A9.8)

Substituting Eqs. (A9.6), (A9.7), and (A9.8) for the three terms in Eqs. (A9.6) and gathering terms leads to the general finite-difference equation for an internal node in a heterogeneous, anisotropic aquifer:

Ah_{i,j}^k = Bh_{i,j-1}^k + Ch_{i+1,j}^k + Dh_{i,j+1}^k + Eh_{i-1,j}^k + F (A9.9)


A = \frac{1}{2\Delta x}^2[(T_x)_{i+1,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}] + \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + 2(T_y)_{i,j} + (T_y)_{i,j-1}] + \frac{S_{i,j}}{\Delta t}

B = \frac{1}{2\Delta y^2}[(T_y)_{i,j} + (T_y)_{i,j-1}]

C = \frac{1}{2\Delta x^2}[(T_x)_{i+1,j} + (T_x)_{i,j}]

D = \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + (T_y)_{i,j}]

E = \frac{1}{2\Delta x^2}[(T_x)_{i,j} + (T_x)_{i-1,j}]

F = \frac{S_{i,j}}{\Delta t} \cdot h_{i,j}^{k-1}

If the aquifer is homogeneous and isotropic, then T_x = T_y = T for all (i, j) and S_{i, j} = S for all (i, j). Under these conditions, and for a square nodal grid with \Delta x = \Delta y, the coefficients of Eq. (A9.9) become

A = \frac{4T}{\Delta x^2} + \frac {S_{i,j}}{\Delta t}

B = C = D = E = \frac{T}{\Delta x^2}

F = \frac{S}{\Delta t} \cdot h_{i,j}^{k-1}

If we divide through by T/\Delta x^2, these coefficients are seen to be the same as those developed in a less rigorous way in Section 8.8 and presented in connection with Eq. (8.57).

Trescott et al. (1976) have suggested that there is some advantage to utilizing the harmonic mean rather than the arithmetic mean in Eq. (A9.5). This approach changes the coefficients in the finite-difference equation but does not alter the concepts underlying the development.

Equation (A9.9) is written in terms of the hydraulic head values at five nodes at time step k and one node at time step k - 1. It is known as a backward-difference approximation. Remson et al. (1971) note that there are some computational advantages to using a central-difference approximation, known as the Crank-Nicholson scheme, that utilizes head values at five nodes at time step k and five nodes at time step k - 1. The alternating-direction implicit procedure (ADIP) utilized by Pinder and Bredehoeft (1968) involves two finite-difference equations, one in the xt plane and one in the yt plane. Each uses head values at three nodes at time step k and three nodes at time step k - 1.