Finite difference method

A finite difference scheme is of order $p$ if its truncation error (the terms dropped from the Taylor series) is proportional to $\Delta x^p$ (or $\Delta t^p$). Higher order means the error falls faster when the grid is refined.

First-order accurate schemes

From Taylor's series,

$$\frac{f_{i+1}-f_i}{\Delta x}=f'_i+\frac{\Delta x}{2}f''_i+\dots\ \Rightarrow\ f'_i=\frac{f_{i+1}-f_i}{\Delta x}+O(\Delta x)$$

The forward difference $(f_{i+1}-f_i)/\Delta x$ and the backward difference $(f_i-f_{i-1})/\Delta x$ are first-order accurate. Applied to $\partial Q/\partial t+c\,\partial Q/\partial x=0$ with forward time and backward space:

$$\frac{Q_i^{n+1}-Q_i^n}{\Delta t}+c\frac{Q_i^n-Q_{i-1}^n}{\Delta x}=0,\qquad \text{error } O(\Delta t,\Delta x)$$

Second-order accurate schemes

Subtract the two Taylor series for $f_{i+1}$ and $f_{i-1}$; the $f''$ terms cancel:

$$\frac{f_{i+1}-f_{i-1}}{2\Delta x}=f'_i+\frac{\Delta x^2}{6}f_i^{(3)}+\dots\ \Rightarrow\ f'_i=\frac{f_{i+1}-f_{i-1}}{2\Delta x}+O(\Delta x^2)$$

This central difference is second-order accurate. Likewise $f''_i=\dfrac{f_{i+1}-2f_i+f_{i-1}}{\Delta x^2}+O(\Delta x^2)$. A second-order scheme in both variables is the four-point (box) scheme centred at $(i+\tfrac12,\ n+\tfrac12)$:

$$\frac{Q_i^{n+1}+Q_{i+1}^{n+1}-Q_i^n-Q_{i+1}^n}{2\Delta t}+c\,\frac{Q_{i+1}^{n+1}+Q_{i+1}^n-Q_i^{n+1}-Q_i^n}{2\Delta x}=0$$

Graphical meaning

 f
 |                     . f(i+1)
 |                 . '
 |            f(i)
 |         .'
 |   f(i-1)
 +----|---------|---------|----> x
     i-1        i        i+1

 forward : chord i  -> i+1
 backward: chord i-1-> i
 central : chord i-1-> i+1 (closest to tangent at i)

The slope of the chord approximates the slope of the tangent at $i$. The forward and backward chords are tilted either side of the tangent (error proportional to $\Delta x$); the central chord is almost parallel to it (error proportional to $\Delta x^2$).

Grid view of the two PDE schemes

 t
 n+1 :   o--------o            o   o
         |        |              \ |
 n   :   o        x         x-----x-----x   
        (i-1)    (i)       (i-1)  (i)  (i+1)
   first order (forward t,      second order (central
   backward x): 2 known pts     in x): 3 known points

Comparison

The finite difference method (FDM) is a numerical technique that replaces the continuous solution domain by a grid of discrete points and replaces the derivatives in a differential equation by differences between the values at neighbouring grid points. The differential equation becomes a set of algebraic equations in the unknown grid values.

Explicit and implicit schemes (example: $\partial Q/\partial t+c\,\partial Q/\partial x=0$, $C=c\Delta t/\Delta x$)

Explicit: the spatial derivative is evaluated at the known time level $n$.

$$\frac{Q_i^{n+1}-Q_i^n}{\Delta t}+c\,\frac{Q_i^n-Q_{i-1}^n}{\Delta x}=0\ \Rightarrow\ Q_i^{n+1}=Q_i^n-C\,(Q_i^n-Q_{i-1}^n)$$

Each unknown is found directly from known values. It is conditionally stable (needs $C\le1$).

Implicit: the spatial derivative is evaluated at the new time level $n+1$.

$$\frac{Q_i^{n+1}-Q_i^n}{\Delta t}+c\,\frac{Q_i^{n+1}-Q_{i-1}^{n+1}}{\Delta x}=0\ \Rightarrow\ Q_i^{n+1}=\frac{Q_i^n+C\,Q_{i-1}^{n+1}}{1+C}$$

