Simulated and defined demand do not match

[Fati-Mon]

What When simulating the demand over time it does not show the same demand that has been given in the code. Example: "The farmer needs to apply irrigation to its field starting from April to September. The irrigated field is 0.32km² and receives 10m³/s and 12m³/s of diverted water during spring and summer, respectively." Thus irrigation demand required in certain periods of time, not the whole year. The rest should remain zero:

# DEMAND (irrigation)
model.user_demand.add(
    Node(6, Point(-1.5, 1.0), name='IrrA'),
    [user_demand.Time(
        demand=[0.0, 10, 12, 0.0],
        return_factor=[0.0, 0.0, 0.0, 0.0],
        min_level=[0.0, 0.0, 0.0, 0.0],
        priority=[5, 1, 1, 5],
        time=[starttime, "2023-04-01", "2023-07-01", "2023-10-01"]
        )
    ]
)

When plotting we get:

It seems to interpolate, but continues to abstract water even though it shouldn't.

How Maybe my understanding of this node is wrong, either way if someone could clarify this to me, that would be helpfull.

Additional context Full code:

import shutil
from pathlib import Path
import pathlib 
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ribasim import Allocation, Model, Node # The main library used for river basin modeling.
from ribasim.nodes import (
    flow_boundary,
    basin,
    tabulated_rating_curve,
    user_demand,
    terminal
)
from shapely.geometry import Point
import subprocess  # For running the model
import plotly.express as px

base_dir = Path("c:/Ribasim9")
model_dir = base_dir / "Crystal_Basin"
data_path = model_dir / "data/input/ACTINFLW.csv"

starttime = "2023-01-01"
endtime = "2024-01-01"
model = Model(
    starttime=starttime,
    endtime=endtime,
    crs="EPSG:4326",
)


#%% FLOW BOUNDARY

# FLOW BOUNDARY
data = pd.read_csv(data_path, sep=";")
data['sum']= data['minor']+data['main']
print('Average inflowQ m3/s:',data['sum'].mean())

model.flow_boundary.add(
    Node(1, Point(0.0, 0.0), name='Main'), 
    [flow_boundary.Time(time=data.time, flow_rate=data.main, #name="Main"
    )]
)

model.flow_boundary.add(
    Node(2, Point(-3.0, 0.0), name='Minor'), 
    [flow_boundary.Time(time=data.time, flow_rate=data.minor, #name="Main"
    )]
)

# BASIN (diversion and confluence)

model.basin.add(
    Node(3, Point(-0.75, -0.5), name='Div'),
    [basin.Profile(area=[950, 40000], level=[0, 15]),
     basin.State(level=[11]),
     basin.Time(time=[starttime, endtime]),
    ],
)

model.basin.add(
    Node(4, Point(-1.5, -1), name='Conf'),
    [basin.Profile(area=[1000, 43000], level=[0, 15]),
     basin.State(level=[11]),
     basin.Time(time=[starttime, endtime]),
    ],
)


model.tabulated_rating_curve.add(
    Node(5, Point(-1.5, -1.5), name='Maingate'),
    [tabulated_rating_curve.Static(
        level=[0.0, 0.1],
        flow_rate=[0.0, 50],
        )
    ]
)

# DEMAND (irrigation)
model.user_demand.add(
    Node(6, Point(-1.5, 1.0), name='IrrA'),
    [user_demand.Time(
        demand=[0.0, 10, 12, 0.0],
        return_factor=[0.0, 0.0, 0.0, 0.0],
        min_level=[0.0, 0.0, 0.0, 0.0],
        priority=[5, 1, 1, 5],
        time=[starttime, "2023-04-01", "2023-07-01", "2023-10-01"]
        )
    ]
)


model.tabulated_rating_curve.add(
    Node(7, Point(-1.125, -0.75), name='Main'),
    [tabulated_rating_curve.Static(
        level=[0.0, 0.1],
        flow_rate=[0.0, 50],
        )
    ]
)


# TERMINAL

model.terminal.add(Node(8, Point(-1.5, -3.0), name="Terminal"))

# EDGES

model.edge.add(model.flow_boundary[1], model.basin[3])
model.edge.add(model.flow_boundary[2], model.basin[4])
model.edge.add(model.basin[3], model.user_demand[6])
model.edge.add(model.user_demand[6], model.basin[4])
model.edge.add(model.basin[3], model.tabulated_rating_curve[7])
model.edge.add(model.tabulated_rating_curve[7], model.basin[4])
model.edge.add(model.basin[4], model.tabulated_rating_curve[5])
model.edge.add(model.tabulated_rating_curve[5], model.terminal[8])

# SCHEMATIZATION

model.plot()

toml_path = model_dir/ "Crystal_1.2/ribasim.toml"
model.write(toml_path)
rib_path = base_dir / "ribasim_windows/ribasim.exe"

subprocess.run([rib_path, toml_path], check=True)

# result = subprocess.run([rib_path, toml_path], capture_output=True)
# result.check_returncode()
# print(result.stderr)

#%%

# Dictionary mapping node_ids to names
edge_names = {
    (1,3): 'Main',
    (2,4): 'Minor',
    (3,6): 'IrrA Demand',
    (6,4): 'IrrA Drain',
    (3,7): 'Div2Main',
    (7,4): 'Main2Conf',
    (4,5): 'Conf2TRC',
    (5,8): 'TRC2Term',
}

