WSIMOD model demonstration - Oxford (.py)¶

Note - this script can also be opened in interactive Python if you wanted to play around. On the GitHub it is in docs/demo/scripts

1. Introduction

2. Data

3. Nodes

   3.1 Freshwater Treatment Works

   3.2 Land

   3.3 Demand

   3.4 Reservoir

   3.5 Distribution

   3.6 Wastewater Treatment Works

   3.7 Sewers

   3.8 Groundwater

   3.9 Node list

4. Arcs

   4.1 Arc parameters

   4.2 Create the arcs

5. Mapping

6. Orchestration

   6.1 Orchestrating an individual timestep

   6.2 Ending a timestep

7. Model object

   7.1 Validation

We will cover a demo WSIMOD case study¶

The glamorous town of Oxford will be our demo case study. Below, we will create these nodes and arcs, orchestrate them into a model, and run simulations.

Although GIS is pretty, the schematic below is a more accurate representation of what will be created. WSIMOD treats everything as a node or an arc.

Imports and forcing data¶

Import packages

In [1]:

import os

import pandas as pd

from wsimod.arcs.arcs import Arc
from wsimod.core import constants
from wsimod.nodes.catchment import Catchment
from wsimod.nodes.demand import ResidentialDemand
from wsimod.nodes.land import Land
from wsimod.nodes.nodes import Node
from wsimod.nodes.sewer import Sewer
from wsimod.nodes.storage import Groundwater, Reservoir
from wsimod.nodes.waste import Waste
from wsimod.nodes.wtw import FWTW, WWTW
from wsimod.orchestration.model import Model

os.environ["USE_PYGEOS"] = "0"
import geopandas as gpd
from matplotlib import pyplot as plt
from shapely.geometry import LineString

import os import pandas as pd from wsimod.arcs.arcs import Arc from wsimod.core import constants from wsimod.nodes.catchment import Catchment from wsimod.nodes.demand import ResidentialDemand from wsimod.nodes.land import Land from wsimod.nodes.nodes import Node from wsimod.nodes.sewer import Sewer from wsimod.nodes.storage import Groundwater, Reservoir from wsimod.nodes.waste import Waste from wsimod.nodes.wtw import FWTW, WWTW from wsimod.orchestration.model import Model os.environ["USE_PYGEOS"] = "0" import geopandas as gpd from matplotlib import pyplot as plt from shapely.geometry import LineString

Load input data

In [2]:

# Select the root path for the data folder. Use the appropriate value for your case.
data_folder = os.path.join(os.path.abspath(""), "docs", "demo", "data")

input_fid = os.path.join(data_folder, "processed", "timeseries_data.csv")
input_data = pd.read_csv(input_fid)
input_data.loc[input_data.variable == "flow", "value"] *= constants.M3_S_TO_M3_DT
input_data.loc[input_data.variable == "precipitation", "value"] *= constants.MM_TO_M
input_data.date = pd.to_datetime(input_data.date)
data_input_dict = input_data.set_index(["variable", "date"]).value.to_dict()
data_input_dict = (
    input_data.groupby("site")
    .apply(lambda x: x.set_index(["variable", "date"]).value.to_dict())
    .to_dict()
)
print(input_data.sample(10))

# Select the root path for the data folder. Use the appropriate value for your case. data_folder = os.path.join(os.path.abspath(""), "docs", "demo", "data") input_fid = os.path.join(data_folder, "processed", "timeseries_data.csv") input_data = pd.read_csv(input_fid) input_data.loc[input_data.variable == "flow", "value"] *= constants.M3_S_TO_M3_DT input_data.loc[input_data.variable == "precipitation", "value"] *= constants.MM_TO_M input_data.date = pd.to_datetime(input_data.date) data_input_dict = input_data.set_index(["variable", "date"]).value.to_dict() data_input_dict = ( input_data.groupby("site") .apply(lambda x: x.set_index(["variable", "date"]).value.to_dict()) .to_dict() ) print(input_data.sample(10))

           site       date     variable          value
58304  cherwell 2009-05-06    magnesium       0.008950
47200  cherwell 2010-11-01    potassium       0.006571
2630   evenlode 2012-05-20  temperature      12.700000
11285    thames 2012-02-29      silicon       0.002878
43542  evenlode 2012-10-11       sodium       0.010512
45967    thames 2011-06-12       sodium       0.031786
50259       ray 2011-03-28    potassium       0.008143
71144  cherwell 2012-08-10         flow  201312.000000
71246  cherwell 2012-11-20         flow  279072.000000
52032    thames 2012-02-08    potassium       0.004900

/tmp/ipykernel_2096/946595242.py:12: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
  .apply(lambda x: x.set_index(["variable", "date"]).value.to_dict())

Input data is stored in dicts

In [3]:

print(data_input_dict["cherwell"][("boron", pd.to_datetime("2010-11-20"))])

print(data_input_dict["cherwell"][("boron", pd.to_datetime("2010-11-20"))])

6.985714285714286e-05

We select dates that are available in the input data

In [4]:

dates = input_data.date.unique()
dates = dates[dates.argsort()]
dates = [pd.Timestamp(x) for x in dates]
print(dates[0:10])

dates = input_data.date.unique() dates = dates[dates.argsort()] dates = [pd.Timestamp(x) for x in dates] print(dates[0:10])

[Timestamp('2009-03-03 00:00:00'), Timestamp('2009-03-04 00:00:00'), Timestamp('2009-03-05 00:00:00'), Timestamp('2009-03-06 00:00:00'), Timestamp('2009-03-07 00:00:00'), Timestamp('2009-03-08 00:00:00'), Timestamp('2009-03-09 00:00:00'), Timestamp('2009-03-10 00:00:00'), Timestamp('2009-03-11 00:00:00'), Timestamp('2009-03-12 00:00:00')]

We can specify the pollutants. In this example we choose based on what pollutants we have input data for.

In [5]:

constants.POLLUTANTS = input_data.variable.unique().tolist()
constants.POLLUTANTS.remove("flow")
constants.POLLUTANTS.remove("precipitation")
constants.POLLUTANTS.remove("et0")
constants.NON_ADDITIVE_POLLUTANTS = ["temperature"]
constants.ADDITIVE_POLLUTANTS = list(
    set(constants.POLLUTANTS).difference(constants.NON_ADDITIVE_POLLUTANTS)
)
constants.FLOAT_ACCURACY = 1e-11
print(constants.POLLUTANTS)

constants.POLLUTANTS = input_data.variable.unique().tolist() constants.POLLUTANTS.remove("flow") constants.POLLUTANTS.remove("precipitation") constants.POLLUTANTS.remove("et0") constants.NON_ADDITIVE_POLLUTANTS = ["temperature"] constants.ADDITIVE_POLLUTANTS = list( set(constants.POLLUTANTS).difference(constants.NON_ADDITIVE_POLLUTANTS) ) constants.FLOAT_ACCURACY = 1e-11 print(constants.POLLUTANTS)