The unknown at point $i$ depends on other unknowns at $n+1$. For problems such as groundwater flow the result is a set of simultaneous equations (tridiagonal matrix). It is unconditionally stable.

Comparison

Example (groundwater, $\lambda=T\Delta t/S\Delta x^2$): explicit $h_i^{n+1}=h_i^n+\lambda(h_{i+1}^n-2h_i^n+h_{i-1}^n)$ needs $\lambda\le\tfrac12$; implicit $-\lambda h_{i-1}^{n+1}+(1+2\lambda)h_i^{n+1}-\lambda h_{i+1}^{n+1}=h_i^n$ is stable for any $\lambda$.

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).

Numerical diffusion

The leading truncation error of an odd-order scheme (e.g. first-order upwind) behaves like a physical diffusion term. For $Q_t+cQ_x=0$ with backward-space, forward-time:

$$Q_t+cQ_x=\frac{c\Delta x}{2}(1-C)\,Q_{xx}+\dots$$

The solution is smeared: peaks are lowered and sharp fronts spread out, although the true equation has no diffusion.

Numerical dispersion

For second-order (central) schemes the leading error term has an odd derivative, $\propto Q_{xxx}$. Different wavelengths then travel at different speeds, so a sharp front produces spurious oscillations (ripples ahead of or behind the front). It is a phase error, with little change in amplitude.

Stability

A scheme is stable if errors (round-off or truncation) do not grow from step to step. For explicit schemes, the von Neumann or Courant-Friedrichs-Lewy condition requires

$$C=\frac{c\,\Delta t}{\Delta x}\le1$$

Implicit schemes (and the four-point scheme with $\theta\ge0.5$) are unconditionally stable, but a large $\Delta t$ still reduces accuracy.

Four-point implicit method

$$\frac{\partial Q}{\partial t}\approx\frac{(Q_i^{n+1}-Q_i^n)+(Q_{i+1}^{n+1}-Q_{i+1}^n)}{2\Delta t}$$ $$\frac{\partial Q}{\partial x}\approx\theta\,\frac{Q_{i+1}^{n+1}-Q_i^{n+1}}{\Delta x}+(1-\theta)\,\frac{Q_{i+1}^n-Q_i^n}{\Delta x}$$

Data: $Q_i^{n+1}=75$, $Q_{i+1}^{n+1}=72$, $Q_i^n=73$, $Q_{i+1}^n=69$; $\Delta t=3600$ s, $\Delta x=500$ m, $\theta=0.55$.

Time derivative

$$\frac{\partial Q}{\partial t}=\frac{(75-73)+(72-69)}{2\times3600}=\frac{5}{7200}=6.9444\times10^{-4}\ \text{m}^3/\text{s}^2$$

Space derivative

$$\frac{\partial Q}{\partial x}=0.55\,\frac{72-75}{500}+0.45\,\frac{69-73}{500}=-0.00330+(-0.00360)=-0.00690\ \text{m}^3/\text{s per m}$$

Answer: $\partial Q/\partial t=6.9444\times10^{-4}$ m$^3$/s$^2$ and $\partial Q/\partial x=-0.00690$ m$^2$/s.

A simple explicit scheme is the Lax diffusive scheme: forward in time, central in space, with the old value at $i$ replaced by the average of its two neighbours. The equations for a rectangular channel are

$$B\frac{\partial y}{\partial t}+\frac{\partial Q}{\partial x}=0,\qquad \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$$

Finite difference form (grid $(i,n+1)$ from $(i-1,n),(i,n),(i+1,n)$):

$$y_i^{n+1}=\frac{y_{i+1}^n+y_{i-1}^n}{2}-\frac{\Delta t}{2B\Delta x}\left(Q_{i+1}^n-Q_{i-1}^n\right)$$ $$Q_i^{n+1}=\frac{Q_{i+1}^n+Q_{i-1}^n}{2}-\frac{\Delta t}{2\Delta x}\left[\frac{(Q^2)_{i+1}}{A_{i+1}}-\frac{(Q^2)_{i-1}}{A_{i-1}}\right]-gA_i\frac{\Delta t}{2\Delta x}\left(y_{i+1}-y_{i-1}\right)+gA_i\Delta t\,(S_0-S_{f,i})$$

Data

$B=50$ m, $S_0=0.0002$, $n=0.035$, $\Delta x=1000$ m, $\Delta t=600$ s, so $\dfrac{\Delta t}{2\Delta x}=0.300$.

