2D groundwater FD scheme, boundary conditions

Derivation of the steady-state finite difference scheme

Consider a small element of a confined aquifer of transmissivity $T$. For steady flow, inflow equals outflow. With Darcy's law $q_x=-T_x\dfrac{\partial\phi}{\partial x}$, $q_y=-T_y\dfrac{\partial\phi}{\partial y}$ (flux per unit width):

$$\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}=0\ \Rightarrow\ \frac{\partial}{\partial x}\left(T_x\frac{\partial\phi}{\partial x}\right)+\frac{\partial}{\partial y}\left(T_y\frac{\partial\phi}{\partial y}\right)=0$$

For a homogeneous, isotropic aquifer ($T_x=T_y=T$ constant) this is the Laplace equation $\dfrac{\partial^2\phi}{\partial x^2}+\dfrac{\partial^2\phi}{\partial y^2}=0$.

Place a grid of spacing $\Delta x=\Delta y$ and write central differences about node $(i,j)$:

$$\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{\Delta x^2}+\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{\Delta y^2}=0$$ $$\boxed{\ \phi_{i,j}=\frac14\left(\phi_{i+1,j}+\phi_{i-1,j}+\phi_{i,j+1}+\phi_{i,j-1}\right)\ }$$

              (i,j+1)
                 |
   (i-1,j) --- (i,j) --- (i+1,j)
                 |
              (i,j-1)

Flow coefficients

Flow between node $(i,j)$ and a neighbour is written as a conductance times the head difference. The coefficients are

$$C_{E}=C_{W}=\frac{T_x\Delta y}{\Delta x},\qquad C_{N}=C_{S}=\frac{T_y\Delta x}{\Delta y}$$

so the flow from neighbour $k$ into the node is $Q_k=C_k(\phi_k-\phi_{i,j})$, and the steady-state balance $\sum Q_k=0$ gives

$$\phi_{i,j}=\frac{C_E\phi_{i+1,j}+C_W\phi_{i-1,j}+C_N\phi_{i,j+1}+C_S\phi_{i,j-1}}{C_E+C_W+C_N+C_S}$$

For an isotropic, homogeneous aquifer on a square grid, all four coefficients are equal ($=T$) and each weight is $\tfrac14$. For a layered aquifer the harmonic mean of the transmissivities of the two grids is used for the face between them.

Boundary conditions

1. Specified head (Dirichlet): the potential is known, e.g. along a river or reservoir, $\phi=h_b$. The node is not solved; its value enters the neighbour's equation.
2. No-flow (Neumann, $\partial\phi/\partial n=0$): along an impervious boundary or flow line. Use a mirror (image) node outside the boundary with the same potential as the inside neighbour. For a boundary on the left of node $(i,j)$:

$$\phi_{i-1,j}=\phi_{i+1,j}\ \Rightarrow\ \phi_{i,j}=\tfrac14\left(2\phi_{i+1,j}+\phi_{i,j+1}+\phi_{i,j-1}\right)$$

3. Specified flux: $q_n$ given (e.g. a pumping well or recharge): the known flow across that boundary is added to the node balance as a known inflow or outflow.
4. Head-dependent (river bed leakage): $Q=C_{riv}(h_{riv}-\phi)$ added to the node balance.

The resulting algebraic equations (one per interior node) are solved iteratively (Gauss-Seidel) or by matrix methods until the heads converge.