ñòð. 3 |

uniform uplift or base level drop for a hillslope 50 m in length. Q3

Q2c

(b) Same as (a) except that results are plotted in terms of

dimensionless length x /L . Fig 2.4 (a) Hanaupah Canyon alluvial fan, showing

mid-Pleistocene (Q2a) to Holocene (Q3) terraces. (b) Aerial

photographs of Q2a, Q2b, Q2c, and Q3 terraces, illustrating

the correlation between age and degree of roundedness of

2.2.5 Evolution of alluvial-fan terraces the terrace slope adjacent to gullies.

following fan-head entrenchment

In this section we consider the evolution of

initially-planar hillslope terraces subject to rapid aggrading mode, the ratio of sediment to water

base-level drop. This example provides an oppor- is high and sediment is deposited across the fan

tunity to use both series and similarity solu- (or the available accommodation space in the en-

tions within the context of a real-world hillslope- trenched zone of an older abandoned terrace) to

evolution application. form a low-relief deposit by repeated episodes of

First some background. In the Basin and channel avulsion (Figure 2.5a). When the ratio of

Range province of the western United States, sediment to water decreases, such as when cooler

many alluvial fans have a suite of terraces that and wetter climate conditions cause revegetation

rise like a ï¬‚ight of stairs from the active channel. of hillslopes and anchoring of the available re-

Figure 2.4 illustrates the classic alluvial fan ter- golith, the sediment deposited during the pre-

races of Hanaupah Canyon in Death Valley, Cali- vious period of high sediment ï¬‚ux is rapidly en-

fornia. These terrace suites are generally consid- trenched (Figure 2.5b). This entrenchment acts as

ered to be caused by variations in the sediment- a rapid base-level drop for the abandoned terrace

to-water ratio in channels draining the mountain and gullies quickly form at the distal end of the

front (Bull, 1979). When alluvial fans are in an terrace and migrate headward (Figure 2.5c). At

2.2 ANALYTIC METHODS AND APPLICATIONS 39

1.0

(a) 1 (b) 2

kt =

0.8

20 m2

A A 0.6

h

A h0 kt =

1280 m2

0.4

3(d) 4

(c) kt =

0.2

2560 m2

0.0

A A

A A 0 50

x(m)

Fig 2.6 Plots of Eq. (2.42) for Îºt values of 50, 100, 200 . . .

B

6400 m2 and a terrace width of 2L = 100 m. The early phase

of the model is characterized by penetration of erosion into

the scarp interior. After Îºt â‰ˆ 1000 m2 no planar remnant

Fig 2.5 Schematic illustration of a cycle of aggradation and remains and the scarp proï¬le can be approximated by a sine

fan-head entrenchment. When alluvial fans are in an function with height decreasing linearly with time.

aggradational mode, the ratio of sediment to water is high

and sediment is deposited across the available

accommodation space to form a low-relief surface by once the gully has formed. This approach pro-

repeated episodes of channel avulsion (a). When the ratio of vides a relative-age dating and correlation tool for

sediment to water decreases, the sediment deposited during

fan terraces (Hsu and Pelletier, 2004). Series solu-

the previous high-sediment-ï¬‚ux period is rapidly entrenched

tions (Eq. (2.20)) are the most comprehensive ap-

leaving an abandoned terrace (b). This entrenchment acts as a

proach to this problem, but similarity solutions

rapid base-level drop for the abandoned terrace and gullies

can be used to model young terraces (those with

quickly form at the distal end of the terrace and migrate

planar remnants) with no loss of accuracy.

headward (c).

First we consider series solutions. Substituting

f (x) = h 0 into Eq. (2.20) gives

âˆž

any point along a line parallel to the mountain (âˆ’1)n

4h 0

h(x, t) =

front, a proï¬le can be plotted that crosses one Ï€ 2n + 1

n=0

or more gullies. In young terraces, gully proï¬les (2n + 1)Ï€ x Îº(2n+1)2 Ï€ 2 t

âˆ’

Ã— cos e (2.42)

4L 2

are characterized by narrow valleys and broad,

2L

planar terrace treads. In older proï¬les, gully ero-

sion penetrates farther into the terrace. Given Equation (2.42) can be used to model terrace gully

a sufï¬ciently old terrace, hillslope erosion will proï¬les through their entire evolution from back-

completely penetrate the terrace tread and no wearing scarp to downwearing gully remnant.

planar remnant will be preserved. The effects of That two-phase evolution is illustrated in Figure

erosional beveling can be seen in the aerial or- 2.6. The early phase of the model is character-

thophotographs shown in Figure 2.4b. As the ter- ized by penetration of the backwearing scarp into

the terrace interior. After Îºt â‰ˆ 250 m2 , no planar

race age decreases from mid-Pleistocene (Q2a) to

latest Pleistocene--early Holocene (Q3), more pla- remnant remains and the scarp proï¬le can be

nar terraces are preserved and the degree of ter- approximated as a sine function (or cosine func-

tion, depending on where x = 0 is deï¬ned) with

race rounding at gullies decreases markedly.

Equations (2.36) and (2.20) can be used to a linearly-decreasing height over time.

model the evolution of gully proï¬les assuming The results of Figure 2.6 nicely illustrate the

that gullies form rapidly following entrenchment ï¬ltering aspects of diffusion. A series solution

and that the gully ï¬‚oor has a constant elevation is composed of sine functions multiplied by an

40 THE DIFFUSION EQUATION

relatively constant over these time periods. If we

1.0

assume a Îº value of 1 m2 /kyr, the estimated age

Q2c

kt = 50 m2

Q2a of the Q2c terrace is 50 kyr and the Q2b terrace is

0.8

150 kyr. These values are reasonably close to the

kt > 240 m2

estimated ages of these terraces based on avail-

0.6

h able geochronology.

Q2b

h0 Normalization is appropriate in comparing

kt = 150 m2

0.4

diffusion-model predictions to natural landforms

because the relative shape of the hillslope does

0.2 not depend on the absolute height of the hill-

slope. In other words, the landform height can

0.0 Hanaupah Canyon, Death Valley, CA be scaled up or down and the relative hillslope

20

10 30

0 40 50

x (m) shape predicted by the diffusion model will re-

main the same.

Fig 2.7 Gully proï¬les of the Q2a, Q2b, and Q2c fan

In the ï¬eld, the length of scarp backwearing

terraces at Hanaupah Canyon, Death Valley, California.

can be used as a tool to correlate terraces. The

Best-ï¬t morphologic ages to the diffusion proï¬le given by

length of erosional penetration into a scarp can

Eq. (2.42) are 240 m2 , 150 m2 , and 50 m2 , respectively. Only a âˆš

be estimated as Î» = Îºt. This is a general result

minimum age can be determined for Q2a because the height

for all diffusion problems. This relationship im-

of the original terrace is unknown.

plies that Q2c- and Q2b-aged terraces wear back

by approximately 7 and 12 m, respectively, assum-

ing Îº = 1 m2 /kyr.

exponential damping term that is a function of

wavelength. Diffusion acts to ï¬lter the small- The similarity solution can be used to esti-

wavelength components of the topography over mate the age of the Q2c and Q2b terraces because

time through the exponential term in the series planar remnants are preserved on these terraces.

solution. For large values of t, all of the small As such, the terrace width can be considered to

wavelength (i.e. large n) terms will be negligible be semi-inï¬nite. To solve the problem with the

and the only signiï¬cant term will be the largest- similarity method, we consider a semi-inï¬nite

wavelength sine function. This is a general re- terrace of height h 0 subject to rapid base-level

drop at h = L . The solution to this problem is

sult: all base-level-controlled hillslopes will ap-

proximate a sine function during their waning

L âˆ’x

h(x, t) = h 0 erf âˆš (2.43)

stages of evolution, regardless of their initial

2 Îºt

conditions.

Equation (2.43) is identical to Eq. (2.42) for small

