Ribasim Delwaq coupling

In order to generate the Delwaq input files, we need a completed Ribasim simulation (typically one with a results folder) that ideally also includes some substances and initial concentrations. Let’s take the basic test model for example, which already has set some initial concentrations.

All testmodels can be downloaded from here.

from pathlib import Path

toml_path = Path("../../generated_testmodels/basic/ribasim.toml")

assert toml_path.is_file()

This Ribasim model already has substance concentrations for Cl and Tracer in the input tables, and we will use these to generate the Delwaq input files.

from ribasim import Model

model = Model.read(toml_path)

display(model.basin.concentration_state)  # basin initial state
display(model.basin.concentration)  # basin boundaries
display(model.flow_boundary.concentration)  # flow boundaries
display(model.level_boundary.concentration)  # level boundaries
model.plot();  # for later comparison
Basin / concentration_state
node_id substance concentration
fid
0 1 Cl 0.0
1 3 Cl 0.0
2 6 Cl 0.0
3 9 Cl 0.0
Basin / concentration
node_id time substance drainage precipitation
fid
0 1 2020-01-01 00:00:00 Cl 0.0 0.0
1 1 2020-01-02 00:00:00 Cl 1.0 1.0
2 1 2020-01-01 00:00:00 Tracer 1.0 1.0
3 3 2020-01-01 00:00:00 Cl 0.0 0.0
4 3 2020-01-02 00:00:00 Cl 1.0 1.0
5 3 2020-01-01 00:00:00 Tracer 1.0 1.0
6 6 2020-01-01 00:00:00 Cl 0.0 0.0
7 6 2020-01-02 00:00:00 Cl 1.0 1.0
8 6 2020-01-01 00:00:00 Tracer 1.0 1.0
9 9 2020-01-01 00:00:00 Cl 0.0 0.0
10 9 2020-01-02 00:00:00 Cl 1.0 1.0
11 9 2020-01-01 00:00:00 Tracer 1.0 1.0
FlowBoundary / concentration
node_id time substance concentration
fid
0 15 2020-01-01 00:00:00 Cl 0.0
1 15 2020-01-01 00:00:00 Tracer 1.0
2 16 2020-01-01 00:00:00 Cl 0.0
3 16 2020-01-01 00:00:00 Tracer 1.0
LevelBoundary / concentration
node_id time substance concentration
fid
0 11 2020-01-01 00:00:00 Cl 34.0
1 17 2020-01-01 00:00:00 Cl 34.0

model.basin.profile
Basin / profile
node_id area level storage
fid
0 1 0.01 0.0 <NA>
1 1 1000.0 1.0 <NA>
2 3 0.01 0.0 <NA>
3 3 1000.0 1.0 <NA>
4 6 0.01 0.0 <NA>
5 6 1000.0 1.0 <NA>
6 9 0.01 0.0 <NA>
7 9 1000.0 1.0 <NA>

Let’s add another tracer to the model, to setup a fraction calculation.

from ribasim.delwaq import add_tracer

add_tracer(model, 11, "Foo")
add_tracer(model, 15, "Bar")
display(model.flow_boundary.concentration)  # flow boundaries
display(model.level_boundary.concentration)  # flow boundaries

model.write(toml_path)
FlowBoundary / concentration
node_id time substance concentration
fid
0 15 2020-01-01 00:00:00 Cl 0.0
1 15 2020-01-01 00:00:00 Tracer 1.0
2 16 2020-01-01 00:00:00 Cl 0.0
3 16 2020-01-01 00:00:00 Tracer 1.0
4 15 2020-01-01 00:00:00 Bar 1.0
LevelBoundary / concentration
node_id time substance concentration
fid
0 11 2020-01-01 00:00:00 Cl 34.0
1 17 2020-01-01 00:00:00 Cl 34.0
2 11 2020-01-01 00:00:00 Foo 1.0
PosixPath('../../generated_testmodels/basic/ribasim.toml')

Given the path to a completed Ribasim simulation, we can call ribasim.delwaq.generate for generating the required input files for Delwaq from scratch.

from ribasim.delwaq import generate

output_path = Path("../../generated_testmodels/basic/delwaq")

graph, substances = generate(toml_path, output_path)
/home/runner/work/Ribasim/Ribasim/python/ribasim/ribasim/delwaq/generate.py:414: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

This call produces a handful of files in the user defined folder. Let’s take a look at them:

list(output_path.iterdir())
[PosixPath('../../generated_testmodels/basic/delwaq/ribasim.nc'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.len'),
 PosixPath('../../generated_testmodels/basic/delwaq/network.csv'),
 PosixPath('../../generated_testmodels/basic/delwaq/delwaq.inp'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.atr'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.vol'),
 PosixPath('../../generated_testmodels/basic/delwaq/B5_bounddata.inc'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim_bndlist.inc'),
 PosixPath('../../generated_testmodels/basic/delwaq/bndlist.csv'),
 PosixPath('../../generated_testmodels/basic/delwaq/dimr_config.xml'),
 PosixPath('../../generated_testmodels/basic/delwaq/flows.csv'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.are'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.poi'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.vel'),
 PosixPath('../../generated_testmodels/basic/delwaq/ribasim.flo'),
 PosixPath('../../generated_testmodels/basic/delwaq/volumes.csv')]

These files form a complete Delwaq simulation, and can be run by either pointing DIMR to the dimr_config.xml file or pointing Delwaq to the delwaq.inp file.

Note that the call to generate produces two output variables; graph and substances that are required for parsing the results of the Delwaq model later on. Nonetheless, we can also inspect them here, and inspect the created Delwaq network.

substances  # list of substances, as will be present in the Delwaq netcdf output
{'Bar',
 'Cl',
 'Continuity',
 'Drainage',
 'FlowBoundary',
 'Foo',
 'Initial',
 'LevelBoundary',
 'Precipitation',
 'Terminal',
 'Tracer',
 'UserDemand'}

As you can see, the complete substances list is a combination of user input (Cl and Tracer in the input tables), a Continuity tracer, and tracers for all nodetypes in the Ribasim model. The latter tracers allow for deeper inspection of the Ribasim model, such as debugging the mass balance by plotting fraction graphs. Let’s inspect the graph next, which is the Delwaq network that was created from the Ribasim model:

import matplotlib.pyplot as plt
import networkx as nx

# Let's draw the graph
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
nx.draw(
    graph,
    pos={k: v["pos"] for k, v in graph.nodes(data=True)},
    with_labels=True,
    labels={k: k for k, v in graph.nodes(data=True)},
    ax=ax[0],
)
ax[0].set_title("Delwaq node IDs")
nx.draw(
    graph,
    pos={k: v["pos"] for k, v in graph.nodes(data=True)},
    with_labels=True,
    labels={k: v["id"] for k, v in graph.nodes(data=True)},
    ax=ax[1],
)
ax[1].set_title("Ribasim node IDs")
fig.suptitle("Delwaq network");

Here we plotted the Delwaq network twice, with the node IDs as used by Delwaq on the left hand side, and the corresponding Ribasim node IDs on the right hand side. As you can see, the Delwaq network is very similar to the Ribasim network, with some notable changes:

1 Parsing the results

With Delwaq having run, we can now parse the results using ribasim.delwaq.parse. This function requires the graph and substances variables that were output by ribasim.delwaq.generate, as well as the path to the results folder of the Delwaq simulation.

from ribasim.delwaq import parse

nmodel = parse(toml_path, graph, substances, output_folder=output_path)

The parsed model is identical to the Ribasim model, with the exception of the added concentration_external table that contains all tracer results from Delwaq.

display(nmodel.basin.concentration_external)
print(substances)
t = nmodel.basin.concentration_external.df
t[t.time == t.time.unique()[2]]
Basin / concentration_external
time node_id concentration substance
fid
0 2020-01-01 00:00:00 1 1.0 Initial
1464 2020-01-01 00:00:00 1 0.0 Foo
2928 2020-01-01 00:00:00 1 0.0 Drainage
4392 2020-01-01 00:00:00 1 0.0 UserDemand
5856 2020-01-01 00:00:00 1 0.0 FlowBoundary
... ... ... ... ...
11711 2020-12-31 00:00:00 9 0.17941 Tracer
13175 2020-12-31 00:00:00 9 0.0 Bar
14639 2020-12-31 00:00:00 9 0.019136 Cl
16103 2020-12-31 00:00:00 9 0.000563 LevelBoundary
17567 2020-12-31 00:00:00 9 0.179993 Continuity

17568 rows × 4 columns

{'Initial', 'Foo', 'Drainage', 'UserDemand', 'FlowBoundary', 'Precipitation', 'Terminal', 'Tracer', 'Bar', 'Cl', 'LevelBoundary', 'Continuity'}
time node_id concentration substance
fid
8 2020-01-03 00:00:00 1 0.035288 Initial
1472 2020-01-03 00:00:00 1 0.0 Foo
2936 2020-01-03 00:00:00 1 0.0 Drainage
4400 2020-01-03 00:00:00 1 0.0 UserDemand
5864 2020-01-03 00:00:00 1 0.984592 FlowBoundary
7328 2020-01-03 00:00:00 1 0.143128 Precipitation
8792 2020-01-03 00:00:00 1 0.0 Terminal
10256 2020-01-03 00:00:00 1 1.127719 Tracer
11720 2020-01-03 00:00:00 1 0.0 Bar
13184 2020-01-03 00:00:00 1 0.0 Cl
14648 2020-01-03 00:00:00 1 0.0 LevelBoundary
16112 2020-01-03 00:00:00 1 1.163008 Continuity
9 2020-01-03 00:00:00 3 0.010151 Initial
1473 2020-01-03 00:00:00 3 0.019046 Foo
2937 2020-01-03 00:00:00 3 0.0 Drainage
4401 2020-01-03 00:00:00 3 0.0 UserDemand
5865 2020-01-03 00:00:00 3 0.006279 FlowBoundary
7329 2020-01-03 00:00:00 3 0.06482 Precipitation
8793 2020-01-03 00:00:00 3 0.0 Terminal
10257 2020-01-03 00:00:00 3 0.071098 Tracer
11721 2020-01-03 00:00:00 3 0.0 Bar
13185 2020-01-03 00:00:00 3 0.647569 Cl
14649 2020-01-03 00:00:00 3 0.019046 LevelBoundary
16113 2020-01-03 00:00:00 3 0.100296 Continuity
10 2020-01-03 00:00:00 6 0.001409 Initial
1474 2020-01-03 00:00:00 6 0.001741 Foo
2938 2020-01-03 00:00:00 6 0.0 Drainage
4402 2020-01-03 00:00:00 6 0.0 UserDemand
5866 2020-01-03 00:00:00 6 0.000573 FlowBoundary
7330 2020-01-03 00:00:00 6 0.005924 Precipitation
8794 2020-01-03 00:00:00 6 0.0 Terminal
10258 2020-01-03 00:00:00 6 0.006497 Tracer
11722 2020-01-03 00:00:00 6 0.0 Bar
13186 2020-01-03 00:00:00 6 0.059185 Cl
14650 2020-01-03 00:00:00 6 0.001741 LevelBoundary
16114 2020-01-03 00:00:00 6 0.009647 Continuity
11 2020-01-03 00:00:00 9 0.003869 Initial
1475 2020-01-03 00:00:00 9 0.00003 Foo
2939 2020-01-03 00:00:00 9 0.0 Drainage
4403 2020-01-03 00:00:00 9 0.0 UserDemand
5867 2020-01-03 00:00:00 9 0.00001 FlowBoundary
7331 2020-01-03 00:00:00 9 0.016115 Precipitation
8795 2020-01-03 00:00:00 9 0.0 Terminal
10259 2020-01-03 00:00:00 9 0.016124 Tracer
11723 2020-01-03 00:00:00 9 0.0 Bar
13187 2020-01-03 00:00:00 9 0.001035 Cl
14651 2020-01-03 00:00:00 9 0.00003 LevelBoundary
16115 2020-01-03 00:00:00 9 0.020023 Continuity

We can use this table to plot the results of the Delwaq model, both spatially as over time.

from ribasim.delwaq import plot_fraction

plot_fraction(nmodel, 1)  # default tracers, should add up to 1
plot_fraction(nmodel, 9, ["Foo", "Bar"])  # custom tracers
plot_fraction(nmodel, 9, ["Continuity"])  # mass balance check

from ribasim.delwaq import plot_spatial

plot_spatial(nmodel, "Bar")
plot_spatial(nmodel, "Foo", versus="Bar")  # ratio of Meuse to Rhine