Using the Landlab Shared Stream Power Model¶

(Ann Thompson, Arizona State University, July 2024)

This notebook is a tutorial of the SharedStreamPower component, which is an extension of the ErosionDeposion component. SharedStreamPower is a landscape evolution model designed to seamlessly transition between bedrock incision (detachment limited erosion) and sediment transport (transport limited erosion). It is based off of the shared stream power model from Hergarten (2021).

Overview¶

Theory¶

Here is the equation for the shared stream power model:

\[E = k_{bedrock}A^mS^n- \frac{k_{bedrock}}{k_{transport}}\frac{Q_s}{A}\]

where \(k_{bedrock}\) is the erodibility with no sediment and \(k_{transport}\) is the ability to transport sediment. The ratio, \(\frac{k_{bedrock}}{k_{transport}}\) determines how much relative bedrock incision and sediment transport is occurring.

For \(\frac{k_{bedrock}}{k_{transport}} = 0\), the model is entirely detachment limited, and behaves as the stream power model:

\[E = k_{bedrock}A^mS^n\]

For \(\frac{k_{bedrock}}{k_{transport}} = \infty\), the model is completely dominated by deposition.

Import Libraries

[ ]:

import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import trange

from landlab import RasterModelGrid
from landlab.components import (
    ChannelProfiler,
    ChiFinder,
    DepressionFinderAndRouter,
    FlowAccumulator,
    SharedStreamPower,
    SteepnessFinder,
)

Make Grid

[ ]:

# Create grid
grid = RasterModelGrid((20, 20), xy_spacing=100.0)

# Leave bottom boundary open
grid.set_closed_boundaries_at_grid_edges(True, True, True, False)


# Create initial topography
np.random.seed(35)  # seed set so our figures are reproducible

# initial noise on elevation grid
# set up the elevation on the grid
grid.at_node["topographic__elevation"] = np.random.rand(grid.number_of_nodes) / 1000.0

Set Parameters

[ ]:

k_bedrock = 0.001  # bedrock erodibility
k_transport = 0.01  # transport coefficient
# For detachment limited behavior, set k_bedrock = 0.001, k_transport = 0.01
# For transport limited behavior, set k_bedrock = 0.01, k_transport = 0.001
# For mixed bedrock-alluvial conditions, set k_bedrock = k_transport = 0.001

F_f = 0.0  # fraction of fines
m_sp = 0.5  # discharge exponent
n_sp = 1.0  # slope exponent
r = 1.0  # m/yr # Define runoff parameter r, where Q=Ar

# time parameters, timestep, run time, print time  (years)
timestep = 10
run_time = 50000
print_time = 10000

# Set elapsed time to zero
elapsed_time = 0

# set uplift rate
rock_uplift_rate = 1e-4  # m/yr

Instantiate Components

[ ]:

# Instantiate flow accumulator and depression finder for each timestep
flow_accumulator = FlowAccumulator(grid, flow_director="D8")
depression_finder = DepressionFinderAndRouter(grid)

# Instantiate components for plotting results
steepness_finder = SteepnessFinder(
    grid, reference_concavity=m_sp / n_sp, min_drainage_area=1000.0
)
chi_finder = ChiFinder(
    grid,
    min_drainage_area=1000.0,
    reference_concavity=m_sp / n_sp,
    use_true_dx=True,
)

# Instantiate the shared stream power component
shared_stream_power = SharedStreamPower(
    grid, k_bedrock=k_bedrock, k_transport=k_transport, m_sp=m_sp, n_sp=n_sp
)

Run

[ ]:

for _ in trange(0, run_time, timestep):
    # Run the flow router
    flow_accumulator.run_one_step()

    # Run the depression finder and router; optional
    depression_finder.map_depressions()

    # Run the SSPM model
    shared_stream_power.run_one_step(dt=timestep)

    # Move  elevation of core nodes upwards relative to base level
    # at the rock uplift rate
    grid.at_node["topographic__elevation"][grid.core_nodes] += (
        rock_uplift_rate * timestep
    )

Plot¶

[ ]:

# Pick individual channels to plot in map view
prf = ChannelProfiler(
    grid,
    number_of_watersheds=2,
    main_channel_only=False,
    minimum_channel_threshold=grid.dx**2,
)
prf.run_one_step()

plt.figure(1)
prf.plot_profiles_in_map_view()
plt.show()

# Plot channel profiles
plt.figure(2)
prf.plot_profiles()
plt.show()

