Masters programme | E-portfolio
Semester I & II
Spatial Analyses

Runoff analyses

Hydrological analyses using the example of an alpine watershed with a reservoir. Spatial delineation of the catchment area, consideration of flow paths with runoff accumulation as well as temporal analyses of runoff volumes.

Table of Contents

Study area & Research question

This analysis focuses on the application of GIS hydrology tools to characterise runoff behaviour in an certain area of interest (AoI). As shown below, the area surrounding the Tauernmoossee – a reservoir in the south of the federal state of Salzburg – was choosen to be the AoI for the purpose of the following analyses. These include three main parts:

  • delineating the catchment draining into the lake and visualising exemplary flow paths
  • calculating flow accumulation across the catchment in order to identify important drainage channels
  • setting up a simple model to simulate the hydrological response of the catchment to precipitation events

Data & Methods

A DEM with 10m resolution, provided as open data by the austrian authorities, forms the basis for all subsequent spatial calculations. Within ArcGISPro the layer was derived by setting a corresponding definition query on the multi-resolution world elevation layer „Terrain“. To speed up the process, the DEM was then clipped to the extent of the AoI. As an exact delineation of the AoI would require knowledge on the catchment area coverage, the AoI was defined somewhat more broadly based on an initial visual delineation of the catchment area.

In order to conduct hydrological analyses, typical preprocessing steps including sink filling and flow direction calculation were performed. Building on this foundation, the following figure shows the subsequent analyses steps as a model builder visualisation grouped into colour-coded compartments. Of course, the different compartments cannot be regarded as completely isolated units since some outputs from one part are required as inputs for a tool in another compartment. The delineated catchment, for example, is needed as an input for extracting flow paths with a minimum through-flow of at least 15% of the catchment (carried out in part III -> Raster Calculator (4)). Regarding the response of the catchment to precipitation, the focus is on analysing the runoff at the lowest point of the catchment on the northern shore of the lake. Hydrographs are simulated for the case that the water level of the reservoir is not actively controlled.  Hence, all runoff into the lake is considered to pass the pour point on the northern shore. Two models are prepared to simulate the time-variant discharge at this point:

a.) The first is simply based on flow lengths, i.e. calculations of downstream distances along flow paths. Apart from the topography-dependent distance, no further factors influencing the runoff behaviour are considered here. The runoff velocity is assumed to be equal to 1 throughout the whole catchment. 

b.) A slightly more elaborated approach creates a velocity field based on the slope and flow accumulation. This method was proposed by Maidment et al. and utilises the following equation:  $$V = V_{m} * \frac{V_{m} * s^{b}A^{c}}{s^{b}A_{m}^{c}}$$ with $V$ the velocity for a single cell, $V_{m}$ the average velocity for all cells, the local slope $s$ and the upstream contributing area $A$. As recommended by Maidment et al. in the absence of any calibration possibilities, the coefficients $b$ and $c$ are set to 0.5 and for $V_{m}$ a value of 0.1 m/s is assumed. As depicted at the bottom-left of the graphical workflow below, the resulting velocity grid is further processed by adjusting all vales for cells located on the lake. Due to the minimal slope, runoff is assumed to happen very slowly here ignoring that runoff propagation in water bodies is not bound to slope at all and follows completely different laws than in terrestrial drainage channels. For the sake of simplicity, a value of 0.5 m/s is assumed to be an appropriate estimate for the velocity with which the runoff travels through the lake. The inverse of the resulting velocity grid is then used as a weight for calculating flow lengths in order to derive a second model for subsequent hydrograph simulations. 

The simulations including their visualisations were then performed in a python environment. Further methodological aspects are explained in the course of the presentation and discussion of the results. The scripted version of the ArcGiS model builder workflow as well as the python analyses can be found in detail at the end of this site.

Results & Discussion

Based on the flow length models, one can create simple unit hydrographs in order to analyse the hypothetical response of the catchment to one unit of spatially homogeneously occurring, constant intensity rainfall. The unit hydrographs are basically a classified representation of the runoff timings, whose spatial variability are shown in the maps above. In addition to the already mentioned differences between the two models, their unit hydrographs also point to differences in the resulting times of highest runoff accumulation at the pour point. For the simple model, the most critical zone corresponds to the interval [0.5h, 1h], whereas for the advanced model the maximum amount of discharge steams from locations that are 2.5 to 3 hours away from the pour point. Taking a look at the location of the corresponding zones in the maps above, there is not only a significant difference in timing but also in the spatial location of the critical zones. Analysing the spatial component of this more closely may be highly relevant in cases where the peak discharge involves a potential hazard and must be minimised as far as possible. In such cases, one may, for example, avoid clear-cutting of forests within the critical zones in order to minimise direct surface runoff through increased interception evaporation. The same logic would apply to restrict surface sealing in favour of increased infiltration capacity and correspondingly delayed groundwater runoff of precipitation. In the present case of a reservoir at the lower end of the catchment, however, these considerations are rather hypothetical and not as significant. It should also be noted that the identification of critical zones as those areas on which precipitation reaches the pour point at a similar time is only valid when considering short precipitation peaks. As soon as a continuous rainfall event occurs, all areas in the catchment contribute to the maximum runoff. This can also been seen from the subsequent simulations.

Simplified unit hydrographs based on the simple (left) and advanced (right) flow length raster

In order to simulate the hydrological response of the catchment to a specific rainfall event, the concept of linear system theory (which also forms the basis of unit hydrographs) was applied. This means that the runoff is simply simulated as superpositions of single runoff curves accounting for the necessary time lags. The response functions for four events, each with a total of 20 mm precipitation, are shown below. The difference along the horizontal axis of the plot grid relates to the precipitation time series being simulated either as following a plateau or gaussian function. Between the top and bottom row of the plots, the time span during which the precipitation event happened was modified. All simulations are based on the advanced flow length model.

As expected, in case of a short rainfall event the response function closely follows the unit hydrograph. Furthermore, the difference between the plateau and gaussian time series are rather small. In both cases the maximum discharge exceeds 50 m3/s and the function follows a similar course. This changes when the time span of the event increases relative to the avergage runoff timings of the catchment. The plots in the bottom row differ significantly both in the shape of the response function and in the maximum discharge.

Simulated hydrographs using different rain series

It should be emphasised that even though the applied model is more elaborated than the simple model, there are still some significant simplifications inherent. First of all, the model assumes that the runoff velocity is discharge invariant, whereas in reality it certainly depends on the amount of water flow. Second, the model does not map the influence of surface roughness. Third, runoff factors are not considered, i.e. it is assumed that precipitation directly translates into surface runoff. Further simplifications have already been mentioned in the methods section above.

Final remarks

Appendix - Scripts

With regard to the function of the reservoir, it would of course be interesting to look at the lake not only as an unchanging surface with a pour point at one end, but to analyse the actual storage function of the lake. Water level changes due to rainfall-induced inflow could be analysed more closely if a terrain model of the reservoir basin or bottom were available. However, since the elevation values of the DEM refer to the water surface, the use of existing tools (e.g. Storage Capacity Tool) makes little sense in the analyses above. However, the total storage volume of the lake of 55 hm3 known from other sources can be put in relation to the precipitation event of 20 mm simulated above. Over the entire catchment area, the runoff into the lake cumulates to about 0.45 hm3, which is less than one hundredth of the maximum storage volume. 

Instead of performing the above presented analyses with ArcGISPro, a perhaps better-suited alternative especially for more in-depth analyses may be GRASS GIS. GRASS provides a dedicated suite of hydrology tools including ready-to-use model implementations such as the physically based TOPMODEL. For a tutorial on how to use GRASS in the field of hydrological analyses including a simulation of flow discharge with consideration of surface roughness, the reader is referred to the following. It may also be worth it to explore this repository, which lists additional python resources for hydrologists. Nevertheless, for the scope of the present task using ArcGIS was sufficient.

ArcGIS Workflow

				
					# -*- coding: utf-8 -*-
"""
Generated automatically by ArcGIS ModelBuilder
"""
import arcpy
from arcpy.ia import *
from arcpy.sa import *
from arcpy.sa import *
from arcpy.sa import *


def workflow():
    # To allow overwriting outputs change overwriteOutput option to True.
    arcpy.env.overwriteOutput = False

    # Check out any necessary licenses.
    arcpy.CheckOutExtension("spatial")
    arcpy.CheckOutExtension("3D")
    arcpy.CheckOutExtension("ImageAnalyst")
    arcpy.CheckOutExtension("ImageExt")

    arcpy.ImportToolbox(
        r"c:\users\felix\appdata\local\programs\arcgis\pro\Resources\ArcToolbox\toolboxes\Conversion Tools.tbx"
    )
    # Model Environment settings
    with arcpy.EnvManager(
        extent="DEFAULT",
        outputCoordinateSystem="""PROJCS['WGS_1984_UTM_Zone_33N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',
        SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],
        PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],
        PARAMETER['Central_Meridian',15.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],
        UNIT['Meter',1.0]]""",
    ):
        dem_snapped = arcpy.Raster("dem10_stubbachtal")
        pour_point = "pour_point"
        assignment_runoff = "C:\\Users\\felix\\Desktop\\assignment_runoff"
        Streams_Salzburg = "Streams_Salzburg"
        Input_false_raster_or_constant_value = 0.5
        Input_true_raster_or_constant_value = 1
        Input_false_raster_or_constant_value_2_ = 0
        Catchments_Salzburg = "Catchments_Salzburg"
        drop_points = "drop_points"

        # Process: Fill (Fill) (sa)
        dem_filled = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\dem_filled"
        )
        Fill = dem_filled
        dem_filled = arcpy.sa.Fill(in_surface_raster=dem_snapped, z_limit=None)
        dem_filled.save(Fill)

        # Process: Flow Direction (Flow Direction) (sa)
        flowdir = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flowdir"
        Flow_Direction = flowdir
        flowdrop = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flowdrop"
        flowdir = arcpy.sa.FlowDirection(
            in_surface_raster=dem_filled,
            force_flow="NORMAL",
            out_drop_raster=flowdrop,
            flow_direction_type="D8",
        )
        flowdir.save(Flow_Direction)

        flowdrop = arcpy.Raster(flowdrop)

        # Process: Slope (Slope) (3d)
        dem_slope = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\dem_slope"
        )
        arcpy.ddd.Slope(
            in_raster=dem_snapped,
            out_raster=dem_slope,
            output_measurement="PERCENT_RISE",
            z_factor=1,
            method="GEODESIC",
            z_unit="METER",
        )
        dem_slope = arcpy.Raster(dem_slope)

        # Process: Watershed (Watershed) (sa)
        catchment = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\catchment"
        )
        Watershed = catchment
        catchment = arcpy.sa.Watershed(
            in_flow_direction_raster=flowdir,
            in_pour_point_data=pour_point,
            pour_point_field="OBJECTID",
        )
        catchment.save(Watershed)

        # Process: Con (2) (Con) (sa)
        flowdir_clipped = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flowdir_clipped"
        )
        Con_2_ = flowdir_clipped
        flowdir_clipped = arcpy.sa.Con(
            in_conditional_raster=catchment,
            in_true_raster_or_constant=flowdir,
            in_false_raster_or_constant="",
            where_clause="",
        )
        flowdir_clipped.save(Con_2_)

        # Process: Flow Length (2) (Flow Length) (sa)
        flow_length_simple = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flow_length_simple"
        Flow_Length_2_ = flow_length_simple
        with arcpy.EnvManager(
            extent="314783.839005856 5218149.23804989 329003.839005856 5228739.23804989"
        ):
            flow_length_simple = arcpy.sa.FlowLength(
                in_flow_direction_raster=flowdir_clipped,
                direction_measurement="DOWNSTREAM",
                in_weight_raster="",
            )
            flow_length_simple.save(Flow_Length_2_)

        # Process: Raster To Other Format (Raster To Other Format) (conversion)
        Updated_Output_Workspace = arcpy.conversion.RasterToOtherFormat(
            Input_Rasters=[flow_length_simple],
            Output_Workspace=assignment_runoff,
            Raster_Format="TIFF",
        )[0]

        # Process: Select (Select) (analysis)
        lake_tauernmoossee = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\lake_tauernmoossee"
        arcpy.analysis.Select(
            in_features=Streams_Salzburg,
            out_feature_class=lake_tauernmoossee,
            where_clause="WIS_NAME = 'Stausee Tauernmoossee'",
        )

        # Process: Feature To Polygon (Feature To Polygon) (management)
        lake_polygon = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\lake_polygon"
        )
        arcpy.management.FeatureToPolygon(
            in_features=[lake_tauernmoossee],
            out_feature_class=lake_polygon,
            cluster_tolerance="",
            attributes="ATTRIBUTES",
            label_features="",
        )

        # Process: Feature to Raster (Feature to Raster) (conversion)
        lake_raster = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\lake_raster"
        )
        arcpy.conversion.FeatureToRaster(
            in_features=lake_polygon,
            field="Shape_length",
            out_raster=lake_raster,
            cell_size="5",
        )
        lake_raster = arcpy.Raster(lake_raster)

        # Process: Is Null (Is Null) (ia)
        lake_boolean = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\lake_boolean"
        )
        Is_Null = lake_boolean
        with arcpy.EnvManager(
            extent="314783.839005856 5218149.23804989 329003.839005856 5228739.23804989"
        ):
            lake_boolean = arcpy.ia.IsNull(in_raster=lake_raster)
            lake_boolean.save(Is_Null)

        # Process: Flow Accumulation (Flow Accumulation) (sa)
        flowacc = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flowacc"
        Flow_Accumulation = flowacc
        flowacc = arcpy.sa.FlowAccumulation(
            in_flow_direction_raster=flowdir,
            in_weight_raster="",
            data_type="FLOAT",
            flow_direction_type="D8",
        )
        flowacc.save(Flow_Accumulation)

        # Process: Con (3) (Con) (sa)
        flowacc_clipped = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flowacc_clipped"
        )
        Con_3_ = flowacc_clipped
        flowacc_clipped = arcpy.sa.Con(
            in_conditional_raster=catchment,
            in_true_raster_or_constant=flowacc,
            in_false_raster_or_constant="",
            where_clause="",
        )
        flowacc_clipped.save(Con_3_)

        # Process: Raster Calculator (Raster Calculator) (sa)
        area_slope_factor = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\area_slope_factor"
        Raster_Calculator = area_slope_factor
        if flowacc_clipped:
            with arcpy.EnvManager(mask=catchment):
                area_slope_factor = SquareRoot("dem_slope") + SquareRoot("flowacc")
                area_slope_factor.save(Raster_Calculator)

        # Process: Get Raster Properties (Get Raster Properties) (management)
        if flowacc_clipped:
            mean_factor = arcpy.management.GetRasterProperties(
                in_raster=area_slope_factor, property_type="MEAN", band_index=""
            )[0]

        # Process: Raster Calculator (2) (Raster Calculator) (sa)
        velocity = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\velocity"
        Raster_Calculator_2_ = velocity
        if flowacc_clipped:
            velocity = 0.1 * ("area_slope_factor" / float(mean_factor))
            velocity.save(Raster_Calculator_2_)

        # Process: Con (Con) (sa)
        velocity_adjusted = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\velocity_adjusted"
        Con = velocity_adjusted
        if flowacc_clipped:
            velocity_adjusted = arcpy.sa.Con(
                in_conditional_raster=lake_boolean,
                in_true_raster_or_constant=velocity,
                in_false_raster_or_constant=Input_false_raster_or_constant_value,
                where_clause="",
            )
            velocity_adjusted.save(Con)

        # Process: Raster Calculator (3) (Raster Calculator) (ia)
        velocity_weight = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\velocity_weight"
        )
        Raster_Calculator_3_ = velocity_weight
        if flowacc_clipped:
            velocity_weight = 1 / velocity_adjusted
            velocity_weight.save(Raster_Calculator_3_)

        # Process: Flow Length (Flow Length) (sa)
        flow_length_advanced = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flow_length_advanced"
        Flow_Length = flow_length_advanced
        if flowacc_clipped:
            flow_length_advanced = arcpy.sa.FlowLength(
                in_flow_direction_raster=flowdir_clipped,
                direction_measurement="DOWNSTREAM",
                in_weight_raster=velocity_weight,
            )
            flow_length_advanced.save(Flow_Length)

        # Process: Raster To Other Format (3) (Raster To Other Format) (conversion)
        if flowacc_clipped:
            Updated_Output_Workspace_3_ = arcpy.conversion.RasterToOtherFormat(
                Input_Rasters=[flow_length_advanced],
                Output_Workspace=assignment_runoff,
                Raster_Format="TIFF",
            )[0]

        # Process: Con (4) (Con) (sa)
        Watershed_binary = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\Watershed_binary"
        Con_4_ = Watershed_binary
        Watershed_binary = arcpy.sa.Con(
            in_conditional_raster=catchment,
            in_true_raster_or_constant=Input_true_raster_or_constant_value,
            in_false_raster_or_constant=Input_false_raster_or_constant_value_2_,
            where_clause="",
        )
        Watershed_binary.save(Con_4_)

        # Process: Zonal Statistics as Table (Zonal Statistics as Table) (ia)
        ZonalSt_Watersh1 = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\ZonalSt_Watersh1"
        arcpy.ia.ZonalStatisticsAsTable(
            in_zone_data=Watershed_binary,
            zone_field="VALUE",
            in_value_raster=Watershed_binary,
            out_table=ZonalSt_Watersh1,
            ignore_nodata="DATA",
            statistics_type="SUM",
            process_as_multidimensional="CURRENT_SLICE",
            percentile_values=90,
            percentile_interpolation_type="AUTO_DETECT",
        )

        # Process: Get_Field_Value (Get Field Value)
        # Get Field Value Utility is not implemented

        # Process: Raster Calculator (4) (Raster Calculator) (sa)
        high_flowacc = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\high_flowacc"
        )
        Raster_Calculator_4_ = high_flowacc
        high_flowacc = Con("flowacc_clipped" > 0.15 * int(Value), "flowacc_clipped")
        high_flowacc.save(Raster_Calculator_4_)

        # Process: Raster to Polygon (Raster to Polygon) (conversion)
        catchment_polygon = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\catchment_polygon"
        with arcpy.EnvManager(outputMFlag="Disabled", outputZFlag="Disabled"):
            arcpy.conversion.RasterToPolygon(
                in_raster=catchment,
                out_polygon_features=catchment_polygon,
                simplify="NO_SIMPLIFY",
                raster_field="Value",
                create_multipart_features="SINGLE_OUTER_PART",
                max_vertices_per_feature=None,
            )

        # Process: Select (2) (Select) (analysis)
        catchment_reference = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\catchment_reference"
        arcpy.analysis.Select(
            in_features=Catchments_Salzburg,
            out_feature_class=catchment_reference,
            where_clause="HZB_CODE IN ('2  8272 77  3 1 0 0 0', '2  8272 77  3 3 0 0 0', '2  8272 77  3 2 0 0 0')",
        )

        # Process: Symmetrical Difference (Symmetrical Difference) (analysis)
        catchment_differences = "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\catchment_differences"
        arcpy.analysis.SymDiff(
            in_features=catchment_polygon,
            update_features=catchment_reference,
            out_feature_class=catchment_differences,
            join_attributes="ALL",
            cluster_tolerance="",
        )

        # Process: Cost Path (Cost Path) (sa)
        flow_paths = (
            "C:\\Users\\felix\\Desktop\\assignment_runoff\\Default.gdb\\flow_paths"
        )
        Cost_Path = flow_paths
        flow_paths = arcpy.sa.CostPath(
            in_destination_data=drop_points,
            in_cost_distance_raster=flowacc,
            in_cost_backlink_raster=flowdir,
            path_type="EACH_CELL",
            destination_field="OBJECTID",
            force_flow_direction_convention="INPUT_RANGE",
        )
        flow_paths.save(Cost_Path)


