HydPy-Musk-MCT (Muskingum-Cunge-Todini routing)

The HydPy-Musk family member musk_mct realises a modification of the variable coefficient Muskingum approach of Cunge (1969) (MC) proposed by Todini (2007) (MCT).

Opposed to the original Muskingum method of (McCarthy, 1940), MC does not rely on predetermined coefficients. Instead, it calculates the coefficients for each simulation step based on a water depth-discharge relationship. Essentially, MC is a kinematic wave approach that does not explicitly account for pressure gradients or other forcings than the channel’s bottom slope explicitly. However, it adjusts its coefficients so that the numerical error of the underlying finite difference discretisation (the Muskingum working formula) mimics the physical diffusion of flood waves due to pressure gradients as well as possible. Hence, it can simulate flood waves even in channels with mild bottom slopes without much calibration (provided that information on the rating curve is available). Compared to more complete solutions of the Saint-Vernant equations, it is less computationally demanding. On the downside, it cannot deal with cases where, for example, effects or rapid changes in channel geometry play a significant role.

musk_mct implements the Muskingum-Cunge modification of Todini (2007), which solves a “mass conservation inconsistency” and a “steady-state inconsistency” of the original formulation. Todini (2007) evaluates the method for simple rectangular, triangular, and trapezoidal cross-section profiles in combination with the Manning-Strickler equation. musk_mct also allows to define more complex profiles by selecting and configuring submodels that follow the CrossSectionModel_V1 interface. However, be aware that the achieved performance or accuracy might be worse than expected if the underlying expectations are violated.

When applying musk_mct, one should be aware that it relies on two nested iterations. First, it uses a root search algorithm to find the current “reference water depth” that results in the previously determined “reference discharge”. While Todini (2007) uses the Newton-Raphson approach, musk_mct applies the Pegasus iteration (Dowell and Jarratt, 1972). By default, musk_mct strives for very high accuracies. There is a good chance you can increase the computation speed significantly in your use case without reducing numerical precision too much (see the documentation on method Calc_ReferenceWaterDepth_V1 for further information). The second iteration addresses the “reference discharge” itself, which’s initial guess can be quite inaccurate (see Calc_ReferenceDischarge_V1). Todini (2007) suggests refining it by starting subsequent runs of the whole MCT procedure (including the Newton-Raphson or Pegasus approach described above). musk_mct follows Todini’s suggestion to run everything twice. You can deviate from this default behaviour by changing the value of parameter NmbRuns.

Integration tests

Note

When new to HydPy, consider reading section Integration Tests first.

All of the following tests are equal to or similar to the examples provided by Todini (2007).

The simulation period covers four days; the simulation step is 30 minutes:

>>> from hydpy import pub, Nodes, Element
>>> pub.timegrids = "2000-01-01", "2000-01-05", "30m"

The model shall retrieve its inflow from two Node objects (input1 and input2) and passe its outflow to node output. We define these nodes and connect them to an Element object that handles a musk_mct instance:

>>> nodes = Nodes("input1", "input2", "output")
>>> stream = Element("stream", inlets=["input1", "input2"], outlets="output")
>>> from hydpy.models.musk_mct import *
>>> parameterstep("1h")
>>> stream.model = model

The last general preparation step is to define a test function object for controlling the following test runs:

>>> from hydpy.core.testtools import IntegrationTest
>>> test = IntegrationTest(stream)
>>> test.plotting_options.axis1 = (states.discharge,)
>>> test.dateformat = "%H:%M"

Base example

This example repeats the base example of Todini (2007) dealing with a trapezoidal profile. A channel with a total length of 100 km is subdivided into 50 segments:

>>> nmbsegments(50)
>>> length(100.0)
>>> bottomslope(0.00025)

We use the wq_trapeze_strickler submodel to configure the shape of a simple trapeze:

>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(1)
...     bottomlevels(0.0)
...     bottomwidths(15.0)
...     sideslopes(5.0)
...     stricklercoefficients(1.0/0.035)

