1d Sediment Transport

Get ready by cloning the exercise repository:

git clone https://github.com/Ecohydraulics/Exercise-SedimentTransport.git

Theory¶

1d Cross-section Averaged Hydrodynamics¶

From the stage-discharge (Manning-Strickler formula) exercise, we recall the formula to calculate the relationship between water depth $h$ (incorporated in the hydraulic radius $R_{h}$) and flow velocity $u$:

where

• $n_m$ is the Manning coefficient in fictional units of (s/m$^{1/3}$).

• $S_{e}$ is the hypothetical energy slope (m/m) and corresponds to the channel slope for steady, uniform flow conditions (non-existing in natural rivers).

• hydraulic radius $R_{h} = A / P$, where (for a trapezoidal cross-section):

  • the wetted (trapezoidal) cross-section area is $A = h \cdot 0.5\cdot (b + B) = h \cdot (b + h\cdot m)$;

  • the wetted perimeter of a trapezoid is $P = b + 2h\cdot(m^2 + 1)^{1/2}$;

  • $b$ (channel base width) and $m$ (bank slope) are illustrated in the figure below to calculate the depth-dependent water surface width $B=b+2\cdot h\cdot m$.

This exercise uses one-dimensional (1d) cross-section averaged hydraulic data produced with the US Army Corps of Engineers HEC-RAS software U.S. Army Corps of Engineeers, 2016, which solves the Manning-Strickler formula numerically for any flow cross-section shape. In this exercise, HEC-RAS provides the hydraulic data needed to determine the sediment transport capacity of a channel cross-section, although no explanations for creating, running, and exporting data from HEC-RAS models are given.

Sediment Transport¶

Fluvial Sediment transport can be distinguished into two modes: (1) suspended load and (2) bedload (see Fig. 3). Finer particles with a weight that can be carried by the fluid (water) are transported as suspended load. Coarser particles rolling, sliding, and jumping on the channel bed are transported as bedload. There is another type of transport, the so-called wash load, which is finer than the coarse bedload, but too heavy (large) to be transported in suspension Einstein, 1950.

In the following, we will look at the bedload transport mode. In this case, a sediment particle located in or on the riverbed is mobilized by shear forces of the water as soon as they exceed a critical value (see figure below). In river hydraulics, the so-called dimensionless bed shear stress or Shields stress Shields, 1936 is often used as the threshold value for the mobilization of sediment from the riverbed (see Fig. 4). This exercise uses one of the dimensionless bed shear stress approaches and the next section provides more explanations.

The Meyer-Peter and Müller (1948) Formula¶