if __name__ == "__main__":
    # Global Environment settings
    with arcpy.EnvManager(
        scratchWorkspace=r"C:\Users\felix\Desktop\assignment_runoff\Default.gdb",
        workspace=r"C:\Users\felix\Desktop\assignment_runoff\Default.gdb",
    ):
        workflow()

				
			

Python - Hydrograph visualisations

				
					import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import numpy as np
import pandas as pd
import os
import rasterio
from datetime import datetime
from datetime import timedelta
from scipy.stats import norm


### part I - basic unit hydrograph ###
work_dir = "C:/Users/felix/Desktop/assignment_runoff/"
exp_raster = ["flow_length_simple_1.tif", "flow_length_advanced_1.tif"]

for i in exp_raster:
    # import raster and ignore na values
    flow_raster = rasterio.open(os.path.join(work_dir, i))
    flow_vals = flow_raster.read(1)
    flow_vals = flow_vals[flow_vals < 10 ** 10]
    # convert time into hours
    flow_vals = flow_vals / 3600
    # plotting of unit hydrograph
    fig, ax1 = plt.subplots()
    bins = np.arange(0, 6, 0.5).tolist()
    ax1.hist(flow_vals, bins=bins, color="darkred", alpha=0.9, rwidth=0.9)
    ax1.set_xlabel("time (h)")
    ax1.set_ylabel("unit discharge")
    ax1.set_xticks(bins)
    ax1.set_yticks([])
    ax1.set_yticklabels([])
    ax1.set_xlim(0, 5.5)
    ax1.set_ylim(0, 70000)
    ax1.spines["top"].set_visible(False)
    ax1.spines["right"].set_visible(False)
    fig.savefig(
        os.path.join(work_dir, f'hydrograph_unit_{i.split("_")[2]}'),
        bbox_inches="tight",
        dpi=300,
    )