musk_mct uses the catchment area for calculating the value of the solver parameter ToleranceDischarge, which determines the accuracy of the iterative estimation of the reference water depth (see methods modify_init() and Calc_ReferenceWaterDepth_V1). We set it 25 000 km², which is approximately the catchment size of the river Main at gauge Frankfurt-Osthafen:

>>> catchmentarea(25000.0)

This setting makes up for an error tolerance of 0.025 m³/s, which is 360 times smaller than the lowest discharge reported for Frankfurt-Ostenhagen (9 m³/s, 1976):

>>> solver.tolerancedischarge.update()
>>> solver.tolerancedischarge
tolerancedischarge(0.025)

Besides “old” discarge values, musk_mct requires the Courant number and the cell Reynolds number as initial conditions. These states must be consistent if you need good simulation results right from the start. If there are none available from a previous simulation run, you can use prepare_states() to calculate them for definable discharge values under the assumption of stationarity:

>>> model.prepare_states(discharge=100.0)
>>> test.inits = ((states.discharge, 100.0),
...               (states.courantnumber, states.courantnumber[0]),
...               (states.reynoldsnumber, states.reynoldsnumber[0]))

Node input1 provides a constant baseflow of 100 m³/s, and node input2 provides a synthetic flood wave with a peak discharge of 800 m³/s, making a total peak discharge of 900 m³/s (which is similar to the average flood discharge at Frankfurt-Ostenhagen) that occurs 24 hours after the start of the simulation period:

>>> import numpy
>>> q_base, q_peak, t_peak, beta = 100.0, 900.0, 24.0, 16.0
>>> ts = pub.timegrids.init.to_timepoints()
>>> nodes.input1.sequences.sim.series = q_base
>>> nodes.input2.sequences.sim.series = (
...     (q_peak - q_base) * ((ts / t_peak) * numpy.exp(1.0 - ts / t_peak)) ** beta)

The following table and figure show the discharge calculated for all 51 channel segment endpoints, including the start and endpoint of the complete channel. Other properties like the reference water depth, the wetted perimeter, or the Courant number represent estimated averages for all 50 segments:

>>> conditions = test("musk_mct_base_example", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

Todini (2007) reports a peak discharge simulated for the channel outlet of 643.72 m³/s, perfectly agreeing with our results:

>>> from hydpy import round_
>>> round_(numpy.max(fluxes.outflow.series))
643.742834

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

Long segments

One advantage of the MC routing method is its relatively low sensitivity to changes in spatial discretisation. However, too long channel segments can cause unrealistic discharge reductions at the start of the rising limb of flood waves. We show this by increasing the segment length to 25 km:

>>> nmbsegments(4)
>>> model.prepare_states(discharge=100.0)
>>> test.inits = ((states.discharge, 100.0),
...               (states.courantnumber, states.courantnumber[0]),
...               (states.reynoldsnumber, states.reynoldsnumber[0]))

>>> conditions = test("musk_mct_long_segments", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

The coarse spatial resolution also leads to a higher peak flow, compensating for the initial discharge reduction:

>>> round_(numpy.max(fluxes.outflow.series))
651.690965
>>> round_(model.check_waterbalance(conditions))
0.0

No baseflow

The accidental discharge reductions discussed in the Long segments example may become more problematic when leading to negative discharge values. musk_mct or, more specifically, method Calc_Discharge_V2 avoids this by simply setting negative discharge values simulated for any segment endpoint to zero. We demonstrate this by removing the base flow from the synthetic flood wave:

>>> nmbsegments(4)
>>> length(40.0)
>>> model.prepare_states(discharge=0.0)
>>> test.inits = ((states.discharge, 0.0),
...               (states.courantnumber, states.courantnumber[0]),
...               (states.reynoldsnumber, states.reynoldsnumber[0]))
>>> nodes.input1.sequences.sim.series = 0.0

>>> conditions = test("musk_mct_no_baseflow", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

Unfortunately, this simple approach for preventing negative discharges introduces an (in this case, relatively small) water balance error:

>>> round_(model.check_waterbalance(conditions))
-0.000149

Little baseflow

Here, we repeat the No baseflow example but assume a small baseflow of 5 m³/s. The trimming of negative discharge values again works and, unfortunately, again introduces some error into the water balance. Besides that, we see a kind of fluctuation progressing from the first to the last segment endpoint with rising amplitude. As a result, the channel outflow first rises to around 8 m³/s and then drops to zero before it finally starts to reflect the expected flood wave shape:

>>> model.prepare_states(discharge=5.0)
>>> test.inits = ((states.discharge, 5.0),
...               (states.courantnumber, states.courantnumber[0]),
...               (states.reynoldsnumber, states.reynoldsnumber[0]))
>>> nodes.input1.sequences.sim.series = 5.0

>>> conditions = test("musk_mct_little_baseflow", get_conditions="2000-01-01")
Click to see the table
Click to see the graph
>>> round_(model.check_waterbalance(conditions))
-0.000044
class hydpy.models.musk_mct.NmbRuns(subvars)[source]

Bases: NmbRuns

The number of (repeated) runs of the RUN_METHODS of the musk_mct model per simulation step [-].

INIT: int | float | bool = 2
name: str = 'nmbruns'

Name of the variable in lowercase letters.

unit: str = '-'

Unit of the variable.

class hydpy.models.musk_mct.Model[source]

Bases: SegmentModel

HydPy-Musk-MCT (Muskingum-Cunge-Todini routing).

The following “inlet update methods” are called in the given sequence at the beginning of each simulation step:

  • Pick_Inflow_V1 Assign the actual value of the inlet sequence to the inflow sequence.

  • Update_Discharge_V1 Assign the inflow to the start point of the first channel segment.

The following “run methods” are called in the given sequence during each simulation step:
The following “outlet update methods” are called in the given sequence at the end of each simulation step:

  • Calc_Outflow_V1 Take the discharge at the last segment endpoint as the channel’s outflow.

  • Pass_Outflow_V1 Pass the channel’s outflow to the outlet sequence.

The following “additional methods” might be called by one or more of the other methods or are meant to be directly called by the user:
The following “submodels” might be called by one or more of the implemented methods or are meant to be directly called by the user:
DOCNAME: DocName = ('Musk-MCT', 'Muskingum-Cunge-Todini routing')
prepare_states(discharge: float, tolerance: float = 1e-06) → None[source]

Determine the Courant and cell Reynolds numbers for the given discharge under the assumption of stationarity.

We prepare a channel consisting of three segments. The first and the last segment are equal to those defined in the Base example:

>>> from hydpy.models.musk_mct import *
>>> simulationstep("30m")
>>> parameterstep()
>>> catchmentarea(25000.0)
>>> nmbsegments(3)
>>> length(6.0)
>>> bottomslope(0.00025)
>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(1)
...     bottomlevels(0.0)
...     bottomwidths(15.0)
...     sideslopes(5.0)
...     stricklercoefficients(1.0/0.035)

We apply prepare_states() for a constant discharge of 100 m³/s:

>>> model.prepare_states(discharge=100.0)
>>> states.courantnumber
courantnumber(0.7207, 0.7207, 0.7207)
>>> states.reynoldsnumber
reynoldsnumber(2.591116, 2.591116, 2.591116)

The smoothness of the discharge simulated for the first timesteps of the Base example indicates that the Courant number and the cell Reynolds number are consistent with the given discharge.

For zero discharge, both the Courant and the cell Reynolds number are zero (even within a triangle profile):

>>> model.prepare_states(discharge=0.0)
>>> states.courantnumber
courantnumber(0.0, 0.0, 0.0)
>>> states.reynoldsnumber
reynoldsnumber(0.0, 0.0, 0.0)
check_waterbalance(initial_conditions: dict[str, dict[str, dict[str, float | ndarray[Any, dtype[float64]]]]]) → float[source]

Determine the water balance error of the previous simulation run in %.

Method check_waterbalance() calculates the balance error as follows:

\(100 \cdot \frac{ Seconds \cdot \sum_{t=t0}^{t1} \Big( Inflow_t - Outflow_t \Big) + 1000 \cdot \sum_{s=1}^{NmbSegments} \Big( Length^s \cdot \big( WettedArea_{t0}^s - WettedArea_{t1}^s \big) \Big) }{ Seconds \cdot \sum_{t=t0}^{t1} \Big( Inflow_t \Big) + 1000 \cdot \sum_{s=1}^{NmbSegments} \Big(Length^s \cdot WettedArea_{t0}^s \Big) }\)

Todini (2007) reports the relative volume error to be less than 0.00 % in all performed experiments. However, we are not sure check_waterbalance() calculates this error in the same way. Please inform us if you encounter larger errors that are not due to preventing negative discharge (see No baseflow).

You can pick the required initial conditions before starting the simulation run via property conditions.

wqmodel

Required submodel that complies with the following interface: CrossSectionModel_V1.

wqmodel_is_mainmodel
wqmodel_typeid
add_wqmodel_v1

Initialise the given submodel that follows the CrossSectionModel_V1 interface and is responsible for calculating discharge and related properties.

>>> from hydpy.models.musk_mct import *
>>> parameterstep()
>>> bottomslope(0.01)
>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(2)
...     bottomlevels(1.0, 3.0)
...     sideslopes(2.0, 4.0)
...     bottomslope
bottomslope(0.01)

>>> model.wqmodel.parameters.control.nmbtrapezes
nmbtrapezes(2)
>>> model.wqmodel.parameters.control.bottomslope
bottomslope(0.01)
REUSABLE_METHODS: ClassVar[tuple[type[ReusableMethod], ...]] = ()
class hydpy.models.musk_mct.Masks[source]

Bases: Masks

Masks applicable to musk_mct.

The following classes are selected:
class hydpy.models.musk_mct.ControlParameters(master: Parameters, cls_fastaccess: type[FastAccessParameter] | None = None, cymodel: CyModelProtocol | None = None)

Bases: SubParameters

Control parameters of model musk_mct.

The following classes are selected:
class hydpy.models.musk_mct.DerivedParameters(master: Parameters, cls_fastaccess: type[FastAccessParameter] | None = None, cymodel: CyModelProtocol | None = None)

Bases: SubParameters

Derived parameters of model musk_mct.

The following classes are selected:
class hydpy.models.musk_mct.FactorSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: FactorSequences

Factor sequences of model musk_mct.

The following classes are selected:
class hydpy.models.musk_mct.FluxSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: FluxSequences

Flux sequences of model musk_mct.

The following classes are selected:
class hydpy.models.musk_mct.InletSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: InletSequences

Inlet sequences of model musk_mct.

The following classes are selected:

  • Q() Runoff [m³/s].

class hydpy.models.musk_mct.OutletSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: OutletSequences

Outlet sequences of model musk_mct.

The following classes are selected:

  • Q() Runoff [m³/s].

class hydpy.models.musk_mct.SolverParameters(master: Parameters, cls_fastaccess: type[FastAccessParameter] | None = None, cymodel: CyModelProtocol | None = None)

Bases: SubParameters

Solver parameters of model musk_mct.

The following classes are selected:

  • ToleranceWaterDepth() Acceptable water depth error for determining the reference water depth [m].

  • ToleranceDischarge() Acceptable discharge error for determining the reference water depth [m³/s].

  • NmbRuns() The number of (repeated) runs of the RUN_METHODS of the musk_mct model per simulation step [-].

class hydpy.models.musk_mct.StateSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: StateSequences

State sequences of model musk_mct.

The following classes are selected: