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