Land nodes - hydrology and agriculture (.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. Basic surface

4. Pervious surface

5. Connecting land nodes in a model

   5.1 Model object

6. Growing surface

Introduction¶

A land node is the object that interfaces with other WSIMOD nodes. Each land node can have one or more surfaces. Each surface can be parameterised differently or follow different equations, similar to the hyrological response unit concept in hydrological models.

The different surfaces within WSIMOD are used to capture hydrological, agricultural, and urban runoff processes. In this demo we will see how they all work.

Imports and forcing data¶

Import packages

In [1]:

import os

import pandas as pd
from matplotlib import pyplot as plt

from wsimod.arcs.arcs import Arc
from wsimod.core import constants
from wsimod.demo.create_oxford import create_timeseries
from wsimod.nodes.land import Land
from wsimod.nodes.nodes import Node
from wsimod.nodes.storage import Groundwater
from wsimod.nodes.waste import Waste
from wsimod.orchestration.model import Model

import os import pandas as pd from matplotlib import pyplot as plt from wsimod.arcs.arcs import Arc from wsimod.core import constants from wsimod.demo.create_oxford import create_timeseries from wsimod.nodes.land import Land from wsimod.nodes.nodes import Node from wsimod.nodes.storage import Groundwater from wsimod.nodes.waste import Waste from wsimod.orchestration.model import Model

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()
)

dates = input_data.date.drop_duplicates()
dates_monthyear = input_data.date.dt.to_period("M").unique()
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() ) dates = input_data.date.drop_duplicates() dates_monthyear = input_data.date.dt.to_period("M").unique() print(input_data.sample(10))

           site       date     variable      value
29618  cherwell 2010-07-14     sulphate   0.073363
1496   evenlode 2009-04-12  temperature  11.625000
54549  evenlode 2011-01-09      calcium   0.099714
61624       ray 2010-06-18    magnesium   0.005529
49015  evenlode 2011-10-26    potassium   0.005743
3921        ray 2011-12-07  temperature   4.885714
19983  evenlode 2012-01-22     chloride   0.029930
62294       ray 2012-04-18    magnesium   0.006214
19396  evenlode 2010-06-14     chloride   0.025746
55664       ray 2010-02-02      calcium   0.115700

C:\Users\bdobson\AppData\Local\Temp\ipykernel_2148\1905528647.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]:

land_inputs = data_input_dict["oxford_land"]

example_date = pd.to_datetime("2009-03-03")

print(land_inputs[("precipitation", example_date)])
print(land_inputs[("et0", example_date)])
print(land_inputs[("temperature", example_date)])

land_inputs = data_input_dict["oxford_land"] example_date = pd.to_datetime("2009-03-03") print(land_inputs[("precipitation", example_date)]) print(land_inputs[("et0", example_date)]) print(land_inputs[("temperature", example_date)])

0.015
0.002
7.625

We are just using some basic pollutants for demonstration

In [4]:

constants.set_simple_pollutants()

print(constants.POLLUTANTS)

constants.set_simple_pollutants() print(constants.POLLUTANTS)

['phosphate', 'temperature']

Basic surface¶

Create a simple Land node with one basic surface. We can pass surfaces as a dictionary when creating a Land node.

In [5]:

surface = {
    "type_": "Surface",
    "surface": "my_surface",
    "area": 10,
    "depth": 1,
    "pollutant_load": {"phosphate": 1e-7},
}

land = Land(name="my_land", data_input_dict=land_inputs, surfaces=[surface])

surface = { "type_": "Surface", "surface": "my_surface", "area": 10, "depth": 1, "pollutant_load": {"phosphate": 1e-7}, } land = Land(name="my_land", data_input_dict=land_inputs, surfaces=[surface])

Each surface object is stored in a list of surfaces within the land node.

In [6]:

print(land.surfaces[0])

print(land.surfaces[0])

