Steady 2D model: seepage under a dam

Physical situation

   h1 (upstream)      dam        h2 (downstream)
 ~~~~~~~~~~~~~~~~ [=========] ~~~~~~~~~~~~~~
 phi=h1  . . . . . no flow . . . . phi=h2
 . . . . . . . . . . . . . . . . . . . . .
 . . . . . . . . . . . . . . . . . . . . .
 ///////// impervious layer (no flow) /////

Water seeps under the dam through the foundation soil from the upstream bed (head $h_1$) to the downstream bed (head $h_2$). The flow is steady and two-dimensional in a vertical section.

Governing equation

Continuity for steady flow, $\dfrac{\partial q_x}{\partial x}+\dfrac{\partial q_y}{\partial y}=0$, with Darcy's law $q_x=-K_x\dfrac{\partial\phi}{\partial x}$, $q_y=-K_y\dfrac{\partial\phi}{\partial y}$ ($\phi$ = potential or hydraulic head). For a homogeneous, isotropic aquifer ($K_x=K_y=K$):

$$\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=0\qquad\text{(Laplace equation)}$$

Finite difference expression

Cover the section by a square grid ($\Delta x=\Delta y$). Central differences from Taylor's series:

$$\frac{\partial^2\phi}{\partial x^2}\approx\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{\Delta x^2},\qquad \frac{\partial^2\phi}{\partial y^2}\approx\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{\Delta y^2}$$

Adding and setting $\Delta x=\Delta y$:

$$\boxed{\ \phi_{i,j}=\frac14\left(\phi_{i+1,j}+\phi_{i-1,j}+\phi_{i,j+1}+\phi_{i,j-1}\right)\ }$$

The potential at each grid point is the average of its four neighbours.

Boundary conditions

• Upstream bed: $\phi=h_1$; downstream bed: $\phi=h_2$ (specified head).
• Base of the dam and the impervious layer: no flow, $\partial\phi/\partial n=0$. Use an image (ghost) point: for a point on a horizontal no-flow boundary, $\phi_{i,j-1}=\phi_{i,j+1}$ so

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

• Far left and right ends: assumed no-flow (or fixed heads equal to $h_1$, $h_2$), placed far enough from the dam.

Iterative procedure (Gauss-Seidel / Liebmann)

1. Set up the grid and number the nodes; give the boundary nodes their values.
2. Guess the potential at all interior nodes (for example, a linear variation from $h_1$ to $h_2$).
3. Sweep the grid row by row. At each interior node compute $\phi_{i,j}^{new}$ with the equation above, using the latest available neighbour values. Optionally over-relax: $\phi^{new}=\phi^{old}+\omega(\phi^{GS}-\phi^{old})$, $1<\omega<2$.
4. Find the largest change $\max|\phi^{new}-\phi^{old}|$ over the grid.
5. If it is greater than the tolerance $\varepsilon$ (e.g. 0.001 m), repeat from step 3; otherwise stop.
6. Seepage rate. Between two adjacent nodes the flow through one cell face (per unit length of dam) is

$$q=K\,\frac{\phi_a-\phi_b}{\Delta x}\,\Delta y=K(\phi_a-\phi_b)\quad(\Delta x=\Delta y)$$

Choose a vertical section between two columns $i$ and $i+1$ below the dam; the total seepage is the sum over all rows:

$$Q=K\sum_{j}\left(\phi_{i,j}-\phi_{i+1,j}\right)$$

(check with a different section; the totals must agree if the solution has converged, and also equal the flow entering at the upstream bed). 7. Print the potentials (to draw equipotential lines and the flow net) and the seepage $Q$ in m$^3$/s per metre length of dam.