['temperature', 'silicon', 'fluoride', 'chloride', 'nitrate', 'sulphate', 'nitrogen', 'sodium', 'potassium', 'calcium', 'magnesium', 'boron']

Create nodes¶

For waste nodes, no parameters are needed, they are just the model outlet

In [6]:

thames_above_abingdon = Waste(name="thames_above_abingdon")

thames_above_abingdon = Waste(name="thames_above_abingdon")

For junctions and abstraction locations, we can simply use the default nodes

In [7]:

farmoor_abstraction = Node(name="farmoor_abstraction")
evenlode_thames = Node(name="evenlode_thames")
cherwell_ray = Node(name="cherwell_ray")
cherwell_thames = Node(name="cherwell_thames")
thames_mixer = Node(name="thames_mixer")

farmoor_abstraction = Node(name="farmoor_abstraction") evenlode_thames = Node(name="evenlode_thames") cherwell_ray = Node(name="cherwell_ray") cherwell_thames = Node(name="cherwell_thames") thames_mixer = Node(name="thames_mixer")

For catchment nodes, we only need to specify the input data (as a dictionary format).

In [8]:

evenlode = Catchment(name="evenlode", data_input_dict=data_input_dict["evenlode"])
thames = Catchment(name="thames", data_input_dict=data_input_dict["thames"])
ray = Catchment(name="ray", data_input_dict=data_input_dict["ray"])
cherwell = Catchment(name="cherwell", data_input_dict=data_input_dict["cherwell"])

evenlode = Catchment(name="evenlode", data_input_dict=data_input_dict["evenlode"]) thames = Catchment(name="thames", data_input_dict=data_input_dict["thames"]) ray = Catchment(name="ray", data_input_dict=data_input_dict["ray"]) cherwell = Catchment(name="cherwell", data_input_dict=data_input_dict["cherwell"])

We can see that, even though we provided mimimal information (name and input data) each node comes with many predefined functions.

In [9]:

print(dir(evenlode))

print(dir(evenlode))

['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'apply_overrides', 'blend_vqip', 'check_basic', 'compare_vqip', 'concentration_to_total', 'copy_vqip', 'data_input_dict', 'ds_vqip', 'ds_vqip_c', 'empty_vqip', 'empty_vqip_predefined', 'end_timestep', 'end_timestep_', 'extract_vqip', 'extract_vqip_c', 'get_avail', 'get_connected', 'get_data_input', 'get_direction_arcs', 'get_flow', 'in_arcs', 'in_arcs_type', 'mass_balance', 'mass_balance_ds', 'mass_balance_in', 'mass_balance_out', 'name', 'node_mass_balance', 'out_arcs', 'out_arcs_type', 'pull_check', 'pull_check_abstraction', 'pull_check_basic', 'pull_check_deny', 'pull_check_handler', 'pull_distributed', 'pull_set', 'pull_set_abstraction', 'pull_set_deny', 'pull_set_handler', 'push_check', 'push_check_accept', 'push_check_basic', 'push_check_deny', 'push_check_handler', 'push_distributed', 'push_set', 'push_set_deny', 'push_set_handler', 'query_handler', 'reinit', 'route', 'sum_vqip', 't', 'total_in', 'total_out', 'total_to_concentration', 'unrouted_water', 'v_change_vqip', 'v_change_vqip_c', 'v_distill_vqip', 'v_distill_vqip_c']

Freshwater treatment works¶

Each type of node uses different parameters (see API reference). Below we create a freshwater treatment works (FWTW)

In [10]:

oxford_fwtw = FWTW(
    service_reservoir_storage_capacity=1e5,
    service_reservoir_storage_area=2e4,
    treatment_throughput_capacity=4.5e4,
    name="oxford_fwtw",
)

oxford_fwtw = FWTW( service_reservoir_storage_capacity=1e5, service_reservoir_storage_area=2e4, treatment_throughput_capacity=4.5e4, name="oxford_fwtw", )

Each node type has different types of functionality available

In [11]:

print(dir(oxford_fwtw))

print(dir(oxford_fwtw))

['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_liquor_multiplier', '_percent_solids', 'apply_overrides', 'blend_vqip', 'calculate_volume', 'check_basic', 'compare_vqip', 'concentration_to_total', 'copy_vqip', 'current_input', 'data_input_dict', 'ds_vqip', 'ds_vqip_c', 'empty_vqip', 'empty_vqip_predefined', 'end_timestep', 'end_timestep_', 'extract_vqip', 'extract_vqip_c', 'get_connected', 'get_data_input', 'get_direction_arcs', 'get_excess_throughput', 'in_arcs', 'in_arcs_type', 'liquor', 'liquor_multiplier', 'mass_balance', 'mass_balance_ds', 'mass_balance_in', 'mass_balance_out', 'name', 'node_mass_balance', 'out_arcs', 'out_arcs_type', 'percent_solids', 'previous_pulled', 'process_parameters', 'pull_check', 'pull_check_basic', 'pull_check_deny', 'pull_check_fwtw', 'pull_check_handler', 'pull_distributed', 'pull_set', 'pull_set_deny', 'pull_set_fwtw', 'pull_set_handler', 'push_check', 'push_check_accept', 'push_check_basic', 'push_check_deny', 'push_check_handler', 'push_distributed', 'push_set', 'push_set_deny', 'push_set_handler', 'query_handler', 'reinit', 'service_reservoir_initial_storage', 'service_reservoir_storage_area', 'service_reservoir_storage_capacity', 'service_reservoir_storage_elevation', 'service_reservoir_tank', 'solids', 'sum_vqip', 't', 'total_deficit', 'total_in', 'total_out', 'total_pulled', 'total_to_concentration', 'treat_current_input', 'treat_water', 'treated', 'treatment_throughput_capacity', 'unpushed_sludge', 'v_change_vqip', 'v_change_vqip_c', 'v_distill_vqip', 'v_distill_vqip_c']

The FWTW node has a tank representing the service reservoirs, we can see that it has been initialised empty.

In [12]:

print(oxford_fwtw.service_reservoir_tank.storage)

print(oxford_fwtw.service_reservoir_tank.storage)

{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}

If we try to pull water from the FWTW, it responds that there is no water to pull.

In [13]:

print(oxford_fwtw.pull_check({"volume": 10}))

print(oxford_fwtw.pull_check({"volume": 10}))

{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}

If we add in some water, we see the pull check responds that water is available.

In [14]:

oxford_fwtw.service_reservoir_tank.storage["volume"] += 25
print(oxford_fwtw.pull_check({"volume": 10}))

oxford_fwtw.service_reservoir_tank.storage["volume"] += 25 print(oxford_fwtw.pull_check({"volume": 10}))

{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 10.0}

When we set a pull request, we see that we successfully receive the water and the tank is updated.

In [15]:

reply = oxford_fwtw.pull_set({"volume": 10})
print(reply)
print(oxford_fwtw.service_reservoir_tank.storage)

reply = oxford_fwtw.pull_set({"volume": 10}) print(reply) print(oxford_fwtw.service_reservoir_tank.storage)

{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 10.0}
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 15.0}