### part II - simulated hydrograph ###
def hydrograph(flow_time_raster, start, end, sim_type, exp_name):
    # import raster, ignore na, get cell_size
    flow_raster = rasterio.open(os.path.join(work_dir, flow_time_raster))
    flow_vals = flow_raster.read(1)
    flow_vals = flow_vals[flow_vals < 10 ** 10]
    flow_cell_size = flow_raster.transform[0] * -flow_raster.transform[4]

    # define time interval
    range_min = (end - start).seconds / 60
    time_steps = [start + timedelta(minutes=1 * x) for x in range(1, int(range_min))]
    time_steps = [start, *time_steps, end]
    n_time_step = list(range(0, len(time_steps)))

    # set/simulate rain values
    total_rain_mm = 20
    if sim_type == "plateau":
        rain_series = [
            *np.linspace(0, 1, int(len(time_steps) / 6)),
            *np.random.uniform(0.9, 1.1, 4 * int(len(time_steps) / 6)),
            *np.linspace(1, 0, int(len(time_steps) / 6)),
        ]
    if sim_type == "gauss":
        rain_series = [
            *[norm.pdf(x) for x in np.linspace(-3, 0, int(len(time_steps) / 2))],
            *[norm.pdf(x) for x in np.linspace(-3, 0, int(len(time_steps) / 2))][::-1],
        ]
    rain_series = [x * total_rain_mm / sum(rain_series) for x in rain_series]

    # construct time series
    cols = {
        "n": n_time_step,
        "time": time_steps,
        "rain": rain_series,
    }
    rain = pd.concat([pd.Series(v, name=k) for k, v in cols.items()], axis=1)
    rain = rain.fillna(0)

    # get the minute-by-minute number of cells accumulating outflow at the pour point
    bins = np.arange(0, flow_vals.max() + 60, 60).tolist()
    flow_binned = np.histogram(flow_vals, bins)[0]

    # calculate single discharge series tied to minutely rain fall
    def single_discharge_series(time_step):
        discharge_isolated = rain["rain"][time_step] * flow_binned
        nulls_upfront = [0] * time_step
        nulls_after = [0] * (max(n_time_step) - time_step)
        discharge_embedded = [*nulls_upfront, *discharge_isolated, *nulls_after]
        return discharge_embedded

    discharge_series = list(map(single_discharge_series, rain["n"]))

    # overlay single series to get overall discharge series
    discharge_series = np.sum(np.array(discharge_series), axis=0)

    # convert discharge series to unit m3/s
    # bind discharge series to rain df
    # convert rain fall to mm/h
    discharge_series = [flow_cell_size * x / 60 / 1000 for x in discharge_series]
    rain_discharge = pd.concat(
        [rain, pd.DataFrame({"discharge": discharge_series})], axis=1
    )
    rain_discharge["n"] = range(0, len(rain_discharge))
    rain_discharge["time"] = list(
        map(lambda x: start + timedelta(minutes=x), rain_discharge["n"])
    )
    rain_discharge["rain"] = [x * 60 for x in rain_discharge["rain"]]

    # plot hydrograph
    fig, ax1 = plt.subplots()
    # plotting discharge on left axis
    ax1.plot("time", "discharge", data=rain_discharge, color="darkred")
    ax1.set_xlabel("time")
    ax1.set_ylabel("discharge (m\u00b3/s)\n")
    ax1.set_xticklabels(ax1.get_xticks(), rotation=45)
    ax1.set_ylim(0)
    # plotting rain fall on right axis
    ax2 = ax1.twinx()
    ax2.bar(
        "time", data=rain_discharge[::3], height="rain", width=0.001, color="darkblue"
    )
    ax2.set_ylabel("\nrainfall (mm/h)")
    ax2.set_ylim(0, 2 * np.max(rain_discharge["rain"]))
    ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
    # set plotting order
    ax1.set_zorder(ax2.get_zorder() + 1)
    ax1.patch.set_visible(False)
    # save fig
    fig.savefig(os.path.join(work_dir, exp_name), bbox_inches="tight", dpi=300)


