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 \times dx \times 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 to build the Cython extension module for the compute kernel.

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 \times dx \times 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 \times dx \times 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). Using the values of $C$ and $D$ computed from the previous time step, the subroutine mesh_diffusive_forward for CNS in diffusive.f90 computes Q at each compute node for the subsequent time step, with the computation progressing from upstream to downstream. Once the Q computation reaches the tailwater of a sub-network at a given time step, the process mesh_diffusive_backward then computes water depth using the computed Q values and the downstream boundary condition for water depth, with the computation now progressing from upstream to downstream. The equation used is

$$ \begin{align} y_i(t+1) = y_{i+1}(t+1) + \left| S_f \right| \times \Delta x \times sign(Q_i^{n+1}). \end{align} $$

Hence, Q is computed from upstream to downstream, followed by computing y from downstream to upstream at each time step, and this process repeats until the end of the simulation period.

2.2 Implementation of the Diffusive Wave Routing within T-Route

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

1) Write a Fortran wrapper

pydiffusive.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. diffusive.pyx is written in Python using Cython syntax to build the Cython extension module for the compute kernel.

4) Generate the Header File for Fortran

fortran_wrappers.pxd and pydiffusive.h 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.

2.3 Channel Connectivity Crosswalk Between Python and Fortran

While the MC method obtains channel network connectivity information, including flow direction, from reach extension module at every time step, the diffusive wave routing retrieves this information from frnw_g in diffusive_utils_v02.py at the start of the simulation. It then maintains this information, along with channel geometry data, within the Fortran compute kernel for the entire duration of the simulation.

In Figure 1, each stream reach is assigned an integer number based on frnw_g, indicating the sequence in which the diffusive wave routing is applied. Each reach is colored differently according to its respective stream junction number. A stream junction is a confluence point where two or more stream reaches meet and combine their flow. The stream junction order reflects the total number of stream junctions downstream of a given reach.

Reaches with indices from 1 to 6 share a common junction order of four and do not transfer water to each other. As a result, channel routing can be applied to these reaches either simultaneously or serially. Once channel routing for the reaches with a junction order of four is complete, routing for the reaches with one less junction order—those with indices from 7 to 11—follows. Similarly, these reaches can also be computed either simultaneously or serially. This process continues until channel routing is applied to reaches with a junction order of zero.

This approach reflects the computation direction for channel routing methods, such as the Muskingum-Cunge (MC) and diffusive wave methods used within T-Route, which moves from upstream to downstream. This direction ensures that flow calculations for upstream reaches are completed before those for downstream reaches.

Conversely, when computation from downstream to upstream is required—such as in calculating water depth based on downstream boundary conditions—the computation direction is reversed, starting from reaches with the least junction order and progressing to those with the highest.

Table 1 provides the contents of frnw_g. The row index corresponds to the sequence of reaches in which diffusive wave routing is applied. The first column displays the number of compute nodes in a stream reach, which is always one more than the number of stream segments within the reach. This is because each compute node represents a terminal point of a segment, with two co-located terminal points—one from the upstream segment and one from the downstream segment—treated as a single terminal point.

The second column indicates the reach index of the reach immediately downstream of the reach associated with the current row index. For example, the stream reach represented by row index 1 corresponds to the reach in the top-left corner with index 1 in Figure 1, and the reach immediately downstream has an index of 7.

The third column indicates the total number of upstream reaches converging at the upstream terminal point of the reach associated with the row index. When this number is greater than one, additional columns are included, each representing the reach index of one of the upstream reaches.

3. Hybrid Routing within T-Route

Currently, the hybrid routing methodology in T-Route allows for the dynamic exchange of output between the Muskingum-Cunge (MC) and reservoir modules during each real-time computation step. Following this, the diffusive wave module is executed separately on the diffusive channel domain, overwriting the previously computed results with the more accurate diffusive wave output. This approach is applied specifically to diffusive domains that are downstream of MC and reservoir domains. However, if a non-diffusive domain is situated downstream of a diffusive domain, it does not received the updated outflow from the upstream diffusive domain. This issue is illustrated in Figure 2. In the diagram illustrating the sequential execution, the reaches marked with red circle do not receive updated flow from the upstream diffusive domain, even though the flow values in the light blue domain has been updated by the diffusive wave output.

However, a prototype that fully integrates the diffusive wave module with the MC and reservoir modules, allowing dynamic exchange of output among all modules, is currently under development. This prototype ensures that the output from the diffusive domain is properly carried downstream to non-diffusive domains, while also enabling the output from non-diffusive domains, such as water depth, to be carried upstream to the diffusive domain during each real-time computation step.