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)
...     calibrationfactors(1.0)

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