<wsimod.nodes.land.Surface object at 0x000001C59BE0DA50>

We can also access a specific surface by name.

In [7]:

print(land.get_surface("my_surface"))

print(land.get_surface("my_surface"))

<wsimod.nodes.land.Surface object at 0x000001C59BE0DA50>

This surface will have all the data we passed to the land node in a dict.

In [8]:

print(land.get_surface("my_surface").area)
print(land.get_surface("my_surface").pollutant_load)

print(land.get_surface("my_surface").area) print(land.get_surface("my_surface").pollutant_load)

10
{'phosphate': 1e-07}

A surface is a generic tank, by default it is initialised as empty.

In [9]:

print(land.get_surface("my_surface").storage)

print(land.get_surface("my_surface").storage)

{'phosphate': 0, 'temperature': 0, 'volume': 0}

Let's try to run a timestep for a land node

In [10]:

land.t = example_date
land.monthyear = land.t.to_period("M")
land.run()

land.t = example_date land.monthyear = land.t.to_period("M") land.run()

We see that a day of phosphate deposition has occured on the surface.

In [11]:

print(land.get_surface("my_surface").storage)

print(land.get_surface("my_surface").storage)

{'phosphate': 1e-06, 'temperature': 0, 'volume': 0}

The run function in land simply calls the run function in all of its surfaces.

In [12]:

land.get_surface("my_surface").run()

land.get_surface("my_surface").run()

After running the surface, we see that another day of phosphate deposition has occured.

In [13]:

print(land.get_surface("my_surface").storage)

print(land.get_surface("my_surface").storage)

{'phosphate': 2e-06, 'temperature': 0, 'volume': 0}

However... no precipitation has occurred (no volume above), despite there being rainfall!

In [14]:

print(land.get_data_input("precipitation"))

print(land.get_data_input("precipitation"))

0.015

To understand why this is, we have to understand how the 'run' function in a surface works. Each surface has a list of functions stored in 'inflows', in 'processes', and in 'outflows'. When 'run' is called, all of the functions in inflows are executed, then processes, then outflows.

If we print the lists of functions, we see only some simple deposition occurs on the basic surface in the inflows...

In [15]:

print(land.get_surface("my_surface").inflows)
print(land.get_surface("my_surface").processes)
print(land.get_surface("my_surface").outflows)

print(land.get_surface("my_surface").inflows) print(land.get_surface("my_surface").processes) print(land.get_surface("my_surface").outflows)

[<bound method Surface.simple_deposition of <wsimod.nodes.land.Surface object at 0x000001C59BE0DA50>>]
[]
[]

We need a specific subclass of surface to implement more sophisticated hydrological processes!

Pervious surface¶

The pervious surface contains a simple lumped hydrological model.

We can create the surface, again via parameters stored in a dictionary that is passed to the Land node. is passed to the Land node. Note that the subclass of surface to be created must be specified by the type_ keyword.

In [16]:

surface = {
    "type_": "PerviousSurface",
    "surface": "my_surface",
    "area": 10,
    "depth": 0.5,
    "pollutant_load": {"phosphate": 1e-7},
    "wilting_point": 0.05,
    "field_capacity": 0.1,
}

land = Land(name="my_land", data_input_dict=land_inputs, surfaces=[surface])

surface = { "type_": "PerviousSurface", "surface": "my_surface", "area": 10, "depth": 0.5, "pollutant_load": {"phosphate": 1e-7}, "wilting_point": 0.05, "field_capacity": 0.1, } land = Land(name="my_land", data_input_dict=land_inputs, surfaces=[surface])

We have lots of functions now! You will have to look at the documentation of PerviousSurface to understand them in detail

We note that the function 'ihacres' in the inflows is the hydrological processes representation, using equations based off the IHACRES model. If a user preferred to use a different hydrological model, they would simply need to substitute out this function in the inflows list.

In [17]:

print(land.get_surface("my_surface").inflows)
print(land.get_surface("my_surface").processes)
print(land.get_surface("my_surface").outflows)

print(land.get_surface("my_surface").inflows) print(land.get_surface("my_surface").processes) print(land.get_surface("my_surface").outflows)

[<bound method Surface.simple_deposition of <wsimod.nodes.land.PerviousSurface object at 0x000001C59BDE28D0>>, <bound method PerviousSurface.ihacres of <wsimod.nodes.land.PerviousSurface object at 0x000001C59BDE28D0>>]
[<bound method PerviousSurface.calculate_soil_temperature of <wsimod.nodes.land.PerviousSurface object at 0x000001C59BDE28D0>>]
[<bound method PerviousSurface.route of <wsimod.nodes.land.PerviousSurface object at 0x000001C59BDE28D0>>]

Now when we run the model, we can see that some rain has happened, because the storage volume has increased (as well as the deposition).

In [18]:

land.t = example_date
land.monthyear = land.t.to_period("M")
land.run()
print(land.get_surface("my_surface").storage)

land.t = example_date land.monthyear = land.t.to_period("M") land.run() print(land.get_surface("my_surface").storage)

{'phosphate': 1e-06, 'temperature': 7.921875000000001, 'volume': 0.14248522706924713}

However, if we look at the land node, which is the parent of the surfaces that interfaces with other WSIMOD components, we see that all of its tanks are still empty.

In [19]:

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

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

{'phosphate': 0.0, 'temperature': 0, 'volume': 0.0}
{'phosphate': 0.0, 'temperature': 0, 'volume': 0.0}
{'phosphate': 0, 'temperature': 0, 'volume': 0.0}

This is because the IHACRES model requires soil moisture to be more than a specified amount (the field capacity) for flows to be generated.

We can run the surface multiple times to fill it up with water.

In [20]:

for i in range(10):
    land.get_surface("my_surface").run()

for i in range(10): land.get_surface("my_surface").run()

There's now a lot of water in the tank

In [21]:

print(land.get_surface("my_surface").storage)

print(land.get_surface("my_surface").storage)

{'phosphate': 5.734530999564807e-06, 'temperature': 7.957258511567751, 'volume': 0.7806952872182779}

Critically, we see that the moisture content is greater than 0.1 (i.e., the field capacity moisture content)

In [22]:

print(land.get_surface("my_surface").get_smc() / land.get_surface("my_surface").depth)

print(land.get_surface("my_surface").get_smc() / land.get_surface("my_surface").depth)

0.39034764360913893

Once soil moisture content is greater than the field capacity, flows will be generated and the land tanks will fill up.

In [23]:

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

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

{'phosphate': 3.949101750326394e-06, 'temperature': 7.903529789012976, 'volume': 0.5435991598154842}
{'phosphate': 1.3163672501087978e-06, 'temperature': 7.903529789012975, 'volume': 0.18119971993849474}
{'phosphate': 0, 'temperature': 7.625, 'volume': 0.03814730946073574}

These tanks represent flow from the soil layer to either rivers or groundwater. However, land nodes expect to be able to route flow onwards to other nodes.

Since the land isn't connected to anything, these won't actually go anywhere if we run it, and they will just build up

In [24]:

for i in range(10):
    land.run()

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

for i in range(10): land.run() print(land.percolation.storage) print(land.subsurface_runoff.storage) print(land.surface_runoff.storage)

{'phosphate': 1.1091929990869039e-05, 'temperature': 7.907975528908189, 'volume': 1.5050766284155752}
{'phosphate': 3.0837386635250816e-06, 'temperature': 7.861013960531203, 'volume': 0.5016922094718587}
{'phosphate': 6.135713334312649e-07, 'temperature': 7.848067449790649, 'volume': 0.10561941252039124}

Connecting land nodes in a model¶

