HydPy-KinW-Implicit-Euler (A-stable routing based on the implicit Euler method)

kinw_impl_euler is a kinematic wave routing approach that does not suffer from stability issues due to high Courant numbers. Its robustness is a result of solving the underlying (non-linear) ordinary differential equations with the implicit Euler method (also known as the Backward Euler method). For the sake of computational efficiency, we implemented this method without any adaptive error control. Hence, it might sometimes be less accurate than routing models as musk_mct (with respect to physical processes) or kinw_williams (with respect to the original differential equations), but it works efficiently when the represented channels are relatively short and the simulation step is relatively long (where the simulation with musk_mct produces unreliable results and the simulation with kinw_williams takes ages). Hence, kinw_impl_euler is best suited for longer (e.g., daily) simulation steps, where the fine details of routing processes are lost anyway.

Note that, in comparison with kinw_williams, kinw_impl_euler is more flexible as it uses submodels for calculating channel geometry properties and discharges (like musk_mct and the members of the more hydrodynamic model family HydPy-SW1D (Shallow Water 1D) do).

Integration tests

Note

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

For good comparability with musk_mct and kinw_williams, we define the following test settings quasi-identical to the ones of the Base example on musk_mct and at least very similar to the ones of the main channel flow on kinw_williams (based on an original test case taken from Todini (2007)):

>>> from hydpy.models.kinw_impl_euler import *
>>> parameterstep()

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

>>> nodes = Nodes("input1", "input2", "output")
>>> stream = Element("stream", inlets=["input1", "input2"], outlets="output")
>>> stream.model = model

>>> length(100.0)
>>> nmbsegments(8)
>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(1)
...     bottomlevels(0.0)
...     bottomwidths(15.0)
...     bottomslope(0.00025)
...     sideslopes(5.0)
...     stricklercoefficients(1.0/0.035)

>>> from hydpy.core.testtools import IntegrationTest
>>> IntegrationTest.plotting_options.axis1 = (
...     fluxes.inflow, fluxes.internalflow, fluxes.outflow
... )
>>> test = IntegrationTest(stream)

>>> import numpy
>>> q_base = 100.0
>>> q_peak = 900.0
>>> t_peak = 24.0
>>> beta = 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)

base example

The simulation results agree well with those of kinw_williams and, in particular, with those of musk_mct. Compared with the latter, the peak flow is underestimated by less than 3 %:

>>> model.prepare_states(100.0)
>>> conditions = test("kinw_impl_euler_base_example", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

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

>>> from hydpy import round_
>>> round_(model.check_waterbalance(conditions))
0.0

overbank flow

Now, we make the channel’s profile more complicated by stacking three trapezes on top of each other, as in the overbank flow example on kinw_williams:

>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(3)
...     bottomlevels(0.0, 6.0, 8.0)
...     bottomwidths(15.0, 200.0, 0.0)
...     bottomslope(0.00025)
...     sideslopes(5.0, 10.0, 100.0)
...     stricklercoefficients(1.0/0.035, 10.0, 10.0)

The simulation results are similar to those of kinw_williams but with more recognisable inflexion points related activation and deactivation of overbank flow, which is due to the “limitation” of kinw_impl_euler not implementing an adaptive numerical integration approach:

>>> model.prepare_states(100.0)
>>> test("kinw_impl_euler_overbank_flow")
Click to see the table
Click to see the graph

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

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

short segments

Short channel segments do neither result in unstable results nor increased computation time:

>>> length(8.0)
>>> model.prepare_states(100.0)
>>> conditions = test("kinw_impl_euler_short_segments", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

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

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

tiny segments

There is no limitation on the channel’s “shortness”; even a length of zero is okay:

>>> length(0.0)
>>> model.prepare_states(100.0)
>>> conditions = test("kinw_impl_euler_tiny_segments", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

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

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

no segments

However, if one wishes to forward any inflow directly as outflow (for testing purposes or so), setting the number of segments to zero is preferable, as it sidesteps unnecessary calculations:

>>> nmbsegments(0)
>>> model.prepare_states(100.0)
>>> conditions = test("kinw_impl_euler_no_segment", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

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

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

numerical accuracy

The numerical parameters WaterVolumeTolerance and WaterDepthTolerance determine the minimum precision of the iterative search for the “right” water depth (“right” from the perspective of the implicit Euler method). This search is the reason for the stability of kinw_impl_euler but causes computational overhead. Hence, setting the values of WaterVolumeTolerance and WaterDepthTolerance is a compromise between numerical accuracy and computational efficiency. We hope we found default values sensible for most applications. However, you might need to refine them in some cases. Here, we set WaterVolumeTolerance to an tremendously large value for illustration, which causes wobbly simulation results but no water balance errors (see method Calc_WaterDepth_V1 for more information):

>>> length(100.0)
>>> nmbsegments(8)
>>> model.prepare_states(100.0)
>>> solver.watervolumetolerance(0.01)
>>> conditions = test("kinw_impl_euler_numerical_accuracy", get_conditions="2000-01-01")
Click to see the table
Click to see the graph
>>> round_(model.check_waterbalance(conditions))
0.0

negative inflow

kinw_impl_euler even sticks to the water balance for negative inflows. It does so by successively removing water from the most upstream to the most downstream channel segment. As soon as the complete channel becomes dry, it forwards the remaining deficit as negative outflow:

>>> solver.watervolumetolerance(1e-10)
>>> model.prepare_states(100.0)
>>> nodes.input1.sequences.sim.series = -200.0
>>> conditions = test("kinw_impl_euler_negative_inflow", get_conditions="2000-01-01")
Click to see the table
Click to see the graph

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

>>> round_(model.check_waterbalance(conditions))
0.0
class hydpy.models.kinw_impl_euler.Model[source]

Bases: SegmentModel

HydPy-KinW-Implicit-Euler (A-stable routing based on the implicit Euler method).

Please note that the following error report appears correct at first glance but is incorrect because the first time when InternalFlow| is required as the flow into a channel’s second segment, it is already calculated as the flow out of its first segment:

>>> from hydpy import prepare_model
>>> from hydpy.core.testtools import check_methodorder
>>> print(check_methodorder(prepare_model("kinw_impl_euler")))
Method Update_WaterVolume_V1 requires the following sequences, which are not among the result sequences of any of its predecessors: InternalFlow
The following “inlet update methods” are called in the given sequence at the beginning of each simulation step:
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:
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:

  • Return_InitialWaterVolume_V1 Calculate and return the initial water volume that agrees with the given final water depth following the implicit Euler method.

  • Return_VolumeError_V1 Calculate and return the difference between the initial water volume that stems from the last simulation step plus the current inflow and the water volume that agrees with the given final water depth following the implicit Euler method.

The following “submodels” might be called by one or more of the implemented methods or are meant to be directly called by the user:

  • PegasusImplicitEuler Pegasus iterator for determining the water level at the end of a simulation step.

DOCNAME: DocName = ('KinW-Implicit-Euler', 'A-stable routing based on the implicit Euler method')
OBSERVER_METHODS: ClassVar[tuple[type[Method], ...]] = ()
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.kinw_impl_euler import *
>>> parameterstep()
>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(2)
...     bottomlevels(1.0, 3.0)
...     bottomslope(0.01)
...     sideslopes(2.0, 4.0)

>>> model.wqmodel.parameters.control.nmbtrapezes
nmbtrapezes(2)
>>> model.wqmodel.parameters.control.bottomslope
bottomslope(0.01)
prepare_states(discharge: float) → None[source]

Determine the initial water volume under the assumption of stationarity.

Due to the assumptions of stationarity and a single profile geometry for the entire channel, the resulting initial volumes are identical for all segments:

>>> from hydpy import pub, Nodes, Element
>>> pub.timegrids = "2000-01-01", "2000-01-05", "30m"
>>> from hydpy.models.kinw_impl_euler import *
>>> parameterstep()
>>> length(25.0)
>>> nmbsegments(2)
>>> with model.add_wqmodel_v1("wq_trapeze_strickler"):
...     nmbtrapezes(1)
...     bottomlevels(0.0)
...     bottomwidths(15.0)
...     bottomslope(0.00025)
...     sideslopes(5.0)
...     stricklercoefficients(1.0/0.035)
>>> model.prepare_states(0.0)
>>> states.watervolume
watervolume(0.0, 0.0)
>>> model.prepare_states(100.0)
>>> states.watervolume
watervolume(1.560986, 1.560986)
>>> model.prepare_states(10000.0)
>>> states.watervolume
watervolume(48.454829, 48.454829)

Calling method prepare_states() can overwrite previously calculated values of some sequences but at least resets the channel length (which it sets to zero for one calculation step) to its original value:

>>> length
length(25.0)
check_waterbalance(initial_conditions: dict[str, dict[str, dict[str, float | ndarray[tuple[Any, ...], dtype[float64]]]]]) → float[source]

Determine the water balance error of the previous simulation run in million m³.

Method check_waterbalance() calculates the balance error as follows:

\[\begin{split}Error = \Sigma InOut - \Delta Vol \\ \\ \Sigma InOut = s \cdot 10^{-6} \cdot \sum_{t=t0}^{t1} I_t - O_t \\ \Delta Vol = \sum_{i=1}^{n} V_{t1}^i - V_{t0}^i \\ \\ s = Seconds \\ I = Inflow \\ O = Outflow \\ n = NmbSegments \\ V = WaterVolume\end{split}\]

The returned error should always be in scale with numerical precision so that it does not affect the simulation results in any relevant manner.

Pick the required initial conditions before starting the simulation run via property conditions. See the integration tests of the application model kinw_impl_euler for some examples.

REUSABLE_METHODS: ClassVar[tuple[type[ReusableMethod], ...]] = ()
class hydpy.models.kinw_impl_euler.ControlParameters(master: Parameters, cls_fastaccess: type[FastAccessParameter] | None = None, cymodel: CyModelProtocol | None = None)

Bases: SubParameters

Control parameters of model kinw_impl_euler.

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

Bases: SubParameters

Derived parameters of model kinw_impl_euler.

The following classes are selected:

  • Seconds() The length of the actual simulation step size in seconds [s].

  • NmbDiscontinuities() Number of points of discontinuity in the rating curve [-].

  • FinalDepth2InitialVolume() A pair of the final water depth and the corresponding initial water volume according to the implicit Euler method for each point of discontinuity in the rating curve [m and million m³].

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

Bases: FactorSequences

Factor sequences of model kinw_impl_euler.

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

Bases: FluxSequences

Flux sequences of model kinw_impl_euler.

The following classes are selected:

  • Inflow() Flow into the first channel segment [m³/s].

  • InternalFlow() Flow between the channel segments [m³/s].

  • Outflow() Flow out of the last channel segment [m³/s].

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

Bases: InletSequences

Inlet sequences of model kinw_impl_euler.

The following classes are selected:

  • Q() Abfluss (runoff) [m³/s].

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

Bases: OutletSequences

Outlet sequences of model kinw_impl_euler.

The following classes are selected:

  • Q() Abfluss (runoff) [m³/s].

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

Bases: SubParameters

Solver parameters of model kinw_impl_euler.

The following classes are selected:

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

  • WaterVolumeTolerance() Targeted accuracy in terms of the relative water volume for the Pegasus search of the final water depth [-].

  • WaterDepthTolerance() Targeted accuracy in terms of the absolute water depth for the Pegasus search of the final water depth [m].

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

Bases: StateSequences

State sequences of model kinw_impl_euler.

The following classes are selected: