Continuity equation; seepage under a dam algorithm

Continuity equation in groundwater flow

Continuity expresses conservation of mass. For a small control volume $\Delta x\,\Delta y\,\Delta z$ of porous medium,

$$\text{inflow}-\text{outflow}=\text{rate of change of water stored}$$

With Darcy fluxes $q_x,q_y,q_z$, specific storage $S_s$ and a source/sink term $W$ (recharge positive):

$$-\left(\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}+\frac{\partial q_z}{\partial z}\right)+W=S_s\frac{\partial h}{\partial t}$$

Substituting Darcy's law ($q_x=-K_x\,\partial h/\partial x$, and similarly for $y$ and $z$) gives the groundwater flow equation

$$\frac{\partial}{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial}{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+\frac{\partial}{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)=S_s\frac{\partial h}{\partial t}-W$$

For steady flow ($\partial h/\partial t=0$) without sources in a homogeneous, isotropic soil, it reduces to the Laplace equation $\nabla^2h=0$, which is used for seepage under a dam.

Algorithm for seepage under a dam

1. Define the problem: draw the vertical section, choose a square grid ($\Delta x=\Delta y$) covering the dam foundation down to the impervious layer, and number the nodes.
2. Read data: upstream head $h_1$, downstream head $h_2$, hydraulic conductivity $K$, grid spacing, tolerance $\varepsilon$ and maximum number of iterations.
3. Boundary conditions: set $\phi=h_1$ on the upstream bed nodes and $\phi=h_2$ on the downstream bed nodes; mark the dam base and the impervious layer as no-flow boundaries (use mirror nodes).
4. Initialise all interior potentials, e.g. by linear interpolation between $h_1$ and $h_2$.
5. Iterate: for each interior node compute

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

(using the no-flow form on impervious boundaries). Record the largest change $\Delta_{max}$. 6. Test: if $\Delta_{max}>\varepsilon$, go to step 5; otherwise stop. 7. Seepage rate: compute $Q=K\sum_j(\phi_{i,j}-\phi_{i+1,j})$ over a vertical section below the dam (per unit length of dam), and check it at another section. 8. Output: the potentials at all nodes (for equipotential lines and the flow net) and the seepage $Q$.