Second-order non-linear kinematic wave; FDM principle

Non-linear kinematic wave model

Continuity: $\dfrac{\partial A}{\partial t}+\dfrac{\partial Q}{\partial x}=q$. Kinematic approximation ($S_f=S_0$) with Manning's equation gives a unique $Q$-$A$ relation, $Q=\dfrac{S_0^{1/2}}{nP^{2/3}}A^{5/3}$, or

$$A=\alpha Q^{\beta},\qquad \alpha=\left(\frac{nP^{2/3}}{S_0^{1/2}}\right)^{0.6},\quad \beta=0.6$$

Therefore

$$\frac{\partial(\alpha Q^\beta)}{\partial t}+\frac{\partial Q}{\partial x}=q$$

Space-time discretisation (second order, four-point/box)

Centre the derivatives at the middle of the cell, $(i+\tfrac12,\ n+\tfrac12)$:

 t
 n+1 :  (i,n+1) o------o (i+1,n+1)   <- Q(i,n+1) known,
                |  *   |                  Q(i+1,n+1) unknown
 n   :  (i,n)   o------o (i+1,n)    <- both known
              i        i+1       --> x

$$\frac{\partial A}{\partial t}\approx\frac{(A_i^{n+1}-A_i^n)+(A_{i+1}^{n+1}-A_{i+1}^n)}{2\Delta t},\qquad \frac{\partial Q}{\partial x}\approx\frac{(Q_{i+1}^{n+1}-Q_i^{n+1})+(Q_{i+1}^n-Q_i^n)}{2\Delta x},\qquad q\approx\bar q$$

Taylor expansion about the cell centre shows that both errors are $O(\Delta t^2)$ and $O(\Delta x^2)$, so the scheme is second-order accurate. Substituting $A=\alpha Q^\beta$:

$$\frac{\alpha\left[(Q_i^{n+1})^\beta+(Q_{i+1}^{n+1})^\beta-(Q_i^n)^\beta-(Q_{i+1}^n)^\beta\right]}{2\Delta t}+\frac{Q_{i+1}^{n+1}+Q_{i+1}^n-Q_i^{n+1}-Q_i^n}{2\Delta x}=\bar q$$

Solution for the one unknown $Q_{i+1}^{n+1}$

Multiply by $2\Delta t$ and collect the unknown $x=Q_{i+1}^{n+1}$ on the left:

$$f(x)=\frac{\Delta t}{\Delta x}\,x+\alpha x^{\beta}-K=0$$ $$K=\alpha\left[(Q_i^n)^\beta+(Q_{i+1}^n)^\beta-(Q_i^{n+1})^\beta\right]-\frac{\Delta t}{\Delta x}\left(Q_{i+1}^n-Q_i^{n+1}-Q_i^n\right)+2\Delta t\,\bar q$$

Because $x$ appears non-linearly, solve by Newton-Raphson:

$$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)},\qquad f'(x)=\frac{\Delta t}{\Delta x}+\alpha\beta x^{\beta-1}$$

starting from $x_0=Q_{i+1}^n$ (or $Q_i^{n+1}$) and repeating until $|x_{k+1}-x_k|$ is small. Marching along $i$ from the upstream boundary and then in time gives the whole hydrograph.

Principle of the finite difference method

The solution domain is covered by a grid; every derivative in the differential equation is replaced by a difference of grid-point values (from the Taylor series), converting the PDE into algebraic equations that are solved for the unknown values at the grid points.