Step 1: Areas and friction slope

$A_{i-1}=93.50$, $A_i=91.00$, $A_{i+1}=98.50$ m$^2$.

$$P_i=B+2y_i=53.64\ \text{m},\quad R_i=\frac{A_i}{P_i}=1.6965\ \text{m},\quad S_{f,i}=\frac{n^2Q_i^2}{A_i^2R_i^{4/3}}=9.475\times10^{-5}$$

Step 2: Flow depth

$$y_i^{n+1}=\frac{1.97+1.87}{2}-0.00600\,(35-38)=1.9200-(-0.0180)=1.9380\ \text{m}$$

Step 3: Discharge

• Average term: $\dfrac{35+38}{2}=36.500$
• Convective term: $\dfrac{Q^2}{A}$ at $i+1$ is 12.437, at $i-1$ is 15.444; difference $=-3.0073$; term $=-0.300(-3.0073)=0.902$
• Pressure term: $-gA_i\dfrac{\Delta t}{2\Delta x}(y_{i+1}-y_{i-1})=-9.81(91.00)(0.300)(0.10)=-26.781$
• Gravity and friction: $gA_i\Delta t(S_0-S_f)=9.81(91.00)(600)(0.0002-9.475\times10^{-5})=56.374$

$$Q_i^{n+1}=36.500+(0.902)+(-26.781)+(56.374)=66.99\ \text{m}^3/\text{s}$$

Answer: $y_i^{n+1}=1.9380$ m and $Q_i^{n+1}=66.99$ m$^3$/s.

Remark: the very large change in $Q$ comes from the time step: $\Delta t=600$ s gives a Courant number $(V+\sqrt{gy})\Delta t/\Delta x\approx2.8>1$, so the explicit scheme is outside its stability limit and a smaller $\Delta t$ would be needed to continue the simulation. One step is computed here as asked.

Numerical diffusion

The leading truncation error of an odd-order scheme (e.g. first-order upwind) behaves like a physical diffusion term. For $Q_t+cQ_x=0$ with backward-space, forward-time:

$$Q_t+cQ_x=\frac{c\Delta x}{2}(1-C)\,Q_{xx}+\dots$$

The solution is smeared: peaks are lowered and sharp fronts spread out, although the true equation has no diffusion.

Numerical dispersion

For second-order (central) schemes the leading error term has an odd derivative, $\propto Q_{xxx}$. Different wavelengths then travel at different speeds, so a sharp front produces spurious oscillations (ripples ahead of or behind the front). It is a phase error, with little change in amplitude.

Stability

A scheme is stable if errors (round-off or truncation) do not grow from step to step. For explicit schemes, the von Neumann or Courant-Friedrichs-Lewy condition requires

$$C=\frac{c\,\Delta t}{\Delta x}\le1$$

Implicit schemes (and the four-point scheme with $\theta\ge0.5$) are unconditionally stable, but a large $\Delta t$ still reduces accuracy.

Four-point implicit method

$$\frac{\partial Q}{\partial t}\approx\frac{(Q_i^{n+1}-Q_i^n)+(Q_{i+1}^{n+1}-Q_{i+1}^n)}{2\Delta t}$$ $$\frac{\partial Q}{\partial x}\approx\theta\,\frac{Q_{i+1}^{n+1}-Q_i^{n+1}}{\Delta x}+(1-\theta)\,\frac{Q_{i+1}^n-Q_i^n}{\Delta x}$$

Data: $Q_i^{n+1}=95$, $Q_{i+1}^{n+1}=92$, $Q_i^n=93$, $Q_{i+1}^n=89$; $\Delta t=3600$ s, $\Delta x=600$ m, $\theta=0.55$.

Time derivative

$$\frac{\partial Q}{\partial t}=\frac{(95-93)+(92-89)}{2\times3600}=\frac{5}{7200}=6.9444\times10^{-4}\ \text{m}^3/\text{s}^2$$

Space derivative

$$\frac{\partial Q}{\partial x}=0.55\,\frac{92-95}{600}+0.45\,\frac{89-93}{600}=-0.00275+(-0.00300)=-0.00575\ \text{m}^3/\text{s per m}$$

Answer: $\partial Q/\partial t=6.9444\times10^{-4}$ m$^3$/s$^2$ and $\partial Q/\partial x=-0.00575$ m$^2$/s.

A simple explicit scheme is the Lax diffusive scheme: forward in time, central in space, with the old value at $i$ replaced by the average of its two neighbours. The equations for a rectangular channel are

$$B\frac{\partial y}{\partial t}+\frac{\partial Q}{\partial x}=0,\qquad \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$$

Finite difference form (grid $(i,n+1)$ from $(i-1,n),(i,n),(i+1,n)$):

$$y_i^{n+1}=\frac{y_{i+1}^n+y_{i-1}^n}{2}-\frac{\Delta t}{2B\Delta x}\left(Q_{i+1}^n-Q_{i-1}^n\right)$$ $$Q_i^{n+1}=\frac{Q_{i+1}^n+Q_{i-1}^n}{2}-\frac{\Delta t}{2\Delta x}\left[\frac{(Q^2)_{i+1}}{A_{i+1}}-\frac{(Q^2)_{i-1}}{A_{i-1}}\right]-gA_i\frac{\Delta t}{2\Delta x}\left(y_{i+1}-y_{i-1}\right)+gA_i\Delta t\,(S_0-S_{f,i})$$

Data

$B=10$ m, $S_0=0.0002$, $n=0.04$, $\Delta x=1000$ m, $\Delta t=300$ s, so $\dfrac{\Delta t}{2\Delta x}=0.150$.

Step 1: Areas and friction slope

$A_{i-1}=19.00$, $A_i=18.50$, $A_{i+1}=20.00$ m$^2$.

$$P_i=B+2y_i=13.70\ \text{m},\quad R_i=\frac{A_i}{P_i}=1.3504\ \text{m},\quad S_{f,i}=\frac{n^2Q_i^2}{A_i^2R_i^{4/3}}=0.00452$$

Step 2: Flow depth

$$y_i^{n+1}=\frac{2.0+1.9}{2}-0.01500\,(37.5-40)=1.9500-(-0.0375)=1.9875\ \text{m}$$

Step 3: Discharge

• Average term: $\dfrac{37.5+40}{2}=38.750$
• Convective term: $\dfrac{Q^2}{A}$ at $i+1$ is 70.312, at $i-1$ is 84.211; difference $=-13.8980$; term $=-0.150(-13.8980)=2.085$
• Pressure term: $-gA_i\dfrac{\Delta t}{2\Delta x}(y_{i+1}-y_{i-1})=-9.81(18.50)(0.150)(0.10)=-2.722$
• Gravity and friction: $gA_i\Delta t(S_0-S_f)=9.81(18.50)(300)(0.0002-0.00452)=-235.358$

$$Q_i^{n+1}=38.750+(2.085)+(-2.722)+(-235.358)=-197.25\ \text{m}^3/\text{s}$$

Answer: $y_i^{n+1}=1.9875$ m and $Q_i^{n+1}=-197.25$ m$^3$/s.

Remark: the very large change in $Q$ comes from the time step: $\Delta t=300$ s gives a Courant number $(V+\sqrt{gy})\Delta t/\Delta x\approx1.9>1$, so the explicit scheme is outside its stability limit and a smaller $\Delta t$ would be needed to continue the simulation. One step is computed here as asked.

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.

The non-linear kinematic wave model is used with $A=\alpha Q^\beta$ and the second-order four-point scheme (centred in space and time), solved by Newton's method.

Step 1: Channel properties at the mean flow

Mean discharge $=(420+580)/2=500$ m$^3$/s. Trapezoid: bottom width $b=350$ m, side slope $z=6$, $S_0=0.015$, $n=0.0359$.

$$Q=\frac{S_0^{1/2}}{n}\frac{A^{5/3}}{P^{2/3}},\quad A=(b+zy)y,\quad P=b+2y\sqrt{1+z^2}$$

Solving for $Q=500$: $y=0.592$ m, $A=209.30$ m$^2$, $P=357.20$ m.

Step 2: Kinematic coefficients

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

Then $A_i^n=\alpha(420)^{0.6}=188.51$ m$^2$ and $A_i^{n+1}=\alpha(580)^{0.6}=228.80$ m$^2$. (Wave celerity $c\approx3.98$ m/s, so $C=c\Delta t/\Delta x\approx4.6>1$: an explicit scheme would be unstable, the four-point scheme is not.)

Step 3: Finite difference equation

Known: $Q_i^n=420$ (initial flow), $Q_{i+1}^n=420$, $Q_i^{n+1}=580$ (upstream flood). $\Delta x=4650$ m, $\Delta t=1.5\times3600=5400$ s, $\Delta t/\Delta x=1.1613$, no lateral inflow.

$$\frac{\Delta t}{\Delta x}x+\alpha x^{0.6}=K,\qquad K=\alpha\left[420^{0.6}+420^{0.6}-580^{0.6}\right]-\frac{\Delta t}{\Delta x}\left(420-580-420\right)=821.779$$

with $x=Q_{i+1}^{n+1}$.

Step 4: Newton iteration

$f(x)=1.1613x+5.0280x^{0.6}-821.779$, $f'(x)=1.1613+3.0168x^{-0.4}$.

Answer: $Q_{i+1}^{n+1}=522.6$ m$^3$/s at 4.65 km downstream after 1.5 h (it rises from 420 towards the upstream value of 580 m$^3$/s).

Governing equations of fluid movement

The unsteady flow of water in an open channel is governed by the Saint-Venant (1D shallow water) equations.

Continuity (conservation of mass):

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q$$

Momentum (conservation of momentum):

$$\underbrace{\frac{\partial Q}{\partial t}}_{\text{local acc.}}+\underbrace{\frac{\partial}{\partial x}\!\left(\frac{Q^2}{A}\right)}_{\text{convective acc.}}+\underbrace{gA\frac{\partial y}{\partial x}}_{\text{pressure force}}-\underbrace{gAS_0}_{\text{gravity}}+\underbrace{gAS_f}_{\text{friction}}=0$$

where $A$ = flow area, $Q$ = discharge, $y$ = flow depth, $q$ = lateral inflow per unit length, $S_0$ = bed slope, $S_f=\dfrac{n^2Q|Q|}{A^2R^{4/3}}$ = friction slope (Manning), $g$ = gravity. These are the Saint-Venant equations (1D unsteady open-channel flow).

In terms of velocity $V$ and depth $y$ for a wide rectangular channel:

$$\frac{\partial y}{\partial t}+V\frac{\partial y}{\partial x}+y\frac{\partial V}{\partial x}=0,\qquad \frac{\partial V}{\partial t}+V\frac{\partial V}{\partial x}+g\frac{\partial y}{\partial x}=g(S_0-S_f)$$

Forward, backward and central differencing

The grid has points $x_i=i\Delta x$. The derivative at point $i$ is approximated using neighbouring values.

Expand $f$ about $x_i$ by Taylor's series:

$$f_{i+1}=f_i+\Delta x f'_i+\frac{\Delta x^2}{2}f''_i+\frac{\Delta x^3}{6}f_i^{(3)}+\dots,\qquad f_{i-1}=f_i-\Delta x f'_i+\frac{\Delta x^2}{2}f''_i-\frac{\Delta x^3}{6}f_i^{(3)}+\dots$$

The same formulas are used for the time derivative with $\Delta t$ (e.g. forward in time: $\partial Q/\partial t\approx (Q_i^{n+1}-Q_i^n)/\Delta t$).

 f
 |                     . f(i+1)
 |                 . '
 |            f(i)
 |         .'
 |   f(i-1)
 +----|---------|---------|----> x
     i-1        i        i+1

 forward : chord i  -> i+1
 backward: chord i-1-> i
 central : chord i-1-> i+1 (closest to tangent at i)

• Forward uses the point ahead; backward uses the point behind (suited to flow from upstream, "upwind"); central uses both neighbours and is more accurate, but is less stable for explicit time marching.

The linear kinematic wave equation is $\dfrac{\partial Q}{\partial t}+c\dfrac{\partial Q}{\partial x}=0$ (no lateral inflow), with a constant celerity $c$ taken from the initial uniform flow. It is solved with the second-order four-point scheme, which uses all four grid points $(i,n)$, $(i+1,n)$, $(i,n+1)$, $(i+1,n+1)$.

Step 1: Initial uniform flow (US units, Manning constant 1.49)

$$Q=\frac{1.49}{n}AR^{2/3}S_0^{1/2},\qquad A=200y,\ P=200+2y,\ R=A/P$$

