Note
This page was generated from a jupyter notebook.
Unit tests and parameterization for AreaSlopeTransporter#
(Greg Tucker, University of Colorado Boulder, October 2022)
This notebook gives background on the default parameters for the AreaSlopeTransporter component, and also provides the calculations for a set of unit tests that are used in the code.
Parameters for the Einstein-Brown transport equation#
The AreaSlopeTransporter component calculates fluvial volumetric sediment flux, \(Q_s\), from a generic area-slope power law:
This section works out the transport parameters \(k_s\), \(m\), and \(n\) for the Einstein-Brown law, based on Willgoose’s 1991 appendix but adding the extra step of integrating across a channel of width \(W = a Q^b\) to get total volumetric transport rate. These parameters are used as the default in the component.
Converting to area-slope form#
(This follows Willgoose et al. (1991a WRR) appendix):
Dimensionless transport rate:
Definition of dimensionless xport rate:
where \(M_s\) sediment discharge in mass/time/width, \(s\) is sed specific grav, \(D\) is grain diameter.
\(\tau_*\) is the usual:
The weirdo factor is:
From plenty of other derivations, we can cast \(\tau_b\) in terms of discharge and slope. Here’s one version using the Manning equation for a wide channel:
Check the units here… we have L\(^{9/5 - 3/5 - 1/5}\) T\(^{-3/5+3/5}\)…ok.
Plug in:
Simplify
So then,
which converts to:
Check units again: L\(^{-3/5 - 3 + 3.6 = 0}\), T\(^{9/5 - 9/5 = 0}\).
Now the annoying \(\phi\):
Now it’s time to bring in the pesky channel-width issue:
Using \(b=1/2\), we have:
Let’s now cast in terms of total flux, again with \(b=1/2\):
That’s instantaneous. Let’s factor in an intermittency factor, \(I\), and combine pieces,
For \(F_1\), we need a grain diameter, densities, and kinematic viscosity. Using 0.05 m for grain diameter,
Next, we’ll convert this into an area-slope formulation. Let’s assume \(Q_{bf} = RA\) (not ideal, given common sub-linear \(A\) scaling, but let’s run with it). Then,
where
Estimating a value for \(k_s\)#
Ok, now for some numbers. Let’s use \(n=0.01\) in SI units. For bankfull width, if we use \(b = 1/2\) and \(a = 5\) s\(^{1/2}\)/m, and a bankfull runoff rate of 10 m/y, what would we have for width?
[ ]:
import numpy as np
a = 5.0
b = 0.5
area = np.array([1.0e6, 10.0e6, 100.0e6, 1000.0e6])
discharge = 1.0e-7 * area
width = a * discharge**b
print("At a bankfull discharge of:", discharge, "cms...")
print("...the corresponding bankfull width would be:", width)
Note that this is just a rough guess at \(a\). One should ideally interrogate hydraulic geometry, but here the purpose is simply to identify a value that is basically reasonable.
Here’s one example. Boulder Creek at Orodell has an annual peak (snow-melt) flow 500 cfs. How wide would it be at a = 1, 5, 10? Drainage area is 102 mi\(^2\), or 264 km\(^2\). The channel is roughly 20 m wide here (though possibly artifically narrowed by the road). If we assume the annual maximum flow is the bankfull event, how wide would the channel be if the coefficient were 5 and the exponent 1/2?
[ ]:
Qbf = 500 * (0.3048**3) # bankfull flow, converted to cms
w = 5 * Qbf**0.5 # predicted width, m
am2 = 102 * (5280 * 0.3048) ** 2 # while we're at it, convert drainage area to m2
print("Peak discharge (cms):", Qbf)
print("Predicted width (m):", w)
print("Drainage area (km2):", am2 / 1e6)
print("Runoff rate (m/s):", Qbf / am2)
print("Runoff rate (m/y):", Qbf * 3600 * 365.25 * 24 / am2)
The above turns out to be quite close to the estimated width of 20 m. Although this is just one example, and not a proper fit by any means, for our purposes it suffices to demonstrate that an \(a\) value of 5 (in SI units) is not completely crazy.
For the sake of estimating a default value for \(k_s\), we’ll choose values for the kinematic viscosity of water (which is temperature-dependent; the value \(10^{-6}\) Pa\(\cdot\)s applies to 20 \(^\circ\)C), grain diameter (5 cm), and sediment immersed specific gravity (1.65, representing quartz). Clearly this means that the coefficient value depends in particular on grain diameter, and to a lesser extent on grain density.
Next, the Einstein-Brown \(F_1\) parameter:
[ ]:
nu = 1.0e-6 # kinematic viscosity of water, m2/s
D = 0.05 # grain diameter, m
sm1 = 1.65 # immersed specific gravity
g = 9.8
term = (36.0 * nu**2) / (g * D**3 * sm1)
F1 = ((2.0 / 3.0) + term) ** 0.5 - term**0.5
print("F1 =", F1)
This allows us to calculate \(k_{sq}\), here with \(n=0.01\) and \(a=5\),
[ ]:
n = 0.01 # Manning roughness coefficient in SI units
a = 5.0 # width-discharge factor, SI
ksq = a * F1 * g**0.5 * ((40.0 * n**4.5) / (sm1**2.5 * D**1.5 * a**4.5))
print("ksq", ksq)
Next, we’ll calculate the corresponding \(k_s\) by applying a bankfull runoff rate (2 m/y) and an intermittency factor (0.01):
[ ]:
R = 2.0 # bankfull runoff rate
intermittency = 0.01 # intermittency factor
ks = ksq * intermittency * R**0.9
print("ks", ks)
For landscape evolution, it makes sense to compute \(Q_s\) in m\(^3\)/y rather than m\(^3\)/s, so we will factor a unit conversion into the coefficient:
[ ]:
# Convert ks to time units of years:
ksy = ks * 365.25 * 24 * 3600.0
print("ksy", ksy)
If we round this off, we have our default value for \(k_s\): 0.0055, with length units of m and time units of years.
Now we’re ready to do a sanity check by examining the slope-area relationship. Let sediment flux be erosion rate times drainage area (ignoring fines and dissolved load), \(Q_s = EA\). Equilibrium means:
[ ]:
print(1 / 2.1)
print(0.4 / 2.1)
What does this imply for gradients given \(E=10^{-4}\) m/y, and $A:nbsphinx-math:in `(10, 100, 1000) :math:`km^2$?
[ ]:
E = 1.0e-4
A = np.array([10.0e6, 100.0e6, 1000.0e6])
print("S=", (E / ksy) ** 0.48 * A ** (-0.19))
These seem like reasonable values for gradient.
Unit tests#
Default parameters#
Based on the derivation above, the default parameters will be \(k_s = 0.0055\), \(m=1.4\), \(n=2.1\).
Slope-area for one grid cell#
Consider a grid cell 1000 x 1000 m with an uplift rate of \(10^{-4}\) m/y. The sediment input rate is therefore 100 m\(^3\)/y. The resulting outflux should equal this:
[ ]:
Seq = (100.0 / (0.0055 * 1.0e6**1.4)) ** (1.0 / 2.1)
elev_eq = Seq * 1000.0
print("Equilibrium slope:", Seq)
print("Equilibrium elevation:", elev_eq)
# test: elevation at the node should be 1,068 cm when rounded to the nearest cm.
print("Rounded elevation in cm:", int(round(elev_eq * 100)))
Transport capacity for a single cell#
Assume an elevation of 1.0 m, \(\Delta x = 100\) m. Drainage area therefore is 10,000 m\(^2\) and slope is 0.01. Therefore with default parameters, the volumetric sediment discharge should be:
[ ]:
Qs = 0.0055 * (100 * 100) ** 1.4 * (0.01) ** 2.1
print("Qs:", Qs)
print("Qs in liters per year (to nearest liter):", int(Qs * 1000))
Sediment rate of change#
Consider two core nodes flowing in one direction toward a single open boundary. Grid is raster with \(\Delta x = 100\) m. Slope is 0.01. Calculate \(dz/dt\) at each grid node.
[ ]:
A = 1.0e4 * np.array([3, 2, 1]) # drainage area (m2); flow is right to left
Qsout = np.zeros(len(A))
Qsin = np.zeros(len(A))
Qsout[1:] = 0.0055 * A[1:] ** 1.4 * 0.01**2.1 # sediment outflux from each node, m3/y
Qsin[:2] = Qsout[1:] # sediment influx to each node, m3/y
dzdt = (Qsin - Qsout) / 1.0e4 # rate of elevation change, m/y
print(Qsin)
print(Qsout)
print(dzdt)
# test: answer should be array([ 0. , 0.365, 0.138])
print(np.round(Qsout, 3))
# test: answer should be array([ 0.365, 0.138, 0. ])
print(np.round(Qsin, 3))
# test: answer should be array([ -2.26400000e-05, -1.38200000e-05])
print(np.round(dzdt[1:], 8))
run_one_step#
Running one time step of 10,000 years with the above case should lower the nodes as follows:
[ ]:
elev = np.array([0.0, 1.0, 2.0])
elev[1:] += dzdt[1:] * 10000.0 # erode for one time step of 10,000 years
print(elev)
# test: should be array([0. , 0.7736, 1.8618])
np.round(elev, 4)
These unit tests are embedded into the code as doctests.