Bidirectional coupling of HydPy-W(ALRUS) and HydPy-Dam

[tyralla]

We recently made some progress in coupling HydPy-W(ALRUS) and SOBEK 3 in a bidirectional manner, allowing us to simulate both effluent in influent flow conditions. The feedback and the possible reverse flow from surface water (SOBEK) to groundwater (HydPy-W) generally make simulations more realistic and help circumvent the "over drainage problem" we and others encountered in previous studies when applying HydPy-W or the original WALRUS model in a standard unidirectional coupling.

This approach (based on BMI) is excellent for detailed studies but requires detailed information on channel geometries for building a SOBEK 3 model and is computationally time-consuming. Hence, we next plan to support a more "lumped" approach via bidirectional coupling of HydPy-W with HydPy-Dam-Pump, HydPy-Dam-Pump-Sluice, and HydPy-Dam-Pump-Sluice.

HydPy-W and SOBEK 3 are coupled via individual Python scripts that equate each trench of the SOBEK model with the surface water storage of a separate HydPy-W instance. For connecting HydPy-W with HydPy-Dam, we strive for a computationally more efficient and general solution that does not require setting up problem-specific scripts, maybe something similar to using a HydPy-Exch-Weir-HBV96 instance for allowing bidirectional flows between two lakes modelled by separate HydPy-Dam-L-Lake instances.

An HydPy-Exch-Weir-HBV96-like approach would introduce some additional configuration burden on the users and introduce some inaccuracies into the results as feedback could only be accounted for with a time delay of on simulation step. Also, HydPy-Dam's water volume to water level relationship and HydPy-W's surface water storage geometry would need to be synchronised, or additional calculation inaccuracies could appear.

Another idea is to allow HydPy-W to use HydPy-Dam as a submodel. This might be the approach with the smallest amount of required stopgaps. However, it would be a very new kind of functionality and the first time we couple two models defined via ordinary differential equations via the submodel approach. Hence, it is hard to foresee the required amount of work.