Land¶

We will now create a land node, it is a bit involved so you might want to skip ahead to demand, or check out the land node tutorial

Data inputs are a single dictionary

In [16]:

land_inputs = data_input_dict["oxford_land"]

print(land_inputs[("precipitation", pd.to_datetime("2010-11-20"))])
print(land_inputs[("et0", pd.to_datetime("2010-11-20"))])
print(land_inputs[("temperature", pd.to_datetime("2009-10-15"))])

land_inputs = data_input_dict["oxford_land"] print(land_inputs[("precipitation", pd.to_datetime("2010-11-20"))]) print(land_inputs[("et0", pd.to_datetime("2010-11-20"))]) print(land_inputs[("temperature", pd.to_datetime("2009-10-15"))])

0.0003
0.002
13.107142857142858

Assign some default pollutant deposition values (kg/m2/d)

In [17]:

pollutant_deposition = {
    "boron": 100e-10,
    "calcium": 70e-7,
    "chloride": 60e-10,
    "fluoride": 0.2e-7,
    "magnesium": 6e-7,
    "nitrate": 2e-9,
    "nitrogen": 4e-7,
    "potassium": 7e-7,
    "silicon": 7e-9,
    "sodium": 30e-9,
    "sulphate": 70e-7,
}

pollutant_deposition = { "boron": 100e-10, "calcium": 70e-7, "chloride": 60e-10, "fluoride": 0.2e-7, "magnesium": 6e-7, "nitrate": 2e-9, "nitrogen": 4e-7, "potassium": 7e-7, "silicon": 7e-9, "sodium": 30e-9, "sulphate": 70e-7, }

Create two surfaces as a list of dicts

In [18]:

surface = [
    {
        "type_": "PerviousSurface",
        "area": 2e7,
        "pollutant_load": pollutant_deposition,
        "surface": "rural",
        "field_capacity": 0.3,
        "depth": 0.5,
        "initial_storage": 2e7 * 0.4 * 0.5,
    },
    {
        "type_": "ImperviousSurface",
        "area": 1e7,
        "pollutant_load": pollutant_deposition,
        "surface": "urban",
        "initial_storage": 5e6,
    },
]

surface = [ { "type_": "PerviousSurface", "area": 2e7, "pollutant_load": pollutant_deposition, "surface": "rural", "field_capacity": 0.3, "depth": 0.5, "initial_storage": 2e7 * 0.4 * 0.5, }, { "type_": "ImperviousSurface", "area": 1e7, "pollutant_load": pollutant_deposition, "surface": "urban", "initial_storage": 5e6, }, ]

Create the land node from these surfaces and the input data

In [19]:

oxford_land = Land(surfaces=surface, name="oxford_land", data_input_dict=land_inputs)

oxford_land = Land(surfaces=surface, name="oxford_land", data_input_dict=land_inputs)

We can see the land node has various tanks that have been initialised empty

In [20]:

print(oxford_land.surface_runoff.storage)
print(oxford_land.subsurface_runoff.storage)
print(oxford_land.percolation.storage)

print(oxford_land.surface_runoff.storage) print(oxford_land.subsurface_runoff.storage) print(oxford_land.percolation.storage)

{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}
{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}
{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 0}

We can see the surfaces have also been initialised, although they are not empty because we provided 'initial_storage' parameters.

In [21]:

rural_surface = oxford_land.get_surface("rural")
urban_surface = oxford_land.get_surface("urban")
print("{0}-{1}".format("rural", rural_surface.storage))
print("{0}-{1}".format("urban", urban_surface.storage))

rural_surface = oxford_land.get_surface("rural") urban_surface = oxford_land.get_surface("urban") print("{0}-{1}".format("rural", rural_surface.storage)) print("{0}-{1}".format("urban", urban_surface.storage))

rural-{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 4000000.0}
urban-{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 5000000.0}

We can run a timestep of the land node with the 'run' command

In [22]:

oxford_land.t = pd.to_datetime("2012-12-22")

oxford_land.run()

oxford_land.t = pd.to_datetime("2012-12-22") oxford_land.run()

We can see that the land and surface tanks have been updated

In [23]:

print(oxford_land.surface_runoff.storage)
print(oxford_land.subsurface_runoff.storage)
print(oxford_land.percolation.storage)

print("{0}-{1}".format("rural", rural_surface.storage))
print("{0}-{1}".format("urban", urban_surface.storage))

print(oxford_land.surface_runoff.storage) print(oxford_land.subsurface_runoff.storage) print(oxford_land.percolation.storage) print("{0}-{1}".format("rural", rural_surface.storage)) print("{0}-{1}".format("urban", urban_surface.storage))

{'temperature': 3.521534008631931, 'silicon': 0.00020062235262918442, 'fluoride': 0.0005732067217976697, 'chloride': 0.0001719620165393009, 'nitrate': 5.732067217976698e-05, 'sulphate': 0.20062235262918443, 'nitrogen': 0.011464134435953396, 'sodium': 0.0008598100826965045, 'potassium': 0.020062235262918445, 'calcium': 0.20062235262918443, 'magnesium': 0.017196201653930095, 'boron': 0.00028660336089883487, 'volume': 12400.000000000002}
{'temperature': 0.9722493642718524, 'silicon': 0.0017554455855053634, 'fluoride': 0.0050155588157296096, 'chloride': 0.0015046676447188828, 'nitrate': 0.000501555881572961, 'sulphate': 1.7554455855053634, 'nitrogen': 0.1003111763145922, 'sodium': 0.007523338223594414, 'potassium': 0.17554455855053633, 'calcium': 1.7554455855053634, 'magnesium': 0.1504667644718883, 'boron': 0.0025077794078648048, 'volume': 58900.0}
{'temperature': 0.3349282031818326, 'silicon': 0.005868203814403644, 'fluoride': 0.016766296612581843, 'chloride': 0.005029888983774552, 'nitrate': 0.001676629661258184, 'sulphate': 5.868203814403645, 'nitrogen': 0.3353259322516368, 'sodium': 0.02514944491887276, 'potassium': 0.5868203814403645, 'calcium': 5.868203814403645, 'magnesium': 0.5029888983774553, 'boron': 0.008383148306290921, 'volume': 176700.00000000003}
rural-{'temperature': 6.203473168254872, 'silicon': 0.1321757282474618, 'fluoride': 0.3776449378498909, 'chloride': 0.11329348135496725, 'nitrate': 0.037764493784989084, 'sulphate': 132.17572824746182, 'nitrogen': 7.552898756997817, 'sodium': 0.5664674067748363, 'potassium': 13.21757282474618, 'calcium': 132.17572824746182, 'magnesium': 11.329348135496726, 'boron': 0.18882246892494545, 'volume': 3980000.0}
urban-{'temperature': 0.13343969249141666, 'silicon': 0.06999999999999999, 'fluoride': 0.2, 'chloride': 0.06, 'nitrate': 0.02, 'sulphate': 70.0, 'nitrogen': 4.0, 'sodium': 0.3, 'potassium': 7.0, 'calcium': 70.0, 'magnesium': 6.0, 'boron': 0.1, 'volume': 5104000.0}

