Tridiagonal matrix, river stage and water table

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$$

Each equation has three unknowns, so a set of $N$ simultaneous equations must be solved at each time step.

Boundary conditions

• River side: $h_0^{j+1}=h_L$ (the river stage, which may change with time). In the first equation the known term moves to the right-hand 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.