<< . .

. 12
( : 51)



. . >>

1 km
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s




q = 10
q = ’10 q = 30
10 20 40
0 80 cm

Fig 2.19 (a) Plot of eolian silt thickness versus downwind distance, with analytic solutions for the two-dimensional model for
representative values of the model parameters. (b) Schematic diagram of model geometry. Depositional topography shown in this
example is an inclined plane located downwind of source (but model can accept any downwind topography). (c) Color maps of
three-dimensional model results, illustrating the role of variable downwind topography. In each case, model parameters are
u = 5 m/s, K = 5 m2 /s, and p = 0.05 m/s. Downwind topography is, from left to right, a ¬‚at plane, an inclined plane, and a
triangular ridge. Width and depth of model domain are both 6 km. (d) Color maps of three-dimensional model results for
deposition downwind of Franklin Lake Playa, illustrating the role of variable wind direction. From left to right, wind direction is
θ = ’10—¦ (0—¦ is due south), 10—¦ , and 30—¦ . (e) Map of three-dimensional model results obtained by integrating model results over a
range of wind conditions, weighted by the wind-rose data in Figure 2.18a. For color version, see plate section. Modi¬ed from
Pelletier and Cook (2005).




topography because mass loss from the plume is the spatially variable deposition that occurs as
still assumed to occur along a horizontal plane in the plume intersects with complex topography.
Eq. (2.68). However, by using the elevated plume Two-dimensional (2D) model solutions iden-
concentration -- i.e., c[x, y, h(x, y)] -- to estimate tify a fundamental length scale for downwind de-
the deposition term pc, the model approximates position. Smith (2003) showed that Eq. (2.68) can
2.3 NUMERICAL TECHNIQUES AND APPLICATIONS 57


be approximated by a power-law function along model corresponding to the wind-rose diagram
the plume centerline: of Figure 2.18a. The thickness values mapped in
Figure 2.19e correspond to a uniform Q value of
3/2
Q x
c(x, 0, 0) = √ 1.0 m/kyr (assuming a Qa3 surface age of 50 ka),
(2.70)
2 π pL p Lp
obtained by scaling the model results for differ-
ent values of Q to match the maximum thickness
where L p is a characteristic length scale given by
L p = uK / p2 . Figure 2.19a, for example, presents values between the model and observed data.
The modeled and observed deposition patterns
two plots of Eq. (2.70) with different sets of model
match each other in most respects. However, the
parameters. These parameter sets were chosen
secondary hot spot located on the western pied-
to have different values of K and p but simi-
mont of the study area is not reproduced by
lar values of L p . The similarity of the resulting
the model. This feature could result from more
plots indicates that the pattern of 2D dust depo-
westerly paleowinds transporting more dust east-
sition is not sensitive to particular parameters;
ward from the Amargosa River relative to modern
it is only sensitive to the scaling parameter L p .
conditions. Alternatively, it could result from en-
Therefore, raising the value of K and p simulta-
hanced hillslope erosion of Qa3 surfaces in this
neously, as in Figure 2.19a, leads to solutions that
area.
are equally consistent with observations down-
wind of Franklin Lake Playa. This nonuniqueness
makes model calibration (in 2D or 3D) dif¬cult.
2.3 Numerical techniques and
The model is best calibrated using wind-speed
data to constrain u and K to the greatest ex-
applications
tent possible, because these values can be in-
dependently determined using readily-available
2.3.1 Forward-Time-Centered-Space
data. The deposition velocity, in contrast, can-
method
not be readily determined. By constraining u and
Analytic solutions are elegant but they are often
K , however, model results can then be run for a
only loosely applicable to real problems. Linear
range of values of p to determine the value most
fault scarps and volcanic cones are special cases
consistent with observations.
where the landform symmetry allows for an ana-
Model results are shown in Figure 2.19c to il-
lytic approach because a 2D landform is reduced
lustrate the role of complex downwind topogra-
to 1D. Numerical solutions are generally needed
phy in controlling deposition patterns. In each
in order to model the evolution of landforms that
case, a circular source with uniform Q was con-
sidered with model parameters u = 5 m/s, K = lack symmetry or are controlled by complex base-
5 m2 /s, and p = 0.05 m/s. From left to right, the level boundary conditions.
In order to solve partial differential equations
downwind topography was considered to be ¬‚at,
on a computer, we must discretize the equa-
an inclined plane, and a triangular ridge. Rela-
tions in both space and time. Discretization in
tive to ¬‚at topography, the inclined plane leads
space means that we represent the spatial do-
to a short plume, and the triangular ridge diverts
main as a discrete set of grid points (e.g. from
and bisects the plume.
position x1 to xn for 1D problems). Discretization
Long-term dust transport does not correspond
in time means that the solution will be solved for-
to a single wind direction as assumed in the sim-
ward in discrete intervals of time. First, we must
pli¬ed cases of Figure 2.19c. The effects of multi-
determine how to best represent ¬rst, second,
ple wind directions can be incorporated into the
and higher derivatives in space within this dis-
model by integrating a series of model runs with
crete framework. The starting point is the Taylor
a range of wind directions, each weighted by the
expansion:
frequency of that wind direction in the wind-
rose diagram. Figure 2.19d illustrates model re-
‚h ( x)2 ‚ 2 h
sults for the Eagle Mountain region with wind h(x + x, t) = h(x, t) + +
x
‚x 2 ‚ x2
directions from θ = ’10—¦ to +30—¦ (θ = 0—¦ is due x,t x,t

+ O (( x) ) 2
South). Figure 2.19e illustrates the integrated (2.71)
58 THE DIFFUSION EQUATION


The simplest formula for the ¬rst derivative is grid-point values at the previous time step. The
obtained by neglecting second-order (( x)2 ) terms FTCS method is limited by the fact that it is sta-
and rearranging Eq. (2.71) to obtain the forward- ble only for small timesteps.
difference approximation: The FTCS method starts with an initial value
0
h i everywhere on the grid. In addition to the ini-
h(x + x, t) ’ h(x, t) ‚h tial condition, boundary conditions on the value
= + O ( x) (2.72)
‚x
x of h or its ¬rst derivative must also be speci¬ed
x,t

at both ends of the grid. These boundary con-
This is called forward difference because we
ditions may be constant or may vary as a func-
look forward to x + x to compute the differ-
tion of time. Equation (2.76) is then applied se-
ence. Equation (2.72) is not unique, however, be-
quentially from one end of the grid to the other
cause we can also calculate the ¬rst derivative
to specify the grid of values at the next time
using a backward difference:
step.
h(x, t) ’ h(x ’ ‚h The FTCS method is generally not useful for
x, t)
= + O ( x) (2.73)
‚x large grids or higher dimensions because very
x x,t
small step sizes must be taken in order to main-
We can also use a centered-difference approxima- tain stability. To see this, consider the sequence of
tion: grid points in Figure 2.20. Initially, the values of
h are zero everywhere on the grid. At time zero,
h(x + x, t) ’ h(x ’ ‚h
x, t)
= + O (( x)2 ) the left-side boundary condition on h is raised to
‚x
2x x,t
a ¬nite value. The right-side boundary condition
(2.74)
is h n = 0. During the ¬rst time step, all of the
L
grid points remain zero except for the second
The second derivative can be calculated by using
grid point from the left. The second derivative
a forward difference to calculate the slope in the
of h, (h 0 ’ 2h 0 + h 0 )/( x)2 is positive, which re-
positive direction and subtracting the backward 0 1 2
sults in a ¬nite value for h 1 . Note that if the time
difference in the negative direction: 1
step is large enough, the value of h 1 can actually
1
h(x + x, t) ’ h(x, t) h(x, t) ’ h(x ’
1 x, t) 0
rise higher than that of h 0 . If this happens, the

x x x model will become unstable and will continue
h(x + x, t) ’ 2h(x, t) + h(x ’ x, t) to diverge from the correct solution in subse-
=
( x)2 quent timesteps. In order to maintain stability,
the timestep t must be less than ( x)2 /2D . The
‚ 2h
= + O (( x)2 ) (2.75)
‚ x x,t
2 time required for diffusion to propagate a dis-
tance L is ≈ L 2 /D . Therefore, the number of time
The discretization of the diffusion equation, steps required to model diffusion is ≈ L 2 /( x)2 .
using Eq. (2.75) and introducing abbreviated dis- To model the diffusion equation in 1D on a grid
crete notation h(i x, n t) = h in , is given by with 1000 grid points, therefore, requires a min-
imum of about one million time steps. Once we
h i+1 ’ 2h in + h i’1
h in+1 ’ h in n n
=D move to 2D and 3D problems with tens or hun-
( x)2
t
dreds of thousands of grid points, one million
Dtn
= h in + (h ’ 2h in + h i’1 )
h in+1 n
(2.76) time steps becomes computationally unfeasible
( x)2 i+1
on even the fastest computers. Another related
This Forward-Time-Centered-Space (FTCS) method disadvantage of the FTCS method is that infor-
is a useful means for solving the 1D diffusion mation travels very slowly from one end of the
equation for small grids. This scheme is numer- grid to the other. Figure 2.20 illustrates that it
ically stable provided that the time step is less would take 1000 time steps for the grid points
than or equal to ( x)2 /2D . Equation (2.76) is also on the right side of the grid to feel any effect
known as an explicit scheme because the value of the boundary condition on the left side of
of each grid point is an explicit function of the the grid. This slow propagation of information
2.3 NUMERICAL TECHNIQUES AND APPLICATIONS 59



(b)
(a)
h h
unstable, ”t too large

h0
h0
h1
h1
stable



x
x
implicit method:
explicit method:



Fig 2.20 Illustration of the difference between (a) explicit
Eq. (2.45)). Substituting Eq. (2.45) into Eq. (2.2)
numerical schemes and (b) implicit schemes. In explicit
gives
schemes, grid values at the previous time steps are used to
solve for values at the next time step. In the example, the
left-side boundary, h 0 , is abruptly raised at t = 0. In the ‚h ‚ 2h ‚h
0
=κ x 2 + (2.77)
explicit scheme, only the ¬rst grid point to the right, h 1 will
‚t ‚x ‚x
1
feel the effects of this change at the boundary in the following
time step. If the value of the time step is chosen to be too
Here we consider the slope-wash-dominated ver-
large, the explicit scheme can overcorrect and make h 1 larger
1
than h 1 . In the implicit scheme, the effect of raising the left sion of instantaneous base-level drop of an ini-
0
boundary is simultaneously felt throughout the grid no matter tially ¬‚at terrace considered in Section 2.2.5.
how large it is, resulting in a more stable solution. Figure 2.21a compares hillslope pro¬les ob-
tained for the nonuniform diffusion equation
(Eq. (2.77)) with those of the classic diffusion
is one reason why the time step must remain equation (Eq. (2.20)). For early times, the two
so low. models show no difference. For later times, how-
So, if the FTCS is so poor, why use it at all? ever, hillslope erosion occurs more rapidly in the
The answer is that it is easy to program and, in nonuniform diffusion case, but otherwise there
some ways, more versatile than other schemes we is little difference in the pro¬le shape.
will discuss. In the next two sections we utilize These results are broadly consistent with
the FTCS method to solve for the evolution of those of Carson and Kirkby (1972), who compared
hillslopes with spatially variable diffusivity and hillslope pro¬les for uniform and nonuniform
nonlinear sediment transport. cases. Figure 2.21b is based on the classic ¬g-
ure from their book, illustrating the signature
hillslope forms for different hillslope processes.
2.3.2 Evolution of hillslopes with In this ¬gure, Carson and Kirkby argued that
spatially-varying diffusivity slopes governed by the uniform diffusion equa-
tion (m = 0 and n = 1 in their terminology) were
On hillslopes where slope wash is the predomi-
nant transport process, the classic diffusion equa- characteristically convex downward. In contrast,
the slope-wash-dominated case (i.e. m = 1 and
tion cannot be used. In such cases it is appropri-
n = 0) was characterized by a straight pro¬le with
ate to model hillslope sediment ¬‚ux as increas-
ing linearly with distance from the divide (i.e. no curvature.
60 THE DIFFUSION EQUATION



1.0
kt =
1.0
20 m2
0.8
0.8


kt = 0.6
0.6 h
h
1280 m2 h0
h0 kt =
1280 m2
0.4
0.4

kt =
kt =
q= k 0.2
0.2 2560 m2
2560 m2 kx
q= L
0.0
0.0
25
0 50
25
0 50
x (m)
x (m)
(a) (b)
kx
slope wash
1.0
q= L
without gulying
l


h
q= k
soil creep
h0
0.5

slope wash
with gullying

0.0
0.5
0 1.0
x/L
(c)
where n is a constant. Substituting Eq. (2.78) into
Fig 2.21 (a) Comparison of hillslope evolution of an initially
Eq. (2.2) gives
horizontal terrace with instantaneous base-level drop for
⎡ ¤
homogeneous diffusion (i.e. ¬‚ux is given by Eq. (2.2)) and
‚h ‚h/‚ x n|‚h/‚ x|n
⎢ ⎥
inhomogeneous diffusion (i.e. ¬‚ux is also proportional to
=κ ⎣1 + ¦
n n
‚t
distance from the divide, or Eq. (2.45)). The two evolution
1 ’ |‚h/‚ x| ‚h/‚ x|
Scn 1 ’
Sc Sc
equations yield almost identical hillslope forms, but the rate
of erosion is higher in the inhomogeneous case if the same κ (2.79)
value is used in both cases. (b) Characteristic slope pro¬les
Andrews and Hanks (1985) proposed using Eq.
based on the characteristic forms of Carson and Kirkby
(2.79) with n = 1, while Andrews and Bucknam
(1972). These authors argued that uniform and nonuniform
(1987), Roering et al. (1999), and Mattson and
diffusion result in distinctly different hillslope forms, in
Bruhn (2001) assumed n = 2. Appendix 1 provides
contrast to the model results in (a).
a program to solve for hillslope evolution gov-
erned by Eq. (2.79).
2.3.3 Evolution of hillslopes
Figure 2.22 presents the nonlinear hillslope
with landsliding
response (Eq. (2.79)) to instantaneous base-level
The literature on scarp analysis has included a
drop with an initially planar terrace. The ini-
vigorous debate on linear versus nonlinear diffu-
tial phase of the model is an instantaneous re-
sion models. For example, Andrews and Bucknam
laxation of the hillslope to the critical angle
(1987) ¬rst proposed that the ¬‚ux of material is
Sc at time t = 0. During the ¬rst 1 kyr of the
proportional to the gradient for small values, but
model (i.e. κt = 1 m2 assuming κ = 1 m2 /s), little
increases towards in¬nity as the gradient reaches
hillslope rounding takes place and the hillslope
a critical value Sc :
wears back to a slightly lower angle. As time pro-
ceeds, downwearing gives way to diffusive round-
‚h/‚ x
q = ’κ ing of the hillslope crest. The rate of hillslope
(2.78)
n
|‚h/‚ x|
1’ evolution decreases rapidly with the overall relief
Sc
2.3 NUMERICAL TECHNIQUES AND APPLICATIONS 61


evolution is very different from that of the linear
1.0
kt = 0
kt = 0 diffusion equation. Appendix 1 provides a pro-
1 m2 gram to solve for hillslope evolution governed by
0.8 10 m2
Eq. (2.79). This program may be used as a tem-
100 m2
plate to solve hillslope evolution problems with a

<< . .

. 12
( : 51)



. . >>