[ ]:

# Plot slope-area plots for the channels
plt.figure(3)
plt.figure(figsize=(6, 2))
for i, outlet_id in enumerate(prf.data_structure):
    for j, segment_id in enumerate(prf.data_structure[outlet_id]):
        if j == 0:
            label = f"channel {i + 1}"
        else:
            label = "_nolegend_"
        segment = prf.data_structure[outlet_id][segment_id]
        profile_ids = segment["ids"]
        color = segment["color"]
        plt.loglog(
            grid.at_node["drainage_area"][profile_ids],
            grid.at_node["topographic__steepest_slope"][profile_ids],
            ".",
            color=color,
            label=label,
        )

plt.legend(loc="lower left")
plt.xlabel("drainage area (m^2)")
plt.ylabel("channel slope [m/m]")
plt.title("Area vs. Slope")


# calculate normalized channel steepness
steepness_finder.calculate_steepnesses()

# plots of normalized channel steepness in the profiled channels
plt.figure(6)
plt.figure(figsize=(6, 2))
for i, outlet_id in enumerate(prf.data_structure):
    for j, segment_id in enumerate(prf.data_structure[outlet_id]):
        if j == 0:
            label = f"channel {i + 1}"
        else:
            label = "_nolegend_"
        segment = prf.data_structure[outlet_id][segment_id]
        profile_ids = segment["ids"]
        distance_upstream = segment["distances"]
        color = segment["color"]
        plt.plot(
            distance_upstream,
            grid.at_node["channel__steepness_index"][profile_ids],
            "x",
            color=color,
            label=label,
        )

plt.xlabel("distance upstream (m)")
plt.ylabel("steepness index")
plt.legend(loc="lower left")
plt.title("Distance Upstream vs. Ksn")

# Plot drainage area vs. sediment flux
plt.figure(7)
plt.figure(figsize=(6, 2))
plt.scatter(
    grid.at_node["drainage_area"],
    grid.at_node["sediment__flux"],
    marker="o",
    color="y",
)

plt.xlabel("drainage area (m^2)")
plt.ylabel("(sediment flux m^3/s)")
plt.title("Area vs. Sediment Flux")

Run with Transient Uplift¶

To observe the landscape response to increased uplift, we set a new uplift rate, run for 10,000 years, and plot the normalized channel steepness at every 1000 years.

Set new parameters

[ ]:

rock_uplift_rate = 0.001  # increased by a factor of     0
time_interval = 1000  # interval between each plot
run_time = 10000  # total run time
elapsed_time = 0  # reset elapsed time to 0

Run loop at plot at each interval

[ ]:

for elapsed_time in trange(0, run_time, timestep):
    if elapsed_time % time_interval == 0:  # if time interval is reached, plot
        prf.run_one_step()
        steepness_finder.calculate_steepnesses()
        plt.figure(6)
        plt.figure(figsize=(6, 2))
        for i, outlet_id in enumerate(prf.data_structure):
            for j, segment_id in enumerate(prf.data_structure[outlet_id]):
                if j == 0:
                    label = f"channel {i + 1}"
                else:
                    label = "_nolegend_"
                segment = prf.data_structure[outlet_id][segment_id]
                profile_ids = segment["ids"]
                distance_upstream = segment["distances"]
                color = segment["color"]
                plt.plot(
                    distance_upstream,
                    grid.at_node["channel__steepness_index"][profile_ids],
                    "x",
                    color=color,
                    label=label,
                )

        plt.xlabel("distance upstream (m)")
        plt.ylabel("steepness index")
        plt.legend(loc="lower left")
        plt.title(f"Steepness index at t = {elapsed_time}")

    # Run the flow router
    flow_accumulator.run_one_step()

    # Run the depression finder and router; optional
    depression_finder.map_depressions()

    # Run the SSPM model for one timestep
    shared_stream_power.run_one_step(dt=timestep)

    # Move  elevation of core nodes upwards relative to base level
    # at the rock uplift rate
    grid.at_node["topographic__elevation"][grid.core_nodes] += (
        rock_uplift_rate * timestep
    )

References¶

Hergarten, S. (2021). The influence of sediment transport on stationary and mobile knickpoints in river profiles. Journal of Geophysical Research: Earth Surface, 126, e2021JF006218. https://doi.org/10.1029/2021JF006218

Generated by nbsphinx from a Jupyter notebook.