Figure 2.7 illustrates normalized surveyed

values of Îºt.

cross sections of the Q2a, Q2b, and Q2c terraces

in Hanaupah Canyon along with their best-ï¬t pro-

2.2.6 Evolution of alpine moraines

ï¬les to Eq. (2.42). The Q2c scarp is best repre-

sented by Eq. (2.42) with Îºt = 50 m2 , while Q2b is Advancing alpine glaciers bulldoze rock and de-

best represented by Îºt = 150 m2 . These Îºt values bris out along their margins. When the glacier

are sometimes referred to as morphologic ages. retreats, a terminal moraine is formed. Terminal

Only a minimum Îºt value can be determined for moraines are initially triangular in cross section.

Q2a because no planar remnant remains and we As such, the solution to the diffusion equation

do not know the height of the original terrace. with a triangular initial condition provides a sim-

Therefore, once the terrace proï¬le is normalized, ple model of moraine evolution.

any Îºt value greater than â‰ˆ 250 m2 matches the In the simplest case of a symmetrical

observed proï¬le equally. These results suggest moraine, only half of the moraine need be con-

that the Q2b terrace on the Hanaupah Canyon sidered. The highest point of the moraine forms a

divide located at x = 0 (Figure 2.8a). At the base of

fan is approximately three times older than the

Q2c terrace, assuming that Îº values have been the moraine two different boundary conditions

2.2 ANALYTIC METHODS AND APPLICATIONS 41

1.0 t=0 (a)

t = 20 kyr (Pinedale)

t = 40 kyr

0.8 t = 60 kyr

...

t = 140 kyr (Bull Lake)

0.6

h

h0

0.4

(c)

moraines

0.2

constant elevation at base of moraine

k = 1m2/kyr

0.0

1.0

(b)

t = 20 kyr (Pinedale) formerly

0.8

glaciated valley

t = 140 kyr (Bull Lake)

0.6

h

h0

0.4

0.2

no local base level control

k = 1 m2/kyr

0.0

0 40

20 30 50

10

x (m)

20 ka) moraine, for example, is rounded at the

Fig 2.8 (a) Diffusion model of moraine evolution plotted

top over a width of approximately 10 m. A Bull

for 20 kyr intervals up to t = 200 kyr, assuming a Îº value of

1 m2 /kyr, with a ï¬xed-elevation boundary condition at the Lake moraine (approx. 140 ka), in contrast, has a

slope base. (b) Same as (a), but with a ï¬xed-elevation rounded crest of approximately 30 m in width as-

boundary condition at x = âˆž. (c) Virtual oblique aerial suming a Îº value of 1 m2 /yr.

photograph of a series of terminal moraines in Owens Valley, Second, we consider the case in which the

Inyo County, California.

moraine deposits sediment at its slope base (such

as when the moraine is surrounded by an allu-

can be considered. First we consider the case of vial piedmont). In this case, sediment shed from

a ï¬xed elevation at a distance L from the divide. the moraine will be deposited at the base of the

As in previous examples, this boundary condition moraine and on the surrounding piedmont. The

is most appropriate for cases in which a stable boundary condition in this case is a ï¬xed eleva-

tion of zero at x = âˆž. For the case of an initially

channel is present at the slope base. Substituting

f (x) = h 0 (L âˆ’ x) into Eq. (2.20) gives the solution triangular hillslope, the solution is

âˆž

8h 0 1

L âˆ’x

h(x, t) = 2 2h 0

h(x, t) = (L âˆ’ x)erf âˆš + (L + x)

Ï€ (2n + 1)2

L

n=0 4Îºt

(2n + 1)Ï€ x Îº(2n+1)2 Ï€ 2 t

L +x

âˆ’ x

Ã— cos e (2.44)

4L 2

Ã— erf âˆš âˆ’ 2xerf âˆš

2L 4Îºt 4Îºt

Îºt (xâˆ’L )2 (x+L )2 x2

Figure 2.8a shows plots of Eq. (2.45) at 20 kyr eâˆ’ + eâˆ’ âˆ’ 2eâˆ’ 4Îºt

+2 4Îºt 4Îºt

Ï€

increments from t = 20 kyr to t = 200 kyr assum-

ing a Îº value of 1 m2 /kyr. A Pinedale-age (approx. (2.45)

42 THE DIFFUSION EQUATION

Figure 2.8b illustrates proï¬les of Eq. (2.45) for and Anderson, 1979; Hanks and Andrews, 1989).

the same ages as in Figure 2.8a. For the case In this approach, the far-ï¬eld gradient (i.e. the

with a depositional slope base the erosion of the gradient of the initially unfaulted surface) must

moraine crest occurs at a slightly lower rate com- ï¬rst be subtracted from the midpoint gradient

pared to the base with a ï¬xed base level. The to take into account the effects of scarps cut into

main effect of the deposition is to widen the sloping surfaces. The plot of midpoint gradient

moraine as sediment progrades out on the ad- minus the far-ï¬eld gradient is compared to char-

jacent slope base. acteristic curves for the diffusion equation with

different parameter values (e.g. Bucknam and

2.2.7 Evolution of pluvial shoreline and Anderson, 1979; Hanks et al., 1984) to determine

fault scarps which parameters best match the observed data.

Both methods are potentially unreliable because

The study of pluvial shoreline and fault scarps

they reduce the information content of the scarp

has had great historical importance in hillslope

to the slope angle or gradient at a single point.

geomorphology. These landforms have been stud-

The scarp midpoint is the least diagnostic point

ied within the context of the diffusion model of

along the entire scarp proï¬le, and it is the point

hillslope evolution for over forty years. In that

that is also most sensitive to uncertainty in the

time, researchers have taken several different ap-

initial scarp angle. As a scarp evolves, the greatest

proaches to analyzing scarps. In this section, we

amount of erosion and deposition occurs where

present the analytic solutions to the diffusion

the magnitude of the hillslope curvature is great-

equation appropriate for pluvial shoreline and

est, near the top and bottom of the scarp. The

fault scarps. In addition, however, we also com-

scarp midpoint is a point of inï¬‚ection with zero

pare and contrast past approaches in order to

curvature, and hence this point changes the least

put the analytic solutions into the historical con-

of all points along the proï¬le for young scarps.

text of other approaches. By comparing different

An alternative approach is to ï¬t the entire scarp

methods for scarp analysis and modeling, we are

(or its derivative, the hillslope gradient) to an-

able to explore some of the subtleties involved in

alytical or numerical solutions of the diffusion

comparing model predictions with natural land-

equation (linear or nonlinear). Avouac and his col-

forms. For example, how do we handle uncer-

leagues (Avouac, 1993; Avouac and Peltzer, 1993)

tainty in the landform initial condition at some

pioneered the use of this technique in fault-scarp

point in the past? How do we handle measure-

studies in Asia. Mattson and Bruhn (2001) ex-

ment uncertainty and random variability in the

tended the technique to include uncertainty in

landscape?

the initial scarp angle.

Pluvial shoreline scarps form by wave-cut ac-

Figure 2.9 presents an oblique aerial view of

tion during a prolonged period of lake-level high

the shoreline in east-central Tule Valley as it cuts

stand. Following lake-level fall, scarps formed in

across several distinct alluvial fan terraces. This

unconsolidated alluvium evolve to an angle of re-

ï¬gure illustrates that scarp height is controlled

pose by mass movements and then further evolve

primarily by the slope of the fan terrace, with

by diffusive hillslope processes.

the tallest scarps formed in steeply dipping al-

Historical scarp-analysis methods come in two

luvial fan terraces and debris cones close to the

basic types: midpoint-slope methods and full-

mountain front.

scarp methods. Mid-point slope methods can be

The geometry of a general scarp is deï¬ned

further divided into two types. In the ï¬rst type,

by the scarp height 2a, initial gradient Î±, and

the midpoint slope of a single scarp is inverted

far-ï¬eld gradient b (Figure 2.10a). The initial and

for diffusion age using the analytic solution

far-ï¬eld gradients can be combined for the pur-

to the diffusion equation (Andrews and Hanks,

poses of scarp analysis into a single variable: the

1985). In the second approach, the midpoint gra-

reduced initial gradient Î± âˆ’ b. The analytic solu-

dient is plotted versus scarp offset for a collection

tion to the diffusion equation for a general scarp

of many scarps with different heights (Bucknam

2.2 ANALYTIC METHODS AND APPLICATIONS 43

House Range

Dome Canyon

Bonneville scarp

lacustrine

deposits

8 m vert.

is mathematically simpler because the derivative

Fig 2.9 Virtual oblique aerial photograph (i.e. aerial photo

draped over digital topography) of pluvial features on the east of Eq. (2.46) is a more compact expression:

side of central Tule Valley. Scarp height is primarily controlled

âˆ‚h(x, t) (Î± âˆ’ b) x + a/(Î± âˆ’ b)

âˆ’b = âˆš

erf

by the slope of the alluvial-fan and debris-cone deposits, i.e.

âˆ‚x 2 4Îºt

wave-cut erosion of steeply dipping deposits tends to

x âˆ’ a/(Î± âˆ’ b)

produce steep, tall scarps. Inset photo shows 8 m tall scarp at

âˆ’ erf âˆš (2.47)

the entrance to Dome Canyon. 4Îºt

The left-hand side of Eq. (2.47) is called the re-

duced gradient proï¬le (i.e. the gradient proï¬le

was given by Hanks and Andrews (1989) as

reduced by b, in order to take the far-ï¬eld gra-

Îºt x+a/(Î±âˆ’b))2 xâˆ’a/(Î±âˆ’b))2 dient into account). Equations (2.46) and (2.47)

eâˆ’ âˆ’ eâˆ’

h(x, t) = (Î± âˆ’ b) 4Îºt 4Îºt

Ï€ are plotted in Figure 2.1 for a range of diffu-

(Î± âˆ’ b) x + a/(Î± âˆ’ b) sion ages Îºt and reduced initial scarp gradients,

a

+ x+ âˆš

erf

Î±âˆ’b Î± âˆ’ b. Figure 2.10b illustrates that the midpoint

2 4Îºt

x âˆ’ a/(Î± âˆ’ b) of the scarp is insensitive to age for young, tall

a

âˆ’ xâˆ’ + bx

âˆš

erf

Î±âˆ’b scarps. Note that in Figure 2.9, the units of the

4Îºt

x-axis are scaled to the scarp width x0 , which, in

(2.46)

turn, is scaled to the scarp height through the re-

lationship x0 = a/(Î± âˆ’ b). Figure 2.10c illustrates

Full-scarp analyses can be performed by compar-

ing the measured elevation proï¬le of a scarp to that, as the reduced initial gradient decreases,

Eq. (2.46), or by taking the derivative of the ele- the age-diagnostic portion of the scarp becomes

vation proï¬le and working with the slope or gra- progressively concentrated away from the mid-

dient proï¬le. Working with the gradient proï¬le point towards the top and bottom of the scarp.

44 THE DIFFUSION EQUATION

These plots suggest that full-scarp methods are

b (a)

a necessary to extract the independent controls

1

of scarp age and initial gradient on scarp mor-

h0

h phology.

h0 Alternatively, a brute-force approach can be

used which computes analytic or numerical so-

x0

0

lutions corresponding to a full range of model

parameters and then plots the error between

the model and measured proï¬les as a function

1 of each parameter value (e.g. diffusion age and

initial gradient). The diffusion age is then con-

1 strained by the model proï¬le that has the lowest

(b)

error relative to the measured data. This brute-

b force approach is most often used, and it has the

advantage that the goodness-of-ï¬t can be evalu-

a-b

ated visually on a proï¬le-by-proï¬le basis. Avouac

0.5 (1993) and Arrowsmith et al. (1998) pioneered this

kt = 1.0 technique, assuming initial scarp angles between

30â—¦ and 35â—¦ . Mattson and Bruhn (2001) improved

= 0.3

kt = 0.01

upon this approach by including a range of initial

= 0.03

= 0.1 gradients in the analysis. Mattson and Bruhnâ€™s er-

0 ror plots are thus two-dimensional, including a

0 1 2 3

range of both diffusion age and reduced initial

angle in the analysis.

0.5 kt = 0.25, 0.5, 1

(c) a âˆ’ b = 0.5 Figures 2.11a and 2.11b illustrate the full-scarp

0.4 method on a short and a tall scarp, respectively,

kt = 0.25, 0.5, 1 each selected from the set of Bonneville scarps.

a - b = 0.35

0.3 The reduced gradient is plotted as a function of

kt = 0.25, 0.5, 1

-b distance along the proï¬le. To the right of each

a âˆ’ b = 0.2

0.2 plot is a grayscale map that illustrates the rel-

ative least-square error between the model and

0.1 observed proï¬les for a range of diffusivity val-

ues and reduced initial gradients. In both maps,

0.0 black colors represent the lowest mismatch and

0 4 6

2

x/x0 X marks the spot of the optimal ï¬t. The loca-

tion of this best-ï¬t solution in the model pa-

Fig 2.10 (a) Analytic solutions of the linear diffusion rameter space simultaneously determines diffu-

equation, starting from an initial gradient Î±, scarp height 2h 0 , sivity and reduced initial gradient on a proï¬le-by-

and far-ï¬eld slope b. The initial scarp width, 2x 0 , is related to

proï¬le basis. In the short-scarp example case, the

height and reduced initial gradient by x 0 = h 0 /(Î±Ã±b). (b)

best-ï¬t solution corresponds to Îº = 0.95 m2 /kyr

Plots of reduced gradient corresponding to the proï¬les in (a).

and Î± âˆ’ b = 0.70. The black band running across

For small Îºt values (< 0.1), no change in the midpoint slope is

the grayscale map illustrates that, although the

observed, illustrating the insensitivity of the midpoint slope to

best-ï¬t solution occurs for Î± âˆ’ b = 0.70, the solu-

diffusion age for young scarps. (c) Plots of reduced gradient

tion is essentially independent of the reduced ini-

versus distance along proï¬le, for families of solutions

corresponding to a range of morphologic ages (Îºt = 0.25, tial gradient. This makes physical sense because,

0.5, 1) and a range of reduced initial gradients (Î± âˆ’ b = 0.2, for a short, narrow scarp, diffusive smoothing

0.35, 0.5). quickly reduces the maximum slope to a value

far below the angle of repose. For a tall, broad

scarp (Figure 2.11b), however, the scarp form is

2.2 ANALYTIC METHODS AND APPLICATIONS 45

of the slope break at the top and bottom of the

0

short scarp (a) scarp. The reduced initial gradient is determined

0.4

by the shape of the ï¬‚at, central portion of the gra-

k

0.3 dient proï¬le shown in Figure 2.11b (also shown

âˆ’b

in the analytic solutions of Figure 2.1f).

0.2 k = 0.95 m2/kyr

a âˆ’ b = 0.70 7

2.2.8 Degradation of archeological ruins

aâˆ’ b 1.4

0.4

0.1

In the summer of 2006 I performed ï¬eld work in

Tibet with Jennifer Boerner. Tibet is full of impor-

0

0

tant archeological ruins, and I wondered whether

tall scarp

0.4

the degradation of stone-wall ruins could be de-

k = 0.78 m2/kyr

a âˆ’ b = 0.48

k

scribed by a diffusion model. Stone walls begin

0.3

as consolidated, vertical structures. Over time

âˆ’b

scales of centuries to millennia following aban-

0.2

donment, however, mortar will degrade, stones

7

aâˆ’ b 1.4

0.4

will topple, and surface processes (creep, biotur-

0.1

bation, eolian deposition, etc.) will act to smooth