The Meyer-Peter & Müller (1948) formula for estimating bedload transport was published by Swiss researchers Eugen Meyer-Peter (founder of the Laboratory of Hydraulics, Hydrology and Glaciology (VAW) and Robert Müller. Their study began one year after the establishment of the VAW in 1931 when Robert Müller was appointed assistant to Eugen Meyer-Peter. The two scientists worked in collaboration with Henry Favre and Albert Einstein’s son Hans Albert. In 1934, the laboratory published for the first time a formula for the calculation of bedload transport and its fundamental relationship between observed $\tau_{x}$ and critical $\tau_{x,cr}$ dimensionless bed shear stresses is used until today. The dimensionless bedload transport rate $\Phi_b$ according to Meyer-Peter & Müller (1948) is:

The other parameters are:

• $s$ $\approx$ 2.68, the dimensionless ratio of sediment grain density $\rho_{s}$ ($\approx$ 2680 kg/m³) and water density $\rho_{w}$ ($\approx$ 1000 kg/m³);

• $D_{char}$, the characteristic grain size in (m). It can be assumed that $D_{char} \approx D_{84}$ (i.e., the grain diameter of which 84% of a sediment mixture is smaller) in line with the scientific literature (e.g., Rickenmann & Recking (2011)).

The Meyer-Peter & Müller formula applies (like any other Sediment transport formula) only to certain rivers that have the following characteristics (range of validity):

• 0.4 $\cdot$ 10$^{-3}$ m $< D_{char}$ < 28.6 $\cdot$ 10$^{-3}$ m

• 10$^{-4}$ m $< Fr <$ 639 ($Fr$ denotes the dimensionless Froude number)

• 0.0004 $< S_{e} <$ 0.02

• 0.0002 m$^3$/(s $\cdot$ m) $< q <$ 2.0 m$^3$/(s $\cdot$ m) ($q$ is the unit discharge, i.e., $q=Q/[0. 5\cdot (b + B)]$)

• 0.25 $< s <$ 3.2

The dimensionless expression for bedload $\Phi_b$ was used to enable information transfer between different channels across scales by preserving geometric, kinematic, and dynamic similarity. The set of dimensionless parameters used results from Buckingham’s $\Pi$ theorem Buckingham, 1915. Therefore, to add dimensions to $\Phi_b$, it needs to be multiplied with the same set of parameters used for deriving the dimensionless expression from Meyer-Peter & Müller. Their set of parameters involves the characteristic grain size $D_{char}$, the grain density $\rho_{s}$, and the gravitational acceleration $g$. Thus, the dimensional unit bedload is (in kg/s and meter width, i.e., kg/(s$\cdot$m):

The cross-section averaged bedload $Q_{b}$ (kg/s) is then:

where $b_{eff}$ is the hydraulically active channel width of the flow cross-section (e.g., for a trapezoid $b_{eff} = 0.5 \cdot (b + B)$).

Code¶

Set the Frame¶

The object-oriented code will use custom classes that we will call within a main.py script. Create the following additional scripts, which will contain the custom classes and functions to control logging.

• fun.py will contain logging functions.

• hec.py will contain a HecSet class to read hydraulic output data from HEC-RAS as structured objects.

• grains.py will contain a GrainReader class to read grain size class information as structured objects.

• bedload.py will contain the class BedCore with basic elements that most bedload formulae have in common.

• mpm.py will contain the class MPM, which inherits from BedCore and calculates bedload as above described (Meyer-Peter & Müller 1948).

We will create the classes and functions in the indicated scripts according to the following flow chart:

To start with the main.py script, add a main function as well as a get_char_grain_size and a calculate_mpm function. Moreover, make the script stand-alone executable:

# This is main.py
import os


def get_char_grain_size(file_name, D_char):
    return None


def calculate_mpm(hec_df, D_char):
    return None


def main():
    pass


if __name__ == '__main__':
    main()

Logging Functions¶

The fun.py script will contain two functions:

1. start_logging to setup logging formats and a log file name as described in the section on logging, and

2. log_actions, which is a function wrapper for the main() (main.py) functions to log script execution messages.

The start_logging function should look like this (change the log file name if desired):

import logging


def start_logging():
    logging.basicConfig(filename="logfile.log", format="[%(asctime)s] %(message)s",
                        filemode="w", level=logging.DEBUG)
    logging.getLogger().addHandler(logging.StreamHandler()

The log_actions wrapper function follows the instructions from the functions theory section:

def log_actions(fun):
    def wrapper(*args, **kwargs):
        start_logging()
        fun(*args, **kwargs)
        logging.shutdown()
    return wrapper

To use the log_actions wrapper throughout the program, we will implement it at the highest level, which is the main() function in main.py:

# main.py
from fun import *

...

@log_actions
def main():
    logging.info("This is a test message (do not keep in the function).")


if __name__ == '__main__':
    main()

Now, we can log messages at different levels (info, warning, error, or others) in all functions called within main() by using for example logging.info("Message"), logging.warning("Message"), or logging.error("Message") rather than the print() function.

Read grain size data¶

Sediment grain size classes (ranging from $D_{16}$ to $D_{max}$) are provided in the file grains.csv (delimiter=",") and can be customized.

Write a GrainReader class that uses the read_csv method from Pandas to read the grain size distribution from grains.csv. Write the class in a separate Python script (e.g., grains.py as indicated in the above figure):

class GrainReader:
    def __init__(self, csv_file_name="grains.csv", delimiter=","):
        self.sep = delimiter
        self.size_classes = pd.DataFrame
        self.get_grain_data(csv_file_name)

The get_grain_data method should look like this for reading the provided grain size classes:

    def get_grain_data(self, csv_file_name):
        self.size_classes = pd.read_csv(csv_file_name,
                                        names=["classes", "size"],
                                        skiprows=[0],
                                        sep=self.sep,
                                        index_col=["classes"])

Implement the instantiation of a GrainReader object in the main.py script in the get_char_grain_size function. The function should receive the string-type arguments file_name (here: "grains.csv") and D_char (i.e., the characteristic grain size to use from grains.csv). The main() function calls the get_char_grain_size function with the arguments file_name=os.path.abspath("..") + "\\grains.csv" and D_char="D84" (corresponds to the first column in grains.csv).

# main.py
import os
from grains import GrainReader

def get_char_grain_size(file_name=str, D_char=str):
    grain_info = GrainReader(file_name)
    return grain_info.size_classes["size"][D_char]

...

@log_actions
def main():
    # get characteristic grain size = D84
    D_char = get_char_grain_size(file_name=os.path.abspath("..") + "\\grains.csv",
                                 D_char="D84")

Read HEC-RAS Input Data¶

The provided HEC-RAS dataset is stored in the xlsx workbook HEC-RAS/output.xlsx and contains the following output:

To load HEC-RAS output data, write a custom class (in a separate script called hec.py) that takes the file name as input argument and reads the HEC-RAS file as pandas data frame:

class HecSet:
    def __init__(self, xlsx_file_name="output.xlsx"):
        self.hec_data = pd.DataFrame
        self.get_hec_data(xlsx_file_name)

The get_hec_data method should look (something) like this:

    def get_hec_data(self, xlsx_file_name):
        self.hec_data = pd.read_excel(xlsx_file_name,
                                      skiprows=[1],
                                      header=[0])

To create a HecSet object in the main() (main.py) function, we need to import and instantiate it for example as hec = HecSet(file_name). In addition, we can already implement passing the pd.DataFrame of the HEC-RAS data to the calculate_mpm function (also in main.py) that we will complete later on.

# main.py
import os
from ...
from hec import HecSet

...

@log_actions
def main():
    D_char = ...

    hec_file = os.path.abspath("..") + "{0}HEC-RAS{0}output.xlsx".format(os.sep)
    hec = HecSet(hec_file)

Create a Bedload Core Class¶

A BedCore class written in the bedload.py script provides variables and methods, which are relevant to many bedload and Sediment transport calculation formulae, such as the Parker-Wong correction Wong & Parker, 2006 or the Smart & Jaeggi (1983) (direct download ). Moreover, the BedCore class contains constants such as the gravitational acceleration $g$ (i.e., self.g=9.81), the ratio of sediment grain and water density $s$ (i.e., self.s=2.68), and the critical dimensionless bed shear stress $\tau_{x,cr}$ (i.e., self.tau_xcr=0.047, which may be re-defined by users). The header of the BedCore class should look (similar) like this:

from fun import *
import numpy as np


class BedCore:
    def __init__(self):
        self.tau_x = np.nan
        self.tau_xcr = 0.047
        self.g = 9.81
        self.s = 2.68
        self.rho_s = 2680.0  # kg/m3 sediment grain density
        self.Se = np.nan  # energy slope (m/m)
        self.D = np.nan  # characteristic grain size
        self.Fr = np.nan  # Froude number
        self.h = np.nan  # water depth (m)
        self.phi = np.nan  # dimensionless bedload
        self.Q = np.nan  # discharge (m3/s)
        self.Rh = np.nan  # hydraulic radius (m)
        self.u = np.nan  # flow velocity (m/s)

Add a method to convert the dimensionless bedload transport $\Phi_b$ into a dimensional value (kg/s). In addition to the variables defined in the __init__ method, the add_dimensions method will require the effective channel width $b_{eff}$ (recall the above calculus):

    def add_dimensions(self, b):
        try:
            return self.phi * b * np.sqrt((self.s - 1) * self.g * self.D ** 3) * self.rho_s
        except ValueError:
            logging.warning("Non-numeric data. Returning Qb=NaN.")
            return np.nan

Many bedload transport formulae involve the dimensionless bed shear stress [$\tau_{x}$ (see above definitions) associated with a set of cross-section averaged hydraulic parameters. Therefore, implement the calculation method compute_tau_x in BedCore:

    def compute_tau_x(self):
        try:
            return self.Se * self.Rh / ((self.s - 1) * self.D)
        except ValueError:
            logging.warning("Non-numeric data. Returning tau_x=NaN.")
            return np.nan

Write a Meyer-Peter & Müller Bedload Assessment Class¶

Create a new script (e.g., mpm.py) and implement an MPM class (Meyer-Peter & Müller) that inherits from the BedCore class. The __init__ method of MPM should initialize BedCore and overwrite (recall Polymorphism) relevant parameters to the calculation of bedload according to Meyer-Peter & Müller (1948). Moreover, the initialization of an MPM object should go along with a check of the validity and the calculation of the dimensionless bedload transport $\Phi_b$ (see above explanations of MPM):

from bedload import *


class MPM(BedCore):
    def __init__(self, grain_size, Froude, water_depth,
                 velocity, Q, hydraulic_radius, slope):
        # initialize parent class
        BedCore.__init__(self)
        # assign parameters from arguments
        self.D = grain_size
        self.h = water_depth
        self.Q = Q
        self.Se = slope
        self.Rh = hydraulic_radius
        self.u = velocity
        self.check_validity(Froude)
        self.compute_phi()

Add the check_validity method to verify if the provided cross-section characteristics fall into the range of validity of the Meyer-Peter & Müller formula (i.e., slope, grain size, ratio of discharge and water depth, and Froude number):

    def check_validity(self, Fr):
        if (self.Se < 0.0004) or (self.Se > 0.02):
            logging.warning('Warning: Slope out of validity range.')
        if (self.D < 0.0004) or (self.D > 0.0286):
            logging.warning('Warning: Grain size out of validity range.')
        if ((self.u * self.h) < 0.002) or ((self.u * self.h) > 2.0):
            logging.warning('Warning: Discharge out of validity range.')
        if (self.s < 0.25) or (self.s > 3.2):
            logging.warning('Warning: Relative grain density (s) out of validity range.')
        if (Fr < 0.0001) or (Fr > 639):
            logging.warning('Warning: Froude number out of validity range.')

To calculate dimensionless bedload transport $\Phi_b$ according to Meyer-Peter & Müller, implement a compute_phi method that uses the compute_tau_x method from BedCore:

   def compute_phi(self):
        tau_x = self.compute_tau_x()
        try:
            if tau_x > self.tau_xcr:
                self.phi = 8 * (0.85 * tau_x - self.tau_xcr) ** (3 / 2)
            else:
                self.phi = 0.0
        except TypeError:
            logging.warning("Could not calculate PHI (result=%s)." % str(tau_x)
            self.phi = np.nan

With the MPM class defined, we can now fill the calculate_mpm function in the main.py script. The function should create a pandas data frame with columns of dimensionless bedload transport $\Phi$ and dimensional bedload transport $Q_{b}$ associated with a channel profile ("River Sta") and flow scenario ("Profile" > "Scenario").

The following code block illustrates an example of the calculate_mpm function that creates the pandas data frame from a Dictionary (mpm_dict). The illustrative function creates the dictionary with void value lists, extracts hydraulic data from the HEC-RAS data frame, and loops over the "River Sta" entries. The loop checks if the "River Sta" entries are valid (i.e., not "Nan") because empty rows that HEC-RAS automatically adds between output profiles should not be analyzed. If the check was successful, the loop appends the profile, scenario, and discharge directly to mpm_dict. The section-wise bedload transport results from MPM objects. After the loop, the function returns mpm_dict as a pd.DataFrame object.

# main.py
from ...
from ...
from mpm import *

...

def calculate_mpm(hec_df, D_char):
    # create dictionary with relevant information about bedload transport with void lists
    mpm_dict = {
            "River Sta": [],
            "Scenario": [],
            "Q (m3/s)": [],
            "Phi (-)": [],
            "Qb (kg/s)": []
    }

    # extract relevant hydraulic data from HEC-RAS output file
    Froude = hec_df["Froude # Chl"]
    h = hec_df["Hydr Depth"]
    Q = hec_df["Q Total"]
    Rh = hec_df["Hydr Radius"]
    Se = hec_df["E.G. Slope"]
    u = hec_df["Vel Chnl"]

    for i, sta in enumerate(list(hec_df["River Sta"]):
        if not str(sta).lower() == "nan":
            logging.info("PROCESSING PROFILE {0} FOR SCENARIO {1}".format(str(hec_df["River Sta"][i]), str(hec_df["Profile"][i]))
            mpm_dict["River Sta"].append(hec_df["River Sta"][i])
            mpm_dict["Scenario"].append(hec_df["Profile"][i])
            section_mpm = MPM(grain_size=D_char,
                              Froude=Froude[i],
                              water_depth=h[i],
                              velocity=u[i],
                              Q=Q[i],
                              hydraulic_radius=Rh[i],
                              slope=Se[i])
            mpm_dict["Q (m3/s)"].append(Q[i])
            mpm_dict["Phi (-)"].append(section_mpm.phi)
            b = hec_df["Flow Area"][i] / h[i]
            mpm_dict["Qb (kg/s)"].append(section_mpm.add_dimensions(b)
    return pd.DataFrame(mpm_dict)

Having defined the calculate_mpm() function, the call to that function from the main() function should now assign a Pandas data frame to the mpm_results variable. To finalize the script, write mpm_results to a workbook (e.g., "bed_load_mpm.xlsx") in the main() function:

# main.py
import os
from ...

...

def calculate_mpm(hec_df, D_char):
    ...

@log_actions
def main():
    ...

    mpm_results = calculate_mpm(hec.hec_data, D_char)
    mpm_results.to_excel(os.path.abspath("..") + os.sep + "bed_load_mpm.xlsx")

Launch and Debug¶

In your IDE, run the script (e.g., in PyCharm, right-click in the main.py script and click > Run 'main'). If the script crashes or raises error messages, trace them back, and fix the issues. Add try - except statements where necessary and recall the debugging instructions.

A successful run of main.py produces a bed_load_mpm.xlsx file that looks like this:

The logfile should look similar to this:

[20XX-XX-XX 14:08:22,900] PROCESSING PROFILE 1970.1 FOR SCENARIO Q mean
[20XX-XX-XX 14:08:22,900] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,901] PROCESSING PROFILE 1970.1 FOR SCENARIO HQ2.33
[20XX-XX-XX 14:08:22,901] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,901] PROCESSING PROFILE 1970.1 FOR SCENARIO HQ5
[20XX-XX-XX 14:08:22,902] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,902] PROCESSING PROFILE 1970.1 FOR SCENARIO HQ10
[20XX-XX-XX 14:08:22,902] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,902] PROCESSING PROFILE 1970.1 FOR SCENARIO HQ100
[20XX-XX-XX 14:08:22,903] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,903] PROCESSING PROFILE 1893.37 FOR SCENARIO Q mean
[20XX-XX-XX 14:08:22,903] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,903] PROCESSING PROFILE 1893.37 FOR SCENARIO HQ2.33
[20XX-XX-XX 14:08:22,903] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,904] PROCESSING PROFILE 1893.37 FOR SCENARIO HQ5
[20XX-XX-XX 14:08:22,904] Warning: Discharge out of validity range.
[20XX-XX-XX 14:08:22,904] PROCESSING PROFILE 1893.37 FOR SCENARIO HQ10
[20XX-XX-XX 14:08:22,904] Warning: Discharge out of validity range.
[...]