# Dictionary mapping basins (node_ids) to names
node_names = {
    3: 'Div',
    4: 'Conf',
}

df_basin = pd.read_feather(model_dir / "Crystal_1.2/results/basin.arrow")
# Adding the 'name' column using the dictionary
df_basin['name'] = df_basin['node_id'].map(node_names)
df_basin['dif_IN_OUT'] = df_basin["inflow_rate"] - df_basin["outflow_rate"]
df_basin_dif = df_basin.pivot_table(
    index="time", columns="name", values=["dif_IN_OUT"]
)
df_basin_dif.plot()

# Create pivot tables and plot for basin data
df_basin_wide = df_basin.pivot_table(
    index="time", columns="name", values=["storage", "level"]
)

# Skip the first timestep
df_basin_wide = df_basin_wide.iloc[1:]



#%% Create subplots for Div and Conf basin data

# Create pivot tables and plot for basin data
df_basin_div = df_basin_wide.xs('Div', axis=1, level=1, drop_level=False)
df_basin_conf = df_basin_wide.xs('Conf', axis=1, level=1, drop_level=False)

def plot_basin_data(ax, ax_twin, df_basin, level_color='b', storage_color='r', title='Basin'):
    # Plot level data
    for idx, column in enumerate(df_basin["level"].columns):
        ax.plot(df_basin.index, df_basin["level"][column], 
                linestyle='-', color=level_color, 
                label=f'Level - {column}')
    
    # Plot storage data
    for idx, column in enumerate(df_basin["storage"].columns):
        ax_twin.plot(df_basin.index, df_basin["storage"][column], 
                     linestyle='--', color=storage_color, 
                     label=f'Storage - {column}')
    
    ax.set_ylabel('Level [m]', color=level_color)
    ax_twin.set_ylabel('Storage [m³]', color=storage_color)
    
    ax.tick_params(axis='y', labelcolor=level_color)
    ax_twin.tick_params(axis='y', labelcolor=storage_color)
    
    ax.set_title(title)

    # Combine legends from both axes
    lines, labels = ax.get_legend_handles_labels()
    lines_twin, labels_twin = ax_twin.get_legend_handles_labels()
    ax.legend(lines + lines_twin, labels + labels_twin, loc='upper left')

# Create subplots
fig, (ax1, ax3) = plt.subplots(2, 1, figsize=(12, 12), sharex=True)

# Plot Div basin data
ax2 = ax1.twinx()  # Secondary y-axis for storage
plot_basin_data(ax1, ax2, df_basin_div, title='Div Basin Level and Storage Over Time')

# Plot Conf basin data
ax4 = ax3.twinx()  # Secondary y-axis for storage
plot_basin_data(ax3, ax4, df_basin_conf, title='Conf Basin Level and Storage Over Time')

# Common X label
ax3.set_xlabel('Time')

fig.tight_layout()  # Adjust layout to fit labels
plt.show()


#%% Plot flow data

# Sample data loading and preparation
df_flow = pd.read_feather(model_dir / "Crystal_1.2/results/flow.arrow")
df_flow["edge"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))
df_flow["name"] = df_flow["edge"].map(edge_names)
df_flow["flow_m3d"] = df_flow.flow_rate * 86400

# Plot the flow data, interactive plot with Plotly
pivot_flow = df_flow.pivot_table(index="time", columns="name", values="flow_rate").reset_index()
fig = px.line(pivot_flow, x="time", y=pivot_flow.columns[1:], title="Flow Over Time [m3/s]")

fig.update_layout(legend_title_text='Edge')
fig.write_html("plot_edges.html")
fig.show()

[gijsber]

Hi @Fati-Mon, Wether it is convenient for the modeller or not is another question, but can you check whether the following provides the result you expect ?: [user_demand.Time( demand=[0.0, 0.0, 10, 10, 12, 12, 0.0], return_factor=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], min_level=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], priority=[5, 5, 1, 1, 1, 1, 5], time=[starttime, "2023-03-31", "2023-04-01", "2023-06-30", "2023-07-01", "2023-09-30", "2023-10-01"]

[Fati-Mon]

Hi @gijsber ,

A slight difference, but unfortunately still the same issue.

[SouthEndMusic]

This problem is a bit subtle, it has to do with all the data that can be packed in the UserDemand / time table in long table format. You have to think about multiple priorities as parallel in time, so you cannot have one demand whose priority changes over time.

What you now specified is this (extrapolation not shown):

[Jingru923]

If trying with one priority, the result gets much more sensible

        user_demand.Time(
            demand=[0.0, 10, 12, 0.0],
            return_factor=[0.0, 0.0, 0.0, 0.0],
            min_level=[0.0, 0.0, 0.0, 0.0],
            priority=[1,1,1,1]#edited
            time=[starttime, "2023-04-01", "2023-07-01", "2023-10-01"],
        )

[Jingru923]

After fix the priority, I ran the models a few times with different demand value. I noticed that flow will approach the first non-zero demand, but once time pass "2023-04-01", the flow start to decrease and never approach the second non-zero demand. This happens even I set

demand=[0.0, 10, 10, 0.0],

The second non-zero demand seems to have no influence on the flow.

[Fati-Mon]

When setting the priority to 1 for every time period, I now get:

[SouthEndMusic]

Your model is for 2023 but your plot shows 2022, what's going on there? Can you show exactly what you changed? It would be nice to have things like this in a PR.