Lumped Parameter Model for Lake Dynamics

In [91]:
# Display graphics inline with the notebook
%matplotlib inline

# Standard Python modules
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import pandas as pd
import os
import datetime

# Module to enhance matplotlib plotting
import seaborn
seaborn.set()

# Modules to display images and data tables
from IPython.display import Image
from IPython.core.display import display

# Data Directory
dir = './data/'

# Styles
from IPython.core.display import HTML
HTML(open("styles/custom.css", "r").read())
Out[91]:

Stage-Volume Relationships

$$ V(h) = a_0 + a_1 h + a_2 h^2 $$$$\begin{align*} A(h) & = \frac{d V}{d h} \\ & = a_1 + 2 a_2 h \\ & = \underbrace{\left(a_1 + 2 a_2 h_0\right)}_{b_0} \; + \; \underbrace{2 a_2}_{b_1} \left(h-h_0\right)\end{align*}$$

Rainy Lake

In [96]:
h = np.array([335.0, 336.0, 336.5, 337.0, 337.5, 338.0, 339.0, 340.0])
v = np.array([112.67, 798.00, 1176.42, 1577.25, 2002.06, 2450.57, 3416.85, 4458.97])

plt.subplot(2,1,1)
plt.scatter(h,v)
plt.xlim(h.min(),h.max())
plt.ylim(0,plt.ylim()[1])
plt.title('Rainy Lake Stage-Volume Relationship')
plt.xlabel('Water Elevation (m)')
plt.ylabel('Volume (million cu. m)')

pf = sp.polyfit(h,v,2)

VRL= sp.poly1d(pf)
plt.hold(True)
plt.plot(h,VRL(h))
plt.hold(False)

ARL = sp.poly1d(np.array([2.0*pf[0],pf[1]]))
plt.subplot(2,1,2)
plt.plot(h,ARL(h))
plt.title('Rainy Lake Stage-Surface Area Relationship')
plt.xlabel('Elevation (m)')
plt.ylabel('Area (sq. km.)')

plt.tight_layout()

df = pd.DataFrame(zip(h,v,VRL(h),ARL(h)),columns=['h','V','Vhat','Ahat'])
print df.to_string(formatters={'V':'  {:.0f}'.format, 
                                'Vhat':'  {:.0f}'.format,
                                'Ahat':'  {:.0f}'.format})
       h      V   Vhat   Ahat
0  335.0    113    111    644
1  336.0    798    799    734
2  336.5   1176   1178    780
3  337.0   1577   1579    825
4  337.5   2002   2003    870
5  338.0   2451   2450    916
6  339.0   3417   3411   1007
7  340.0   4459   4463   1097
In [98]:
h = np.array([337.0, 338.0, 338.5, 339.0, 339.5, 340.0,
     340.5, 341.0, 341.5, 342.0, 343.0])
v = np.array([65.33, 259.95, 364.20, 475.58, 592.46, 712.28,
     836.56, 966.17, 1099.79, 1239.68, 1540.75])

plt.subplot(2,1,1)
plt.scatter(h,v)
plt.xlim(h.min(),h.max())
plt.ylim(0,plt.ylim()[1])
plt.title('Namakan Lake Stage-Volume Relationship')
plt.xlabel('Water Elevation (m)')
plt.ylabel('Volume (million cu. m)')

pf = sp.polyfit(h,v,2)

VNL= sp.poly1d(pf)
plt.hold(True)
plt.plot(h,VNL(h))
plt.hold(False)

ANL = sp.poly1d(np.array([2.0*pf[0],pf[1]]))
plt.subplot(2,1,2)
plt.plot(h,ANL(h))
plt.title('Namakan Lake Stage-Surface Area Relationship')
plt.xlabel('Elevation (m)')
plt.ylabel('Area (sq. km.)')

plt.tight_layout()

df = pd.DataFrame(zip(h,v,VNL(h),ANL(h)),columns=['h','V','Vhat','Ahat'])
print df.to_string(formatters={'V':'  {:.0f}'.format, 
                                'Vhat':'  {:.0f}'.format,
                                'Ahat':'  {:.0f}'.format})
        h      V   Vhat  Ahat
0   337.0     65     66   185
1   338.0    260    260   205
2   338.5    364    365   215
3   339.0    476    475   225
4   339.5    592    591   235
5   340.0    712    711   246
6   340.5    837    836   256
7   341.0    966    966   266
8   341.5   1100   1102   276
9   342.0   1240   1242   286
10  343.0   1541   1539   306

Stage-Discharge Relationships

Rainy Lake

In [108]:
h = np.array([335.4, 336.0, 336.5, 336.75, 337.0, 337.25, 337.5, 
              337.75, 338.0, 338.5, 339, 339.5, 340.0])
d = np.array([0.0, 399., 425., 443., 589., 704., 792., 909.,
              1014., 1156., 1324., 1550., 1778.])

#plt.subplot(2,1,1)

plt.hold(True)
plt.scatter(d,h)
plt.plot(d,h)

plt.ylim(h.min(),h.max())
plt.xlim(0,plt.xlim()[1])
plt.title('Rainy Lake Stage-Discharge Relationship')
plt.ylabel('Water Elevation (m)')
plt.xlabel('Discharge (cu. per sec.)')


# Get historical flowrates and levels on RR and RL
RR = pd.read_pickle(dir+'RR.pkl')['1970':]
RL = pd.read_pickle(dir+'RL.pkl')['1970':]

Q = pd.concat([RR,RL],axis=1)
Q.columns = ['RR','RL']

A = Q['1970':'2000']

plt.scatter(A['RR'],A['RL'],marker='+',color='b')
B = Q['2000':]
plt.scatter(B['RR'],B['RL'],marker='+',color='r')

ax = plt.axis()
plt.plot([ax[0],ax[1]],[336.70,336.70],'y--')
plt.plot([ax[0],ax[1]],[337.20,337.20],'y--')
plt.plot([ax[0],ax[1]],[337.75,337.75],'m--')
plt.plot([ax[0],ax[1]],[337.90,337.90],'c--')

plt.xlabel('Upper Rainy River Flow [cubic meters/sec]')
plt.ylabel('Rainy Lake Level [meters]')
ax = plt.axis()
plt.axis([0,ax[1],ax[2],ax[3]])
plt.legend(['Winter Drought Line','Summer Drought Line','URC_max','All Gates Open','1970-2000','2000-2010'],loc="upper left")

plt.hold(False)
#plt.savefig('./images/RainyRiverDischarge.png')
In [112]:
Q = lambda x: np.interp(x,h,d)
x = np.linspace(336.,339.)
plt.plot(x,Q(x))
Out[112]:
[<matplotlib.lines.Line2D at 0x10f27e4d0>]
In [ ]: