This notebook works out the expected hillslope sediment flux, topography, and soil thickness for steady state on a 4x7 grid. This provides "ground truth" values for tests.
Let the hillslope erosion rate be $E$, the flux coefficient $D$, critical gradient $S_c$, and slope gradient $S$. The regolith thickness is $H$, with bare-bedrock production rate $P_0$ and depth-decay $H_*$. Finally, we set the transport decay scale the same as the production depth-decay scale. Then we have the hillslope flux as a function of distance from ridgetop, $x$, as
$q_s = E x = \left( DS + \frac{D}{S_c^2} S^3 \right) \left(1 - e^{ -H/H_*} \right)$
Parameter values: let $D = 0.01 m^2 y^{-1}$, $S_c = 0.8$, $H_* = 0.5 m$, $P_0 = 0.0002$, and $E = 0.0001 m y^{-1}$:
In [7]:
D = 0.01
Sc = 0.8
Hstar = 0.5
E = 0.0001
P0 = 0.0002
With that, calculate the expected equilibrium $H$:
$E = P_0 e^{-H/H_*}$
$H = -H_* \ln (E/P_0)$
Plugging in the numbers:
In [8]:
import math
H = -Hstar * math.log(E / P0)
H
Out[8]:
Double check: if we plug this $H$ back in, do we recover $E$?
In [9]:
P0 * math.exp(-H / Hstar)
Out[9]:
Yes, good.
Now, our geometry consists of a hillslope discretized into seven nodes. The two on either end are zero-elevation fixed boundaries, so we have to find the elevations of the five interior ones. But the hillslope should be symmetrical, so we really only have to find 1, 2, and 3 as in
0 --- 1 --- 2 --- 3 --- etc.
where node 3 is the top of the hill.
The slope between nodes 1 and 0 must be positive (uphill to right). It must be just steep enough to carry all the sediment from its own cell plus the sediment from node 2's cell, plus half the sediment from node 3's cell. We'll assume all cells have width $dx = 10 m$. Therefore, we have to transport sediment produced in strip 25 m x 1 m, or 25 m2. Our expected flux is then:
In [33]:
qs = 25 * E
qs
Out[33]:
In fact, for each interface between cells, the slope at that interface is given by the following polynomial:
$f\frac{D}{S_c^2} S^3 + 0 S^2 + fDS - qs = 0$
Here the $f$ is shorthand for $1 - \exp (-H/H_*)$. I've included the zero in front of the $S^2$ term just to make it explicit.
So, for the slope between nodes 0 and 1, we need first to define our polynomial coefficients, $p$. Then we'll invoke numpy's roots function to solve for $S$. To be consistent with roots usage, we'll call the coefficient of the highest (cubic) term $p_0$, the next highest (square) $p_1$, etc. So:
$p_0 S^3 + p_1 S^2 + p_2 S + p_3 = 0$
Clearly, we'll need $f$, so let's calculate that first:
In [34]:
f = 1.0 - math.exp(-H / Hstar)
f
Out[34]:
Now, let's calculate the coefficients:
$p_0 = f D / S_c^2$
$p_1 = 0$
$p_2 = f D$
$p_3 = -q_s$
Clearly, only $p_3$ will vary from node to node. Here are the numbers:
In [35]:
import numpy as np
p = np.zeros(4)
p[0] = (f * D) / (Sc ** 2)
p[1] = 0.0
p[2] = f * D
p[3] = -qs
p
Out[35]:
Now let's find the roots of this cubic polynomial:
In [36]:
my_roots = np.roots(p)
my_roots
Out[36]:
There's just one real root here: $S \approx 1.33$. Let's plug that back in and see if we recapture the correct $qs$:
In [37]:
Spred = 0.4
qspred = (D * Spred + (D / (Sc * Sc)) * (Spred ** 3)) * (1.0 - np.exp(-H / Hstar))
qspred
Out[37]:
Great! That's extremely close. Let's try with the slope between nodes 1 and 2. The only difference here is that the flux $qs$ now derives from just $15 m^2$, so $qs = 0.0015:
In [38]:
p[3] = -0.0015
my_roots = np.roots(p)
my_roots
Out[38]:
Once again, let's test:
In [39]:
Spred = 0.269437
qspred = (D * Spred + (D / (Sc * Sc)) * (Spred ** 3)) * (1.0 - np.exp(-H / Hstar))
qspred
Out[39]:
Finally, the slope between 2 and 3, which needs to carry half a cell's worth of sediment, or $qs = 0.0005$:
In [40]:
p[3] = -0.0005
my_roots = np.roots(p)
my_roots
Out[40]:
And check this:
In [42]:
Spred = 0.0985
qspred = (D * Spred + (D / (Sc * Sc)) * (Spred ** 3)) * (1.0 - np.exp(-H / Hstar))
qspred
Out[42]:
Fabulous. Now to find the predicted elevations: just add up slope x distance for each node, going inward from the boundaries:
In [43]:
elev = np.zeros(7)
elev[1] = 0.4 * 10.0
elev[5] = elev[1]
elev[2] = elev[1] + 0.269437 * 10.0
elev[4] = elev[2]
elev[3] = elev[2] + 0.0985 * 10.0
elev
Out[43]:
So, at equilibrium, our model should create a symmetrical hill with a peak elevation a little under 8 m and a soil thickness of 0.347 m.
What time step size would be reasonable? Start by defining an "effective D" parameter, which is the linearized coefficient in front of the cubic term:
$D_{eff} = D (S / S_c)^2$
Then take the steepest steady state slope:
In [44]:
S = 0.4
Deff = D * ((S / Sc) ** 2)
Deff
Out[44]:
Now, maximum time step size should be $\Delta x^2 / 2 D_{eff}$:
In [45]:
10.0*10.0/(2.0*Deff)
Out[45]:
There's also a constraint for the weathering piece. The characteristic time scale is $T = H_* / P_0$, which in this case is:
In [48]:
Hstar / P0
Out[48]:
So, this calculation suggests that weathering is the limiting factor on time-step size. We might choose 250 years for a reasonably smooth solution.
The time it would take for baselevel fall to bring the crest of the hill up to its ten times its equilibrium elevation of 8 m:
In [54]:
80.0 / E
Out[54]:
So let's say we run for 800,000 years at 250 year time steps:
In [55]:
8.0e5/250.
Out[55]:
So, make it 3200 iterations of 250 years each.
In [ ]: