<<

. 3
( 14)



>>

Fig 2.3 (a) Approach to steady state in a hillslope subject to
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-

’ = (2.65)
gosa Valley, California, where Franklin Lake Playa
d·2

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
( 14)



>>