Trial for $Q=2000$ cfs gives $y=1.680$ ft, so $A=336.09$ ft$^2$, $P=203.36$ ft, $R=1.653$ ft and $Q=2000$ cfs (check).

$$V=\frac{2000}{336.09}=5.951\ \text{ft/s}$$

Step 2: Kinematic wave celerity

With Manning's equation $Q=\alpha A^{5/3}$, $c=\dfrac{dQ}{dA}=\dfrac53V$:

$$c=\tfrac53(5.951)=9.918\ \text{ft/s},\qquad C=\frac{c\,\Delta t}{\Delta x}=\frac{9.918\times180}{3000}=0.5951$$

Step 3: Finite difference equation

$$\frac{Q_i^{n+1}+Q_{i+1}^{n+1}-Q_i^n-Q_{i+1}^n}{2\Delta t}+c\,\frac{Q_{i+1}^{n+1}+Q_{i+1}^n-Q_i^{n+1}-Q_i^n}{2\Delta x}=0$$

Solving for the unknown:

$$Q_{i+1}^{n+1}=Q_i^n+\frac{1-C}{1+C}\left(Q_{i+1}^n-Q_i^{n+1}\right)$$

Step 4: Substitution

$Q_i^n=2000$, $Q_{i+1}^n=2000$ (uniform flow at $t=0$), $Q_i^{n+1}=2250$ (upstream, $t=3$ min); the point at 3000 ft is $i+1$.

$$Q_{i+1}^{n+1}=2000+\frac{1-0.5951}{1+0.5951}(2000-2250)=2000+(0.2539)(-250)=1936.5\ \text{cfs}$$

Answer: the discharge at 3000 ft downstream after 3 min is $Q=1936.5$ cfs (about 1937 cfs); the flood wave has only started to arrive.

Governing equations

For 1D unsteady open-channel flow (Saint-Venant equations):

Continuity (conservation of mass):

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q$$

Momentum (conservation of momentum):

$$\underbrace{\frac{\partial Q}{\partial t}}_{\text{local acc.}}+\underbrace{\frac{\partial}{\partial x}\!\left(\frac{Q^2}{A}\right)}_{\text{convective acc.}}+\underbrace{gA\frac{\partial y}{\partial x}}_{\text{pressure force}}-\underbrace{gAS_0}_{\text{gravity}}+\underbrace{gAS_f}_{\text{friction}}=0$$

where $A$ = flow area, $Q$ = discharge, $y$ = flow depth, $q$ = lateral inflow per unit length, $S_0$ = bed slope, $S_f=\dfrac{n^2Q|Q|}{A^2R^{4/3}}$ = friction slope (Manning), $g$ = gravity. These are the Saint-Venant equations (1D unsteady open-channel flow).

In terms of velocity $V$ and depth $y$ for a wide rectangular channel:

$$\frac{\partial y}{\partial t}+V\frac{\partial y}{\partial x}+y\frac{\partial V}{\partial x}=0,\qquad \frac{\partial V}{\partial t}+V\frac{\partial V}{\partial x}+g\frac{\partial y}{\partial x}=g(S_0-S_f)$$

Kinematic wave approximation

In many flood-routing problems on steep slopes the local acceleration, convective acceleration and pressure terms are small compared with the gravity and friction terms. Dropping them, the momentum equation becomes

$$S_0-S_f=0\quad\Rightarrow\quad S_f=S_0$$

so the flow is always at its normal-depth value for the local area and a unique relation exists:

$$Q=\alpha_kA^{\beta},\qquad \beta=\frac53\ \text{(Manning)}$$

Together with continuity, $\dfrac{\partial A}{\partial t}+\dfrac{\partial Q}{\partial x}=q$, this gives the kinematic wave equation

$$\frac{\partial Q}{\partial t}+c\,\frac{\partial Q}{\partial x}=c\,q,\qquad c=\frac{dQ}{dA}=\beta V$$

The wave moves downstream at speed $c$ without attenuation (no backwater effects).

Courant number

The Courant number is the ratio of the distance a wave travels in one time step to the grid spacing:

$$C=\frac{c\,\Delta t}{\Delta x}$$

For explicit schemes, stability requires $C\le1$. For the dynamic wave, $c=V\pm\sqrt{gy}$.