Saint-Venant, second-order explicit scheme

The dynamic-wave model uses the full Saint-Venant equations. A second-order accurate explicit scheme is obtained by replacing every derivative with a central difference in both space and time (the leapfrog scheme).

Governing equations

For a channel of top width $B$ (so $\partial A/\partial t=B\,\partial y/\partial t$):

$$B\frac{\partial y}{\partial t}+\frac{\partial Q}{\partial x}=q$$ $$\frac{\partial Q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{Q^2}{A}\right)+gA\frac{\partial y}{\partial x}-gA(S_0-S_f)=0$$

Grid and central differences

 t
 n+1 :        o (i, n+1)   unknown
 n   :  x-----x-----x      (i-1, i, i+1) known
 n-1 :        x (i, n-1)   known
       i-1    i    i+1   --> x

For any variable $f$ at $(i,n)$:

$$\frac{\partial f}{\partial t}\approx\frac{f_i^{n+1}-f_i^{n-1}}{2\Delta t},\qquad \frac{\partial f}{\partial x}\approx\frac{f_{i+1}^{n}-f_{i-1}^{n}}{2\Delta x}$$

Both are accurate to $O(\Delta t^2)$ and $O(\Delta x^2)$ (the $f''$ terms cancel in the Taylor series).

Finite difference equations

Continuity:

$$B\,\frac{y_i^{n+1}-y_i^{n-1}}{2\Delta t}+\frac{Q_{i+1}^n-Q_{i-1}^n}{2\Delta x}=q_i^n$$ $$y_i^{n+1}=y_i^{n-1}-\frac{\Delta t}{B\,\Delta x}\left(Q_{i+1}^n-Q_{i-1}^n\right)+\frac{2\Delta t}{B}q_i^n$$

Momentum:

$$\frac{Q_i^{n+1}-Q_i^{n-1}}{2\Delta t}+\frac{1}{2\Delta x}\left[\left(\frac{Q^2}{A}\right)_{i+1}^n-\left(\frac{Q^2}{A}\right)_{i-1}^n\right]+gA_i^n\frac{y_{i+1}^n-y_{i-1}^n}{2\Delta x}-gA_i^n\left(S_0-S_{f,i}\right)=0$$ $$Q_i^{n+1}=Q_i^{n-1}-\frac{\Delta t}{\Delta x}\left[\left(\frac{Q^2}{A}\right)_{i+1}^n-\left(\frac{Q^2}{A}\right)_{i-1}^n\right]-gA_i^n\frac{\Delta t}{\Delta x}\left(y_{i+1}^n-y_{i-1}^n\right)+2g\Delta t\,A_i^n\left(S_0-S_{f,i}\right)$$

with $S_{f,i}=\dfrac{n^2Q_i^{n-1}|Q_i^{n-1}|}{(A_i^{n-1})^2(R_i^{n-1})^{4/3}}$. The friction term is lagged to level $n-1$ (a central evaluation at $n$ would be unstable for this explicit scheme).

Solution procedure

1. Initial condition gives $Q_i^0,y_i^0$; the first step to level 1 is made with a one-step scheme (e.g. forward in time) because level $-1$ does not exist.
2. For each later step, compute $y_i^{n+1}$ and $Q_i^{n+1}$ at all interior points from the equations above.
3. Boundaries: upstream $Q(t)$ from the inflow hydrograph; downstream a stage or rating curve $Q=f(y)$.
4. Time step limit (Courant condition): $\Delta t\le\dfrac{\Delta x}{|V|+\sqrt{gA/B}}$.

Note: an equally common second-order explicit alternative is the two-step Lax-Wendroff/MacCormack scheme (predictor-corrector).