(b)

0

the surface.

âˆ’30 âˆ’20 âˆ’10 0 10 20 30

x (m) Figure 2.12 illustrates the application of the

diffusion equation to stone ruins in western Ti-

Fig 2.11 Full-scarp curve-ï¬tting method illustrated for (a)

bet. Figure 2.12a illustrates the solution to the

short and (b) tall example Bonneville scarps. The reduced

diffusion equation with an initially vertical struc-

gradient is plotted as a function of distance from the scarp

midpoint. Grayscale maps to the right of each plot illustrate ture. This equation is identical to the solution for

the relative error between the observed and modeled scarp the gradient proï¬le of a scarp, Eq. (2.47):

as a function of diffusivity Îº and reduced initial gradient Î± âˆ’ b.

x + x0 x âˆ’ x0

h0

In both maps, black colors represent the lowest mismatch h(x, t) = âˆ’ erf

âˆš âˆš

erf (2.48)

2x0 4Îºt 4Îºt

and X marks the spot of the optimal ï¬t. The location of this

best-ï¬t solution in the model parameter space simultaneously

where h 0 is the initial height of the wall and

determines diffusivity and reduced initial gradient on a

x0 is the half-width. Figure 2.12b presents an

proï¬le-by-proï¬le basis. In the short-scarp example case in (a),

oblique view of one example of stone ruins in

the best-ï¬t solution corresponds to Îº = 0.95 m2 /kyr and

Tibet. These degraded walls are up to 0.5 m tall

Î± âˆ’ b = 0.70. The black band running across the grayscale

and have spread out to cover a width of several

map illustrates that, although the best-ï¬t solution occurs for

Î± âˆ’ b = 0.70, the solution is essentially independent of the meters. Figure 2.12d presents a contour plot of

reduced initial gradient. This makes physical sense because, a detailed topographic survey of a small section

for a short, narrow scarp, diffusive smoothing quickly reduces of one of these ruins. Figure 2.12c plots two to-

the maximum slope to a value far below the angle of repose.

pographic proï¬les extracted from Figure 2.12d.

For a tall, broad scarp (b), however, the solution is sensitive

These plots are compared with Eq. (2.48) for Îºt

to the initial angle. The best-ï¬t model parameters in this case

values of 2.0 m2 and 1.3 m2 . Radiocarbon ages in-

have well-deï¬ned values for both diffusion age and reduced

dicate a site abandonment age of 1.7 ka at this

initial gradient (Îº = 0.78 m2 /kyr and Î± âˆ’ b = 0.48).

location. Diffusivity values for these walls, there-

fore, are 1.2 and 0.8 m2 /kyr, respectively. These

values are broadly comparable to the represen-

sensitive to the initial angle. The best-ï¬t model

tative value for the Basin and Range.

parameters in this case have well-deï¬ned values

for both diffusion age and reduced initial gradi-

ent (Îº = 0.78 m2 /kyr and Î± âˆ’ b = 0.48). Physically, 2.2.9 Radionuclide dispersion in soils

ï¬tting the full scarp works best because the scarp Thus far we have only considered applications

proï¬le contains independent information about of the diffusion equation to landform evolution.

the diffusion age and the initial gradient. The dif- Many applications of the diffusion equation are

fusion age is determined primarily by the width applicable in other aspects of surface processes,

46 THE DIFFUSION EQUATION

(a) 1 (b)

h(x)

0.5

= 1.0

= 0.3

kt/x02 = 0.01

= 0.03

degraded walls

= 0.1

0

âˆ’3 âˆ’2 âˆ’1 0 1 2 3

x/x0

(d) 5cm contours

(c)

0.6

5m

h(x)

kt = 1.3 m2

(m)

kt = 2.0 m2

0.4

0.2

profiles

0

2 4

0 8

6

x (m)

and 95%, suggesting that very little 137 Cs has dif-

Fig 2.12 Diffusion modeling of the stone-wall degradation

fused below 6 cm in these proï¬les. The radionu-

in western Tibet. (a) Solutions to the diffusion equation with

an initially vertical structure. (b) Oblique view of degraded clide proï¬le shape may reï¬‚ect the inï¬‚uence of

walls in Tibet. (c) Plots of topographic proï¬les perpendicular both surface erosion/deposition and redistribu-

to the wall, with best-ï¬t solutions corresponding to tion within the soil column. However, in the case

Îºt = 2.0 m2 and 1.3 m2 . (d) Contour map of a

of Fortymile Wash it is reasonable to assume that

microtopographic survey of degraded walls in western Tibet

surface erosion and deposition were negligible

(near Kyunglung), with proï¬le locations noted.

during the past 50 yr. Except for the active chan-

nel, all of the older surfaces have not experienced

however. Here we consider the migration of signiï¬cant ï¬‚ooding for several thousand years.

radionuclides into the subsurface within a soil Eolian deposition rates inferred for the Fortymile

proï¬le. Wash fan by Reheis et al. (1995) also suggest that

Beginning in the mid-1950s and continuing eolian erosion/deposition is less than or equal to

through the mid-1960s, above-ground nuclear several millimeters over the 50-yr time scale since

tests introduced radioactive fallout 137 Cs into the the introduction of fallout 137 Cs. Due to the rela-

soil (He and Walling, 1997). Fourteen 137 Cs pro- tive stability of the Fortymile Wash fan surfaces

ï¬les from the mapped area of Fortymile Wash fan to both ï¬‚uvial and eolian erosion/deposition, we

expect the shape of the 137 Cs proï¬le to predom-

were collected and analyzed (Table 2.1) in order to

determine the rate of Cs migration into the soil inantly reï¬‚ect redistribution processes. The ex-

over the past 50 years. The Fortymile Wash fan is ception is the active channel. In this case the Cs

of particular interest because it is the depozone proï¬le may be strongly inï¬‚uenced by ï¬‚uvial mix-

for sediment eroded from the Nevada Test Site. ing in addition to inï¬ltration and other mixing

Bulk samples were collected at 0--3 cm, 3--6 cm, processes.

and 6--9 cm. The fraction of total 137 Cs in the up- The equation describing the evolution of ra-

per half of the proï¬le typically varied between 80 dionuclide concentration C by diffusive processes

2.2 ANALYTIC METHODS AND APPLICATIONS 47

âˆ«C(z,t)/(C0 dw )

C(z,t)/(C0 dw )

0.2 0.4 0.6 0.8 1

0

0.4

0.2

0 0.6

0

0

2

t = 30 5

t = 20

4

z (cm) sandyD = 6.1D = 5.1

peat,

t = 10 yr

z (cm) loam,

calc. sandy loam, D = 3.5

10 clayey loam, D = 2.0

6

8

cumulative

Harwell Airfield

t = 4.8 yr concentration vs. depth

concentration vs. depth 15

10

(a) (b)

diffusivity vs. surface age

Qa7

0.15

D

(cm2/yr)

0.10

Qa6

0.05 Qa3

Qa5 Qa4

0.00

1 10 100

age (kyr)

(c)

Fig 2.13 (a) Diffusion-model results for normalized where C 0 is the initial concentration within a

concentration vs. depth from t = 10 to 100 yr with thin surface layer of thickness dw , z is depth

D = 0.1 cm2 /yr. (b) Cumulative 137 Cs concentration as a

in the soil proï¬le, D is the diffusivity (called

function of depth in four different soil types at Harwell,

D in this section to emphasize that it is not a

England (data from Gale (1964)), along with best-ï¬t results

hillslope diffusivity), and t is time following de-

for the diffusion model. (c) Plot of diffusivity values vs. surface

position. Equation (2.49) also approximates the

age for the 137 Cs proï¬les on Fortymile Wash alluvial fan.

permeable soil layer as semi-inï¬nite. This is an

Horizontal error bars represent the range of age estimates

accurate approximation for man-made fallout

for each surface based on Whitney et al. (2004). Vertical bars

represent the standard deviation of diffusivity values for all of proï¬les, even in a soil with a calcic horizon

the samples collected on each surface. at depth, because radionuclides do not pene-

trate more than several cm into desert soils over

decadal time scales. It should also be noted that

is given by the solution to the diffusion equation

radioactive decay need not be considered explic-

in a semi-inï¬nite medium with no-ï¬‚ux bound-

itly in this analysis because decay does not af-

ary condition at the surface and a fallout mass

fect the relative radionuclide concentration at

C 0 dw input at z = 0 at t = 0 (Carslaw and Jaeger,

different depth intervals, only their absolute val-

1959):

ues. Figure 2.13a gives an example of the diffu-

1 sion model (Eq. (2.49) using D = 0.1 cm2 /yr and

2 /4D t

eâˆ’z

C (z, t) = C 0 dw âˆš (2.49)

Ï€Dt t = 10--100 yr).

48 THE DIFFUSION EQUATION

137 Cs

Table 2.1 data and inferred diffusivity values for geomorphic surfaces on the Fortymile Wash alluvial

fan

total 137 Cs (pCi/g) D (cm2 /yr)

Sample ID fraction at 3 cm Geomorphic Unit

Cs-071802-A 0.313 0.8274 0.047 Qa6

Cs-071802-B 0.271 0.5387 0.165 Qa7

Cs-071802-C 0.258 0.8100 0.052 Qa3

Cs-071802-E 0.208 0.7644 0.063 Qa3

Cs-071802-G 0.118 0.9593 0.021 Qa3

Cs-071802-I 0.388 0.9614 0.021 Qa4

Cs-071802-J 0.155 0.6387 0.108 Qa6

Cs-071802-K 0.340 0.9558 0.022 Qa5

Cs-071802-N 0.218 0.9082 0.031 Qa5

Cs-071802-P 0.245 0.9428 0.024 Qa4

Cs-071802-Q 0.205 0.9951 0.011 Qa5

Cs-071802-R 0.237 0.9578 0.021 Qa3

Cs-071802-S 0.285 0.8807 0.037 Qa3

Cs-071802-BB 0.321 0.7943 0.056 Qa5

Radionuclide concentrations measured in the Table 2.1 summarizes the data and results

ï¬eld are bulk measurements rather than point from Fortymile Wash. To calculate D , the fraction

measurements because the concentration within of total activity at 3 cm was ï¬rst computed by di-

different depth intervals is measured. For the viding the activity from 0--3 cm by the total activ-

purposes of extracting model parameters from ity from 0--6 cm. Equation (2.50) was then used to

bulk measurements, it is most accurate to rep- infer the value of the error function argument,

âˆš

resent measured data cumulatively as the frac- equal to z/ 4D t, corresponding to the fraction of

tion of total concentration to a given depth. To activity at 3 cm after 50 yr of diffusion following

compare the diffusion model predictions to this nuclear testing. A table of calculated error func-

normalized cumulative curve, it is necessary to tion values was used for this purpose. The value

integrate Eq. (2.49) to give of the error function argument was then used to

solve for D (column 3 in Table 2.1). Figure 2.13c il-

z

C (Î·, t) z

dÎ· = erf âˆš lustrates the relationship between diffusivity val-

(2.50)

C 0 dw 2 Dt

0

ues and geomorphic-surface age on the Fortymile

Wash fan. The vertical bars are the standard de-

Figure 2.13b compares the diffusion model pre-

viation of D values obtained on different proï¬les

diction to 137 Cs proï¬les measured during a 4.7-yr

of the same soil-geomorphic unit. The horizontal

ï¬eld experiment at Harwell Airï¬eld in Oxford-

bars correspond to the age range for each unit

shire, England (Gale et al., 1964). Although these

based on the available age control (Whitney et al.,

data are from a very different climatic regime

2004). The plot exhibits an inverse relationship

than Fortymile Wash, the results illustrate the

between diffusivity and age, with values increas-

precision of the diffusion model in a way that

ing slightly on Qa3 (the oldest surface) relative to

proï¬les in desert soils cannot show because shal-

those on Qa4.

low penetration limits the spatial resolution of

Radionuclides may be redistributed within

desert-soil proï¬les. To obtain the results in Fig-

the soil proï¬le through entrainment and re-

ure 2.13b, the measured 137 Cs proï¬les were ï¬t to

deposition by inï¬ltrating runoff and/or by

Eq. (2.50). Diffusivity values are lowest for clayey

bulk-mixing processes including freeze/thaw,

and calcareous soils at Harwell, increasing in

wetting/drying, and bioturbation. The rate of

value in sandy soils and peat.

2.2 ANALYTIC METHODS AND APPLICATIONS 49

transport by inï¬ltration can be expected to corre- The diffusion equation with radial symmetry

late positively with hydrologic conductivity, and is given by

hence negatively with surface age based on the

âˆ‚h 1âˆ‚ âˆ‚h

= Îºr

results of Young et al. (2004). As such, we can con- (2.52)

âˆ‚t r âˆ‚r âˆ‚r

clude that transport of radionuclides during in-

If Îº is a constant, Eq. (2.52) becomes

ï¬ltration is a signiï¬cant component of transport.

For long time scales, the ï¬nite depth of soil

âˆ‚h âˆ‚ 2h 1 âˆ‚h

=Îº + (2.53)

proï¬les must be considered. In such cases, the

âˆ‚t âˆ‚r 2 r âˆ‚r

concentration of radionuclides in the soil proï¬le

The solution to Eq. (2.53) with an initial radi-

will be given by the series solution

ally symmetric topography f (r ) and a constant-

âˆž

elevation boundary condition h(a, t) = 0 at r = a

dw 2 1 nÏ€dw

C (z, t) = C 0 + sin

Ï€

L n L is given by

n=1

âˆž

2 J 0 (Î±n r )

nÏ€ z âˆ’n2 Ï€ 2 D t/L 2 2

eâˆ’ÎºÎ±n t

h(r, t) = 2

Ã— cos e (2.51) 2

a J 1 (Î±n a)

L n=1

a

Ã— r f (r )J 0 (Î±n r )dr (2.54)

where L is the depth of the soil proï¬le. In some 0

cases, L may be the depth to an impermeable soil

where Î±n , n = 1, 2, . . . are the positive roots of

layer in the subsurface. In other cases, it may be

J 0 (Î±a) = 0 and J 0 (r ) and J 1 (r ) are Bessel functions

more appropriate to take L to be the depth of

of the ï¬rst and second kind (Culling, 1963). This

the wetting front (i.e. the maximum distance to

solution utilizes two boundary conditions. First,

which radionuclides are transported by inï¬ltra-

at r = 0 Eq. (2.54) implicitly assumes

tive transport).

âˆ‚h

=0 (2.55)

âˆ‚r r =0

2.2.10 Evolution of cinder cones

Volcanic cinder cones often have a high degree because Eq. (2.54) is a series comprised of even

of radial symmetry and are composed largely of functions of radius (i.e. J 0 (r )). The Bessel func-

tion J 0 (r ) has a derivative of zero at r = 0, so

unconsolidated fallout tephra. As such, solutions

to the diffusion equation in polar coordinates a series comprised of a sum of these functions

can be used to model their evolution. Cinder must also obey that condition. The second bound-

cones are particularly important landforms for ary condition is a constant elevation of zero at

r = a. For simplicity we assumed the elevation

this type of analysis because their initial con-

ditions can often be precisely inferred by mea- of the surrounding alluvial ï¬‚at equal to zero,

suring the dip of subsurface bedding planes on but any constant value could be used. In hill-

opposite sides of the crater rim. The initial mor- slope geomorphology, a ï¬xed-elevation boundary

phology of pluvial shoreline and fault scarps, in condition applies to a channel that is capable of

contrast, cannot be inferred with comparable pre- transporting all of the sediment delivered by the

cision because the initial mass-wasting phase of hillslope, resulting in neither deposition nor ero-

the scarp typically erases any memory of the ini- sion at the channel location. In the cinder cone

case, a ï¬xed elevation at r = a could correspond

tial shape. In addition, 40 Ar/39 Ar dating enables

the eruption age of cinder cones to be precisely to a channel that wraps around the base of the

measured, and hundreds of cones in the west- cone, but this is not a common occurrence. In

ern United States have been dated in this way. most cases, volcanic cones are surrounded by al-

The framework for analyzing volcanic cones with luvial ï¬‚ats that do not readily transport mate-

analytic solutions to the diffusion equation was rial from the base of the cone. Equation (2.54)

described by Culling (1963) and Hanks (2000). To can still be used for these cases, but a must be

date, however, no analytic solution has been com- chosen to be much larger than the radius of the

pared to speciï¬c natural cones. cone. This way, debris will be removed from the

50 THE DIFFUSION EQUATION

The three integrals that appear in Eq. (2.57)

1

kt/a2 = 0

kt/a2 = 0.0001 can also be written as inï¬nite series. In practice,

kt/a2 = 0.0002 however, it is more accurate to evaluate those in-

0.8 kt/a2 = 0.0004

tegrals numerically because the series converge

kt/a2 = 0.0128 very slowly. Therefore, although closed-form ana-

0.6

h lytic solutions for volcanic cones can be written

h0 down, calculating and plotting the solutions re-

0.4

quires some numerical work.

First, we must solve for the roots of J 0 (Î±a) = 0.

0.2 This is done numerically using a root-ï¬nding

technique (Chapter 9 of Press et al. (1992) pro-

0 vides many techniques for doing this). The Bessel

0 0.33 0.66

function has an inï¬nite number of roots, so how

r/a

many do we need to compute? For our purposes,

Fig 2.14 Plots of analytic solutions to the 3D diffusion the ï¬rst one hundred roots are adequate, but

equation for a volcanic cinder cone for a range of times more roots may be required for high-resolution

following eruption.

proï¬les or very young cones. Second, we must

evaluate the integrals in Eq. (2.57) numerically.

To do this, the algorithms in Chapter 4 of Press et

cone and deposited on the surrounding ï¬‚at. In

al. (1992) can be used. Appendix 1 provides a pro-

these cases, the boundary condition at r = a will

gram for computing the integrals and perform-

simply serve to maintain the alluvial ï¬‚at or pied-

ing the sums in Eq. (2.57).

mont at a constant elevation very far from the

Figure 2.14 illustrates radial proï¬les of Eq.

base of the cone.

(2.57) for the initial cone and at eight sub-

For a volcanic cone of radius rc , crater-rim ra-

sequent times from Îºt/a 2 = 0.0001 to Îºt/a 2 =

dius rr , colluvial ï¬ll radius of r f , and maximum

0.0128. Both axes are normalized. The y-axis is

height of h 0 (Figure 2.14), f (r ) is given by

normalized to the initial cone height h 0 . As in

âŽ§ âŽ« all diffusion problems, the solution can be scaled

r f âˆ’rr

âŽª h0 1 + âŽª

if r < r f

âŽª âŽª up or down in height with no change in the rela-

âŽª âŽª

rc âˆ’rr

âŽª âŽª

âŽª âŽª

âŽ¨h 1+ if r f â‰¥ r < rr âŽ¬ tive cone shape. The r -axis is scaled to the model

r âˆ’rr

0

f (r ) = rc âˆ’rr (2.56) domain length a.

âŽªh 1+ if rr â‰¥ r < rc âŽª

âŽª0 âŽª

rr âˆ’r

âŽª âŽª

âŽª âŽª

rc âˆ’rr

âŽª âŽª In the coneâ€™s early-stage evolution, the great-

âŽ© âŽ

if r â‰¥ rc

0 est change occurs at the crater rim. This is not

surprising since this is where the proï¬le curva-

ture is greatest. The crater rim is locally trian-

As a technical point, it should be noted that

gular and for early times the crater evolution

the colluvial ï¬ll radius r f must be ï¬nite in or-

is broadly similar to moraine evolution consid-

der for Eq. (2.54) to apply because otherwise the

ered in Section 2.2.6. For intermediate times,

boundary condition given by Eq. (2.55) is violated.

the position of the crater rim migrates inward.

Physically, r f corresponds to the width of the

This migration is associated with the additional

colluvium that ï¬lls the crater shortly following

advective term in Eq. (2.52) that is not present

eruption. Substituting Eq. (2.56) into Eq. (2.54)

gives

âˆž

r2

2h 0 J 0 (Î±n r ) 2rr rr

2

eâˆ’ÎºÎ±n t

f

h(r, t) = J 1 (Î±n r f ) âˆ’ J 1 (Î±n rr ) + 1 + rc J 1 (Î±n rc )

Î±n J 1 (Î±n a) r c âˆ’ rr r c âˆ’ rr r c âˆ’ rr

2

a2 n=1

Î±n rr rf rc

+ r 2 J 0 (Î±n rr )dr + r 2 J 0 (Î±n r f )dr + r 2 J 0 (Î±n rc )dr

2

r c âˆ’ rr 0 0 0

(2.57)

2.2 ANALYTIC METHODS AND APPLICATIONS 51

tine et al., 2005). The model solutions compare

(a) 117Â° W 116Â° W

well to the observed proï¬les, and best-ï¬t morpho-

Yucca

Mtn.

logic ages are approximately 400 and 4000 m2 ,

Nevada

Black Cone

respectively. The ratios of the morphologic age

N Ca

ev if Amargosa

ad or

to the radiometric age provide estimates for the

Valley

Funeral a nia

36.5Â° N Mtns.

time-averaged Îº values for these cones: 5.2 and

3.6 m2 /kyr. These values are at the upper range of

Death

Valley

Îº values inferred from pluvial shoreline and fault

Black

Mtns.

Lathrop Wells

N 10 km

36.0Â° N

scarps in the southwestern US (Hanks, 2000).

This approach enables us to estimate how

much height the cone has lost since its erup-

tion. Figure 2.15 suggests that the crater rim of

(b) 1 Lathrop Wells has been eroded by approximately

Lathrop Wells 8% of its original height, or 9--10 m. Black Cone,

0.8 in contrast, has been degraded by approximately

50% according to this approach.

A ï¬nal word about Lathrop Wells should be

0.6 kt = 400 m2

h noted. The dating of this cone has been the

h0 Black Cone

subject of considerable debate because it is the

kt = 400 m2

0.4

youngest volcano in the Yucca Mountain region,

the site of the proposed high-level nuclear waste

0.2

repository. Initial attempts to date the cone ge-

omorphically concluded that the cone was no

0

older than approximately 10 ka based on its com-

0.3 0.4

0.2

0 0.1

r (km) plete lack of drainage channels. However, cinder

cones are comprised of coarse (pebble-to-gravel-

Fig 2.15 (b) Observed (thick line) and best-ï¬t model

sized) tephra. As such, it is likely that little or no

proï¬les (thin line) for a relatively young (Lathrop Wells Cone)

surface runoff can occur until sufï¬cient time has

and a relatively old (Black Cone) cone in the Crater Flat

passed for eolian dust to clog surface pore spaces.

volcanic ï¬eld, Nevada (location map in (a)).

Therefore, runoff generation on these cones may

be impossible for tens of thousands of years fol-

in 2D problems. Eventually, the crater is ï¬lled lowing eruption. In the absence of any radio-

by diffusion with debris and the late-stage cone metric dates for this cone, diffusion-morphologic

morphology is described by J 0 (ar ) with decreas- dating would suggest that the cone is 20 ka to

80 ka assuming a conservative range of Îº values

ing amplitude over time. At this point, the expo-

from 0.5 to 2.0 m2 /kyr. This range includes the

nential time-dependent term in Eq. (2.54) has ï¬l-

tered all of the high-frequency components in the now-accepted radiometric age but does not in-

topography and only the lowest-frequency term is clude the original estimate of 10 ka. Of course,

signiï¬cant in the series. The inset photo in Fig- all morphologic-age estimates suffer from poor

constraints on Îº values. However, since Îº val-

ure 2.14 shows the proï¬le of an early Quaternary

volcanic cone in the Cima volcanic ï¬eld, Califor- ues can be calibrated on a cone-by-cone basis as

nia. The proï¬le shape is similar to the late-stage shown here, they offer great potential for bet-

ter regional calibration of Îº values and improved

proï¬les plotted in Figure 2.14.

Figure 2.15 compares the observed topogra- understanding of their underlying climatic and

phy of two cones in the Crater Flat volcanic ï¬eld, textural controls.

Amargosa Valley, Nevada, to best-ï¬t model solu-

2.2.11 Delta progradation

tions. Lathrop Wells Cone has been radiometri-

cally dated to be 77 ka (Heitzler et al., 1999), The application of the diffusion equation to

while Black Cone has an age of 1.1 Ma (Valen- large-scale depositional landforms (i.e. deltas and

52 THE DIFFUSION EQUATION

alluvial fans) is problematic for two reasons. First, U0

sediment

input

depositional landforms have complex internal

(a)

processes that are not captured by the diffusion

model. For example, deposition in deltas and al-

luvial fans takes place by a series of avulsion

h

events. At short time scales, deposition is local-

ized along one or more lobes that are active un- t2

t1

til the lobe becomes superelevated with respect

to the surrounding terrain. A large ï¬‚ood can

cause the lobe to shift rapidly to another portion

of the fan or delta. If the diffusion equation ap-

plies at all to the large-scale evolution of these

x

systems, it applies only at a time scale that av-

1.0

erages the effects of many avulsion events. Sec- (b)

ond, variations in sediment texture can produce

0.8

spatial variations in sediment mobility. The sedi-

ment texture in the proximal portion of an allu-

vial fan, for example, can be one to two orders of 0.6

h

magnitude larger than sediment in the distal re- h0

gion (i.e. boulders compared to sand). Variations 0.4

in channel morphology (and hence ï¬‚ow depth)

and sediment texture both combine to create spa-

0.2

tial variations in the effective diffusivity that are

not well constrained. 0.0

2.0

1.0 3.0

0.0

Kenyon and Turcotte (1985) applied the 1D 4.0

U0

(x âˆ’ U0 t)

diffusion equation to model the progradation of

k

a single delta lobe. Their model does not ap-

ply to the long-term 3D morphology of deltas. Fig 2.16 (a) Schematic diagram of the 2D model for a

Here we review their approach and propose a prograding delta lobe in the Kenyon and Turcotte (1985)

model for the 3D evolution of deltas over geologic model. (b) Solution for the proï¬le shape of a prograding delta

time scales. Both models are limited in that they (Eq. (2.61)).

neglect the ï¬‚exural--isostatic subsidence of the

basin in response to sediment loading. This lim-

diffusion equation gives

itation will be corrected in Chapter 5 where we

consider the geometry of foreland basins within d2 h

dh

âˆ’U 0 =Îº 2 (2.59)

the context of a model that couples the diffusion dÎ¾ dÎ¾

equation to ï¬‚exural subsidence.

Equation (2.59) must be solved with boundary

The Kenyon and Turcotte model assumes that

conditions h(0) = h 0 and h(âˆž) = 0. The solution is

sediment is introduced to the basin at an eleva-

tion h 0 above the basin ï¬‚oor. As the basin ï¬lls, the U0Î¾

h(Î¾ ) = h 0 eâˆ’ (2.60)

Îº

position of the river mouth migrates towards the

basin at a rate U 0 . Generally speaking, solving dif- Substituting Eq. (2.58) into Eq. (2.60) gives

fusion problems with moving boundaries poses a

U0

considerable challenge. In this case, however, the h(x, t) = h 0 eâˆ’ (xâˆ’U 0 t)

(2.61)

Îº

moving boundary can be handled by introducing

a new variable Î¾ given by The solution to this model is given in Figure 2.16.

The solution is an exponential function with a

Î¾ = x âˆ’ U 0t (2.58)

width controlled by a balance between diffusion

The sediment source is ï¬xed in this new coor- (which transport sediment away from the source)

dinate system. Substitution of Eq. (2.58) into the and the movement of the source. This solution

2.2 ANALYTIC METHODS AND APPLICATIONS 53

does a good job of matching the observed stratig- 1.0 kt =

raphy of the southwest pass segment of the Mis-

104 km2

sissippi River delta (Turcotte and Schubert, 2002). 2 Ã— 104 km2

0.8

More generally, however, the exponential form of 4 Ã— 104 km2

Eq. (2.61) is inconsistent with the generally sig- 8 Ã— 104 km2

0.6

h

moidal nature of deltaic deposits.

h0

Over geologic time scales delta lobes will

0.4

avulse and a radial, 3D landform will be cre-

ated. In order to model this 3D development, we

consider the diffusion equation in polar coordi- 0.2

nates (Eq. (2.52)). We also assume a ï¬xed sedi-

0.0

ment source (in contrast to the moving sedi-

0 100 200

r (km)

ment source of Kenyon and Turcotte) and a spa-

tially variable diffusivity that is greatest near the

Fig 2.17 Proï¬les of the 3D radially symmetric model of

river mouth and decreases with distance from

delta formation (Equation (2.67)).

the mouth as 1/r :

Îº

Îº(r ) = (2.62)

r deposits, which is more consistent than the

This spatially variable diffusivity represents the exponential form of the Kenyon and Turcotte

loss of transport capacity as rivers reach the model.

ocean. The speciï¬c 1/r dependence we have

assumed is not well constrained, but it pro-

2.2.12 Dust deposition downwind

vides a useful starting point for representing the

of playas

decrease in sediment mobility as distance from

The dust cycle in arid environments is character-

the shoreline increases.

ized by a net eolian transfer of dust from playas

The 3D diffusion equation with radial symme-

to piedmonts (Pye, 1987). Understanding the pro-

try is given by Eq. (2.52). Substituting Eq. (2.62)

cesses and rates of this transfer is important for

into Eq. (2.52) gives

many basic and applied geologic problems. Dust

âˆ‚h Îº âˆ‚ 2h transport in the atmosphere also provides a nice

= (2.63)

âˆ‚t r âˆ‚r 2 illustration of diffusion acting in concert with

advection. In this section we describe the phe-

Equation (2.63) can be reduced to an ordinary dif-

nomenology of dust deposition in a well-studied

ferential equation by introducing the similarity

ï¬eld site in California, and then apply an analytic

variable

solution for the advection, diffusion, and gravi-

r3

Î·= (2.64) tational settling of dust from complex, spatially-

9Îºt

distributed sources, following the work of

In terms of Î·, Eq. (2.63) becomes Pelletier and Cook (2005).

d2 Î¸ The study area is located in southern Amar-

dÎ¸

âˆ’ = (2.65)

gosa Valley, California, where Franklin Lake Playa

dÎ·2

dÎ·

abuts the Eagle Mountain piedmont (Figure

where Î¸ = h/ h 0 . The solution to this equation is

2.18a). The water table in Franklin Lake Playa is

Î¸ = eâˆ’Î· less than 3 m below the surface (Czarnecki, 1997).

(2.66)

Czarnecki identiï¬ed several distinct geomorphic

Transforming back to the original coordinates: surfaces on Franklin Lake Playa that can be read-

ily identiï¬ed in LANDSAT imagery (Figure 2.18a).

r3

h(r, t) = h 0 eâˆ’ 9Îºt (2.67)

For the purposes of numerical modeling, the ac-

Figure 2.17 presents plots of Eq. (2.67) for tive portion of the playa was mapped based on

different values of Îºt. The model predicts the Czarneckiâ€™s map of playa surfaces with signiï¬-

classic sigmoidal (S-shaped) geometry of deltaic cant dust-emitting potential (Figure 2.18a).

54 THE DIFFUSION EQUATION

(a) (b)

116Â°W

117Â°W

Yucca

Mtn.

Eagle

Mtn.

Nevada

Ne Amargosa

va Valley

da

Ca

Funeral

36.5Â°N

lif

Mtns.

or

ni Eagle

a

Mtn.

Death

Qa2

Valley

Qa3

Qa4

Black

Qa5â€“Qa7

Mtns.

N 10 km

36.0Â°N predominant

Soil-geomorphic map and playa

N

sample-pit locations wind direction pits

(c)

outline of active

wind

modern playa

direction

Eagle

3â€“6 m/s

calm

Mtn.

6â€“9 m/s

90%

Franklin Lake

Playa

1%

2%

3%

Amargosa Eagle

River Mtn.

secondary hotspot

1 km

10 20 40

0 80 cm

Silt thickness on Qa3

Fig 2.18 (a) Location map and LANDSAT image of Eagle conditions (see wind-rose diagram in Figure

Mountain piedmont and adjacent Franklin Lake Playa,

2.18a). Silt-rich eolian deposits occur directly un-

southern Amargosa Valley, California. Predominant wind

derneath the desert pavements of Eagle Moun-

direction is SSE, as shown by the wind-rose diagram (adapted

tain piedmont, varying in thickness from 0 to

from January 2003â€“January 2005 data from Western Regional

80 cm based on soil-pit measurements. The tech-

Climate Center, 2005). Calm winds are deï¬ned to be those

nical term for these deposits is cumulic eolian

less than 3 m/s. (b) Soil-geomorphic map and oblique aerial

epipedons (McFadden et al., 2005), but here we

perspective of Eagle Mountain piedmont, looking southeast.

Terrace map units are based on the regional classiï¬cation by refer to them as eolian or silt layers for sim-

Whitney et al. (2004). Approximate ages: Qa2 â€“ middle plicity. These layers are predominantly composed

Pleistocene, Qa3 â€“ middle to late Pleistocene, Qa4 â€“ late of silt but also include some ï¬ne sand and sol-

Pleistocene, Qa5â€“Qa7 â€“ latest Pleistocene to active. (c) Map

uble salts. The homogeneity of these deposits,

of eolian silt thickness on Qa3 (middle to late Pleistocene)

combined with their rapid transition to gravelly

surface, showing maximum thicknesses of 80 cm close to the

alluvial-fan deposits below, makes the layer thick-

playa source, decreasing by approximately a factor of 2 for

ness a reasonable proxy for the total dust content

each 1 km downwind. Far from the playa, background values

of the soil.

of approximately 20 cm were observed. For color version,

Alluvial-fan terraces on Eagle Mountain pied-

see plate section. Modiï¬ed from Pelletier and Cook (2005).

mont have a range of ages. Eagle Mountain pied-

mont exhibits the classic sequence of Quaternary

alluvial-fan terraces widely recognized in the

The Eagle Mountain piedmont acts as the

southwestern United States (Bull, 1991). Mapping

depositional substrate for dust emitted from

and correlation of the units was made based on

Franklin Lake Playa under northerly wind

2.2 ANALYTIC METHODS AND APPLICATIONS 55

the regional chronology of Whitney et al. (2004). material deposited over time) and that Qa4 sur-

