Case studies¶

In this page we propose the model configurations used in publications.

Dal Molin et al., 2020, HESS¶

This section contains instructions for the implementation of the semi-distributed hydrological model M02 presented in the article:

Dal Molin, M., Schirmer, M., Zappa, M., and Fenicia, F.: Understanding dominant controls on streamflow spatial variability to set up a semi-distributed hydrological model: the case study of the Thur catchment, Hydrol. Earth Syst. Sci., 24, 1319–1345, https://doi.org/10.5194/hess-24-1319-2020, 2020.

In this application, the Thur catchment has been divided in 10 subcatchments and 2 hydrological response units (HRUs). Please refer to the article for the details; here we propose only the code needed to setup a model that corresponds to the one used in the publication.

Model structure¶

The two HRUs are represented using the same model structure represented in the figure.

The structure is similar to the one of HYMOD; its conversion to SuperflexPy is presented here

Note that also the temperature has been threated as a flux: this choice is not forced by the framework but, in this particular case, where it is the first element that needs it, this is particularly convenient. An alternative solution would have been designing the snow reservoir in such a way that the temperature becomes a state of the reservoir; this solution would have been preferable in the case where the element that needed the flux was not at the beginning of the structure.

Defining the elements¶

We here assume that all the elements are already existing; therefore they just need to be imported.

from superflexpy.implementation.elements.thur_model_hess import SnowReservoir, UnsaturatedReservoir, HalfTriangularLag, FastReservoir
from superflexpy.implementation.elements.structure_elements import Transparent, Junction, Splitter
from superflexpy.implementation.computation.pegasus_root_finding import PegasusPython
from superflexpy.implementation.computation.implicit_euler import ImplicitEulerPython

After this, all the elements must be initialized, defining the initial state and the parameters.

solver = PegasusPython()
approximator = ImplicitEulerPython(root_finder=solver)

upper_splitter = Splitter(
    direction=[
        [0, 1, None],    # P and T go to the snow reservoir
        [2, None, None]  # PET goes to the transparent element
    ],
    weight=[
        [1.0, 1.0, 0.0],
        [0.0, 0.0, 1.0]
    ],
    id='upper-splitter'
)

snow = SnowReservoir(
    parameters={'t0': 0.0, 'k': 0.01, 'm': 2.0},
    states={'S0': 0.0},
    approximation=approximator,
    id='snow'
)

upper_transparent = Transparent(
    id='upper-transparent'
)

upper_junction = Junction(
    direction=[
        [0, None],
        [None, 0]
    ],
    id='upper-junction'
)

unsaturated = UnsaturatedReservoir(
    parameters={'Smax': 50.0, 'Ce': 1.0, 'm': 0.01, 'beta': 2.0},
    states={'S0': 10.0},
    approximation=approximator,
    id='unsaturated'
)

lower_splitter = Splitter(
    direction=[
        [0],
        [0]
    ],
    weight=[
        [0.3],   # Portion to slow reservoir
        [0.7]    # Portion to fast reservoir
    ],
    id='lower-splitter'
)

lag_fun = HalfTriangularLag(
    parameters={'lag-time': 2.0},
    states={'lag': None},
    id='lag-fun'
)

fast = FastReservoir(
    parameters={'k': 0.01, 'alpha': 3.0},
    states={'S0': 0.0},
    approximation=approximator,
    id='fast'
)

slow = FastReservoir(
    parameters={'k': 1e-4, 'alpha': 1.0},
    states={'S0': 0.0},
    approximation=approximator,
    id='slow'
)

lower_transparent = Transparent(
    id='lower-transparent'
)

lower_junction = Junction(
    direction=[
        [0, 0]
    ],
    id='lower-junction'
)

Defining the HRUs structure¶

Once all the elements have been created we can connect them composing the two HRUs.

consolidated = Unit(
    layers=[
        [upper_splitter],
        [snow, upper_transparent],
        [upper_junction],
        [unsaturated],
        [lower_splitter],
        [slow, lag_fun],
        [lower_transparent, fast],
        [lower_junction],
    ],
    id='consolidated'
)

unconsolidated = Unit(
    layers=[
        [upper_splitter],
        [snow, upper_transparent],
        [upper_junction],
        [unsaturated],
        [lower_splitter],
        [slow, lag_fun],
        [lower_transparent, fast],
        [lower_junction],
    ],
    id='unconsolidated'
)

Defining the catchments¶

Now that the HRUs have been created, we need to assign them to the catchments

andelfingen = Node(
    units=[consolidated, unconsolidated],
    weights=[0.24, 0.76],
    area=403.3,
    id='andelfingen'
)

appenzell = Node(
    units=[consolidated, unconsolidated],
    weights=[0.92, 0.08],
    area=74.4,
    id='appenzell'
)

frauenfeld = Node(
    units=[consolidated, unconsolidated],
    weights=[0.49, 0.51],
    area=134.4,
    id='frauenfeld'
)

halden = Node(
    units=[consolidated, unconsolidated],
    weights=[0.34, 0.66],
    area=314.3,
    id='halden'
)

herisau = Node(
    units=[consolidated, unconsolidated],
    weights=[0.88, 0.12],
    area=16.7,
    id='herisau'
)

jonschwil = Node(
    units=[consolidated, unconsolidated],
    weights=[0.9, 0.1],
    area=401.6,
    id='jonschwil'
)

mogelsberg = Node(
    units=[consolidated, unconsolidated],
    weights=[0.92, 0.08],
    area=88.1,
    id='mogelsberg'
)

mosnang = Node(
    units=[consolidated],
    weights=[1.0],
    area=3.1,
    id='mosnang'
)

stgallen = Node(
    units=[consolidated, unconsolidated],
    weights=[0.87, 0.13],
    area=186.6,
    id='stgallen'
)

waengi = Node(
    units=[consolidated, unconsolidated],
    weights=[0.63, 0.37],
    area=78.9,
    id='waengi'
)

Note that all the catchments incorporate the information about their area that will then be used by the network.

Not all the catchment must have all the HRUs; if an HRU is not present in a catchment (e.g. unconsolidated in Mosnang, line 50) it can be simply omitted.

Defining the network¶

The last step consists in creating the network that connects all the catchments previously declared.

model = Network(
    nodes=[
        andelfingen,
        appenzell,
        frauenfeld,
        halden,
        herisau,
        jonschwil,
        mogelsberg,
        mosnang,
        stgallen,
        waengi,
    ],
    topography={
        'andelfingen': None,
        'appenzell': 'stgallen',
        'frauenfeld': 'andelfingen',
        'halden': 'andelfingen',
        'herisau': 'halden',
        'jonschwil': 'halden',
        'mogelsberg': 'jonschwil',
        'mosnang': 'jonschwil',
        'stgallen': 'halden',
        'waengi': 'frauenfeld',
    }
)

Running the model¶

Now that all the components have been initialized, we can run the model.

The first step is to assign the input fluxes to the single elements. For this we assume that the data is available as a pandas DataFrame and that the columns are named P_name_of_the_catchment, T_name_of_the_catchment, and PET_name_of_the_catchment.

The inputs can be set using a for loop

for cat, cat_name in zip(catchments, catchments_names):
    cat.set_input([
        df['P_{}'.format(cat_name)].values,
        df['T_{}'.format(cat_name)].values,
        df['PET_{}'.format(cat_name)].values,
    ])

After this, the last thing to be done before actually running the model is setting the time step used in the simulations. This can be done directly at the network level and it will be set to all the components.

thur_catchment.set_timestep(1.0)

The only thing left to do is running the model!

thur_catchment.get_output()