Linear reservoir

This is a simple rainfall runoff model. Precipitation is randomly generated from an exponential distribution. Streamflow is governed by $\frac{dS}{dt} = -kS$, where $S$ is the storage in the reservoir.

In [1]:
import sys
In [2]:
from stockflow import simulation
import numpy as np
from matplotlib import pyplot as plt

Set up the timestep vector and sample $P$ values. t should contain all the points where you want to evaluate the state variables; it is not necessarily where the ODE solver will sample. More on that in a bit.

In [7]:
tmax = 100
dt = 1
t = np.arange(0,tmax,dt)
P = 5*np.random.standard_exponential(len(t),)

Obviously you'll want to pay attention to units here. This is just an example.

In [8]:
k = 0.3

Now create the simulation object, create stocks with their initial conditions (in a dictionary), and create flows one by one with functions describing their behavior as a function of t.

In [9]:
s = simulation(t)
s.stocks({'S': 0})

# discrete forcing - use discrete=True keyword when calling
s.flow('P', start=None, end='S', f=lambda t: P[t])
s.flow('Q', start='S', end=None, f=lambda t: k*s.S)

The functions can either be lambdas, or defined in the usual way with def my_func(t): .... If you're using discrete forcing data, note that P[t] may be given index values t that are not integers. By default NumPy rounds these down to the nearest integer. This won't matter if your forcing is continuous in t.

Now that the simulation is done, we can access all stock/flow timeseries using s.Q or s.S, for example. Let's plot streamflow.

In [10]:
plt.plot(t, s.Q)

Pretty easy. What does the output look like if we change $k$?

In [11]:
k = 0.95
In [13]:
plt.plot(t, s.Q)

You can see here the response is much more "flashy" with $k=0.95$ than with $k=0.3$ above. If we had the real $Q$ data we could compare them, calculate error, etc.