Two alternative solutions that come into my mind but do not seem worth further contemplation are to implement specific models that combine the required HydPy-W and HydPy-Dam functionalities (very ugly) or somehow insert user-made functions that do the coupling (we started discussing such a "hook mechanism" in issue #94 but still made no progress there; so, even if we consider this approach as feasible, it would be a long way to go).

Other ideas are welcome, of course. If no one objects, I would first experiment with the submodel idea to get a first impression of the expected workload.

[sivogel]

Wouldn't it be nice to exchange between dam and walrus via remote nodes as it is already done in the dam models?

[tyralla]

Yes, but one instance must control the exchange logic (the synchronisation between a SOBEK trench and a HydPy-W surface water storage).

[tyralla]

One disadvantage of the submodel approach is that it is likely not applicable to LIAS-like models. It would be great to support using LIAS instead of SOBEK without the additional external scripting efforts required for the SOBEK coupling. And one extensive network model instance cannot be the submodel of multiple HydPy-W instances. If any, the HydPy-W instances would have to be the submodels of individual channel model instances, but that seems like a too-hard break with the standard HydPy project structure. Additionally, HydPy-W and LIAS rely on different numerical integration strategies, so coupling them very tightly is impossible.

So, this consideration speaks more for an HydPy-Exch-Weir-HBV96-like approach (where the communication between models takes place via "remote nodes").

[tyralla]

We had some in-depth discussions and decided to implement low-weighted submodels that introduce the necessary flexibility to use HydPy-W either in the standard way (with weir-like outflow) or the proposed way (with bidirectional coupling).

Only the "bidirectional coupling" submodel has a receiver node for making the water level of a downstream HydPy-Dam or HydPy-LIAS model instance available upstream to affect the HydPy-W calculations.

We prefer to send information on the water level instead of the water volume upstream, as this seems more robust regarding different water volume-to-water level assumptions of the involved models. This approach should never introduce water balance errors despite exchanging levels instead of volumes because WALRUS' surface water storage degrades to an "auxiliary storage" just used to help calculate the actual exchange (bidirectional flow). When using a "bidirectional coupling" submodel, WALRUS' outflow is no longer defined by RH but by all other flows entering or exiting the surface water storage. We will add another sequence for clarification. First suggestion ($NIS$ for "net input into the surface water storage"; can we find a better wording?):

$$ \frac{dHS}{dt} = PS - ETS + FXS + \frac{ALR \cdot \left(AGR \cdot FGS + FQS \right) - RH}{ASR} = NIS- \frac{RH}{ASR} $$

Another outcome of our discussions is to add the control parameters GL (ground level) and BL (channel bottom level) and turn CD (channel depth) into a derived parameter. This change aims to make the parametrisation of HydPy-W more descriptive when coupling it to other models in a manner where differences in elevation are relevant.

[tyralla]

HydPy-W is defined via ordinary differential equations solved via an adaptive Runge-Kutta method. It seems favourable to make the submodels as simple as possible so that they do not need to know anything about this and could be coupled with other kinds of models, too.

I suggest to introduce "rating curve" interfaces. In the simplest case, such an interface could describe submodels that serve as plain functions to convert water levels into discharges:

class RatingCurveModel_V1(modeltools.SubmodelInterface):

    typeid: ClassVar[Literal[1]] = 1

    @modeltools.abstractmodelmethod
    def calc_discharge(self, waterlevel: float) -> float:
        """Calculate the discharge based on the water level given in m and return it in 
        m³/s."""

Currently, we would only need to implement one model following this interface, maybe something like "HydPy-RC-WALRUS" (WALRUS-like weir overflow based on the Manning equation). However, based on this idea, we could later improve HydPy-LIAS and HydPy-Musk-MCT as well, which currently only work with trapezoidal profiles, by introducing models like "HydPy-RC-Manning-Single-Trapeze" or "HydPy-RC-Manning-Triple-Trapeze".

Note that HydPy-RC-WALRUS would not need to own sequences. Instead, its main models would use it to calculate the values for their sequences (in the case of HydPy-W, RH). However, submodels that comply with RatingCurveModel_V1 would usually be more than plain functions, as they would need configuration (e.g. with the weir's crest height).

The second interface would target the bidirectional exchange. A first draft that keeps the interface extremely small:

class RatingCurveModel_V2(modeltools.SubmodelInterface):

    typeid: ClassVar[Literal[2]] = 2

    @modeltools.abstractmodelmethod
    def get_initialwaterlevel(self) -> float:
        """Determine the initial water level and return it in m."""

The possible workflow from HydPy-W's perspective:

1. Query the water level from the downstream HydPy-Dam model (via method get_initialwaterlevel).
2. Calculate all fluxes and states as usual, but do not use the submodel to calculate RH. Instead, always set it to zero.
3. At the end of a simulation step, use the initial water level (returned by the submodel) and the final water level to calculate the net outflow.

Technically, all of this should work (at least, I hope so), and I like the interfaces' simplicity. However, some points might need improvement:

• The term "rating curve" does not really fit, especially not for the second interface.
• One would still have to define a value for the state sequence HS in a condition file, but when using the second interface, this value would be overwritten by the downstream model's water level. during runtime So, it seems reasonable to turn the "initial condition" HS into a property of the submodel but to leave it in HydPy-W for the purpose of numerical integration. I am unsure if a simple and transparent way exists to implement this.
• RH would always be zero when using a submodel that complies with the second interface. However, I think this is unavoidable when working with such "function-like" submodels. Also, it is less annoying than the case of HS because users do not need to specify values for RH and should seldom work directly with it.

[tyralla]

• The term "rating curve" does not really fit, especially not for the second interface.

It seems more correct to name it WaterLevelModel_V1 and add it to stateinterfaces.py:

class WaterLevelModel_V1(modeltools.SubmodelInterface):

    typeid: ClassVar[Literal[1]] = 1

    @modeltools.abstractmodelmethod
    def get_initialwaterlevel(self) -> float:
        """Determine the initial water level and return it in m."""

However, this brings the new situation that a main model (HydPy-W) can either use a water level or a rating curve submodel, and we need to handle cases where users define none or both.

The water level submodel does not necessarily require further configuration, so we could let HydPy-W use it by default.

Also, we should probably consider implementing separate attributes for each (model.waterlevelmodel and model.ratingcurvemodel, for example). This causes more synchronisation efforts but clarifies the different purposes (compared to model.waterlevel_or_ratingcurve_model) and keeps both models distinguishable (which would not be possible based on typeid, which is one in both cases).

[tyralla]

• The term "rating curve" does not really fit, especially not for the second interface.

Regarding the first interface, maybe just DischargeModel_V1? Or is this name too general?

[tyralla]

Suggested differences between the (new) submodel and the old HydPy-W components.

As done for other submodels, I suggest to introduce more descriptive names:

• CD -> ChannelDepth
• HSMin -> CrestHeight
• SH -> CrestHeightTolerance
• CS -> BankfullDischarge
• XS -> DischargeExponent

Regarding the units, first, I suggest switching from mm to m. Second, I suggest m³/s for the bankfull discharge and the returned actual discharge instead of mm/T. This forces HydPy-W to perform internal unit conversions, but it is what routing models expect. On the downside, BankfullDischarge would depend on the subcatchment size.

[tyralla]

2. Calculate all fluxes and states as usual, but do not use the submodel to calculate RH. Instead, always set it to zero.

We decided on another option. If there is no submodel, method Calc_RH_V1 now equates the runoff height with all other flows in and out of the surface water storage:

$$ RH = ASR \cdot (PS - ES + FXS) + ALR \cdot (AGR \cdot FGS + FQS) $$

• RH would always be zero when using a submodel that complies with the second interface. However, I think this is unavoidable when working with such "function-like" submodels. Also, it is less annoying than the case of HS because users do not need to specify values for RH and should seldom work directly with it.

With the above change, we avoid this problem. Also, the coupling with a downstream HydPy-Dam model should become easier to implement and comprehend. While the first suggestion implied a possibly infinite increase in surface water storage, the new suggestion means that the surface water storage's content only changes if a second submodel uses downstream water levels to update it (see below).

[tyralla]

1. Query the water level from the downstream HydPy-Dam model (via method get_initialwaterlevel). 3. At the end of a simulation step, use the initial water level (returned by the submodel) and the final water level to calculate the net outflow.

I found a better solution. Now, HydPy-W uses the WaterLevelModel_V1 submodel to query the downstream water level at the end of the simulation step via a receiver node.

• One would still have to define a value for the state sequence HS in a condition file, but when using the second interface, this value would be overwritten by the downstream model's water level during runtime. So, it seems reasonable to turn the "initial condition" HS into a property of the submodel but to leave it in HydPy-W for the purpose of numerical integration. A simple and transparent way exists to implement this.

Due to this change, HS is always required, as it receives data at the end of a simulation step (or complete simulation) and provides this data at the beginning of the next simulation step (or complete simulation).

So, this approach is more consistent with the current model implementation standards (and it also allows inserting water level measurements as observed time series via remote nodes).

[tyralla]

Due to the abovementioned change, the naming get_initialwaterlevel no longer seems right. So, the interface is now:

class WaterLevelModel_V1(modeltools.SubmodelInterface):
    """Pure getter interface for querying the current water level."""

    typeid: ClassVar[Literal[1]] = 1

    @modeltools.abstractmodelmethod
    def get_waterlevel(self) -> float:
        """Determine the water level and return it in m."""

[tyralla]

Everything has been implemented and pushed to the master branch, so I will close this issue. However, we can re-open it if experiences in a first real-world application (January 2024) suggest the implemented functionalities require improvement.