1d Sediment Transport#
Goals
This exercise features the application of the Meyer-Peter & Müller (1948) bedload transport formulae to a valid application: 1d, cross-section averaged hydraulics. Write object-oriented code with custom classes for tailored interactions with xlsx workbooks. The homework involves built-in methods of Pandas DataFrames and plotting.
Requirements
Read and understand data handling with NumPy and Pandas as well as Object Orientation and Classes.
Get ready by cloning the exercise repository:
git clone https://github.com/Ecohydraulics/Exercise-SedimentTransport.git
Fig. 26 The Arbogne River in Switzerland (source: Sebastian Schwindt 2013).#
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 [USACoEngineeers16], 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. 27). 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 [Ein50].
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 [Shi36] is often used as the threshold value for the mobilization of sediment from the riverbed (see Fig. 28). This exercise uses one of the dimensionless bed shear stress approaches and the next section provides more explanations.
Fig. 28 The principle of sediment mobilization.#
The Meyer-Peter and Müller (1948) Formula#
The Meyer-Peter and Müller [MPM48] 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 and Müller [MPM48] is:
Bed shear stress
\(\tau_{x,cr}\) \(\approx\) 0.047 (up to 0.07 in mountain rivers), and
\(\tau_{x}\) = \(R_{h} \cdot S_{e} / [(s - 1) \cdot D_{char}]\)
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 and Recking [RR11]).
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 [Buc15]. 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):
Dimensional unit bedload
The cross-section averaged bedload \(Q_{b}\) (kg/s) is then:
Dimensionless cross-section-averaged bedload
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.pywill contain logging functions.hec.pywill contain aHecSetclass to read hydraulic output data from HEC-RAS as structured objects.grains.pywill contain aGrainReaderclass to read grain size class information as structured objects.bedload.pywill contain the classBedCorewith basic elements that most bedload formulae have in common.mpm.pywill contain the classMPM, which inherits fromBedCoreand 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:
start_loggingto setup logging formats and a log file name as described in the section on logging, andlog_actions, which is a function wrapper for themain()(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"])
Challenge
Add a __call__() method to the GrainReader class.
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:
Alphabetic Col. |
Variable |
Type/Unit |
Description |
|
|---|---|---|---|---|
Col. 01 |
A |
Reach |
string |
River (reach) name |
Col. 02 |
B |
River Sta |
[m] |
Position on the longitudinal river axis |
Col. 03 |
C |
Profile |
string |
Name of flow scenario profile (e.g., HQ2.33) |
Col. 04 |
D |
Q Total |
[m³/s] |
River discharge |
Col. 05 |
E |
Min Ch El |
[m a.s.l.] |
Minimum elevation (level) of channel cross-section |
Col. 06 |
F |
W.S. Elev |
[m a.s.l.] |
Water surface elevation (level) |
Col. 07 |
G |
Vel Chnl |
[m] |
Flow velocity main channel |
Col. 08 |
H |
Flow Area |
[m²] |
Wetted cross-section area A (see above) |
Col. 09 |
I |
Froude# Chl |
[-] |
Froude number of the channel (if 1, computation error - do not use!) |
Col. 10 |
J |
Hydr Radius |
[m] |
Hydraulic radius |
Col. 11 |
K |
Hydr Depth |
[m] |
Water depth (active cross-section average) |
Col. 12 |
L |
E.G. Slope |
[m/m] |
Energy Gradeline slope |
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 [WP06] or the Smart and Jaeggi [SJ83] (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)
Note
Import fun (the script with logging functions) to enable the usage of logging.warning(...) messages in the methods of BedCore and its child classes.
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.')
Note
The here shown check_validity method takes the Froude number as input argument. Alternatively, assign the Froude number already in __init__ and use self.Fr.
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.
Note
The program intentionally produces warning messages because some of the profile characteristics do not fulfill the Meyer-Peter & Müller formula’s validity range.
A successful run of main.py produces a bed_load_mpm.xlsx file that looks like this:
River Sta |
Scenario |
Q (m3/s) |
Phi (-) |
Qb (kg/s) |
|
|---|---|---|---|---|---|
0 |
1970.1 |
Q mean |
1 |
||
1 |
1970.1 |
HQ2.33 |
13 |
0.548377243 |
42.72291418 |
2 |
1970.1 |
HQ5 |
17 |
0.682792055 |
54.58338633 |
3 |
1970.1 |
HQ10 |
19 |
0.765834516 |
62.56010505 |
4 |
1970.1 |
HQ100 |
25 |
0.905542967 |
77.92848176 |
5 |
1893.37 |
Q mean |
1 |
0.193642263 |
5.075423967 |
6 |
1893.37 |
HQ2.33 |
13 |
0.144406226 |
14.00424884 |
7 |
1893.37 |
HQ5 |
17 |
0.203854633 |
20.40484039 |
8 |
1893.37 |
HQ10 |
19 |
0.229078172 |
23.1352098 |
9 |
1893.37 |
HQ100 |
25 |
0.297767546 |
31.25225316 |
… |
… |
… |
… |
… |
… |
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.
[...]
Closing remarks
Sediment transport calculations based on local cross-section averaged hydraulics have extremely high error rates because of higher-order controls [SNMudiagaOjemuH23].
Finally, there are many possible solutions to this exercise and any solution that results in the same outcome (workbook and logfile) is valid. The key challenge is to use an object-oriented approach with at least one class inheriting from another class.
Homeworks
Implement the Parker-Wong correction [WP06] for the Meyer-Peter and Müller [MPM48] formula:\(\Phi_{b,pw} \approx 4.93 \cdot (\tau_{x} - \tau_{x,cr})^{1.6}\). Implement the formula in the
MPMclass either use an optional keyword argument incompute_phior a new method.Use the
openpyxllibrary to add a background color to the headers of output tables.Choose and extract 3 profiles from
mpm_resultsand plot the dimensional bedload transport \(Q_{b}\) (y-axis) against the discharge \(q\) (x-axis).