Implement the "local inertial approximation of the shallow water equations"

[tyralla]

So far, HydPy provides only "kinematic" routing approaches, which approximate the diffusion of flood waves by fine-tuning numerical diffusion but cannot account for backwater effects. Currently, this limits our capabilities in modelling the anthropogenically controlled flow regimes in northern Germany (Schleswig-Holstein), where terrains are very flat and sometimes even below sea level. Therefore, we decided to close this gap by introducing a new model family responsible for calculating 1-dimensional water movement in channels in a "hydrodynamic" manner by approximating the shallow water equations.

The first decision is already settled. Our first application model (or selectable mathematical approximation) will be the "local inertial approximation of the shallow water equations" introduced by Bates et al. (2010) and "stabilised" by Almeida et al. (2012). The method's original purpose is modelling inundations in two dimensions, but we see no reason not to use it for 1-dimensional routing. The authors provide 1-dimensional test cases (Almeida and Bates, 2013), wflow uses the original method after Bates et al. (2010) for river routing, and comparisons of a prototype with Sobek under simplified conditions were successful.

We are not aware of a catchy name for the method. We agreed spontaneously on Lias-1D (Local Inertial Approximation of the Shallow water equations). ("1D", in case we add a 2-dimensional hydrodynamic method to HydPy later).

In the following, we speak of Lias-1D for simplicity when referring to all potential members (or methods) of the new model family.

