%load_ext autoreload

Examples of using GRLP

GRLP is a model that solves for the long-profile evolution of a fluvial channel using a built-in semi-implicit solver. In particular, it simulates the evolution of a gravel-bed river that shapes its bed by eroding, transporting, and depositing sediments. The model takes a few basic inputs, and since the numerical solver has been set up implicitly, the model can take long timesteps while remaining stable. As its basic inputs, the model takes a sediment and water input, a spatial domain, an uplift or subsidence rate (or equivalently, other local source or sink of material), and a floodplain width, plus a slew of appropriate constants. The model can also replace some of these with a set of power-law scaling relationships between downstream distance and (1) valley width, (2) drainage area, and (3) discharge. If this approach is taken, it allows an analytical solution, assuming that these variables actually follow a power-law downstream. However, for the numerical solution, these can be an arbitrary functions up and downstream, instead of power laws.

Wickert & Schildgen (2019, ESurf) describe the theory behind this model and the equations that are implemented numerically here.

First, we will present an example where we compare the numerical solution with power-law functions for important variables to the corresponding analytical solution, which will demonstrate some basic outcomes and methods of interacting with the model. The very first thing to do is load the relevant libraries.

# Import numerical and plotting libraries
import numpy as np
from matplotlib import pyplot as plt
from copy import deepcopy

# Import the GRLP module
import grlp

%load_ext autoreload
%autoreload

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

GRLP is an object-oriented model, so next, we instantiate grlp model object,(’lp’) which creates the data structure, and associates methods with lp.

lp = grlp.LongProfile()

To provide lp with the necessary starting values, we set some constants and create the domain, using some handy built-in methods that set a suite of reasonable constants for Earth rivers.

# Set up the x domain
lp.set_x(dx=1000, nx=90, x0=10000)

# S0 is the upstream-end slope that determines the sediment input to the catchment
S0 = 1.5E-2

# Set up a starting set of channel-bed elevations (z) on a uniform slope (S0)
lp.set_z(S0=-S0, z1=0)

# Intermittency: What fraction of the total time is the river experiencing a
# geomorphically-effective flood? This assumes a binary on--off state, common 
# for gravel-bed rivers with floodplains (see Blom et al., 2017:
# https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017JF004213).
# For an evaluation of the intermittency in general and its variability across
# a range of fluvial systems, see Hayden et al. (2021):
# https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL092598
lp.set_intermittency(0.80)

# Set up the number of iterations in semi-implicit solver; 3 is a good default.
lp.set_niter(3)

# Utility functions to assign reasonable constants 
# (defined by Wickert & Schildgen, 2019)
lp.basic_constants()
lp.bedload_lumped_constants()
lp.set_hydrologic_constants()

As noted in the comments, the “intermittency” is a factor that simplifies a river channel’s behavior into times of “shaping its channel” and times of performing no morphological work.

• For a discussion of the intermittency and its appropriateness as an approximation to the full hydrograph for the purposes of fluvial geomorphology, see Blom et al., 2017, GRL

• For a catalog of intermittency values across a wide range of fluvial systems and a broader discussion of its significance, see Hayden et al., 2021, JGR - Earth Surface

Pursuant to our assumption that valley width, discharge, and drainage area vary as power-law functions downstream, we will use handy setting functions that provide these values for the whole x domain.

• Distance downstream is known to scale with drainage area through Hack’s Law (countless publications)

• The drainage-area–discharge relationship is well established, see:

  • 50 streams in PA and NJ, USA: Aron and Miller, 1978

  • All USGS gages, including PR: O’Connor and Costa, 2004

  • Theoretical approach: Sólyom and Tucker, 2004

• The valley width – drainage-area scaling is described in detail by Langston & Temme (2019, GRL), and they note a common exponent of ~0.4 in the power-law relation. Other relevant references include:

  • Tomkin et al. (2003): exponent 0.76

  • May et al. (2013): exponent 0.574

  • Schanz & Montgomery (2016): exponents 0.22, 0.34

# Set up transfer functions between drainage area (A), discharge (Q), and
# valley width (B).

# drainage area: A = k_xA * x**P_xA
lp.set_A(k_xA=1., P_xA = 7/4)

# discharge: Q = k_xQ * x**P_xQ
lp.set_Q(k_xQ = 1.43e-5, P_xQ = 7/4*0.7)

# Valley width: B = k_xB * x**P_xB
lp.set_B(k_xB = 4, P_xB = 0.4)

# Set the uplift rate [m/s]; positive upwards
lp.set_uplift_rate(0)

# Set the base level 
# This is currently redundant, but will be useful later. 
# the set_z function already assumes that z_bl starts at 0.
lp.set_z_bl(0.)

