Note
This page was generated from a jupyter notebook.
Modeling Hillslopes and Channels with Landlab¶
The original version of this exercise was donated by Andy Wickert at the University of Minnesota. This notebook was created by Nicole Gasparini at Tulane University.
For tutorials on learning Landlab, click here: https://github.com/landlab/landlab/wiki/Tutorials
What is this notebook?
This notebook illustrates a landscape evolution model in which the landscape evolves according to the equation:
Here, \(K\) are coefficients on the fluvial (\(sp\)) and hillslope (\(hs\)) parts of the equation, and \(m_{sp}\) and \(n_{sp}\) are positive exponents, usually thought to have a ratio, \(m_{sp}/n_{sp} \approx 0.5\). \(A\) is drainage area and \(S\) is the slope of steepest descent (\(-\frac{dz}{dx}\)) where \(x\) is horizontal distance (positive in the downslope direction) and \(z\) is elevation. (If slope is negative there is no fluvial erosion.) \(U\) is an externally-applied uplift field.
The first term on the right hand side of the equation is the fluvial erosion term, which is also known as the stream power equation. The second term on the right hand side of the equation is elevation changes via linear diffusion, and linear diffusion is one way in which to describe hillslope sediment transport.
For more information on the fluvial erosion term, please see:
Whipple, K.X. and Tucker, G.E., 1999. Dynamics of the stream‐power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs. Journal of Geophysical Research: Solid Earth.
For more information on linear diffusion applied to hillslopes (and other fun hillslope models) see:
Roering, J.J., 2008. How well can hillslope evolution models “explain” topography? Simulating soil transport and production with high-resolution topographic data. Geological Society of America Bulletin.
The ideas behind what this notebook does are presented nicely in the two papers below. Neither of them is exactly the same as this notebook, but they discuss drainage density and transitions from hillslope to fluviall processes.
Tucker, G.E. and Bras, R.L., 1998. Hillslope processes, drainage density, and landscape morphology. Water Resources Research.
Perron, J.T., Kirchner, J.W. and Dietrich, W.E., 2009. Formation of evenly spaced ridges and valleys. Nature.
What will you do?
In this exercise you will modify the code to get a better understanding of how different processes and forces control landscape evolution, landscape form and drainage density (as interpreted from slope-area data). It is expected that you have already learned the basics about fluvial bedrock incision (and the stream power equation) and sediment transport through creep on hillslopes (and the diffusion equation). (See references above.)
Start by sequentially running each code block without changing anything. To run an individual code cell, put the cursor in the cell and type shift-enter, or got to the Cell pulldown menu at the top and choose Run Cells. At the end of the notebook you will see the questions that you need to answer by changing parts of the code and rerunning it.
Remember that you can always go to the Kernel pulldown menu at the top and choose Restart & Clear Output or Restart & Run All if you change things and want to start afresh. If you just change one code block and rerun only that code block, only the parts of the code in that code block will be updated. (E.g. if you change parameters but don’t reset the code blocks that initialize run time or topography, then these values will not be reset.)
Questions to answer before starting this assignment.
If hillslope diffusivity (\(K_{hs}\)) is fixed, but fluvial erodibility (\(K_{sp}\)) increases, what do you think will happen to the total relief and drainage density of a landscape?
If fluvial erodibility (\(K_{sp}\)) is fixed but hillslope diffusivity (\(K_{hs}\)) increases, what do you think will happen to the total relief and drainage density of a landscape?
If parameters stay fixed (\(K_{hs}\) and \(K_{sp}\)), but the uplift rate increases, what do you think will happen to the total relief and drainage density of a landscape?
Now on to the code…
First we have to import the parts of Python and Landlab that are needed to run this code. You should not have to change this first code block.
[ ]:
# below is to make plots show up in the notebook
%matplotlib inline
[ ]:
# Code Block 1
import numpy as np
from matplotlib import pyplot as plt
from landlab import RasterModelGrid, imshow_grid
from landlab.components import FlowAccumulator, LinearDiffuser, StreamPowerEroder
Now we set parameters
This part you will need to change for the different exercises.
Note that Landlab does not impose units, but it assumes that all units are consistent. We will assume that everything is given in meters (m) and years (yr).
[ ]:
# Code Block 2
uplift_rate = 0.001 # [m/yr], initially set at 0.001
K_sp = 1.0e-5 # units vary depending on m_sp and n_sp, initially set at 1e-5
m_sp = 0.5 # exponent on drainage area in stream power equation, initially 0.5
n_sp = 1.0 # exponent on slope in stream power equation, initially 1.
K_hs = 0.05 # [m^2/yr], initially 0.05
Make a grid.
This part you may want to change.
[ ]:
# Code Block 3
ncells_side = 150 # number of raster cells on each side, initially 150
dxy = 50 # side length of a raster model cell, or resolution [m], initially 50
# Below is a raster (square cells) grid, with equal width and height
mg = RasterModelGrid((ncells_side, ncells_side), dxy)
Set some variables related to time.
[ ]:
# Code Block 4
dt = 1000 # time step [yr], initially 5000
total_time = 0 # amount of time the landscape has evolved [yr]
tmax = 1e6 # time for the model loop to run [yr], initially 1e6
t = np.arange(0, tmax, dt) # each of the time steps that the code will run
Here we make the initial grid of elevation.
[ ]:
# Code Block 5
np.random.seed(0) # seed set to zero so our figures are reproducible
mg_noise = np.random.rand(mg.number_of_nodes) / 1000.0 # intial noise on elevation grid
# set up the elevation on the grid
zr = mg.add_zeros("topographic__elevation", at="node")
zr += mg_noise
Intializing all of the process components that do the work.
[ ]:
# Code Block 6
# intialize flow routing
frr = FlowAccumulator(mg)
# initialize stream power incision
spr = StreamPowerEroder(mg, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0)
# initialize linear diffusion
dfn = LinearDiffuser(mg, linear_diffusivity=K_hs, deposit=False)
Now for the code loop. Note that you can rerun this code block many times, and as long as you don’t rerun any of the code boxes above, it will take the already evolved landscape and evolve it even more.
[ ]:
# Code Block 7
for ti in t:
zr[mg.core_nodes] += uplift_rate * dt # uplift the landscape
dfn.run_one_step(dt) # diffuse the landscape
frr.run_one_step() # route flow
# df.map_depressions()
spr.run_one_step(dt) # fluvial incision
total_time += dt # update time keeper
print(total_time)
Plot the topography.
[ ]:
# Code Block 8
plt.figure(1)
imshow_grid(
mg, "topographic__elevation", grid_units=("m", "m"), var_name="Elevation (m)"
)
title_text = f"$K_{{sp}}$={K_sp}; $K_{{hs}}$={K_hs}; $time$={total_time}; $dx$={dxy}"
plt.title(title_text)
max_elev = np.max(zr)
suptitle_text = "Maximum elevation is " + str(max_elev)
plt.suptitle(suptitle_text)
print("Maximum elevation is ", np.max(zr))
Plot the slope and area data at each point (in log-log space).
[ ]:
# Code Block 9
plt.figure(2)
indices = np.where(mg.status_at_node[mg.at_node["flow__receiver_node"]] == 0)
plt.loglog(
mg.at_node["drainage_area"][indices],
mg.at_node["topographic__steepest_slope"][indices],
"b.",
)
plt.ylabel("Topographic slope")
plt.xlabel("Drainage area (m^2)")
plt.title(title_text)
Has the landscape reached steady state yet?
Answer: Not quite. At perfect steady state, there will be no scatter in the fluvial part of the slope-area relationship (given this model set-up).
What to do.
Answer the following questions using the code above. All solutions should be typed, and supporting figures (produced using the code) should be embedded in your final document. (Download or screenshoot the figures.) You may want to make other plots with the data you collect using this model. You are free to make those plots using whatever software you choose.
In parts of this exercise you need to work with your classmates. You are encouraged to discuss how to use the model and model results with your classmates, however the write-up that you hand in must be your own work.
Hillslope vs. Fluvial Processes. Using the parameters provided in the initial notebook, run the landscape to steady state, or the point at which the topography and the slope-area relationship stop changing (i.e. erosion equals rock uplift). (You can keep rerunning Code Block 7 until steady state is reached. Steady state is reached asymptotically, so exact steady state is less important than very close.) Use the plots of slope and area to estimate where the hillslope–fluvial transition is (or in otherwords, the threshold drainage area for channel heads. There is usually a range of values. You should be consistent in your method to determine drainage density and describe how you determined it in your write-up). Also record the maximum elevation. Now try keeping \(K_{sp}\) the same but increase and decrease \(K_{hs}\) (change and run Code Block 2, then rerun Code Blocks 6 and 7). How do the maximum elevation and the transition from hillslope to channel change with changes in \(K_{hs}\)? Now try keeping \(K_{hs}\) the same but increase and decrease \(K_{sp}\) (change and run Code Block 2, then rerun Code Blocks 6 and 7). How do the maximum elevation and transition from hillslope to channel change with changes in \(K_{sp}\)? You can work in teams with your classmates so that you can explore more parameter combinations. Produce a relationship between the different values of \(K_{sp}\), \(K_{hs}\) and the threshold drainage area and maximum elevation. Remember to run all of your experiments with different parameter values to steady state before estimating the threshold drainage area and maximum elevation. Describe how the different parameters affect drainage density in words, and how this is seen in the relationship that you generate. You do not have to include plots of every single run, but include plots to illustrate at least three landscapes with different drainage densities.
Uplift and erosion. Now, perform a similar set of exercises as you did in exercise 1, but also systematically vary uplift rate (Code Block 2). Work in teams, and each person should choose two combinations of \(K_{sp}\) and \(K_{hs}\) and three uplift rates (for a total of 6 runs). Make sure the parameter values that you choose do not overlap with your group members. Make sure you document the transition from hillslope to fluvial process (make sure all of the team members are using the same method to determine threshold area for drainage density), and also note the maximum steady-state elevation for each combination of uplift, \(K_{sp}\) and \(K_{hs}\). Produce relationships to show how the area threshold and maximum elevation change with the different variables. Describe how uplift rate affects drainage density in words, and how this is seen in the relationship that you generate. You do not have to include plots of every single run, but include some plots to illustrate the changes that you describe. (Note whom your group members were in your write-up.)
Free-form exploration. (Optional) Try changing the grid type (Code Block 3), grid size (Code Block 3), stream power exponents (Code Block 2), distribution of uplift rate (e.g., what happens if you have just part of the landscape experience uplift, a bit trickier, as uplift in Code Block 2 will need to be changed to an array), etc. Based on what you observe, create a consistent geomorphic history of the system. Creativity is expected here!
Final reflection. Was your initial insight into how parameters would affect the landscape correct? Discuss in less than 5 sentences.