Channel Routing

1. Muskingum-Cunge Method in T-Route

1.1 Governing Equations of the Muskingum-Cunge Method in T-Route

T-Route has currently implemented two channel routing methods: the Muskingum Cunge (MC) and a diffusive wave. The MC is applied to only a synthetic channel cross section that is a compound channel cross section with a trapezoidal main and rectangular floodplain. The flow in each segment at each time step is updated using a finite difference equation that is derived using a continuity equation and a storage equation (The NCAR WRF-Hydro Modeling System Technical Description).

$$ \begin{align} Q_d^c = C_1 Q_u^p + C_2 Q_u^c + C_3 Q_d^p + (dt \quad dx\quad ql)/D \end{align} $$

, where

$$ \begin{align} C_1 = \frac{c(dt/dx) - 2X}{2(1-x) + c(dt/dx)} \\ C_2 = \frac{c(dt/dx) + 2X}{2(1-x) + c(dt/dx)} \\ C_3 = \frac{2(1-x) - c(dt/dx)}{2(1-x) + c(dt/dx)} \\ D = K(1-X) + dt/2 \\ c &= \text{wave celerity} \\ K &= \text{parameter for the time for a wave to travel} \\ X &= \text{weighting coefficient related to the effects of inflow and outflow on storage} \end{align} $$

The wave celerity for a trapezoidal main channel is

$$ \begin{align} c_m = \sqrt{\frac{S_o}{n}}\left( -\frac{2}{3}R^\frac{5}{3}\frac{2\sqrt{1 + z^2}}{B_w + 2 z h} + \frac{5}{3}R^\frac{2}{3} \right) \end{align} $$

, where

$$ \begin{align} S_o &= \text{channel slope} \\ n &= \text{Manning's roughness coefficient for a main channel} \\ R &= \text{hydraulic radius} \\ z &= \text{channel side slope} \\ B_w &= \text{channel bottom width} \\ h &= \text{channel depth from channel bottom until bank-full depth} \end{align} $$

Assuming $\frac{dP}{dA} \approx 0$, The wave celerity for a wide rectangular floodplain is

$$ \begin{align} c_f = \sqrt{\frac{S_o}{n_f}}\left( \frac{5}{3}h_f^\frac{2}{3} \right) \end{align} $$

, where

$$ \begin{align} n_f &= \text{Manning's roughness coefficient for a floodplain} \\ h_f &= \text{water depth after bank-full depth} \\ \end{align} $$

For overbank flow, the MC algorithm implemented within T-Route computes an weighted wave celerity by considering the celerities of the main channel and floodplain according to their relative areas. The weighting is based on the ratio of the main channel flow area to the total flow area and the ratio of the floodplain flow area to the total flow area.

Applying the aforementioned equations, the Muskingum-Cunge (MC) algorithm searches for the water depth that minimizes the difference between $Q_d^c$ and the uniform flow calculated using Manning's formula. This is achieved by repeatedly updating the water depth until the difference between successive approximations is within a desired tolerance or the maximum number of iterations is reached, employing the Secant method.

1.2 Implementation of the Muskingum-Cunge Method within T-Route

The Fortran compute kernel of the MC is MCsingeSegStime_f2py_NOLOOP.f90. Integration of the Fortran kernel with Python using Cython involves several steps:

1) Write a Fortran wrapper

pyMCsingleSegStime_NoLoop.f90 is the Fortran code that is used to create a C-compatible wrapper for the kernel. This wrapper helps to interface the Fortran code with Cython by providing C-callable functions.

2) Compile Fortran kernel and Build Shared Library

Use compiler.sh to run makefile to compile the kernel and Fortran wrapper into shared libraries.

3) Create a Cython Wrapper

Write a Cython wrapper file to interface with the compiled Fortran codes. reach.pyx is written in Python using Cython syntax.

4) Generate the Header File for Fortran

fortran_wrappers.pxd integrates the Fortran codes with C/C++ or Cython by providing the declarations for the Fortran functions that will be used in the C/C++ or Cython code.

5) Build Cython Extension Module

setup.py compiles the Cython module and link it with the Fortran object files.

6) Use Cython Extension Module

mc_reach.pyx can now import and use the compiled Fortran compute kernel for the MC method in Cython.

In mc_reach.pyx, the MC extension module troute.routing.fast_reach.reach is executed within the compute_reach_kernel Cython function. To facilitate parallel computing, T-Route originally made the assumption that $Q_u^c$ remain the same as $Q_u^p$ over a short time interval. As a result, the original finite difference equation

$$ \begin{align} Q_d^c = C_1 Q_u^p + C_2 Q_u^c + C_3 Q_d^p + (dt \quad dx\quad ql)/D \end{align} $$

becomes

$$ \begin{align} Q_d^c = C_1 Q_u^p + C_2 Q_u^p + C_3 Q_d^p + (dt \quad dx\quad ql)/D. \end{align} $$

, allowing, in theory, each stream segment within a given sub-network domain to simultaneously perform the MC computation without waiting for the completions of $Q_d^c$ computations from all upstream segments.

2. Diffusive Wave Routing in T-Route

2.1 Governing Equations of the Diffusive Wave in T-Route

The governing equation of the diffusive wave routing within T-Route is

$$ \begin{align} \frac{∂Q}{∂t} + C\frac{∂Q}{∂x} = D\frac{∂²Q}{∂x²} + Cq_l \end{align} $$

, where

$$ \begin{align} C = \frac{5/3 \left| S_f \right|^{0.3} \left| Q \right|^{0.4}}{B^{0.4}{\frac{1}{S_k}}^{0.6}} \\ D = \frac{Q}{2 S_f B} \end{align} $$

with

$$ \begin{align} S_f &= \text{energy slope} \\ B &= \text{water top width} \\ S_k &= \text{inverse of Manning's n}\\ C &= \text{wave celerity} \\ D &= \text{diffusivity coefficient} \end{align} $$

The partial differential equation is solved using the Crank-Nicolson scheme to non-uniform space grids, combined with an adapted Talyor series and Hermite Interpolation method. This approach is referred to as Crank-Nicolson over Space (CNS).

$∂u/∂t = c² ∂²u/∂x²$

$$ ( x = {a + b} ) $$