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.