Simulation of Ground water flow

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.

A tridiagonal matrix arises when the 1D groundwater equation is written implicitly for every grid along a line from the river to the aquifer boundary.

 river                                     barrier (no flow)
  hL  |  grid1  grid2  grid3  ...  gridN  |
  ~~~~|---o------o------o-------o------o--|
      0   1      2      3      N-1    N  N+1
      |<-dx->|

Governing equation

For a one-dimensional unconfined/confined aquifer next to a river (x measured from the river), with transmissivity $T$, storage coefficient $S$ and head $h$:

$$T\frac{\partial^2h}{\partial x^2}=S\frac{\partial h}{\partial t}$$

Implicit finite difference form

Divide the aquifer into $N$ grids of width $\Delta x$ (grid 1 next to the river). Use a backward difference in time and a central difference in space, with the space terms at the new time level $j+1$:

$$T\,\frac{h_{i-1}^{j+1}-2h_i^{j+1}+h_{i+1}^{j+1}}{\Delta x^2}=S\,\frac{h_i^{j+1}-h_i^j}{\Delta t}$$

With $\lambda=\dfrac{T\Delta t}{S\Delta x^2}$:

$$-\lambda\,h_{i-1}^{j+1}+(1+2\lambda)\,h_i^{j+1}-\lambda\,h_{i+1}^{j+1}=h_i^j,\qquad i=1,\dots,N$$

This has three unknowns, so a set of $N$ simultaneous equations is needed at each time step.

Boundary conditions

• River side: $h_0^{j+1}=h_L$ (river stage, may change with time). In the first equation the known term moves to the right side: $(1+2\lambda)h_1-\lambda h_2=h_1^j+\lambda h_L$.
• Barrier (no flow) at the far end: $\partial h/\partial x=0\Rightarrow h_{N+1}=h_N$, so the last equation is $-\lambda h_{N-1}+(1+\lambda)h_N=h_N^j$. (If the far end has a specified head $h_R$, the last row is instead $-\lambda h_{N-1}+(1+2\lambda)h_N=h_N^j+\lambda h_R$.)

Tridiagonal coefficient matrix

$$\begin{bmatrix} 1+2\lambda&-\lambda&&&\\ -\lambda&1+2\lambda&-\lambda&&\\ &\ddots&\ddots&\ddots&\\ &&-\lambda&1+2\lambda&-\lambda\\ &&&-\lambda&1+\lambda \end{bmatrix} \begin{bmatrix}h_1\\h_2\\ \vdots\\h_{N-1}\\h_N\end{bmatrix}^{j+1} = \begin{bmatrix}h_1^j+\lambda h_L\\h_2^j\\ \vdots\\h_{N-1}^j\\h_N^j\end{bmatrix}$$

The matrix is symmetric, tridiagonal and diagonally dominant, so the system has a unique solution and can be solved efficiently by the Thomas algorithm.

Thomas algorithm for the solution

Write the rows as $a_ih_{i-1}+b_ih_i+c_ih_{i+1}=d_i$ with $a_i=c_i=-\lambda$. Forward elimination:

$$c_1'=\frac{c_1}{b_1},\ d_1'=\frac{d_1}{b_1};\qquad c_i'=\frac{c_i}{b_i-a_ic_{i-1}'},\ \ d_i'=\frac{d_i-a_id_{i-1}'}{b_i-a_ic_{i-1}'}$$

Back substitution: $h_N=d_N'$, $h_i=d_i'-c_i'h_{i+1}$. The new heads become the right-hand side for the next time step, with $h_L$ updated to the new river stage.

The 1D implicit model predicts how the water table in an aquifer beside a river responds to changes in river stage. It treats the aquifer as a line of grids perpendicular to the river.

 river stage hL(t)           aquifer                barrier
   ~~~~~~~~~~|====================================|
             h1     h2     h3     ...     hN

Assumptions

Homogeneous, isotropic aquifer; one-dimensional horizontal flow; constant $T$ and $S$ (Dupuit assumptions); no recharge or pumping; river fully penetrating and connected, so the head at the river boundary equals the river stage.

Governing equation

For a one-dimensional unconfined/confined aquifer next to a river (x measured from the river), with transmissivity $T$, storage coefficient $S$ and head $h$:

$$T\frac{\partial^2h}{\partial x^2}=S\frac{\partial h}{\partial t}$$

Implicit finite difference form

Divide the aquifer into $N$ grids of width $\Delta x$ (grid 1 next to the river). Use a backward difference in time and a central difference in space, with the space terms at the new time level $j+1$:

$$T\,\frac{h_{i-1}^{j+1}-2h_i^{j+1}+h_{i+1}^{j+1}}{\Delta x^2}=S\,\frac{h_i^{j+1}-h_i^j}{\Delta t}$$

With $\lambda=\dfrac{T\Delta t}{S\Delta x^2}$:

$$-\lambda\,h_{i-1}^{j+1}+(1+2\lambda)\,h_i^{j+1}-\lambda\,h_{i+1}^{j+1}=h_i^j,\qquad i=1,\dots,N$$

This has three unknowns, so a set of $N$ simultaneous equations is needed at each time step.

Boundary conditions

• River side: $h_0^{j+1}=h_L$ (river stage, may change with time). In the first equation the known term moves to the right side: $(1+2\lambda)h_1-\lambda h_2=h_1^j+\lambda h_L$.
• Barrier (no flow) at the far end: $\partial h/\partial x=0\Rightarrow h_{N+1}=h_N$, so the last equation is $-\lambda h_{N-1}+(1+\lambda)h_N=h_N^j$. (If the far end has a specified head $h_R$, the last row is instead $-\lambda h_{N-1}+(1+2\lambda)h_N=h_N^j+\lambda h_R$.)

Tridiagonal coefficient matrix

$$\begin{bmatrix} 1+2\lambda&-\lambda&&&\\ -\lambda&1+2\lambda&-\lambda&&\\ &\ddots&\ddots&\ddots&\\ &&-\lambda&1+2\lambda&-\lambda\\ &&&-\lambda&1+\lambda \end{bmatrix} \begin{bmatrix}h_1\\h_2\\ \vdots\\h_{N-1}\\h_N\end{bmatrix}^{j+1} = \begin{bmatrix}h_1^j+\lambda h_L\\h_2^j\\ \vdots\\h_{N-1}^j\\h_N^j\end{bmatrix}$$

The matrix is symmetric, tridiagonal and diagonally dominant, so the system has a unique solution and can be solved efficiently by the Thomas algorithm.

Procedure for the river-stage problem

1. Compute $\lambda=\dfrac{T\Delta t}{S\Delta x^2}$.
2. Give initial water-table heads $h_i^0$ and the river stage $h_L$ at the new time.
3. Form the right-hand side vector and solve the tridiagonal system (Thomas algorithm) for $h_i^{j+1}$.
4. Use these as initial heads for the next step; repeat for all time steps with the new river stage.
5. Flux exchanged with the river at any time: $q=-T\dfrac{h_1-h_L}{\Delta x}$ per unit length of river (positive = flow towards the aquifer when $h_L>h_1$).

Stability and remarks

• The explicit form of the same equation requires $\lambda\le\tfrac12$, which forces very small $\Delta t$. The implicit form is unconditionally stable, so large time steps can be used.
• Accuracy still falls if $\Delta t$ is very large (first-order in time).
• A rising river stage ($h_L>h$) pushes the water table up from the river side (bank storage); a falling stage drains the aquifer towards the river.

The aquifer is divided into three grids ($1200/400=3$). The river head $h_L=365$ m is applied at the river boundary next to grid 1, and the barrier at the right is a no-flow boundary ($h_4=h_3$). Recharge is zero.

Step 1: Coefficient

$$\lambda=\frac{T\Delta t}{S\Delta x^2}=\frac{600\times1}{0.15\times400^2}=0.025$$

Step 2: Implicit equations

$$-\lambda h_{i-1}^{j+1}+(1+2\lambda)h_i^{j+1}-\lambda h_{i+1}^{j+1}=h_i^j$$

• Grid 1: $(1+2\lambda)h_1-\lambda h_2=h_1^0+\lambda h_L=349.13+0.025(365)=358.255$
• Grid 2: $-\lambda h_1+(1+2\lambda)h_2-\lambda h_3=347.97$
• Grid 3 (barrier): $-\lambda h_2+(1+\lambda)h_3=339.36$

$$\begin{bmatrix} 1.050 & -0.025 & 0 \\ -0.025 & 1.050 & -0.025 \\ 0 & -0.025 & 1.025 \end{bmatrix}\begin{bmatrix}h_1\\h_2\\h_3\end{bmatrix}=\begin{bmatrix}358.255\\347.97\\339.36\end{bmatrix}$$

Step 3: Solution (Thomas algorithm)

Forward elimination: $c_1'=-0.0238$, $d_1'=341.1952$; $c_2'=-0.0238$, $d_2'=339.7163$; $d_3'=339.5660$. Back substitution: $h_3=339.57$; $h_2=d_2'-c_2'h_3=347.81$; $h_1=d_1'-c_1'h_2=349.48$.

Check in grid 2: $-0.025(349.48)+1.05(347.81)-0.025(339.57)=347.97\approx347.97$.

Answer: $h_1=349.48$ m, $h_2=347.81$ m, $h_3=339.57$ m.

The subscripts in the data are partly unclear. They are read as the four neighbours of the central grid $(i,j)$: $h_{i+1,j}=50$ m, $h_{i-1,j}=75$ m (along $x$) and $h_{i,j+1}=95$ m, $h_{i,j-1}=10$ m (along $y$).

Step 1: Finite difference equation

For steady flow without withdrawal, $K_{xx}\dfrac{\partial^2h}{\partial x^2}+K_{yy}\dfrac{\partial^2h}{\partial y^2}=0$:

$$K_{xx}\frac{h_{i+1,j}-2h_{i,j}+h_{i-1,j}}{\Delta x^2}+K_{yy}\frac{h_{i,j+1}-2h_{i,j}+h_{i,j-1}}{\Delta y^2}=0$$ $$h_{i,j}=\frac{\dfrac{K_{xx}}{\Delta x^2}\left(h_{i+1,j}+h_{i-1,j}\right)+\dfrac{K_{yy}}{\Delta y^2}\left(h_{i,j+1}+h_{i,j-1}\right)}{2\left(\dfrac{K_{xx}}{\Delta x^2}+\dfrac{K_{yy}}{\Delta y^2}\right)}$$

Step 2: Coefficients

$$\frac{K_{xx}}{\Delta x^2}=\frac{0.95\times10^{-5}}{100^2}=9.500\times10^{-10},\qquad \frac{K_{yy}}{\Delta y^2}=\frac{0.99\times10^{-5}}{80^2}=1.547\times10^{-9}$$

Step 3: Hydraulic head at the central grid

$$h_{i,j}=\frac{9.500\times10^{-10}(50+75)+1.547\times10^{-9}(95+10)}{2(9.500\times10^{-10}+1.547\times10^{-9})}=\frac{2.812\times10^{-7}}{4.994\times10^{-9}}=56.30\ \text{m}$$

Step 4: Discharges into the central grid (per unit thickness)

Flow from a neighbour into the central grid (Darcy's law, per unit thickness): $Q=K\dfrac{h_{nb}-h_{i,j}}{\Delta}\times(\text{face width})$. The face width is $\Delta y=80$ m for the $x$-neighbours and $\Delta x=100$ m for the $y$-neighbours.

A negative value means flow out of the central grid towards that neighbour. Check: the four discharges add up to zero (continuity satisfied). Multiply by the saturated thickness to get m$^3$/s.

Answer: $h_{i,j}=56.30$ m. Inflows: from $(i-1,j)$ $1.421\times10^{-4}$, from $(i,j+1)$ $4.789\times10^{-4}$ m$^3$/s per m; outflows to $(i+1,j)$ $4.792\times10^{-5}$ and to $(i,j-1)$ $5.730\times10^{-4}$ m$^3$/s per m.

Courant condition

The Courant (Courant-Friedrichs-Lewy) condition is the stability limit for an explicit finite difference scheme: the numerical information must not travel more than one grid in one time step. For wave-type equations, $C=\dfrac{c\Delta t}{\Delta x}\le1$. For the groundwater (diffusion) equation the corresponding number is

$$\lambda=\frac{T\,\Delta t}{S\,\Delta x^2}\le\frac12\quad\text{(explicit scheme)}$$

The implicit scheme used below has no such limit.

Coefficient

"Transmissivity $=0.024$ m/s" is read as $0.024$ m$^2$/s (a transmissivity has units of m$^2$/s).

$$\lambda=\frac{T\Delta t}{S\Delta x^2}=\frac{0.024\times86400}{0.015\times100^2}=\frac{2073.6}{150}=13.824$$

Since $\lambda\gg\frac12$, an explicit scheme would be unstable; the implicit scheme is used.

Implicit equations for the 5 grids

$$-\lambda h_{i-1}^{j+1}+(1+2\lambda)h_i^{j+1}-\lambda h_{i+1}^{j+1}=h_i^j$$

Coefficients: sub- and super-diagonal $=-\lambda=-13.824$; diagonal $=1+2\lambda=28.648$; the last diagonal (no-flow barrier, $h_6=h_5$) $=1+\lambda=14.824$.

Right-hand side: grid 1 includes the river head, $190.10+\lambda(195)=190.10+2695.68=2885.78$; other grids use the initial heads.

Matrix form

$$\begin{bmatrix} 28.648 & -13.824 & 0 & 0 & 0 \\ -13.824 & 28.648 & -13.824 & 0 & 0 \\ 0 & -13.824 & 28.648 & -13.824 & 0 \\ 0 & 0 & -13.824 & 28.648 & -13.824 \\ 0 & 0 & 0 & -13.824 & 14.824 \end{bmatrix}\begin{bmatrix}h_1\\h_2\\h_3\\h_4\\h_5\end{bmatrix}=\begin{bmatrix}2885.78\\190.20\\190.30\\190.40\\190.50\end{bmatrix}$$

Solution (Thomas algorithm)

$$h_1=194.02,\ h_2=193.31,\ h_3=192.84,\ h_4=192.54,\ h_5=192.41\ \text{m}$$

(the water table rises most near the river, as expected for a river stage of 195 m).

Answer: coefficients $--13.824$, $28.648$ (and $14.824$ in the last row); the heads after one day are as above.

Darcy's law for a confined (transmissivity) aquifer gives the flux per unit width across the face between two grids as

$$q=T\,\frac{\phi_{\text{from}}-\phi_{\text{to}}}{\Delta}$$

where $\Delta$ is the grid spacing in the direction of flow. Here each flux is measured out of the central grid in the direction of its arrow (positive = leaving the central grid, negative = entering it). The face width is $\Delta Y$ for $x$-direction fluxes and $\Delta X$ for $y$-direction fluxes.

              j-1  (qB, up)
               ^
  i-1 (qC) <-- (i,j) --> (qA) i+1
               v
              j+1  (qD, down)

Fluxes ($T_x=0.013$, $T_y=0.015$ m$^2$/s; $\Delta X=20$ m, $\Delta Y=25$ m)

$$q_A=T_x\frac{\phi_{i,j}-\phi_{i+1,j}}{\Delta X}=0.013\,\frac{13-14}{20}=-6.500\times10^{-4}\ \text{m}^2/\text{s}$$ $$q_B=T_y\frac{\phi_{i,j}-\phi_{i,j-1}}{\Delta Y}=0.015\,\frac{13-13.5}{25}=-3.000\times10^{-4}\ \text{m}^2/\text{s}$$ $$q_C=T_x\frac{\phi_{i,j}-\phi_{i-1,j}}{\Delta X}=0.013\,\frac{13-12}{20}=6.500\times10^{-4}\ \text{m}^2/\text{s}$$ $$q_D=T_y\frac{\phi_{i,j}-\phi_{i,j+1}}{\Delta Y}=0.015\,\frac{13-11}{25}=1.200\times10^{-3}\ \text{m}^2/\text{s}$$

Flow through each face (flux $\times$ face width)

Answer: $q_A=-6.500\times10^{-4}$, $q_B=-3.000\times10^{-4}$, $q_C=6.500\times10^{-4}$, $q_D=1.200\times10^{-3}$ m$^2$/s. Negative signs mean water flows into the central grid from the A and B sides; water leaves through C and D.

Check: the net outflow is $\Sigma Q=0.018$ m$^3$/s $\ne0$, so the given potentials are not a steady-state set. For exact steady state the central potential would be $\phi_{i,j}=12.681$ m (from $\sum$ fluxes $=0$).

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.

Continuity equation in groundwater flow

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

$$\text{net inflow}-\text{net 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=-K\,\partial h/\partial x$ 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$, 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 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$, maximum 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 impervious layer as no-flow (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: potentials at all nodes (equipotential lines, flow net) and the seepage $Q$.