In this example, we focus on silt thicknesses of faces have received a higher average dust ï¬‚ux

the Qa3 terrace unit (middle to late Pleistocene) because they have existed primarily during the

because of its extensive preservation and lim- warm, dry Holocene. The Qa5 unit (latest Pleis-

ited evidence of hillslope erosion. Silt thicknesses tocene) exhibited uniformly thin deposits under-

were measured at locations with undisturbed, lying a weak pavement, independent of distance

planar terrace remnants to the greatest extent from the playa, suggesting that the trapping abil-

possible. A color map of silt layer thickness for ity of young surfaces is limited by weak pavement

this surface is shown in Plate 2.18b, draped over development.

the US Geological Survey (USGS) 30 m digital el- As discussed in Chapter 1, atmospheric trans-

evation model and an orthophotograph of the port of particulate matter can be modeled as a

area. Silt thickness is observed to decrease by combination of turbulent diffusion, downwind

roughly a factor of 2 for every 1 km of distance advection, and gravitational settling. The 3D con-

from the playa. Several kilometers downwind, a centration ï¬eld for a point source located at

(x , y , 0), obtained by solving Eqs. (1.17) and (1.18),

background thickness of 15--20 cm was observed

on Qa3 surfaces. The strongly localized nature of is given by

âŽ¡

exp âˆ’ 4K (xâˆ’x )

uz

u(y âˆ’ y ) p2 (x âˆ’ x )

2

Q p pz

âŽ£

c p (x, y, z, x , y ) = exp âˆ’ âˆ’ +

exp

4K (x âˆ’ x ) Ï€uK (x âˆ’ x ) uK K Ku

4Ï€ K (xâˆ’x )

u

âŽ¤

xâˆ’x

u âŽ¦

Ã—erfc z+ p (2.68)

4K (x âˆ’ x ) Ku

where Q is the source emission rate and erfc

downwind deposition in this area suggests that

is the complementary error function. Equation

Franklin Lake Playa is the source for nearly all

(2.68) combines the two-dimensional solution of

of the eolian deposition on Eagle Mountain pied-

Smith (1962) with an additional term required to

mont. Localized deposition also implies that dust

describe crosswind transport (Huang, 1999). Equa-

deposition rates may vary regionally by an order

tion (2.68) assumes that the settling velocity q

of magnitude or more, down to spatial scales of

is small compared with the deposition velocity

1 km or less.

p. Stokesâ€™ Law (Allen, 1997) implies a settling ve-

Figure 2.19a illustrates plots of silt thickness

locity of less than 1 cm/s for silt particles (i.e.

as a function of downwind distance, for compar-

particles less than 0.05 mm) in air. Deposition

ison with two-dimensional model results. This

velocities consistent with the spatial distribu-

plot includes all measurements collected from

tion of deposition on Eagle Mountain piedmont,

terrace units late Pleistocene in age or older

however, are approximately 5 cm/s, or at least

(Qa4--Qa2) within a 1 km wide swath of western

ï¬ve times greater than q. For nonpoint sources,

Eagle Mountain piedmont. Data from the Qa4

Eq. (2.68) can be integrated to give

and Qa2 units show a similarly rapid downwind

decrease in silt thickness, illustrating that the âˆž

c(x, y, z) = c p (x , y )dx dy

Qa3 pattern is robust. Silt thicknesses were rela- (2.69)

âˆ’âˆž

tively similar on the three terrace units at com-

parable distances from the playa. This similar- Using Eqs. (2.68) and (2.69), the deposition rate

ity was not expected, given the great differences on a ï¬‚at surface is given by pc(x, y, 0). Deposition

in age between the three surfaces. This may be on a complex downwind surface can be estimated

partly explained by the fact that Qa2 surfaces as pc[x, y, h(x, y)], where h(x, y) is the elevation

have undergone extensive hillslope erosion (and of the downwind topography. This approach is

hence preserve only a portion of the total eolian only an approximation of the effects of complex

56 THE DIFFUSION EQUATION

80

(a) (b)

Qa2

u = 5 m/s

Qa3

K = 5 m2/s

Qa4

p = 0.05 m/s

60

z

analytic

silt z depositional

thickness

topography

(cm)

40

(x', y')

cross-wind

direction

y

20

a

u = 5 m/s

K = 10 m2/s source

p = 0.075 m/s

wind x

0

2 3 4

1

0 direction

distance downwind from playa (km)

(c) 2 km

0.1 0.2 0.4 1.0

0

wind direction

ñòð. 3 |