# simulate different scenerios
hydrograph(
    flow_time_raster="flow_length_advanced_1.tif",
    start=datetime.strptime("01/07 10:00:00", "%d/%m %H:%M:%S"),
    end=datetime.strptime("01/07 11:00:00", "%d/%m %H:%M:%S"),
    sim_type="gauss",
    exp_name="hydrograph_sim_1",
)

hydrograph(
    flow_time_raster="flow_length_advanced_1.tif",
    start=datetime.strptime("01/07 10:00:00", "%d/%m %H:%M:%S"),
    end=datetime.strptime("01/07 11:00:00", "%d/%m %H:%M:%S"),
    sim_type="plateau",
    exp_name="hydrograph_sim_2",
)

hydrograph(
    flow_time_raster="flow_length_advanced_1.tif",
    start=datetime.strptime("01/07 10:00:00", "%d/%m %H:%M:%S"),
    end=datetime.strptime("01/07 18:00:00", "%d/%m %H:%M:%S"),
    sim_type="gauss",
    exp_name="hydrograph_sim_3",
)

hydrograph(
    flow_time_raster="flow_length_advanced_1.tif",
    start=datetime.strptime("01/07 10:00:00", "%d/%m %H:%M:%S"),
    end=datetime.strptime("01/07 18:00:00", "%d/%m %H:%M:%S"),
    sim_type="plateau",
    exp_name="hydrograph_sim_4",
)