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

(A9.1)

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

(A9.2)

where Tx and Ty are the principal components of the transmissivity tensor defined by T = Kb, and S is the storativity defined by S = Ssb. 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:

(A9.3)

where the subscript means that the bracketed quantity is evaluated at the midpoint between nodes (i, j) and \left(T_x\frac{\partial h}{\partial 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}\left(T_x\frac{\partial h}{\partial 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}(T_x)_{i+1/2,j}(T_x)_{i-1/2,j}(T_x)_{i+1/2,j} \simeq \frac{(T_x)_{i+1,j} + (T_x)_{i,j}}{2}\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial 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\}\frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial 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\}S\frac{\partial h}{\partial t} \simeq S_{i,j}\left(\frac{h_{i,j}^k – h_{i,j}^{k-1}}{\Delta t}\right)Ah_{i,j}^k = Bh_{i,j-1}^k + Ch_{i+1,j}^k + Dh_{i,j+1}^k + Eh_{i-1,j}^k + FA = \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}\Delta x = \Delta yA = \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}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.