Base level is the downstream boundary: the point (x, z) at the river mouth. set_z_bl sets its elevation (used here); set_x_bl sets its position along the valley; and set_bl(x, z) sets both at once. In v3, base level can also migrate horizontally — for example a shoreline tracking a continental shelf as sea level rises and falls — see the Río Santa Cruz example.

For the last part, we will set the input sediment discharge, which arrives in the upstream (in our case left-most) model node.

# Input sediment discharge: this is set based on your defined S0, above.
# (this ficticious boundary-condition slope is the transport slope for the
#  amount of sediment being supplied)

Qs0 = lp.k_Qs * lp.Q[0] * (S0)**(7/6.)

lp.set_Qs_input_upstream(Qs0)

GRLP v3 also lets you force this upstream boundary directly by slope: lp.set_S0(S0) performs exactly the Qs0 = k_Qs * Q[0] * S0**(7/6) conversion we just did by hand and applies it through set_Qs_input_upstream. Use whichever is more natural for your problem — a known sediment supply (set_Qs_input_upstream) or a known slope (set_S0), e.g. an equilibrium slope measured from a DEM.

Let’s see what the initial model domain looks like, just making a plot of the long-profile. You’ll notice that it’s just a flat line, because of lp.set_z(S0=-S0, z1 = 0) above.

# Plot
plt.figure(figsize=(12,6))
plt.plot(lp.x/1000., lp.z, '0.6', linewidth=6, label='Numerical')
plt.xlim(lp.x.min()/1000, lp.x.max()/1000)
plt.ylim(0, lp.z.max())
plt.xlabel('Downstream distance [km]', fontsize=26)
plt.ylabel('Elevation [m]', fontsize=26)
plt.tick_params(axis='both', which='major', labelsize=16)
plt.legend()
plt.tight_layout()
plt.show()

Now, we can evolve the river profile, and calculate the analytical and numerical solutions.

# Numerical: (number of time steps, length of time step [s])
# note that this is a very long time step (31 Ma)
lp.evolve_threshold_width_river(1, 1E15)
# Analytical, no time steps, just solve for equilibrium profile.
zanalytical = lp.analytical_threshold_width() # suppress printing returned data.

We can now plot the analytical solution, and the numerical solution together. They match!

# Plot
plt.figure(figsize=(12,6))
plt.plot(lp.x/1000., lp.z, '0.6', linewidth=6, label='Numerical')
plt.plot(lp.x/1000., lp.zanalytical, 'k', linewidth=2, label='Analytical')
plt.xlim(lp.x.min()/1000, lp.x.max()/1000)
plt.ylim(0, lp.z.max())
plt.xlabel('Downstream distance [km]', fontsize=26)
plt.ylabel('Elevation [m]', fontsize=26)
plt.tick_params(axis='both', which='major', labelsize=16)
plt.legend()
plt.tight_layout()
plt.show()

So what is the geometry of this channel? The model object lp has a lot of attributes and methods. You can query them (like any other python object), using dir().

dir(lp)

['__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattr__',
 '__getattribute__',
 '__getstate__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_net',
 '_network',
 'analytical_threshold_width',
 'analytical_threshold_width_perturbation',
 'analytical_threshold_width_uplift',
 'compute_Q_s',
 'compute_Qs_gain',
 'compute_Qs_lag',
 'compute_Qs_series_terms',
 'compute_diffusivity',
 'compute_e_folding_time',
 'compute_equilibration_time',
 'compute_length',
 'compute_series_coefficient',
 'compute_wavenumber',
 'compute_z_gain',
 'compute_z_lag',
 'compute_z_series_terms',
 'evolve_threshold_width_river',
 'evolve_threshold_width_river_adaptive',
 'segment',
 'set_adaptive_timestep',
 'set_iteration_tolerance',
 'set_niter',
 'set_time_integration',
 'slope_area']

And we can see that the object has an attribute B for valley width, and b for river width. If we query B…

print(lp.B)

[159.24286822 165.43106491 171.29018947 176.86314045 182.18440794
 187.2821968  192.17990944 196.8972089  201.4507986  205.85500554
 210.12222435 214.26326041 218.28759884 222.20361786 226.01876004
 229.739671   233.37231295 236.9220583  240.3937675  243.79185421
 247.12034023 250.38290209 253.58291079 256.72346587 259.80742475
 262.8374282  265.81592238 268.74517822 271.62730829 274.4642817
 277.25793726 280.00999512 282.7220671  285.39566593 288.03221349
 290.63304818 293.19943157 295.73255438 298.23354186 300.70345866
 303.1433133  305.55406215 307.93661309 310.29182888 312.62053016
 314.9234983  317.20147792 319.45517922 321.68528021 323.89242863
 326.07724384 328.24031847 330.38222003 332.50349232 334.60465683
 336.68621393 338.7486441  340.79240894 342.81795224 344.82570089
 346.81606575 348.78944251 350.74621239 352.68674294 354.61138864
 356.52049159 358.41438207 360.29337908 362.15779093 364.00791565
 365.84404154 367.66644754 369.47540366 371.2711714  373.05400409
 374.82414725 376.5818389  378.32730991 380.06078426 381.78247934
 383.49260621 385.19136984 386.87896939 388.55559839 390.22144498
 391.87669211 393.52151774 395.15609502 396.7805925  398.39517423]

but if we query b we get…

print(lp.b)

None

because this information about channels (width, depth, etc.) is computed after the fact, and it requires knowledge of the grain size and sediment discharge everywhere, which either cancels out of our is calculated within the equations for long-profile evolution. The LongProfile class contains methods to calculate these, but because we follow the near-threshold channel-geometry closure of Parker (1978) (Gravel on banks exactly at the threshold of motion at a bankfull, or equivalently in this framework, geomorphologically effective flow), computing the channel width and flow depth (= channel depth) also requires us to know the grain size. Let’s say that this input grain size is 10 cm.

lp.compute_Q_s()
lp.D = 0.1
lp.compute_channel_width()
lp.compute_flow_depth()
print(lp.b)

[0.67600781 0.69714545 0.69714498 0.69721462 0.69733664 0.69749883
 0.69769256 0.6979115  0.69815093 0.69840724 0.6986776  0.69895976
 0.69925191 0.69955254 0.69986039 0.7001744  0.70049367 0.70081738
 0.70114486 0.70147547 0.70180866 0.70214394 0.70248085 0.70281896
 0.70315789 0.70349729 0.70383682 0.70417617 0.70451504 0.70485315
 0.70519024 0.70552606 0.70586037 0.70619293 0.70652353 0.70685196
 0.707178   0.70750147 0.70782217 0.70813992 0.70845453 0.70876584
 0.70907367 0.70937786 0.70967825 0.70997468 0.710267   0.71055506
 0.71083871 0.71111781 0.71139223 0.71166181 0.71192643 0.71218596
 0.71244026 0.71268921 0.71293268 0.71317056 0.71340271 0.71362902
 0.71384938 0.71406366 0.71427175 0.71447355 0.71466894 0.71485781
 0.71504006 0.71521557 0.71538425 0.715546   0.7157007  0.71584826
 0.71598858 0.71612156 0.71624711 0.71636512 0.71647551 0.71657819
 0.71667305 0.71676    0.71683896 0.71690984 0.71697255 0.71702699
 0.71707309 0.71711075 0.7171399  0.71716044 0.7171723  0.71717539]

And we can plot these things too.

fig, axes = plt.subplots(3, 1, figsize=(12,10))
axes[0].plot(lp.x/1000., lp.b, 'xkcd:blue', linewidth=4)
axes[1].plot(lp.x/1000., lp.B, 'xkcd:red', linewidth=4)
axes[2].plot(lp.x/1000., lp.h, 'xkcd:orange', linewidth=4)
axes[2].set_xlabel('Downstream distance [km]', fontsize=16)
axes[0].set_ylabel('Channel Width [m]', fontsize=16)
axes[1].set_ylabel('Valley Width [m]', fontsize=16)
axes[2].set_ylabel('Flow Depth [m]', fontsize=16)
plt.tick_params(axis='both', which='major', labelsize=16)
plt.tight_layout()

This is the basic usage of the model, which is for a given long-profile of a single stream, to calculate the equilibrium profile. Now, there are plenty of situations where the analytical profile would not be expected to match the numerical profile, because there are several simplifying assumptions applied to the analytical solution to render the math tractable.

1. Uplift is assumed to match river incision exactly for the analytical solution.

2. Discharge, Drainage Area, and Valley Width must be specified functions of downstream distance. However, the model does not specify why that should be, only that it is. So if the downstream function of these variables is anything but a power law, the solutions should be different.

So, for fun, let’s violate #1 and see how it goes. First, we will make a copy of the existing object using deepcopy, then we will change the boundary conditions, and run each of them forward.

lp_uplift = deepcopy(lp)
lp_uplift.set_uplift_rate(1e-10) # in [m/s]

lp_uplift.evolve_threshold_width_river(1, 1E15) # evolve for the same amount of time as lp
_ = lp_uplift.analytical_threshold_width()

lp.evolve_threshold_width_river(1, 1E15) # this is unnecessary, but it's not expensive
_ = lp.analytical_threshold_width()

lp.compute_Q_s()
lp.compute_channel_width()
lp.compute_flow_depth()

lp_uplift.compute_Q_s()
lp_uplift.compute_channel_width()
lp_uplift.compute_flow_depth()

fig, axes = plt.subplots(2, 1, figsize=(12,10))
axes[0].plot(lp.x/1000., lp.z, 'xkcd:green', linewidth=8, label = 'Numerical', alpha = 0.3)
axes[0].plot(lp.x/1000., lp.zanalytical, 'xkcd:green', linewidth=2, label = 'Analytical', alpha = 1)
axes[1].plot(lp.x/1000., lp_uplift.z, 'xkcd:purple', linewidth=8, label = 'Numerical', alpha = 0.3)
axes[1].plot(lp.x/1000., lp_uplift.zanalytical, 'xkcd:purple', linewidth=2, label = 'Analytical', alpha = 1)
axes[1].set_xlabel('Downstream distance [km]', fontsize=16)
axes[0].set_title('Uplift = Incision', fontsize = 20)
axes[1].set_title('Uplift ≠ Incision', fontsize = 20)
for a in axes:
    a.set_ylabel('Elevation [m]', fontsize=16)
    a.legend()
    a.tick_params(axis='both', which='major', labelsize=10)
plt.tight_layout()

Violating #2 does other interesting stuff.

# Screw with transfer functions between drainage area (A), discharge (Q), and valley width (B).
# Instead of power laws, these are now either constants, or a linear function 
# (drainage area has to increase downstream --- that's just geometry)

lp_constantQBA = deepcopy(lp)

lp_constantQBA.set_Q(Q = np.repeat(lp.Q[0], len(lp.x)))
lp_constantQBA.set_B(B = np.repeat(lp.B[0], len(lp.x)))
lp_constantQBA.set_A(A = lp.x * np.repeat(lp.A[0], len(lp.x)))

lp_constantQBA.evolve_threshold_width_river(1, 1E15) # evolve for the same amount of time as lp
_ = lp_constantQBA.analytical_threshold_width()

lp_constantQBA.compute_Q_s()
lp_constantQBA.compute_channel_width()
lp_constantQBA.compute_flow_depth()

lp.evolve_threshold_width_river(1, 1E15) # this is unnecessary, but it's not expensive
_ = lp.analytical_threshold_width()

lp.compute_Q_s()
lp.compute_channel_width()
lp.compute_flow_depth()

fig, axes = plt.subplots(2, 1, figsize=(12,10))
axes[0].plot(lp.x/1000., lp.z, 'xkcd:gold', linewidth=8, label = 'Numerical', alpha = 0.3)
axes[0].plot(lp.x/1000., lp.zanalytical, 'xkcd:gold', linewidth=2, label = 'Analytical', alpha = 1)
axes[1].plot(lp.x/1000., lp_constantQBA.z, 'xkcd:dark grey', linewidth=8, label = 'Numerical', alpha = 0.3)
axes[1].plot(lp.x/1000., lp_constantQBA.zanalytical, 'xkcd:dark grey', linewidth=2, label = 'Analytical', alpha = 1)
axes[1].set_xlabel('Downstream distance [km]', fontsize=16)
axes[0].set_title('QBA vary as power-law', fontsize = 20)
axes[1].set_title('QBA are all constants', fontsize = 20)
for a in axes:
    a.set_ylabel('Elevation [m]', fontsize=16)
    a.legend()
    a.tick_params(axis='both', which='major', labelsize=10)
plt.tight_layout()

The river can also evolve through time! Till now, this tutorial uses really big timesteps (300 My), but the model can also resolve shorter time. Here’s a quick example, where an sudden increase in uplift rate is applied to the upstream end of our original channel profile lp.

lp_increase = deepcopy(lp)
fig = plt.figure(figsize=(12,5))
ax1 = fig.add_subplot(1,1,1)
plt.xlabel('Downstream distance [km]', fontsize=14, fontweight='bold')
plt.ylabel('Elevation [m]', fontsize=14, fontweight='bold')
plt.tight_layout()

# Initial condition
ax1.plot(lp.x/1000., lp_increase.z, color='.5', linewidth=3, label = 'Start')

# Transient
U = 1E-3
lp_increase.set_uplift_rate(U/3.15E7)
for i in range(5):
    lp_increase.evolve_threshold_width_river(1, 1E12) # five timesteps of 31 ka
    ax1.plot(lp_increase.x/1000., lp_increase.z, color='.5', linewidth=1)

# New equilibrium
lp_increase.evolve_threshold_width_river(1, 1E14) # final state 3 Ma later
ax1.plot(lp.x/1000., lp_increase.z, color='0', linewidth=3, label = 'Finish')

ax1.legend()

<matplotlib.legend.Legend at 0x71e0c098e870>

Next steps

This tutorial covered a single river long profile. From here:

• River networks — every GRLP solution is a network solution, and a single profile is just the one-edge case. Build tributaries and confluences by hand (NewNetwork_5segments.py), or generate and analyze whole random networks (shreve_random_network.py).

• Full documentation and API reference: grlp.readthedocs.io.