Muskingum routing in the IJssel - verification scriptΒΆ

Code adjusted from https://github.com/quaquel/epa1361_open. The output of this model is equal to my own model, as can be seen in this script.

Source: https://onlinelibrary.wiley.com/doi/abs/10.1111/jfr3.12532

Ciullo, A., de Bruijn, K. M., Kwakkel, J. H., & Klijn, F. (2019). Accounting for the uncertain effects of hydraulic interactions in optimising embankments heights: Proof of principle for the IJssel River. Journal of Flood Risk Management, e12532.

[1]:

import sys
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

[2]:

import generate_network
from functions_ijssel_muskingum import Muskingum

[3]:

G, dike_list = generate_network.get_network()

[4]:

class DikeNetwork(object):
    def __init__(self):
        # planning steps
        self.num_events = 30

        # load network
        G, dike_list = generate_network.get_network()

        self.Qpeaks = 2000 #np.random.uniform(1000,16000,100)

        self.G = G
        self.dikelist = dike_list

    def printG(self):
        print(G.nodes.data())

    def getG(self):
        return G


    def init_node(self,value, time):
        init = np.repeat(value, len(time)).tolist()
        return init

    def _initialize_hydroloads(self, node, time, Q_0):
        #node['cumVol'], node['wl'], node['Qpol'], node['hbas'] = (
        #    self.init_node(0, time) for _ in range(4))
        node['Qin'], node['Qout'] = (self.init_node(Q_0, time) for _ in range(2))
        #node['status'] = self.init_node(False, time)
        #node['tbreach'] = np.nan
        return node

    def calc_wave(self,timestep=1):
        startnode = G.node['A.0']
        waveshape_id = 0
        Qpeak = self.Qpeaks#[0]
        dikelist = self.dikelist
        time = np.arange(0, startnode['Qevents_shape'].loc[waveshape_id].shape[0],
                             timestep)
        startnode['Qout'] = Qpeak * startnode['Qevents_shape'].loc[waveshape_id]

        # Initialize hydrological event:
        for key in dikelist:
            node = G.node[key]
            #Q_0 = int(G.node['A.0']['Qout'][0])
            Q_0 = G.node['A.0']['Qout'][0]
            self._initialize_hydroloads(node, time, Q_0)

        # Run the simulation:
        # Run over the discharge wave:
        for t in range(1, len(time)):
            # Run over each node of the branch:
            for n in range(0, len(dikelist)):
                # Select current node:
                node = G.node[dikelist[n]]
                if node['type'] == 'dike':
                    # Muskingum parameters:
                    C1 = node['C1']
                    C2 = node['C2']
                    C3 = node['C3']

                    prev_node = G.node[node['pnode']]
                    # Evaluate Q coming in a given node at time t:
                    node['Qin'][t] = Muskingum(C1, C2, C3,
                                                   prev_node['Qout'][t],
                                                   prev_node['Qout'][t - 1],
                                                   node['Qin'][t - 1])

                    node['Qout'][t] = node['Qin'][t]

    def __call__(self, timestep=1, **kwargs):
        G = self.G
        Qpeaks = self.Qpeaks
        dikelist = self.dikelist

[5]:

dikeNetwork = DikeNetwork()

[6]:

dikeNetwork.calc_wave()

[7]:

G = dikeNetwork.getG()

[8]:

#G.nodes['A.1']

[9]:

plt.figure(figsize=(20,10))
plt.plot(G.node['A.0']['Qout'])
df = pd.DataFrame({'Qin':G.node['A.0']['Qout']})

dikelist = dikeNetwork.dikelist
for n in range(0, len(dikelist)):
    node = G.node[dikelist[n]]
    plt.plot(node['Qin'])
    df[dikelist[n]] = node['Qin']
../_images/A-ijssel_ijssel_muskingum_verification_10_0.svg
[10]:

df

[10]:
Qin A.1 A.2 A.3 A.4 A.5
0 449.789315 449.789315 449.789315 449.789315 449.789315 449.789315
1 428.608569 434.884329 436.410432 438.873036 442.977241 445.538383
2 448.566119 444.021351 443.074220 441.913907 439.794119 440.381511
3 460.705731 456.117734 454.886401 452.680952 449.137213 445.484723
4 438.468008 444.056434 445.300738 447.005789 449.982643 450.533818
5 407.024353 417.559643 420.263158 424.917618 432.513852 438.950471
6 382.125718 391.800518 394.417410 399.299993 407.124116 415.138465
7 374.944239 379.181877 380.453989 383.155573 387.363824 392.886278
8 448.863039 427.885339 422.888978 415.150818 402.120269 395.258323
9 650.612166 586.260166 570.082344 542.780756 497.895938 463.519715
10 967.408445 859.509791 831.655087 782.773655 703.222663 634.204853
11 1247.376736 1140.893765 1112.294971 1059.323219 974.276058 888.789152
12 1562.398600 1445.838133 1414.832676 1357.720201 1265.766741 1176.514477
13 1841.356795 1733.284216 1704.088685 1649.309501 1561.558739 1471.585192
14 1972.873757 1910.338320 1892.431930 1856.290187 1799.349871 1731.354583
15 2000.000000 1978.325419 1971.500439 1955.974030 1932.043529 1898.332588
16 1952.257016 1961.676439 1963.432997 1964.500137 1966.988321 1961.870892
17 1836.171355 1872.621123 1881.726035 1896.794261 1921.657699 1939.918283
18 1635.558763 1702.948030 1720.251431 1750.378964 1799.511149 1841.089996
19 1424.241434 1501.549490 1522.038890 1559.333880 1619.471340 1677.268291
20 1196.013888 1280.495371 1302.973060 1344.299312 1410.837652 1475.526321
21 975.827326 1059.490621 1081.948139 1123.743472 1190.834698 1258.158144
22 832.736179 893.378025 910.215355 942.950282 994.960356 1052.585254
23 704.739634 755.888695 769.838052 796.557678 839.222367 884.064722
24 597.517524 640.441107 652.155593 674.535992 710.261457 748.077740
25 530.017304 559.377727 567.588646 583.756825 609.387049 638.315325
26 463.353857 489.508609 496.599998 510.101701 531.695728 554.022547
27 439.228066 452.080036 455.858388 463.720873 476.020882 491.640820
28 418.041637 427.121729 429.648454 434.684398 442.674746 451.495496
29 398.997576 406.620356 408.700256 412.691631 419.059141 425.839215
30 381.892728 388.623129 390.450030 393.916690 399.460510 405.215853