As mentioned, the land node expects to be able to discharge to groundwater, and rivers (where a river could be the River node or just a generic Node). We also provide a Waste node, which is just a model outlet.

In [25]:

node = Node(name="my_river")
gw = Groundwater(name="my_groundwater", area=10, capacity=100)
outlet = Waste(name="my_outlet")

node = Node(name="my_river") gw = Groundwater(name="my_groundwater", area=10, capacity=100) outlet = Waste(name="my_outlet")

We use arcs to join up all of the different nodes according to a standard hydrological representation.

In [26]:

arc1 = Arc(in_port=land, out_port=node, name="quickflow")
arc2 = Arc(in_port=land, out_port=gw, name="percolation")
arc3 = Arc(in_port=gw, out_port=node, name="baseflow")
arc4 = Arc(in_port=node, out_port=outlet, name="outflow")

arc1 = Arc(in_port=land, out_port=node, name="quickflow") arc2 = Arc(in_port=land, out_port=gw, name="percolation") arc3 = Arc(in_port=gw, out_port=node, name="baseflow") arc4 = Arc(in_port=node, out_port=outlet, name="outflow")

If we run the land a few more times, we see that these tanks start to empty (though percolation by nature empties rather slowly!!!)

In [27]:

for i in range(10):
    land.run()

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

for i in range(10): land.run() print(land.percolation.storage) print(land.subsurface_runoff.storage) print(land.surface_runoff.storage)

{'phosphate': 1.569101705155093e-05, 'temperature': 7.909762826961828, 'volume': 2.1142184152368344}
{'phosphate': 1.2131965997487093e-06, 'temperature': 7.896243556724506, 'volume': 0.1715050264585061}
{'phosphate': 0.0, 'temperature': 7.625, 'volume': 0.0}

Model object¶

We can put these nodes and arcs into the Model object to have a functioning hydrological model.

We start by creating a model object.

In [28]:

my_model = Model()

my_model = Model()

Since we have already created our nodes/arcs, we use the add_instantiated functions

In [29]:

my_model.add_instantiated_nodes([land, node, gw, outlet])
my_model.add_instantiated_arcs([arc1, arc2, arc3, arc4])

my_model.add_instantiated_nodes([land, node, gw, outlet]) my_model.add_instantiated_arcs([arc1, arc2, arc3, arc4])

Store dates

In [30]:

my_model.dates = dates

my_model.dates = dates

We have run the surfaces a few times, so will just set all of the model tanks to empty to give a clean start for the model.

In [31]:

my_model.reinit()

my_model.reinit()

We can run the model with the run function

In [32]:

results = my_model.run()

results = my_model.run()

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

 25%|██████████████████████████████████████████████████▉                                                                                                                                                       | 367/1456 [00:00<00:00, 3643.03it/s]

 50%|█████████████████████████████████████████████████████████████████████████████████████████████████████▌                                                                                                    | 732/1456 [00:00<00:00, 3530.23it/s]

 75%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉                                                   | 1086/1456 [00:00<00:00, 3500.54it/s]

 99%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍  | 1437/1456 [00:00<00:00, 3494.33it/s]

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

.. and plot the results!

In [33]:

flows = pd.DataFrame(results[0])

f, axs = plt.subplots(2, 1)
flows.groupby("arc").get_group("outflow").set_index("time").flow.plot(ax=axs[0])
flows.groupby("arc").get_group("outflow").set_index("time").phosphate.plot(ax=axs[1])

flows = pd.DataFrame(results[0]) f, axs = plt.subplots(2, 1) flows.groupby("arc").get_group("outflow").set_index("time").flow.plot(ax=axs[0]) flows.groupby("arc").get_group("outflow").set_index("time").phosphate.plot(ax=axs[1])

Out[33]:

<Axes: xlabel='time'>

Growing surface¶

Hydrology is nice, but anyone using WSIMOD probably isn't interested in hydrology only! The GrowingSurface adds a lot of sophisticated behaviour for agriculture and water quality.

Our GrowingSurface needs a bit more data than other surfaces, for fertiliser, manure and atmospheric deposition of ammonia, nitrate and phosphate. We will make up this data.

Surface pollution data varies at a monthly timestep rather than daily, though it is applied each day.

In [34]:

surface_input_data = {}
for pollutant in ["srp", "nhx", "noy"]:
    for source in ["manure", "fertiliser", "dry", "wet"]:
        amount = 1e-7  # kg/m2/timestep
        ts = create_timeseries(
            amount, dates_monthyear, "{0}-{1}".format(pollutant, source)
        )
        ts = ts.set_index(["variable", "date"]).value.to_dict()
        surface_input_data = {**surface_input_data, **ts}

print(surface_input_data[("nhx-manure", example_date.to_period("M"))])

surface_input_data = {} for pollutant in ["srp", "nhx", "noy"]: for source in ["manure", "fertiliser", "dry", "wet"]: amount = 1e-7 # kg/m2/timestep ts = create_timeseries( amount, dates_monthyear, "{0}-{1}".format(pollutant, source) ) ts = ts.set_index(["variable", "date"]).value.to_dict() surface_input_data = {**surface_input_data, **ts} print(surface_input_data[("nhx-manure", example_date.to_period("M"))])

1e-07

Because the surface represents the nitrogen and phosphorus cycles, we need to simulate a greater number of pollutants

In [35]:

constants.set_default_pollutants()

constants.set_default_pollutants()

As with the other surfaces, we create the surface by passing it as a dictionary to the created Land node. I will use the parameters for Maize for this growing surface

In [36]:

crop_factor_stages = [0.0, 0.0, 0.3, 0.3, 1.2, 1.2, 0.325, 0.0, 0.0]
crop_factor_stage_dates = [0, 90, 91, 121, 161, 213, 244, 245, 366]
sowing_day = 91
harvest_day = 244
ET_depletion_factor = 0.55
rooting_depth = 0.5

surface = {
    "type_": "GrowingSurface",
    "surface": "my_growing_surface",
    "area": 10,
    "rooting_depth": rooting_depth,
    "crop_factor_stage_dates": crop_factor_stage_dates,
    "crop_factor_stages": crop_factor_stages,
    "sowing_day": sowing_day,
    "harvest_day": harvest_day,
    "ET_depletion_factor": ET_depletion_factor,
    "data_input_dict": surface_input_data,
    "wilting_point": 0.05,
    "field_capacity": 0.1,
}

land = Land(name="my_land", data_input_dict=land_inputs, surfaces=[surface])

crop_factor_stages = [0.0, 0.0, 0.3, 0.3, 1.2, 1.2, 0.325, 0.0, 0.0] crop_factor_stage_dates = [0, 90, 91, 121, 161, 213, 244, 245, 366] sowing_day = 91 harvest_day = 244 ET_depletion_factor = 0.55 rooting_depth = 0.5 surface = { "type_": "GrowingSurface", "surface": "my_growing_surface", "area": 10, "rooting_depth": rooting_depth, "crop_factor_stage_dates": crop_factor_stage_dates, "crop_factor_stages": crop_factor_stages, "sowing_day": sowing_day, "harvest_day": harvest_day, "ET_depletion_factor": ET_depletion_factor, "data_input_dict": surface_input_data, "wilting_point": 0.05, "field_capacity": 0.1, } land = Land(name="my_land", data_input_dict=land_inputs, surfaces=[surface])

We see that the inflows includes IHACRES from the pervious surface. But has also added a range of other functions related to deposition.

In [37]:

print(land.get_surface("my_growing_surface").inflows)

print(land.get_surface("my_growing_surface").inflows)

[<bound method GrowingSurface.calc_crop_cover of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method Surface.atmospheric_deposition of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method Surface.precipitation_deposition of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method PerviousSurface.ihacres of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.effective_precipitation_flushing of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.fertiliser of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.manure of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>]

And in particular with the bio-chemical processes occurring within nutrient pools

In [38]:

print(land.get_surface("my_growing_surface").processes)

print(land.get_surface("my_growing_surface").processes)

[<bound method PerviousSurface.calculate_soil_temperature of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.calc_temperature_dependence_factor of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.calc_soil_moisture_dependence_factor of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.soil_pool_transformation of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.calc_crop_uptake of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.erosion of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.denitrification of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>, <bound method GrowingSurface.adsorption of <wsimod.nodes.land.GrowingSurface object at 0x000001C59DB6DB90>>]

Let's recreate our model with this Maize surface and run the model

In [39]:

node = Node(name="my_river")
gw = Groundwater(name="my_groundwater", area=10, capacity=100)
outlet = Waste(name="my_outlet")

arc1 = Arc(in_port=land, out_port=node, name="quickflow")
arc2 = Arc(in_port=land, out_port=gw, name="percolation")
arc3 = Arc(in_port=gw, out_port=node, name="baseflow")
arc4 = Arc(in_port=node, out_port=outlet, name="outflow")
my_model = Model()

my_model.add_instantiated_nodes([land, node, gw, outlet])
my_model.add_instantiated_arcs([arc1, arc2, arc3, arc4])

my_model.dates = dates

results = my_model.run()

node = Node(name="my_river") gw = Groundwater(name="my_groundwater", area=10, capacity=100) outlet = Waste(name="my_outlet") arc1 = Arc(in_port=land, out_port=node, name="quickflow") arc2 = Arc(in_port=land, out_port=gw, name="percolation") arc3 = Arc(in_port=gw, out_port=node, name="baseflow") arc4 = Arc(in_port=node, out_port=outlet, name="outflow") my_model = Model() my_model.add_instantiated_nodes([land, node, gw, outlet]) my_model.add_instantiated_arcs([arc1, arc2, arc3, arc4]) my_model.dates = dates results = my_model.run()

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

 13%|█████████████████████████▉                                                                                                                                                                                | 187/1456 [00:00<00:00, 1867.58it/s]

 26%|███████████████████████████████████████████████████▉                                                                                                                                                      | 374/1456 [00:00<00:00, 1829.74it/s]

 38%|█████████████████████████████████████████████████████████████████████████████▍                                                                                                                            | 558/1456 [00:00<00:00, 1804.61it/s]

 51%|██████████████████████████████████████████████████████████████████████████████████████████████████████▌                                                                                                   | 739/1456 [00:00<00:00, 1802.41it/s]

 63%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋                                                                          | 920/1456 [00:00<00:00, 1800.55it/s]

 76%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉                                                 | 1101/1456 [00:00<00:00, 1793.92it/s]

 88%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊                        | 1281/1456 [00:00<00:00, 1791.88it/s]

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

We can now plot the results

In [40]:

flows = pd.DataFrame(results[0])

f, axs = plt.subplots(2, 1)
flows.groupby("arc").get_group("outflow").set_index("time").flow.plot(ax=axs[0])
flows.groupby("arc").get_group("outflow").set_index("time").phosphate.plot(ax=axs[1])

flows = pd.DataFrame(results[0]) f, axs = plt.subplots(2, 1) flows.groupby("arc").get_group("outflow").set_index("time").flow.plot(ax=axs[0]) flows.groupby("arc").get_group("outflow").set_index("time").phosphate.plot(ax=axs[1])

Out[40]:

<Axes: xlabel='time'>

Observe the differences between the two sets of timeseries: Flows look more or less the same (dynamically), which makes sense since they both use IHACRES for hydrology. Only small differences will arise because the crops change the evapotranspiration coefficient

Meanwhile, phosphate levels look much more interesting with the GrowingSurface, and are not solely dependent on the hydrology.