Residential demand¶

The residential demand node requires population, per capita demand and a pollutant_load dictionary that defines how much (weight in kg) pollution is generated per person per day.

In [24]:

oxford = ResidentialDemand(
    name="oxford",
    population=2e5,
    per_capita=0.15,
    pollutant_load={
        "boron": 500 * constants.UG_L_TO_KG_M3 * 0.15,
        "calcium": 150 * constants.MG_L_TO_KG_M3 * 0.15,
        "chloride": 180 * constants.MG_L_TO_KG_M3 * 0.15,
        "fluoride": 0.4 * constants.MG_L_TO_KG_M3 * 0.15,
        "magnesium": 30 * constants.MG_L_TO_KG_M3 * 0.15,
        "nitrate": 60 * constants.MG_L_TO_KG_M3 * 0.15,
        "nitrogen": 50 * constants.MG_L_TO_KG_M3 * 0.15,
        "potassium": 30 * constants.MG_L_TO_KG_M3 * 0.15,
        "silicon": 20 * constants.MG_L_TO_KG_M3 * 0.15,
        "sodium": 200 * constants.MG_L_TO_KG_M3 * 0.15,
        "sulphate": 250 * constants.MG_L_TO_KG_M3 * 0.15,
        "temperature": 14,
    },
    data_input_dict=land_inputs,
)
# pollutant_load calculated based on expected effluent at WWTW

oxford = ResidentialDemand( name="oxford", population=2e5, per_capita=0.15, pollutant_load={ "boron": 500 * constants.UG_L_TO_KG_M3 * 0.15, "calcium": 150 * constants.MG_L_TO_KG_M3 * 0.15, "chloride": 180 * constants.MG_L_TO_KG_M3 * 0.15, "fluoride": 0.4 * constants.MG_L_TO_KG_M3 * 0.15, "magnesium": 30 * constants.MG_L_TO_KG_M3 * 0.15, "nitrate": 60 * constants.MG_L_TO_KG_M3 * 0.15, "nitrogen": 50 * constants.MG_L_TO_KG_M3 * 0.15, "potassium": 30 * constants.MG_L_TO_KG_M3 * 0.15, "silicon": 20 * constants.MG_L_TO_KG_M3 * 0.15, "sodium": 200 * constants.MG_L_TO_KG_M3 * 0.15, "sulphate": 250 * constants.MG_L_TO_KG_M3 * 0.15, "temperature": 14, }, data_input_dict=land_inputs, ) # pollutant_load calculated based on expected effluent at WWTW

Reservoir¶

A reservoir node is used to make abstractions from rivers and supply FWTWs

In [25]:

farmoor = Reservoir(
    name="farmoor", capacity=1e7, initial_storage=1e7, area=1.5e6, datum=62
)

farmoor = Reservoir( name="farmoor", capacity=1e7, initial_storage=1e7, area=1.5e6, datum=62 )

Distribution¶

We use a generic Node as a junction to represent the distribution network between the FWTW and households

In [26]:

distribution = Node(name="oxford_distribution")

distribution = Node(name="oxford_distribution")

Wastewater treatment works¶

Wastewater treatment works (WWTW) are nodes that can store sewage water temporarily in storm tanks, and reduce the pollution amounts in water before releasing them onwards to rivers.

In [27]:

oxford_wwtw = WWTW(
    stormwater_storage_capacity=2e4,
    stormwater_storage_area=2e4,
    treatment_throughput_capacity=5e4,
    name="oxford_wwtw",
)

oxford_wwtw = WWTW( stormwater_storage_capacity=2e4, stormwater_storage_area=2e4, treatment_throughput_capacity=5e4, name="oxford_wwtw", )

Sewers¶

Sewer nodes enable water to transition between households and WWTWs, and between impervious surfaces and rivers or WWTWs. They use a timearea diagram to represent travel time, which assigns a specified percentage of water to take a specified duration to pass through the sewer node.

In [28]:

combined_sewer = Sewer(
    capacity=4e6, pipe_timearea={0: 0.8, 1: 0.15, 2: 0.05}, name="combined_sewer"
)

combined_sewer = Sewer( capacity=4e6, pipe_timearea={0: 0.8, 1: 0.15, 2: 0.05}, name="combined_sewer" )

Groundwater¶

Groundwater nodes implement a simple residence time to determine baseflow

In [29]:

gw = Groundwater(capacity=3.2e9, area=3.2e8, name="gw", residence_time=20)

gw = Groundwater(capacity=3.2e9, area=3.2e8, name="gw", residence_time=20)

Create a nodelist¶

To keep all the nodes in one place, we put them into a list

In [30]:

nodelist = [
    thames_above_abingdon,
    evenlode,
    thames,
    ray,
    cherwell,
    oxford,
    distribution,
    farmoor,
    oxford_fwtw,
    oxford_wwtw,
    combined_sewer,
    oxford_land,
    gw,
    farmoor_abstraction,
    evenlode_thames,
    cherwell_ray,
    cherwell_thames,
    thames_mixer,
]

print(nodelist)

nodelist = [ thames_above_abingdon, evenlode, thames, ray, cherwell, oxford, distribution, farmoor, oxford_fwtw, oxford_wwtw, combined_sewer, oxford_land, gw, farmoor_abstraction, evenlode_thames, cherwell_ray, cherwell_thames, thames_mixer, ] print(nodelist)

[<wsimod.nodes.waste.Waste object at 0x7f97577b54d0>, <wsimod.nodes.catchment.Catchment object at 0x7f97577a36d0>, <wsimod.nodes.catchment.Catchment object at 0x7f97577a3590>, <wsimod.nodes.catchment.Catchment object at 0x7f97577a2750>, <wsimod.nodes.catchment.Catchment object at 0x7f97577a1790>, <wsimod.nodes.demand.ResidentialDemand object at 0x7f975778bd90>, <wsimod.nodes.nodes.Node object at 0x7f975773fd50>, <wsimod.nodes.storage.Reservoir object at 0x7f975773df90>, <wsimod.nodes.wtw.FWTW object at 0x7f975804fd90>, <wsimod.nodes.wtw.WWTW object at 0x7f9757732f10>, <wsimod.nodes.sewer.Sewer object at 0x7f97577311d0>, <wsimod.nodes.land.Land object at 0x7f9757744750>, <wsimod.nodes.storage.Groundwater object at 0x7f97578e7110>, <wsimod.nodes.nodes.Node object at 0x7f97577b4890>, <wsimod.nodes.nodes.Node object at 0x7f97577b42d0>, <wsimod.nodes.nodes.Node object at 0x7f97577aa710>, <wsimod.nodes.nodes.Node object at 0x7f97577a9950>, <wsimod.nodes.nodes.Node object at 0x7f97577a8b90>]

Arcs¶

Arcs link nodes. An example arc is the link between a FWTW and the distribution node

In [31]:

# Standard simple arcs
fwtw_to_distribution = Arc(
    in_port=oxford_fwtw, out_port=distribution, name="fwtw_to_distribution"
)
print(fwtw_to_distribution)

# Standard simple arcs fwtw_to_distribution = Arc( in_port=oxford_fwtw, out_port=distribution, name="fwtw_to_distribution" ) print(fwtw_to_distribution)

<wsimod.arcs.arcs.Arc object at 0x7f97578e5610>

As with nodes, even though we only gave it a few parameters, the arc comes with a lot built in

In [32]:

print(dir(fwtw_to_distribution))

print(dir(fwtw_to_distribution))

['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'apply_overrides', 'arc_mass_balance', 'blend_vqip', 'capacity', 'compare_vqip', 'concentration_to_total', 'copy_vqip', 'ds_vqip', 'ds_vqip_c', 'empty_vqip', 'empty_vqip_predefined', 'end_timestep', 'extract_vqip', 'extract_vqip_c', 'flow_in', 'flow_out', 'get_excess', 'in_port', 'mass_balance', 'mass_balance_ds', 'mass_balance_in', 'mass_balance_out', 'name', 'out_port', 'preference', 'reinit', 'send_pull_check', 'send_pull_request', 'send_push_check', 'send_push_request', 'sum_vqip', 'total_to_concentration', 'v_change_vqip', 'v_change_vqip_c', 'v_distill_vqip', 'v_distill_vqip_c', 'vqip_in', 'vqip_out']

We can see that the arc links the two nodes

In [33]:

print(fwtw_to_distribution.in_port)

print(fwtw_to_distribution.in_port)

<wsimod.nodes.wtw.FWTW object at 0x7f975804fd90>

In [34]:

print(fwtw_to_distribution.out_port)

print(fwtw_to_distribution.out_port)

<wsimod.nodes.nodes.Node object at 0x7f975773fd50>

And that it has updated the nodes that it is connecting.

In [35]:

print(oxford_fwtw.out_arcs)

print(oxford_fwtw.out_arcs)

{'fwtw_to_distribution': <wsimod.arcs.arcs.Arc object at 0x7f97578e5610>}

In [36]:

print(distribution.in_arcs)

print(distribution.in_arcs)

{'fwtw_to_distribution': <wsimod.arcs.arcs.Arc object at 0x7f97578e5610>}

We use arcs to send checks and requests..

In [37]:

print(fwtw_to_distribution.send_pull_check({"volume": 20}))

print(fwtw_to_distribution.send_pull_check({"volume": 20}))

{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 15.0}

In [38]:

reply = fwtw_to_distribution.send_pull_request({"volume": 20})
print(reply)

reply = fwtw_to_distribution.send_pull_request({"volume": 20}) print(reply)

{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 15.0}

They convey this information to the nodes that they connect to, which update their state variables

In [39]:

print(oxford_fwtw.service_reservoir_tank.storage)

print(oxford_fwtw.service_reservoir_tank.storage)

{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 0.0}

In turn, the arcs update their own state variables.

In [40]:

print(fwtw_to_distribution.flow_in)
print(fwtw_to_distribution.flow_out)

print(fwtw_to_distribution.flow_in) print(fwtw_to_distribution.flow_out)

15.0
15.0

Arc parameters¶

Besides the in/out ports and names, arcs can have a parameter for their capacity, to limit the flow that may pass through it each timestep. A typical example would be on river abstractions to a reservoir

In [41]:

abstraction_to_farmoor = Arc(
    in_port=farmoor_abstraction,
    out_port=farmoor,
    name="abstraction_to_farmoor",
    capacity=5e4,
)

abstraction_to_farmoor = Arc( in_port=farmoor_abstraction, out_port=farmoor, name="abstraction_to_farmoor", capacity=5e4, )

A bit more sophisticated is the 'preference' parameter. We use preference to express where we would prefer the model to send water. In this example, the sewer can send water to both the treatment plant and directly into the river. Of course we would always to prefer to send water to the plant, so we give it a very high preference. Discharging into the river should only be done if there is no capacity left at the WWTW, so we give the arc a very low preference.

In [42]:

sewer_to_wwtw = Arc(
    in_port=combined_sewer, out_port=oxford_wwtw, preference=1e10, name="sewer_to_wwtw"
)
sewer_overflow = Arc(
    in_port=combined_sewer,
    out_port=thames_mixer,
    preference=1e-10,
    name="sewer_overflow",
)

sewer_to_wwtw = Arc( in_port=combined_sewer, out_port=oxford_wwtw, preference=1e10, name="sewer_to_wwtw" ) sewer_overflow = Arc( in_port=combined_sewer, out_port=thames_mixer, preference=1e-10, name="sewer_overflow", )

Create arcs¶

Arcs are a bit less interesting than nodes because they generally don't capture complicated physical behaviours. So we just create all of them below.

In [43]:

evenlode_to_thames = Arc(
    in_port=evenlode, out_port=evenlode_thames, name="evenlode_to_thames"
)

thames_to_thames = Arc(
    in_port=thames, out_port=evenlode_thames, name="thames_to_thames"
)

ray_to_cherwell = Arc(in_port=ray, out_port=cherwell_ray, name="ray_to_cherwell")

cherwell_to_cherwell = Arc(
    in_port=cherwell, out_port=cherwell_ray, name="cherwell_to_cherwell"
)

thames_to_farmoor = Arc(
    in_port=evenlode_thames, out_port=farmoor_abstraction, name="thames_to_farmoor"
)

farmoor_to_mixer = Arc(
    in_port=farmoor_abstraction, out_port=thames_mixer, name="farmoor_to_mixer"
)

cherwell_to_mixer = Arc(
    in_port=cherwell_ray, out_port=thames_mixer, name="cherwell_to_mixer"
)

wwtw_to_mixer = Arc(in_port=oxford_wwtw, out_port=thames_mixer, name="wwtw_to_mixer")

mixer_to_waste = Arc(
    in_port=thames_mixer, out_port=thames_above_abingdon, name="mixer_to_waste"
)

distribution_to_demand = Arc(
    in_port=distribution, out_port=oxford, name="distribution_to_demand"
)

reservoir_to_fwtw = Arc(in_port=farmoor, out_port=oxford_fwtw, name="reservoir_to_fwtw")

fwtw_to_sewer = Arc(in_port=oxford_fwtw, out_port=combined_sewer, name="fwtw_to_sewer")

demand_to_sewer = Arc(in_port=oxford, out_port=combined_sewer, name="demand_to_sewer")

land_to_sewer = Arc(in_port=oxford_land, out_port=combined_sewer, name="land_to_sewer")

land_to_gw = Arc(in_port=oxford_land, out_port=gw, name="land_to_gw")

garden_to_gw = Arc(in_port=oxford, out_port=gw, name="garden_to_gw")

gw_to_mixer = Arc(in_port=gw, out_port=thames_mixer, name="gw_to_mixer")

evenlode_to_thames = Arc( in_port=evenlode, out_port=evenlode_thames, name="evenlode_to_thames" ) thames_to_thames = Arc( in_port=thames, out_port=evenlode_thames, name="thames_to_thames" ) ray_to_cherwell = Arc(in_port=ray, out_port=cherwell_ray, name="ray_to_cherwell") cherwell_to_cherwell = Arc( in_port=cherwell, out_port=cherwell_ray, name="cherwell_to_cherwell" ) thames_to_farmoor = Arc( in_port=evenlode_thames, out_port=farmoor_abstraction, name="thames_to_farmoor" ) farmoor_to_mixer = Arc( in_port=farmoor_abstraction, out_port=thames_mixer, name="farmoor_to_mixer" ) cherwell_to_mixer = Arc( in_port=cherwell_ray, out_port=thames_mixer, name="cherwell_to_mixer" ) wwtw_to_mixer = Arc(in_port=oxford_wwtw, out_port=thames_mixer, name="wwtw_to_mixer") mixer_to_waste = Arc( in_port=thames_mixer, out_port=thames_above_abingdon, name="mixer_to_waste" ) distribution_to_demand = Arc( in_port=distribution, out_port=oxford, name="distribution_to_demand" ) reservoir_to_fwtw = Arc(in_port=farmoor, out_port=oxford_fwtw, name="reservoir_to_fwtw") fwtw_to_sewer = Arc(in_port=oxford_fwtw, out_port=combined_sewer, name="fwtw_to_sewer") demand_to_sewer = Arc(in_port=oxford, out_port=combined_sewer, name="demand_to_sewer") land_to_sewer = Arc(in_port=oxford_land, out_port=combined_sewer, name="land_to_sewer") land_to_gw = Arc(in_port=oxford_land, out_port=gw, name="land_to_gw") garden_to_gw = Arc(in_port=oxford, out_port=gw, name="garden_to_gw") gw_to_mixer = Arc(in_port=gw, out_port=thames_mixer, name="gw_to_mixer")

Again, we keep all the arcs in a tidy list together.

In [44]:

arclist = [
    evenlode_to_thames,
    thames_to_thames,
    ray_to_cherwell,
    cherwell_to_cherwell,
    thames_to_farmoor,
    farmoor_to_mixer,
    cherwell_to_mixer,
    wwtw_to_mixer,
    sewer_overflow,
    mixer_to_waste,
    abstraction_to_farmoor,
    distribution_to_demand,
    demand_to_sewer,
    land_to_sewer,
    sewer_to_wwtw,
    fwtw_to_sewer,
    fwtw_to_distribution,
    reservoir_to_fwtw,
    land_to_gw,
    garden_to_gw,
    gw_to_mixer,
]

arclist = [ evenlode_to_thames, thames_to_thames, ray_to_cherwell, cherwell_to_cherwell, thames_to_farmoor, farmoor_to_mixer, cherwell_to_mixer, wwtw_to_mixer, sewer_overflow, mixer_to_waste, abstraction_to_farmoor, distribution_to_demand, demand_to_sewer, land_to_sewer, sewer_to_wwtw, fwtw_to_sewer, fwtw_to_distribution, reservoir_to_fwtw, land_to_gw, garden_to_gw, gw_to_mixer, ]

Mapping¶

Remember, WSIMOD is an integrated model. Because it covers so many different things, it is very easy to make mistakes. Thus it is always good practice to plot your data!

Below we load the node location data and create arcs from the information in the arclist.

In [45]:

location_fn = os.path.join(data_folder, "raw", "points_locations.geojson")
nodes_gdf = gpd.read_file(location_fn).set_index("name")
arcs_gdf = []

for arc in arclist:
    arcs_gdf.append(
        {
            "name": arc.name,
            "geometry": LineString(
                [
                    nodes_gdf.loc[arc.in_port.name, "geometry"],
                    nodes_gdf.loc[arc.out_port.name, "geometry"],
                ]
            ),
        }
    )

arcs_gdf = gpd.GeoDataFrame(arcs_gdf, crs=nodes_gdf.crs)

location_fn = os.path.join(data_folder, "raw", "points_locations.geojson") nodes_gdf = gpd.read_file(location_fn).set_index("name") arcs_gdf = [] for arc in arclist: arcs_gdf.append( { "name": arc.name, "geometry": LineString( [ nodes_gdf.loc[arc.in_port.name, "geometry"], nodes_gdf.loc[arc.out_port.name, "geometry"], ] ), } ) arcs_gdf = gpd.GeoDataFrame(arcs_gdf, crs=nodes_gdf.crs)

Because we converted the information as GeoDataFrames, we can simply plot them below

In [46]:

f, ax = plt.subplots()
arcs_gdf.plot(ax=ax)
nodes_gdf.plot(color="r", ax=ax, zorder=10)

f, ax = plt.subplots() arcs_gdf.plot(ax=ax) nodes_gdf.plot(color="r", ax=ax, zorder=10)

Out[46]:

<Axes: >

Orchestration¶

Orchestration is making the simulation happen by calling functions in the nodes. These functions simulate physical behaviour within the node, and cause pulls/pushes to happen which in turn triggers physical behaviour in other nodes.

Orchestrating an individual timestep¶

We will start below by manually orchestrating a single timestep.

We start by setting the date, so that every node knows what forcing data to read for this timestep.

In [47]:

date = dates[0]

for node in nodelist:
    node.t = date

print(date)
print(oxford_fwtw.t)

date = dates[0] for node in nodelist: node.t = date print(date) print(oxford_fwtw.t)

2009-03-03 00:00:00
2009-03-03 00:00:00

We can see the service reservoirs are empty but the supply reservoir is not!

In [48]:

print(oxford_fwtw.service_reservoir_tank.storage)
print(farmoor.tank.storage)

print(oxford_fwtw.service_reservoir_tank.storage) print(farmoor.tank.storage)

{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 0.0}
{'temperature': 0, 'silicon': 0, 'fluoride': 0, 'chloride': 0, 'nitrate': 0, 'sulphate': 0, 'nitrogen': 0, 'sodium': 0, 'potassium': 0, 'calcium': 0, 'magnesium': 0, 'boron': 0, 'volume': 10000000.0}

If we call the FWTW's treat_water function it will pull water from the supply reservoir and update its service reservoirs

In [49]:

oxford_fwtw.treat_water()

print(oxford_fwtw.service_reservoir_tank.storage)
print(farmoor.tank.storage)

oxford_fwtw.treat_water() print(oxford_fwtw.service_reservoir_tank.storage) print(farmoor.tank.storage)

{'temperature': 0.0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 43641.0}
{'temperature': 0, 'silicon': 0.0, 'fluoride': 0.0, 'chloride': 0.0, 'nitrate': 0.0, 'sulphate': 0.0, 'nitrogen': 0.0, 'sodium': 0.0, 'potassium': 0.0, 'calcium': 0.0, 'magnesium': 0.0, 'boron': 0.0, 'volume': 9955000.0}

This information is tracked in the arcs that enter the FWTW

In [50]:

print(oxford_fwtw.in_arcs)

print(oxford_fwtw.in_arcs)

{'reservoir_to_fwtw': <wsimod.arcs.arcs.Arc object at 0x7f97578d1450>}

In [51]:

print(reservoir_to_fwtw.flow_in)
print(reservoir_to_fwtw.flow_out)

print(reservoir_to_fwtw.flow_in) print(reservoir_to_fwtw.flow_out)

45000.0
45000.0

Although none of that water has yet entered the distribution network (only some small flow from the earlier demonstration)

In [52]:

print(fwtw_to_distribution.flow_in)

print(fwtw_to_distribution.flow_in)

15.0

That is because no water consumption demand had yet been generated.

If we call the demand node's create_demand function we see that the distribution arc becomes utilised.

In [53]:

oxford.create_demand()
print(fwtw_to_distribution.flow_in)

oxford.create_demand() print(fwtw_to_distribution.flow_in)

30015.0

We also see that this gets pushed onwards into the sewer system

In [54]:

print(demand_to_sewer.flow_in)

print(demand_to_sewer.flow_in)

30000.0

Many nodes have functions intended to be called during orchestration. These functions are described in the documentation. For example, we see in the Land node API reference that the 'run' function is intended to be called from orchestration.

In [55]:

oxford_land.run()

oxford_land.run()

Below we call the functions for other nodes

In [56]:

# Discharge GW
gw.distribute()

# Discharge sewers (pushed to other sewers or WWTW)
combined_sewer.make_discharge()

# Run WWTW model
oxford_wwtw.calculate_discharge()

# Make abstractions
farmoor.make_abstractions()

# Discharge WW
oxford_wwtw.make_discharge()

# Route catchments
evenlode.route()
thames.route()
ray.route()
cherwell.route()

# Discharge GW gw.distribute() # Discharge sewers (pushed to other sewers or WWTW) combined_sewer.make_discharge() # Run WWTW model oxford_wwtw.calculate_discharge() # Make abstractions farmoor.make_abstractions() # Discharge WW oxford_wwtw.make_discharge() # Route catchments evenlode.route() thames.route() ray.route() cherwell.route()

Ending a timestep¶

Because mistakes happen, it is essential to carry out mass balance testing. Each node has a mass balance function that can be called. We see a mass balance violation resulting from the demonstration with the FWTW earlier.

In [57]:

for node in nodelist:
    in_, ds_, out_ = node.node_mass_balance()

for node in nodelist: in_, ds_, out_ = node.node_mass_balance()

mass balance error for volume of 15.0 in oxford_distribution
mass balance error for volume of -15.0 in oxford_fwtw

We should also call the end_timestep function in nodes and arcs. This is important for mass balance testing and capturing the behaviour of some dynamic processes in nodes.

In [58]:

for node in nodelist:
    node.end_timestep()

for arc in arclist:
    arc.end_timestep()

for node in nodelist: node.end_timestep() for arc in arclist: arc.end_timestep()

Model object¶

Of course it would be a massive pain to manually orchestrate every timestep. So instead we store node and arc information in a model object that will do the orchestration for us.

Because we have already created the nodes/arcs above, we simply need to add the instantiated lists above.

In [59]:

my_model = Model()
my_model.add_instantiated_nodes(nodelist)
my_model.add_instantiated_arcs(arclist)
my_model.dates = dates

my_model = Model() my_model.add_instantiated_nodes(nodelist) my_model.add_instantiated_arcs(arclist) my_model.dates = dates

The model object lets us reinitialise the nodes/arcs, and run all of the orchestration with a 'run' function.

In [60]:

my_model.reinit()

flows, _, _, _ = my_model.run()

my_model.reinit() flows, _, _, _ = my_model.run()

  0%|          | 0/1456 [00:00<?, ?it/s]

  2%|▏         | 34/1456 [00:00<00:04, 334.21it/s]

  5%|▍         | 68/1456 [00:00<00:04, 332.50it/s]

  7%|▋         | 102/1456 [00:00<00:04, 328.87it/s]

  9%|▉         | 136/1456 [00:00<00:03, 330.27it/s]

 12%|█▏        | 170/1456 [00:00<00:03, 329.55it/s]

 14%|█▍        | 203/1456 [00:00<00:03, 328.25it/s]

 16%|█▌        | 236/1456 [00:00<00:03, 327.16it/s]

 18%|█▊        | 269/1456 [00:00<00:03, 325.20it/s]

 21%|██        | 302/1456 [00:00<00:03, 324.05it/s]

 23%|██▎       | 335/1456 [00:01<00:03, 323.37it/s]

 25%|██▌       | 368/1456 [00:01<00:03, 321.71it/s]

 28%|██▊       | 401/1456 [00:01<00:03, 319.41it/s]

 30%|██▉       | 433/1456 [00:01<00:03, 319.06it/s]

 32%|███▏      | 465/1456 [00:01<00:03, 316.97it/s]

 34%|███▍      | 498/1456 [00:01<00:02, 320.00it/s]

 36%|███▋      | 531/1456 [00:01<00:02, 320.95it/s]

 39%|███▊      | 564/1456 [00:01<00:02, 321.60it/s]

 41%|████      | 597/1456 [00:01<00:02, 322.85it/s]

 43%|████▎     | 630/1456 [00:01<00:02, 320.32it/s]

 46%|████▌     | 663/1456 [00:02<00:02, 322.03it/s]

 48%|████▊     | 696/1456 [00:02<00:02, 321.77it/s]

 50%|█████     | 729/1456 [00:02<00:02, 322.05it/s]

 52%|█████▏    | 762/1456 [00:02<00:02, 322.78it/s]

 55%|█████▍    | 795/1456 [00:02<00:02, 323.48it/s]

 57%|█████▋    | 828/1456 [00:02<00:01, 324.25it/s]

 59%|█████▉    | 861/1456 [00:02<00:01, 324.29it/s]

 61%|██████▏   | 894/1456 [00:02<00:01, 322.93it/s]

 64%|██████▎   | 927/1456 [00:02<00:01, 323.24it/s]

 66%|██████▌   | 960/1456 [00:02<00:01, 324.11it/s]

 68%|██████▊   | 993/1456 [00:03<00:01, 324.69it/s]

 70%|███████   | 1026/1456 [00:03<00:01, 324.24it/s]

 73%|███████▎  | 1059/1456 [00:03<00:01, 322.12it/s]

 75%|███████▌  | 1092/1456 [00:03<00:01, 323.63it/s]

 77%|███████▋  | 1125/1456 [00:03<00:01, 323.64it/s]

 80%|███████▉  | 1158/1456 [00:03<00:00, 320.70it/s]

 82%|████████▏ | 1191/1456 [00:03<00:00, 321.03it/s]

 84%|████████▍ | 1224/1456 [00:03<00:00, 321.01it/s]

 86%|████████▋ | 1257/1456 [00:03<00:00, 321.43it/s]

 89%|████████▊ | 1290/1456 [00:03<00:00, 322.19it/s]

 91%|█████████ | 1323/1456 [00:04<00:00, 322.61it/s]

 93%|█████████▎| 1356/1456 [00:04<00:00, 322.46it/s]

 95%|█████████▌| 1389/1456 [00:04<00:00, 322.56it/s]

 98%|█████████▊| 1422/1456 [00:04<00:00, 320.74it/s]

100%|█████████▉| 1455/1456 [00:04<00:00, 321.08it/s]

100%|██████████| 1456/1456 [00:04<00:00, 322.86it/s]

The model outputs flows as a dictionary which can be converted to a dataframe

In [61]:

flows = pd.DataFrame(flows)

print(flows.sample(10))

flows = pd.DataFrame(flows) print(flows.sample(10))

                        arc           flow       time  temperature  \
22454      farmoor_to_mixer  831337.786760 2012-02-05     3.322159   
4411       thames_to_thames  145152.000000 2009-09-29    18.200000   
17322            land_to_gw    1537.626321 2011-06-05     7.448602   
12008     reservoir_to_fwtw   30934.213240 2010-09-25    11.406995   
4307        ray_to_cherwell   11232.000000 2009-09-24    17.142857   
28488       demand_to_sewer   30000.000000 2012-11-18    12.202857   
5781      cherwell_to_mixer  946080.000000 2009-12-03     7.150372   
19975     thames_to_farmoor  121305.600000 2011-10-10    15.075275   
22851  cherwell_to_cherwell  158112.000000 2012-02-24     7.250000   
2054      reservoir_to_fwtw   30934.213240 2009-06-08     3.355968   

           silicon    fluoride      chloride      nitrate      sulphate  \
22454  4732.286844  125.928754  31717.696425  5552.174444  44575.165921   
4411   2236.502016   14.515200   6314.112000   819.142212   8418.816000   
17322     0.014657    0.041876      0.012563     0.004188     14.656652   
12008   155.848387    3.373687    812.274922   147.329528   1232.149042   
4307     98.549568    1.925486   1143.257143    82.945114   1107.956571   
28488   600.000000   12.000000   5400.000000  1800.000000   7500.000000   
5781   8276.341166  222.760183  32238.000000  7911.726849  65255.162606   
19975   877.184660   17.797536   5896.717426   572.873570   8286.039195   
22851   724.089715   33.203520  12490.848000   993.464451  13344.652800   
2054     29.794381    0.819281    236.483733    51.542925    357.986986   

          nitrogen        sodium    potassium       calcium    magnesium  \
22454  6257.871985  18937.103659  3753.586498  89938.503372  4091.512337   
4411    856.396800   4325.529600   812.851200  14747.443200   798.336000   
17322     0.837523      0.062814     1.465665     14.656652     1.256284   
12008   154.740755    513.069458   104.715785   2562.969352   121.451412   
4307     92.744229    800.039314   180.032914   1111.326171    69.157029   
28488  1500.000000   6000.000000   900.000000   4500.000000   900.000000   
5781   7773.393189  20698.823314  5695.734857  96635.314286  6062.613943   
19975   609.008544   4352.406171   873.791589  12165.649509   678.122743   
22851  1079.114400   7170.379200   972.388800  17336.980800  1225.368000   
2054     53.508436    152.736044    30.240173    870.617651    40.692423   

           boron  
22454  37.804009  
4411    8.563968  
17322   0.020938  
12008   1.247807  
4307    1.612594  
28488  15.000000  
5781   67.537152  
19975   7.665652  
22851  11.265480  
2054    0.374769  

Validation plots¶

Of course we shouldn't trust a model without proof, and this is doubly true for an integrated model whose errors may easily propagate.

Thankfully we have rich in-river water quality sampling in Oxford that we can use for validation

We load and format this data for dates and pollutants that overlap with what we have simulated

In [62]:

mixer_val_df = pd.read_csv(os.path.join(data_folder, "raw", "mixer_wims_val.csv"))
mixer_val_df.date = pd.to_datetime(mixer_val_df.date)
mixer_val_df = mixer_val_df.loc[mixer_val_df.date.isin(dates)]
val_pollutants = mixer_val_df.variable.unique()
mixer_val_df = mixer_val_df.pivot(index="date", columns="variable", values="value")

mixer_val_df = pd.read_csv(os.path.join(data_folder, "raw", "mixer_wims_val.csv")) mixer_val_df.date = pd.to_datetime(mixer_val_df.date) mixer_val_df = mixer_val_df.loc[mixer_val_df.date.isin(dates)] val_pollutants = mixer_val_df.variable.unique() mixer_val_df = mixer_val_df.pivot(index="date", columns="variable", values="value")

We get rid of the first day of flows because our tanks were initialised empty and will not be informative

In [63]:

flows = flows.loc[flows.time != dates[0]]

flows = flows.loc[flows.time != dates[0]]

We convert our flows, which are simulated in kg/d into a concentration value (kg/m3/d).

In [64]:

flows_plot = flows.copy()

# Convert to KG/M3
for pol in set(val_pollutants):
    if pol != "temperature":
        flows_plot[pol] /= flows_plot.flow

flows_plot = flows.copy() # Convert to KG/M3 for pol in set(val_pollutants): if pol != "temperature": flows_plot[pol] /= flows_plot.flow

We pick the model outlet as a validation location and make a pretty plot - fantastic!

In [65]:

plot_arc = "mixer_to_waste"
f, axs = plt.subplots(val_pollutants.size, 1)
for pol, ax in zip(val_pollutants, axs):
    ax.plot(
        flows_plot.loc[flows_plot.arc == plot_arc, [pol, "time"]].set_index("time"),
        color="b",
        label="simulation",
    )
    ax.plot(mixer_val_df[pol], ls="", marker="o", color="r", label="spot sample")
    ax.set_ylabel(pol)
plt.legend()

plot_arc = "mixer_to_waste" f, axs = plt.subplots(val_pollutants.size, 1) for pol, ax in zip(val_pollutants, axs): ax.plot( flows_plot.loc[flows_plot.arc == plot_arc, [pol, "time"]].set_index("time"), color="b", label="simulation", ) ax.plot(mixer_val_df[pol], ls="", marker="o", color="r", label="spot sample") ax.set_ylabel(pol) plt.legend()

Out[65]:

<matplotlib.legend.Legend at 0x7f9757b8eb90>