Note

This page was generated from a jupyter notebook.

Using the Landlab BedParcelInitializer components to initialize river bed sediment parcels for the NetworkSedimentTransporter component

This tutorial illustrates how to initialize bed sediment to model the transport of coarse sediment through a river network using the NetworkSedimentTransporter Landlab component.

In this example we will: - load a river network shapefile to create a Landlab grid to represent a river network - initialize sediment ‘parcels’ on that river network using the four available BedParcelInitializer components - conduct a simple exploration of the initial parcel attributes

Import the necessary libraries, plus a bit of magic so that we can plot within this notebook:

[ ]:

# import warnings
# warnings.filterwarnings('ignore')

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

from landlab import ExampleData
from landlab.components import (
    BedParcelInitializerArea,
    BedParcelInitializerDepth,
    BedParcelInitializerDischarge,
    BedParcelInitializerUserD50,
)
from landlab.io import read_shapefile
from landlab.plot import graph, plot_network_and_parcels

%matplotlib inline

1. Load a shapefile that represents the river network

First, we need to create a Landlab NetworkModelGrid to represent the river network. Each link on the grid represents a reach of river. Each node represents a break between reaches. All tributary junctions must be associated with grid nodes.

[ ]:

datadir = ExampleData("io/shapefile", case="redb").base

shp_file = datadir / "a001_network.shp"
points_shapefile = datadir / "a001_nodes_att.shp"

grid = read_shapefile(
    shp_file,
    points_shapefile=points_shapefile,
    node_fields=["usarea_m2", "elev_m"],
    link_fields=["usarea_km2", "Length_m", "Slope"],
    link_field_conversion={
        "usarea_km2": "drainage_area",
        "Slope": "channel_slope",
        "Length_m": "reach_length",
    },
    node_field_conversion={
        "usarea_m2": "drainage_area_m2",
        "elev_m": "topographic__elevation",
    },
    threshold=0.01,
)

Let’s plot our simple network:

[ ]:

graph.plot_graph(grid, at="node,link")

Our network consists of 12 links (grid.number_of_links) between 13 nodes (grid.number_of_nodes). In the plot above, X and Y represent the plan-view coordinates of the node locations.

Next, we need to populate the grid with the relevant topographic and hydrologic information:

[ ]:

grid.at_node["bedrock__elevation"] = grid.at_node["topographic__elevation"].copy()

grid.at_link["channel_width"] = 1 * np.ones(grid.number_of_links)  # m

grid.at_link["flow_depth"] = 0.5 * np.ones(grid.number_of_links)  # m

We must distinguish between topographic elevation (the top surface of the bed sediment) and bedrock elevation (the surface of the river in the absence of modeled sediment).

2. Create initial sediment ‘parcels’ in a DataRecord

We represent sediment in the network as discrete parcels (or packages) of grains of uniform size and characteristics. Each parcel is tracked through the network grid according to sediment transport capacity and stratigraphic constraints.

Parcels are tracked using the Landlab DataRecord.

There are four available components for initiating bed sediment parcels, each of which take different parameters to determine the initial median grain size (d50) of the bed sediment. All initializers return a DataRecord of parcels with a set of parcels with a lognormal distribution of parcel grain sizes. By default, the d84 is 2.1 times d50, though the user can change the with of the grain size distribution.

3. BedParcelInitializerUserD50:

The user can specify either a d50 for each link in the network, or a single d50 for all links.

[ ]:

d50_each_link = np.full(grid.number_of_links, 80.0)
# or
single_d50 = 0.05

initialize_parcels = BedParcelInitializerUserD50(
    grid,
    # user_d50 = d50_each_link,
    user_d50=single_d50,
)

parcels = initialize_parcels()

Many options available for specifying :math:`D_{50}`

[ ]:

date = "2017-11-08 17:50:00"
starting_time = pd.to_datetime(date)

time = [starting_time]

depth = np.full(grid.number_of_links, 1.2)
slope_depend_Shields = 0.15 * grid.at_link["channel_slope"] ** 0.25
tau_c_multiplier = 2.3

initialize_parcels = BedParcelInitializerDepth(
    grid,
    flow_depth_at_link=depth,
    tau_c_50=slope_depend_Shields,
    tau_c_multiplier=tau_c_multiplier,
    median_number_of_starting_parcels=100,
    extra_parcel_attributes=["lithology", "Be10_c"],
)

parcels = initialize_parcels()

parcels.dataset["lithology"].values = np.full(parcels.number_of_items, "quartzite")
parcels.dataset["Be10_c"].values = np.full(parcels.number_of_items, "222.2")

In order to track sediment motion, we classify parcels as either active (representing mobile surface sediment) or inactive (immobile subsurface) during each timestep. The active parcels are the most recent parcels to arrive in the link. During a timestep, active parcels are transported downstream (increasing their location_in_link, which is a normalized value ranging from 0 to 1) according to a sediment transport formula.

We begin by assigning each parcel an arbitrary (and small) arrival time and location in the link.

In addition to the required parcel attributes listed above, you can designate optional parcel characteristics, depending on your needs. We now collect the arrays into a dictionary of variables, some of which will be tracked through time (["item_id", "time"]), and others of which will remain constant through time. With all of the required attributes collected, we can create the parcels DataRecord. Often, parcels will eventually transport off of the downstream-most link. To track these parcels, we have designated a “dummy_element” here, which has index value -2.

4. Plot to explore starting parcels

What grain size did we end up with? How many parcels per link (check)? Etc.

First, how many parcels are on the network?

[ ]:

number_of_pacels = parcels.number_of_items
number_of_pacels

What are the parcel volumes?

[ ]:

print(
    "The largest parcel is",
    np.round(np.max(parcels.dataset["volume"].values), 3),
    "m^3",
)
print(
    "The median parcel is",
    np.round(np.median(parcels.dataset["volume"].values), 3),
    "m^3",
)
print(
    "The smallest parcel is",
    np.round(np.min(parcels.dataset["volume"].values), 3),
    "m^3",
)

Next, what is the median grain size at each link?

Note: we’re not weighting the calculation by parcel volume, which is a trivial simplification since nearly all parcels are the same volume

[ ]:

D50 = parcels.calc_aggregate_value(xr.Dataset.median, "D", at="link", fill_value=0.0)

fig = plot_network_and_parcels(
    grid,
    parcels,
    link_attribute=D50,
    link_attribute_title="Bulk parcel D50 (m)",
    network_linewidth=5,
    parcel_alpha=0,
)

What is the volume of sed originating on each link?

[ ]:

vol_on_link = parcels.calc_aggregate_value(
    xr.Dataset.sum,
    "volume",
    at="link",
    # filter_array=parcelfilter, # no need to filter to a specific timestep - we just have one
    fill_value=0.0,
)

fig = plot_network_and_parcels(
    grid,
    parcels,
    link_attribute=vol_on_link,
    link_attribute_title="Vol of sed originating on link x",
    network_linewidth=5,
    parcel_alpha=0,
)

What is the relationship between D50 and channel slope?

[ ]:

plt.plot(grid.at_link["channel_slope"], D50, ".")
plt.xlabel("Channel slope")
plt.ylabel("D50 (m)")

Generated by nbsphinx from a Jupyter notebook.