It is foreseeable that the new model family will heavily use the new submodel concept (#90), which still needs to be fully developed. However, for the same reasons discussed in #95, we prefer to add another example before finalising the implementation of all submodel functionalities. In cases like those covered by #95, each main model will likely always require at most one instance of a submodel of a specific type (e.g. to calculate potential evapotranspiration after Turc-Wendling). In contrast, for LIAS-1D, a main model (representing a channel network) might reference a higher number of submodels of the same type, for example, to model the effects of multiple weirs within the network.

LIAS-1D differs structurally from the already implemented routing methods because it needs to simulate the whole network of connected channels in parallel, so that water can flow in any direction. Hence, a first attempt would be to define one model instance of LIAS-1D per basin model. However, that could result in inconveniences. Usually, we differentiate a basin model into multiple selections, each reflecting the subareas and substreams between two gauges or above the most upstream gauge. This makes the model configuration files more manageable. And it allows modelling only specific subareas of a basin, which is very important for saving computation time.

Maybe we should allow defining multiple parts of the channel network in different control files. Each control file would correspond to a single Element instance and a single Model instance. Before starting a simulation, the separate LIAS-1D instances would have to connect. Unconnectable instances would instead require alternative boundary conditions. Would it make sense to define default boundary conditions in the control files that are only applied when no connections to other LIAS-1D instances are possible?

When allowing separating channel networks on different model instances, we need to find a way to trigger their parallel execution. A very first idea: We could introduce two new base classes, like ParallelController and ParallelModel. Method HydPy.update_devices searches for ParallelModel instances and hands them to a suitable ParallelController instance. When determining the execution order, it does not take the ParallelModel instances but only the ParallelController instances into account. During simulation, ParallelController would only have to control the execution of the ParallelModel instances. So, they also could be considered submodels.

Following this idea, it seems most likely that ParallelController and ParallelModel need to share much logic. So maybe it is preferable to derive LIAS-1D from ParallelModel and let it return an appropriate controller upon request. For the methodology at hand, the controller, for example, would have to determine the suitable numerical step size required to stabilise all ParallelModel instances.

Maybe, the same type of controller could work for different types of ParallelModel subtypes so that selecting different approximations for different parts of a channel network would work?

Many additional aspects require discussion, for example:

• What types of boundary conditions do we want to support?
• Which ways to define channel geometries should we support?
• Must the supplied geometries be ready for simulation, or do we implement some preprocessing (to, for example, allow changing the length of the individual river sections)?
• How can we structure such complex data in the usual control file format (or do we need another data format)?
• In the long run: Should we create a small GUI that helps structure and maintain the channel data (or is such a GUI already available that we could use)?

We must remember that HydPy focuses on catchment hydrology, not river hydraulics. So it's preferable to keep things simple from the hydrologist's perspective.

[tyralla]

Control file structure/ submodels

After some fruitful discussions, our picture of the overall structure becomes clearer. We might heavily rely on the new submodel concept and need to extend it a little. This is the very first draft:

from hydpy.models.sw1d_channel import *
from hydpy.models import sw1d_inflow, sw1d_lias, sw1d_storage, sw1d_trapeze, sw1d_weir

parameterstep()

nmbsegments(3)

with model.add_storagemodel_v1(sw1d_storage, cell=(0, 1, 2)):
    length(2.0)
    bottomlevel(100.0)
    with model.add_profilemodel_v1(sw1d_trapeze):
        bottomwidth(3.0)
        sideslope(2.0)

with model.add_routingmodel_v1(sw1d_inflow, edge=0):
   pass

with model.add_routingmodel_v1(sw1d_lias, edge=1):
    timestepfactor(0.7)
    diffusionfactor(0.2)
    stricklercoefficient(20.0)
    with model.add_profilemodel_v1(sw1d_trapeze):
        bottomwidth(3.0)
        sideslope(2.0)

with model.add_routingmodel_v1(sw1d_lias, edge=2):
    timestepfactor(0.7)
    diffusionfactor(0.2)
    stricklercoefficient(20.0)
    with model.add_profilemodel_v1(sw1d_trapeze):
        bottomwidth(3.0)
        sideslope(2.0)

with model.add_routingmodel_v1(sw1d_weir, edge=3):
    crestlevel(101.0)
    crestwidth(3.0)

The main ideas:

• The primary purpose of the sw1d_channel model is to handle the different submodels which do the actual calculations.
• Following the above terminology, sw1d_channel is the ParallelModel. sw1d_network is the corresponding ControllerModel. (These names suggest that all future "Shallow Water 1D" models rely on staggered grids. Alternatively, we could name both models sw1d_channel_staggeredgrid and sw1d_network_staggeredgrid. However, we will probably not soon add models based on other spatial discretisations. So I would prefer the shorter, less technically sounding names.)
• sw1d_channel allows specifying the number of channel segments and adding an individual storage model for each cell and an individual routing model for each segment edge. (Others use the term "interface" instead of "edge", but I tried to come up with another word to avoid confusing "cell interfaces" with "submodel interfaces".)
• Storage models are responsible for updating the actual water volume and calculating the water level in the segments' centre.
• Routing models are responsible for calculating the discharge.
• sw1d_channel could either control the submodel's data exchange or connect the submodels with each other to let them communicate directly. The first option seems like some overhead (sw1d_lias asks sw1d_channel to ask the neighbour storage models for the current water level) and requires additional methods to let sw1d_channel communicate with its neighbour sw1d_channel instances. The second option requires many model connections (sw1d_lias would be connected to one sw1d_channel, two storage, and two routing model instances) that may differ between model types (sw1d_weir probably requires no connections to other routing models). Currently, I slightly favour the second option.
• Factoring out profile models (sw1d_trapeze) now is a spontaneous idea. It would be nice to have in a current project as it would allow us to switch to double trapeze profiles where the inundation of river banks plays a significant role. However, models like musk_mct would also profit from such flexibility. So we should think about general compatibility and possibly creating a separate profile model family.
• Two sw1d_channel instances are connected by one routing model. We will check if exactly one control file defines such a model at the relevant position and raise an error for missing and double definitions (without checking if the doubles are identical).

[tyralla]

How to apply the q-centred scheme at junctions?

We decided to ignore different angles at channel junctions. That means if A and B flow into C, all are assumed to be aligned in the same direction:

Hence, A and B are principally handled identically regarding the q-centred scheme that introduces numerical diffusion for stability reasons. However, we need an additional assumption here because the discharge at the end of channel C's first segment is much higher than the discharge in the last segments of channels A and B.

Let us take the sw1d_lias submodel between A and B as an example. Taking the original equation without modification, the q-centred term would be:

$$ (1 - \theta) \cdot Q_{A-C,i} + \theta / 2 \cdot ( Q_{A, i-1} + Q_{C, i+1}) $$

We suggest to adjust $Q_{C, i+1}$ based on the relation of $Q_{A-C,i}$ and $Q_{B-C,i}\ $:

$$ (1 - \theta) \cdot Q_{A-C,i} + \theta / 2 \cdot \left( Q_{A, i-1} + \frac{Q_{A-C,i}}{Q_{A-C,i} + Q_{B-C,i}} \cdot Q_{C, i+1} \right) $$

[tyralla]

Preventing too small (and possibly negative) water depths

We will introduce a minimum water depth for the individual cells for stability. (I think we agreed on 1 cm. Or was it 1 mm?).

A local solution would be to prevent any outflow of a cell that would cause the water depth to fall below this threshold without considering the possible inflow from the other edge (or lateral inflow). But I am afraid this could result in undesired on-off-switching situations.

A global solution would be to sum the deficits of all too-dry cells, set them to the minimum depth, and reduce the water volume of all remaining cells by a factor which ensures keeping the water balance. This is stable but would result in infinitely fast water redistribution in large areas. And, theoretically, to cycles where water flows from nearly dry channels upstream to wet channels downstream, which is then "beamed" back.

If no one objects, I try the first solution. We could exchange the discontinuous maximum function with a smoother alternative later if on-off-switching introduces instability.

[tyralla]

Infinite roughness does not completely suppress discharge

This note might help us avoid running into the same problem later.

I tried to define a test with nearly zero Strickler coefficients and expected this would suppress discharge (so that water level gradients would never change). To my surprise, this did not work because the channel roughness (the Strickler coefficient) does not matter at all when the previous discharge value is zero:

$$ Q_{i-1/2}^{n+1} = \frac{\left[ (1-\theta) \cdot Q_{i-1/2}^n + \frac{\theta}{2} \cdot \left( Q_{i-3/2}^n + Q_{i+1/2}^n\right) \right] + g \cdot A_f \cdot \Delta t^n \cdot (y_{i-1}^n - y_i^n) / \Delta x} {1 + g \cdot \Delta t \cdot k_{st}^{-2} \cdot \left| Q_{i-1/2}^n \right| \cdot {P_f^n}^{\frac{4}{3}} \cdot {A_f^n}^{-\frac{7}{3}}} $$

So, this unexpected behaviour is due to the underlying method rather than our current implementation. I do not know if there is any "real-life" setting where it would result in relevant problems.

[tyralla]

Uneven spatial discretisation

The original approach is to be applied on equidistant rectangular grids. For 1-dimensional routing problems, it is seldom possible to work with even discretisations, as "special points" like confluences require special care. This might pose some problems.

Consider a 20 km long channel with a constant inflow of 1 m³/s at the upper and a weir at the lower end. The weir's crest height, crest width, and flow coefficient are 7 m, 5 m, and 0.58, respectively. The channel is divided into eight segments with identical rectangular profiles of 5 m width. Their bottom level is 5 m, their Manning coefficient is 0.03, and the initial water level is 7 m.

I ran an experiment where the first, third, fifth, and seventh segments were 2 km and all others 3 km long. With a "normal" time step safety factor of 0.7 and a "normal" diffusion factor of 0.2., I ran into stability issues (the yellow line exceeds 1 m³/s):

So I reduced the time step safety factor to the extremely small value of 0.1. This stabilised the calculation but introduced an unexpected "up and down" pattern:

Increasing the diffusion coefficient magnifies the problem:

If I set all segments to the same length of 2.5 km, the problem vanishes:

On first thought, the "up and down" problem seems to be due to weighting the upstream and downstream discharges with the same factor without considering potential differences between the distances to the upstream and downstream stations.

[tyralla]

Another diffusion problem

I encountered another unexpected behaviour. The following test case includes neither inflow nor outflow. Instead, I set an initial depth of 3 m for the first four and 1 m for the last four segments to enforce water movement. With settings otherwise identical to those explained above, I get the following eight water level trajectories:

There are some "edges" at the start of the simulation period. Again, I tried to smooth them by decreasing the time safety factor, which seemed to slow the water level equalisation. Repeating this experiment with the tiny safety factor of 0.01 gives the following results:

Again, this seems to be due to the weighting approach introduced by Almeida et al. (2012). After setting the diffusion factor to zero, the results look sufficiently similar for the safety factors 0.7 (first graph) and 0.01 (second graph):

I checked the code thoroughly but could not find a bug. The only "solution" I could come up with is to reduce the diffusion factor for small internal numerical step sizes ($TimeStep$) compared to the external simulation step size ($Seconds$) :

$$ DiffusionFactor_{modified} = DiffusionFactor_{original} \cdot \frac{TimeStep}{Seconds} $$

Then, the results for a time safety factor of 0.01 and the original diffusion coefficient of 0.2 look like this:

Of course, this solution would be neither theoretically grounded nor based on much empirical evidence. At least, it would reduce the problem discussed above...

...to the following:

Ignoring the diffusion term seems like the preferable option for applications to me for now.

[tyralla]

LIAS is now in the master branch. Its documentation is online: HydPy-SW1D.

Unfortunately, the documentation on function combine_channels of the sw1d_network model is missing for technical reasons. But we should discuss this in a separate issue on further developing the automatic documentation generation.

For further information on how to include LIAS into an HydPy project, search for collective, unite_collectives.

Note that we had to modify our idea on applying the q-centred scheme at junctions, as it could result in infinitely high discharge estimates for situations where water flows in different directions. See, for example, Get_PartialDischargeDownstream_V1.