. 2
( 14)


The distal portion of the Whiskey Pete™s allu-
1.2.2 Sand-dominated eolian systems vial fan appears lighter in color than the prox-
Many playas are dominated by sand rather than imal portion of the fan. In the ¬eld, this dif-
silt or clay. Roach Playa, located along the south- ference in surface re¬‚ectivity corresponds to an
ern Nevada--California border, is one example of abrupt increase in the sand content of soils be-
a sand-dominated playa. Sand from Roach Playa low a threshold elevation on the fan. Above this
is transported off the playa and onto the nearby elevation, sand sourced from the playa and from
Whiskey Pete™s alluvial fan (Figure 1.21). Because channels draining the mountain accumulates lo-
transport takes place by saltation rather than cally and is not transported across the fan. Be-
large-scale atmospheric mixing, the spatial pat- low this elevation, both coarse and ¬ne eolian
tern of sand deposition downwind from Roach sand from nearby playa and channel sources is

readily transported across the fan. These distal-
fan eolian corridors are controlled by a thresh-
old fan-surface relief (Cook and Pelletier, 2007).
(mm) ripples dunes
When along-strike relief falls below a threshold mega
value, an eolian transportational surface devel- dunes
ops. When the along-strike relief rises above the
threshold value, sand is trapped locally in low 102
10’1 101 103
10’2 104
l (m)
spots and a continuous surface of transporta-
tion is prevented from developing. Along-strike
Fig 1.22 Relationship between coarsest fraction (20th
relief promotes the storage of windblown sand in
percentile) of sand grains and bedform wavelength, », based
two ways. First, topographic obstacles create aero-
on measurements of ripples, dunes, and megadunes in the
dynamic recirculation zones on their lee sides.
Sahara Desert (after Wilson, 1972).
These recirculation zones are characterized by
very low bed shear stresses. Sand that is deposited
behind those obstacles may be stored inde¬nitely.
Second, the presence of topographic obstacles in-
creases the force necessary to move the sand. It
is not only suf¬cient for sand particles to be en-
trained from the surface, they must be picked up
with suf¬cient force to be transported over large
steps downwind. The transport of sand over a
complex surface, therefore, involves topographic
controls on both the shear stress exerted on the
particle and on the value of threshold shear
(b) top view side view
stress necessary to move the particle past down- VP
wind topographic obstacles.
Eolian bedforms can be classi¬ed into three
basic types depending on their spatial scale and wind
position within the bedform hierarchy: ripples,
dunes, and megadunes (Wilson, 1972). Ripples
are the smallest of the three bedforms, and are
typically spaced by 0.1--1 m and have heights of
1--10 cm with higher, more widely spaced rip-
ples forming in areas with stronger winds and/or
coarser sand. Dunes are the next level in the Fig 1.23 (a) Barchan dunes of the Salton Sea area,
California. (b) Schematic diagram illustrating how a barchan
bedform hierarchy and form only when ripples
gets its shape. The migration velocity of a dune is proportional
are present. Dunes commonly have spacings of
to the ratio of the perimeter to the area of the cross section
10--100 m and heights of 1--10 m. Megadunes form
perimeter. Therefore, the central core of the dune will move
when dunes are present, and may attain spac-
slower than the sides, causing arms to form over time (c).
ings of several kilometers and heights of sev-
eral hundred meters. All three bedform types ex-
hibit a correlation between spacing and grain trols the air ¬‚ow above it, and the pattern of ero-
size in the Sahara Desert (Wilson, 1972) (Figure sion and deposition and subsequent modi¬cation
1.22). This correlation is not observed everywhere, of the bedform size and shape through time.
however, partly because geographic variations in Eolian dunes come in a variety of types de-
wind speed also control bedform spacing. The pending primarily on sand supply and wind-
precise mechanisms responsible for bedform gen- direction variability. Barchan dunes (Figure 1.23)
esis are still a subject of active research. All three provide a nice example of both of these controls.
bedform types, however, likely involve a positive Barchan dunes have a central core ¬‚anked by two
feedback between the bedform shape, which con- arms oriented downwind. Barchan dunes form

in a steady wind direction when the sand sup- of deposition also depends on whether the sand
ply is relatively low (leaving some areas of bare had been transported into the lee side of an incip-
ground exposed). Dune formation under such ient dune. In the model, lee sides are mapped by
conditions requires that sand coalesce and ac- de¬ning a shadow zone (i.e. all areas that were lo-
cated in shadow given a sun angle of 15—¦ oriented
cumulate in patches separated by bare ground.
How does this ˜˜coalescence™™ happen? The reason parallel to the wind direction) with a high proba-
why sand tends to collect in patches rather than bility of deposition. Finally, sand deposited on the
spread out uniformly on the valley ¬‚oor has to bed rolls down the direction of steepest descent if
do with the higher coef¬cient of restitution of deposition causes the surface to be steeper than
soft sandy surfaces relative to hard bare ground. the angle of repose.
Sand in saltation will tend to bounce off of bare Werner™s model is capable of reproducing the
ground and land more softly on patches of sand. four principal dune types (transverse, barchan,
As a result, sand saltating across a valley ¬‚oor has star, and longitudinal) by varying the sand supply
a higher probability of deposition in sand-rich ar- and wind direction variability (e.g. Figure 1.24b).
eas than sand-free areas. This success is remarkable considering that his
Next consider how different portions of an in- model does not include any details on the micro-
cipient sand dune will migrate downwind. Fig- scopic physics of grain-to-grain interactions that
ure 1.23b illustrates an idealized incipient dune had long been assumed to be essential for under-
with a semi-circular cross-sectional shape. Dunes standing the formation of eolian bedforms. Al-
migrate by moving sand along their surfaces, and though Werner™s model does not successfully re-
the rate of sand movement is proportional to the produce all aspects of dunes (e.g. topographic pro-
perimeter, P , of the dune cross section. In order ¬les), his model includes the essential feedbacks
to move the dune a given distance downwind, of dune formation necessary for understanding
however, all of the sand in the cross sectional how different types of dunes are formed. First,
area, A, has to be moved. The rate of dune mi- Werner™s model illustrates how sand blown from
gration, therefore, is proportional to the ratio of a river bed tends to coalesce into discrete mounds
the perimeter to area. In the center of an in- rather than spreading out. Because the probabil-
cipient dune, the ratio of perimeter to area is ity of sand deposition is higher on an area al-
lower than it is on the sides of the dune. As ready covered by sand, a steady supply of sand
such, the sides will migrate faster, causing the particles saltating across an area will be prefer-
development of ˜˜arms™™ oriented downwind (Fig- entially deposited on areas already covered by
ure 1.23c). The wind direction must be steady for sand. In this way, discrete patches of sand will
barchan dunes to form because a slight change form and increase their heights relative to the
in wind direction is suf¬cient for dunes of this surrounding topography. Second, as the heights
type to become unstable. A shift in wind direc- of incipient bedforms increase, the length of
tion can cause sand from one arm of the barchan their shadow zones (which mimics the recircula-
to be transferred to the other arm, lengthening tion zone downwind of the dune crest) will also
one arm at the expense of the other over time. lengthen. This lengthening allows dunes to grow
Eventually, a longitudinal dune is formed with a higher since their heights are limited only by the
crestline oriented parallel to the wind. angle of repose and the lengths of their shadow
Werner (1995) developed a numerical model zones. Werner™s model is further discussed in
for the formation of eolian dunes using sim- Chapter 7.
ple geometric rules. Discrete units of sand in
Werner™s model are picked up from a sand bed
at random and transported a characteristic dis-
1.3 A tour of the glacial system
tance downwind. Sand units are deposited back
down on the bed with a probability that is rela-
tively low on bedrock surfaces compared to sandy Today, ice sheets and glaciers are restricted to
surfaces (re¬‚ecting the higher coef¬cient of resti- a relatively small portion of the Earth™s surface.
tution of a hard versus a soft bed). The probability Large ice sheets are found only on Greenland and

5 km

(a) (b)
elevations. This positive feedback towards thicker
Fig 1.24 (a) Shaded-relief image of the Algodones Dunes
and thicker ice is limited by the shear strength
located along the California“Arizona border near Yuma,
Arizona. (b) Example output of Werner™s model for eolian of the ice and friction at the base of the ice. Even-
dunes for the transverse-dune case with large sand supply and tually, the thickness and slope of the ice become
constant wind direction.
great enough to overcome the shear strength of
ice or the basal friction (whichever is lower), caus-
ing the ice to surge outward by internal defor-
Antarctica, and alpine glaciers are restricted to el-
mation or basal sliding. Zones of higher eleva-
evations above approximately 5 km a.s.l. in tem-
tion have positive mass balance (i.e. more snow
perate latitudes. In contrast, ice covered nearly
is precipitated on the landscape than melts away)
all of Canada, the northernmost United States,
and areas of lower elevation have negative mass
and most of northern Europe during the LGM.
balance. Ice ¬‚ows from areas of positive mass bal-
Alpine glaciers expanded to cover areas with el-
ance (i.e. accumulation) to areas of negative mass
evations approximately 1 km lower than today
balance (i.e. ablation) when suf¬cient thickness
(Porter et al., 1983). While geomorphology in the
and slope have been built up. Although the me-
southernmost United States means wind and wa-
chanics of ice sheets and glaciers is very complex
ter, geomorphology in Canada and the northern
in detail, at large scales their shapes are often
United States is primarily a story of ice (Flint,
well adjusted to the threshold condition for shear
strength or basal friction. Increases in snow ac-
Ice sheets and glaciers can be thought of as
cumulation cause shear stresses to increase above
conveyor belts of ice. Some of the snow that
the threshold condition, causing rapid surging of
falls at high elevations sticks around through
the ice sheet or glacier in response. Therefore, the
the following summer. Successive years produce
shapes of ice sheets and glaciers on Earth are pri-
more snow that acts as a weight on underlying
marily a function of their shear strength or basal
snow, turning the snow into ¬rn (a state inter-
friction. The speed of these ice bodies is primar-
mediate between snow and ice) and eventually
ily a function of their accumulation and ablation
ice as more snow is deposited above it. As snow
and ice accumulate and the ground surface in-
Ice sheets and glaciers move by a complex
creases in elevation, more snow and ice accumu-
combination of internal ice ¬‚ow, sliding of ice
late because of the cooler temperatures at higher

ice margin

basal velocity
till plane

predominant erosion little or no erosion or deposition
(frozen bed or low velocity)
The erosional potential of ice sheets and
Fig 1.25 The subglacial environment can be broadly divided
into three zones: (1) a zone of little or no erosion or glaciers is not fully understood, but most re-
deposition, located near the divide, where the ice sheet is searchers agree with Hallet (1979; 1996) that the
frozen to the bed or where velocities are too small to erosion rate is proportional to a power-law func-
signi¬cantly modify the surface, (2) a zone of predominant
tion of the basal sliding velocity V :
erosion, where basal velocities and meltwater concentrations
= ’aV b
increase toward the margin, transporting all of the bedrock (1.19)
debris removed from the bed, and (3) a zone of predominant
where a is a coef¬cient that depends primarily
deposition near the ice margin, where basal velocities wane
on the bedrock erodibility, and b is an empirical
and moraines and till plains are commonly formed. Also
exponent equal to ≈ 2 (Figure 1.25).
shown are the zones of accumulation, ablation, internal ¬‚ow
lines, and their relationship to the major subglacial zones. The ability of a warm-based ice sheet or
glacier to erode its bed varies as a function of
on bedrock, and sliding of ice on glacial till. distance along the longitudinal pro¬le. Ice ¬‚ow
Whether an ice sheet or glacier deforms inter- velocities increase with distance downslope in a
nally or slides along its base depends primarily manner very similar to that of ¬‚uvial drainage
on the thermal conditions at the base. In ˜˜cold- networks. In a ¬‚uvial drainage network, water
based™™ ice sheets and glaciers the ice is frozen to ¬‚ows in the direction of the hillslope or chan-
its bed, which precludes basal sliding. In ˜˜warm- nel orientation or aspect. Fluvial channels acco-
based™™ ice sheets and glaciers, motion occurs by modate this continually increasing discharge pri-
a combination of basal sliding and internal ice marily by increasing in width and depth down-
¬‚ow. The frictional force between an ice body and slope. Channel ¬‚ow velocities also increase down-
its bed depends on many factors, but sliding of slope, but at rates typically much lower than
ice over glacial till is generally thought to require changes in width and depth. Flow in an ice sheet
lower basal shear stresses than sliding of ice over or glacier also occurs in the direction of the ice-
bedrock. The distinction between cold and warm- surface gradient. As discharge accumulates with
based ice sheets is very signi¬cant for geomor- distance downslope, the ice velocity is the vari-
phology because only warm-based ice sheets per- able that changes most greatly to accommodate
form signi¬cant glacial erosion. Millions of years the increasing ¬‚ow. The shape of the ice sheet or
of glacial activity have produced little measur- glacier does not vary as greatly as the velocity be-
able glacial erosion in areas dominated by cold- cause a small change in thickness or slope causes
based glaciers (e.g. Dry Valleys of Antarctica). a large change in discharge due to the threshold

a Newtonian ¬‚uid, strain rate is proportional to
n = 10
n= shear stress. The coef¬cient of proportionality is
the viscosity, a measure of stiffness of the ¬‚uid.
3 In a non-Newtonian ¬‚uid, strain rate and shear
du a
n= 3 stress are nonlinearly related:
2 n

=a (1.20)
n= 1 dz
1 where u is the ¬‚uid velocity, z is the distance
along the pro¬le, „ is the shear stress, „0 is a ref-
erence or yield stress, and a and n are rheological
parameters. The values of a and n can be consid-
1.0 1.5 2.0
t/t0 ered to be constant in some cases, but more of-
ten they depend on the temperature or density of
Fig 1.26 Plot of normalized velocity gradient versus the ¬‚ow. Experiments indicate that ice deforms
normalized shear stress for a Newtonian ¬‚uid (n = 1), a according to Eq. (1.20) with n ≈ 3 (i.e. Glen™s Flow
non-Newtonian ¬‚uid (e.g. n = 3 and n = 10), and a perfectly
Law (Glen, 1955)). Some materials, including rock
plastic material (n = ∞).
in Earth™s upper crust, have highly nonlinear rhe-
ologies characterized by Eq. (1.20) with n = 10 or
nature of ice ¬‚ow. Close to an ice divide, there- greater. Such materials deform very slowly un-
fore, ice-¬‚ow velocities are quite low because the til brittle fracture occurs, after which deforma-
contributing area is small. Ice ¬‚ow velocities are tion can be accommodated relatively easily by
a combination of internal ¬‚ow and basal slid- slippage along fractures. Equation (1.20) does not
ing, so not all of the increased ¬‚ow downstream resolve the detailed structure of fracture zones
occurs as basal sliding. Nevertheless, increasing in such materials, but it can be used to study
ice discharge is usually associated with increased the large-scale behavior of deformation in such
basal sliding velocities. Erosion rates close to an materials.
In the limit n = ∞, the material is called per-
ice divide are, therefore, relatively low according
to Eq. (1.19). Further downslope, ¬‚ow velocities fectly plastic. Perfectly plastic materials are com-
pletely rigid below their yield stress, „0 , but ¬‚ow
and erosion rates increase to accommodate ice
delivered from upslope. Basal sliding velocities entirely without resistance above it. The perfectly
increase with distance downice, usually reach- plastic model is an idealization, but in some
ing a maximum beneath the equilibrium line. cases it is very useful within certain limitations.
Past the equilibrium line, ice-¬‚ow velocities de- Glaciers, ice sheets, lava ¬‚ows, and debris ¬‚ows,
crease due to ablation. This decrease, combined are all examples of ¬‚ows driven by gravity. Flow
with the accumulation of erosional debris from in these systems occurs when the material thick-
further upice, usually causes a zone of sediment ens locally, increasing the shear stress at the base,
deposition beneath the ablation zone of an ice causing the ¬‚ow to surge forward, thereby low-
sheet or glacier. ering the shear stress back below a threshold
value. The dynamics of this type of surge be-
1.3.1 How ice deforms havior depends sensitively on the value of n in
The ¬‚ow of liquid water and air over the Earth™s Eq. (1.20). If one is not interested in resolving
surface are both examples of Newtonian ¬‚ows. the detailed dynamics of individual surge events,
Understanding glacial erosion, however, requires however, then the perfectly plastic model can be
quantifying the behavior of a non-Newtonian very useful. In nature, there is always some ¬-
¬‚uid. The difference between Newtonian and nite time lag between changes in accumulation
non-Newtonian ¬‚uids is illustrated in Figure 1.26, and the associated glacier response. In the per-
where graphs of velocity gradient versus shear fectly plastic model, however, glacier surges start
stress are illustrated for three types of ¬‚uids. In and stop instantly, because any increase in ¬‚ow

mer ice sheets. These mantle-¬‚ow patterns are
dif¬cult to interpret uniquely in terms of ice-
sheet size and shape, however, because the re-
sponse is a function of lithospheric rigidity and
2 km
the mantle-viscosity pro¬le in addition to the size
and shape of the ice sheet. Therefore, many dif-
3 km
ferent ice-sheet sizes and shapes may be equally
consistent with the observed data pattern of sea-
level change. Landforms such as drumlins, ¬‚utes,
and mega-scale lineations can also be used to re-
construct ice-¬‚ow directions and provide infor-
mation on the locations of ice domes (Clarke
et al., 2000). While these geomorphic markers pro-
vide solid evidence that ice once ¬‚owed in a cer-
tain direction, it is often dif¬cult to reconstruct
Fig 1.27 Model reconstruction of the Laurentide Ice Sheet
a spatially complete ¬‚ow ¬eld or to determine
at 18 ka constructed using a perfectly plastic ice sheet model
the time that the observed ¬‚ow directions were
and assuming isostatic balance.
imprinted on the landscape.

1.3.2 Glacial landforms
thickness that causes the shear stress to increase
above a threshold value instantaneously triggers Upstate New York provides some of the most
¬‚ow to correct the imbalance. Perfectly plastic striking examples of erosional and depositional
glaciers and ice sheets, therefore, are constantly glacial landforms on Earth. Figure 1.28a illus-
in balance with the pattern of accumulation and trates the topography of a portion of upstate New
ablation. York immediately south of Lake Ontario. This On-
Figure 1.27 illustrates a model reconstruction tario Lowlands region has the largest drumlin
of the Laurentide Ice Sheet over North America ¬eld in North America, with nearly 10 000 drum-
at 18 ka. This reconstruction was developed us- lins. Drumlins are subglacial bedforms elongated
ing the perfectly plastic model with basal shear parallel to the ice-¬‚ow direction and composed
stresses similar to modern ice sheets (e.g. Green- primarily of subglacial till or sediment. Drum-
land, East Antarctica) and assuming isostatic bal- lins in the New York State drumlin ¬eld come in
ance. Isostatic balance in this context means that a wide variety of shapes and sizes. The controls
the weight of the ice sheet causes the crust to on drumlin size and shape are not well under-
subside into the mantle until the point where stood, but we will explore a possible model for
the additional load is balanced by the buoyancy drumlin formation in Chapter 7. Subsurface bed-
associated with the displaced mantle material. ding in drumlins often parallels the topographic
Despite decades of research, the shape and thick- surface, suggesting that they form by localized,
ness of the Laurentide Ice Sheet (LIS) is still a upward ¬‚ow of sediment into low pressure zones
subject of active research and debate. Sophis- beneath the ice sheet.
ticated mass-balance thermomechanical models Along the southern portion of Figure 1.28a,
are used to reconstruct the ice sheet, but un- the elevation rises up to the Allegheny Plateau
certainty in paleoclimatic, rheological, and basal- and the Finger Lakes Region. The Finger Lakes
¬‚ow parameters results in uncertainty in the re- Region is characterized by ¬ve major elongated
sulting reconstructions (Marshall et al., 2000). In- lakes that fan out like the ¬ngers of a hand.
terpretation of past sea-level changes is another More generally, ¬nger lakes are any kind of elon-
method by which the shape of the former ice gated glacial lake formed near the margin of a
sheet can be inferred. In this method, the past former ice sheet. Figure 1.28b illustrates the to-
≈100 yr of sea-level records are used to determine pography in the vicinity of the town of Ithaca,
the present pattern of mantle ¬‚ow beneath for- located near the southern tip of Cayuga Lake (the


(b) resolution limit

N (>A) N (>A) ∝ A

102 uncompensated,
"great" lakes
compensated lakes

104 105
100 101
A (km 2)
Fig 1.28 Shaded relief images of two portions of upstate
New York. (a) Drumlin ¬eld of the Ontario Lowlands. (b) Fig 1.29 Analysis of lakes in central Canada. (a) Radarsat
Southern Cayuga Lake (one of the ¬ve Finger Lakes) in the mosaic image of central Canada (resolution 250 m), including
vicinity of the town of Ithaca. a closeup of Dubawnt Lake region to illustrate the full
resolution of the image. (b) Cumulative frequency-size
distribution of lakes extracted from (a). Also included in this
largest of the Finger Lakes). The topography of the plot are the Great Lakes, even though they are not included
in (a). A power-law relationship N(>A ) ∝ A ’1 is observed
Ithaca area is characterized by wide (>1 km), deep
from A = 100 to A ≈ 104 km2 . Above A ≈ 104 km2 , the
(>100 m), glacially scoured valleys cut into the
largest lakes do not ¬t the power-law relationship. The larger,
broad, low-relief Allegheny Plateau (von Engeln,
more uniform size of the great lakes is consistent with the
1961). Fluvial channels in the region are said to
in¬‚uence of lithospheric ¬‚exure on their formation.
be under¬t because channels occupy only a small
portion of the broad valley ¬‚oors that were once
occupied by ice up to 1 km in thickness. Flu-
vial drainage density is low, re¬‚ecting the short cumulative frequency-size distribution of lakes is
period of geologic time since small-scale land- the number of lakes larger than a given area. Con-
forms were last smoothed by glacial processes. sider Figure 1.29b, which illustrates the cumula-
The deeply scoured valleys and smooth uplands tive frequency-size distribution of all the lakes in
illustrate that glacial processes smooth the land- Canada resolvable in a 250 m-resolution Radarsat
scape at small (<5 km) scales, while increasing mosaic of the country. Certainly not all Cana-
relief through trough formation at larger scales. dian lakes are erosional in origin; some of the
In addition to elongated ¬nger lakes, the smaller lakes occur in kettles and hummocky
northern US and Canada contain innumerable moraine, for example. Nonetheless, the majority
glacially carved lakes that come in a wide range of Canadian lakes, particularly the largest lakes,
of sizes (Figure 1.29a). The cumulative frequency- are erosional in origin. Figure 1.29a shows a sub-
size distribution is one tool for quantifying set of the Radarsat mosaic, illustrating the mo-
the population of lakes in a given region. The saic image data at both small and large scales.




largest lakes in North America may be controlled
Fig 1.30 Sur¬cial geologic map of eastern Canada,
by some additional process not involved in the
illustrating moraines (thick curves oriented perpendicular to
erosion of smaller lakes. In Chapter 6, we will
ice ¬‚ow), eskers (long, thin curves oriented parallel to ice
propose that the bending of the lithosphere be-
¬‚ow) and drumlins (short, thin curves parallel to ice ¬‚ow).
neath the ice is that additional process that con-
trols the formation of the largest North American
glacial lakes.
Lakes were extracted from the image by scanning
Depositional landforms beneath ice sheets in-
the grid for pixels with water coverage. Once a
clude moraines, eskers, drumlins, kettles, and
water-covered pixel was identi¬ed, all neighbor-
kames. Figure 1.30 includes a map of deposi-
ing pixels were searched recursively, calculating
tional landforms in eastern Canada. As ice sheets
the number of pixels in the connected lake ˜˜do-
have waxed and waned, nearly all glaciated ar-
main.™™ In this way, the size of each lake in the im-
eas have experienced a combination of erosion
age was determined. We have also included the
and deposition. As such, glacial landscapes are
¬ve Great Lakes in this population, even though
classic palimpsests. This poses a problem for the
they are not included in the mosaic image. The
interpretation of formerly glaciated topography.
distribution of Figure 1.29b includes two dis-
tinct lake types. Lakes smaller than ≈ 104 km2 in For example, predominantly erosional topogra-
phy may be covered with till of varying thickness,
area follow a simple power-law distribution char-
acterized by N (> A) ∝ A ’1 . Above 104 km2 , the obscuring the bedrock erosional surface.
Alpine glaciers differ from ice sheets primar-
largest lakes do not follow the trend of smaller
ily in their size and the extent to which ice ¬‚ow
lakes. Instead, these lakes are larger than the size
is controlled by subglacial topography. In an ice
predicted by an extrapolation of the trend for
sheet, the thickness of the ice is often several
smaller lakes. This distribution suggests that the

high peaks

plateau surface
U-shaped valley

hanging valley glacial

Fig 1.32 Oblique virtual aerial photograph of the Wind
River Range, illustrating the three distinct topographic levels
Fig 1.31 U-shaped and hanging valley of the Beartooth
characteristic of many glaciated mountain ranges of the
Mountains, Montana.
western US.

times greater than the relief of the bed topog-
raphy (except near the margins). As a result, the a main valley takes place much faster than ero-
pattern of ice ¬‚ow is controlled principally by sion of a side tributary, leaving the side trib-
the planform shape of the ice sheet with valleys utary ˜˜hung™™ above the ¬‚oor of the main val-
and ridges beneath the ice sheet playing a rel- ley. Hanging valleys are, therefore, a consequence
atively minor role. In contrast, ice thickness in of the more spatially discontinuous nature of
a glacier is typically much less than the relief of glacial erosion compared with ¬‚uvial erosion. It
the underlying bed topography (e.g. ≈ 100 m for a should be noted, however, that U-shaped valleys
glacier spanning ≈ 1 km in bed elevation). Alpine and hanging valleys are not unique to glacial sys-
glaciers are strongly controlled by bed topogra- tems. Fluvial channels that incise very rapidly
phy, with ice thickness greatest in areas of low- can oversteepen adjacent hillslopes to the point
est bed slope (i.e. valley bottoms). The tendency where water and sediment draining through side
of ice thickness and ¬‚ow velocity to be greatest tributaries will separate from the channel bed
in valley bottoms suggests that glaciers may in- (Wobus et al., 2006). As such, ¬‚uvial systems in tec-
crease subglacial landscape relief more readily tonically very active areas or in areas with strong
than ice sheets. Alpine glaciers are also strongly structural control can also form hanging valleys.
controlled by elevation. The relationship between Figure 1.32 illustrates the large-scale geomor-
the ELA (Equilibrium Line Altitude) and glacial phology of alpine glacial terrain using the Wind
erosion is much tighter in alpine glaciers than River Range as a type example. The Wind River
in ice sheets. Range has three distinct topographic levels. The
Figure 1.31 illustrates examples of two clas- high peaks of the Wind River Range (≈ 4 km a.s.l.)
stand ≈ 1 km above the surrounding low-relief
sic alpine glacial landforms: U-shaped valleys
and hanging valleys. Fluvial valleys are classi- plateau surface. This large-scale pattern of nar-
cally V-shaped, re¬‚ecting the contrasting ero- row high peaks surrounded by a broad plateau
sional power of ¬‚uvial hillslopes process regimes. is common in the glaciated mountain ranges of
In a V-shaped valley, hillslope erosion requires the western US. The formation of these plateau
steep slopes in order to ˜˜keep up™™ with the con- surfaces is not well understood. Did glaciers em-
centrated power of rivers to incise narrowly into anating from the high peaks ˜˜plane™™ the topog-
the valley bottom. In alpine glacial systems, in raphy to a relatively smooth plateau, or was the
contrast, ice ¬lls the valley ¬‚oor, resulting in a plateau exhumed by isostatic rebound respond-
˜˜plug™™ of ice that erodes the valley ¬‚oor and side ing to erosion concentrated in the high peaks?
slopes more uniformly than in the ¬‚uvial case. The most intense glacial erosion occurs near
A hanging valley forms when glacial erosion of the equilibrium line in most alpine glacial

systems. Figure 1.33a illustrates the topography of
(a) a portion of the Uinta Mountains using a shaded-
relief image of a US Geological Survey Digital Ele-
vation Model (DEM). The topography of the Uinta
Mountains is characterized by relatively ¬‚at di-
vide surfaces, steep cirque walls, and overdeep-
ened cirque ¬‚oors. The divides are relatively ¬‚at
because of the limited erosion that takes place
at elevations far above the equilibrium line. Near
the equilibrium line, glacial erosion is concen-
trated spatially, decreasing both up and down-
slope. This local maximum in erosion causes lo-
calized glacial scour. In river systems, bed scour
decreases the local slope, which acts as a negative
2 km
feedback to limit scour. In glacial systems, how-
ever, scour continues as long as the ice surface
(b) cirques
points downslope. Alpine glaciers can maintain
sloping ice surfaces even when the bed beneath
the ice has been scoured to form a closed de-
pression. When the ice retreats, a cirque lake is
formed. Figure 1.33b illustrates the cirque form
at the margin of the Beartooth Mountains in
Fig 1.33 Alpine glacial landforms near the equilibrium line:
(a) Shaded relief image of glacial topography in the Uinta Montana.
Mountains, Utah, illustrating steep cirque walls and The formation of overdeepenings beneath
overdeepened cirque lakes. (b) Cirques of the Beartooth alpine glaciers is illustrated schematically in Fig-
Mountains, Montana.
ure 1.34. Figure 1.34 shows a longitudinal pro¬le


evolves to



of an alpine valley glacier. The basal sliding ve-
Fig 1.34 Conceptual model of ¬‚ow in an alpine valley
locity beneath the glacier is a maximum beneath
glacier. Basal sliding velocity and glacial erosion are
the ELA. Over time, glacial erosion concentrated
concentrated beneath the equilibrium line. Over time,
beneath the ELA will cause the formation of an
concentrated erosion leads to the formation of an
overdeepened bed that forms a cirque lake when the ice overdeepened bed. Numerical models exist that
retreats. form realistic cirques (MacGregor et al., 2000), but

P-T Granite

30 km

Huayna Potosi

P-T Granite


(a) (b)
modern ELA

exhumation begins at 10 Ma,
accelerates towards present
(Benjamin, 1986)

Fig 1.35 Topography, ice cover, and structural geology of
may actually increase the elevations of the high
the Cordillera Real, central Andes. (a) View of the glaciated
peaks of the range. In this ˜˜relief-production™™ hy-
peak of Huayna Potosi, one of the high peaks of the
pothesis, glacial erosion causes isostatic imbal-
Cordillera Real, with the Altiplano in the foreground. (b)
ance. Because the erosion is concentrated at val-
Structural geology of a portion of the Cordillera Real,
ley bottoms with minimal erosion of the high
showing the intrusion of Permo-Triassic granite of the high
peaks, however, the presence of isostatic rebound
peaks into Paleozoic metasediments of the Eastern
together with the absence of signi¬cant erosion
Cordillera. Apatite ¬ssion track (AFT) data indicate that the
high peaks of the Cordillera Real were exhumed rapidly causes the peaks to grow even higher in eleva-
beginning at 10 Ma, accelerating towards the present day. tion. The difference between the glacial-buzzsaw
and relief-production hypotheses is partly an is-
sue of scale. At scales of 100 km or larger, glacial
many questions remain. For example, do most
erosion most likely does limit the relief of moun-
cirques primarily re¬‚ect stable ELAs during LGM
tain belts. At scales of a single valley and peak,
conditions, or are their morphologies controlled
however, the relief-production model appears to
by ¬‚uctuations in ELAs through time?
be most consistent with the results of glacial
The tight coupling between glacial ice cover,
landform evolution models. Feedbacks between
glacial erosion, and elevation suggests that
glaciers and topography may also involve climate
some interesting feedback relationships may ex-
through enhanced orographic precipitation in ar-
ist in glaciated mountain belts. For example, the
eas of higher/steeper topography.
˜˜glacial-buzzsaw™™ hypothesis states that glaciers
The Cordillera Real of the central Andes is one
limit the elevation of mountain belts because
potential example of the role of glacial erosion
any uplift that raises the height of mountain
in enhancing topographic relief. The Cordillera
belts above the ELA will quickly be met with en-
Real comprises the boundary between the Alti-
hanced ice cover and glacial erosion that will
plano to the west and the Eastern Cordillera to
lower elevations back below the ELA (Whipple
the east (Figure 1.35). This region is topograph-
et al., 2004). This model predicts a tight corre-
ically and structurally distinct from both of its
lation between the peak elevation of mountain
surrounding regions. Topographically, the range
ranges and the local ELA. On the other hand, the
is a very narrow (i.e. ≈ 15--20 km) series of peaks
fact that glaciers erode most vigorously in valley
up to ≈ 2 km above the surrounding topography.
bottoms and cirque ¬‚oors suggests that glaciers

Fig 1.36 Schematic model of
ELA isostatic uplift and erosion feedbacks
evolves to
(5 km)
in alpine glacial systems.

isostatic rock uplift

evolves to
(5 km)

broad uplift zone

Structurally, the Cordillera Real is comprised of This rebound would have uplifted the ¬‚anks of
granite exhumed between the Paleozoic sedimen- the Cordillera Real (in addition to its glaciated
tary rocks that make up the Eastern Cordillera. high peaks), increasing the area of the Cordillera
If the Cordillera Real displayed the usual signs Real that stands above the ELA. This increase
of active uplift (modern seismicity, active Qua- in glaciated area would further increase glacial
ternary faults), then its steep, high topography erosion and isostatic uplift in a positive feed-
could simply be attributed to active tectonics. back, causing localized exhumation of igneous
However, active tectonic uplift in the central An- and metamorphic rocks from the subsurface.
des occurs only in the Subandean Zone today,
which is the lowest part of the Andes (i.e. ≈
1.4 Conclusions
1--2 km a.s.l.). Exhumation rates in the Cordillera
Real measured thermochronologically (e.g. Ben-
The remainder of the book is organized accord-
jamin et al., 1987; Gillis et al., 2006), however,
ing to mathematical themes. Starting with Chap-
show high rates from 10 Ma to the present.
ter 2, each chapter focuses on a particular type
Figure 1.36 illustrates a schematic model for
of equation or algorithm and includes example
the exhumation and peak uplift of the Cordillera
applications that span Earth surface processes.
Real. According to this model, the Cordillera Real
This structure serves to emphasize the common
was a broad anticline with peak elevations near
5 km before 10 Ma. Global cooling of ≈ 5—¦ at mathematical language that underlies many of
the disparate process zones of Earth™s surface.
10 Ma lowered the ELA from above 6 km to ap-
Chapter 2, for example, focuses on techniques
proximately 5 km. This ELA lowering would have
used to solve the classic diffusion equation and
initiated glacial erosion in the highest portions
the general class of equations that are similar
of the Cordillera Real. Because glacial erosion in-
to it. Since the diffusion equation can be used
creases as a function of distance from the di-
to describe hillslope evolution, channel-bed evo-
vide, the high peaks of the Cordillera Real would
lution, delta progradation, hydrodynamic disper-
not have been subject to signi¬cant erosion,
sion in groundwater aquifers, turbulent disper-
however, while lower areas within the glaciated
sion in the atmosphere, and heat conduction
zone would have been more intensely eroded.
in soils and the Earth™s crust, Chapter 2 will
As such, an initially broad anticline would have
be relevant to many aspects of Earth surface
been sculpted to a narrow peak. Glacial ero-
sion would also have triggered isostatic rebound.
Chapter 2

The diffusion equation

where ρ is the bulk density of sediment on the
2.1 Introduction hillslope, h is the elevation, κ is the diffusivity
in units of L 2 /T (where L is length and T is
time) and x is the distance from the divide. In
The diffusion equation is perhaps the most
Eq. (2.1), q is a mass ¬‚ux with units of M/L /T
widely used differential equation in science.
(where M is mass, L is length, and T is time).
It is the equation that describes heat trans-
Mass must also be conserved. In this context, con-
port in solids and in many turbulent ¬‚uids.
servation of mass means that any change in sed-
In geomorphology, the diffusion equation quan-
iment ¬‚ux along a hillslope results in either an
ti¬es how landforms, especially hillslopes, are
increase or decrease in the elevation. One way
smoothed over time. Diffusional processes act to
to think about conservation of mass is to con-
smooth elevation in the same way that conduc-
sider a small segment of a hillslope pro¬le (e.g.
tion smooths temperature in solids. The diffusion
the section between x3 and x4 in Figure 2.1).
equation has played a particularly important role
If more sediment enters the segment from up-
in the study of well-dated landforms, such as plu-
slope than leaves the segment downslope, that
vial shoreline scarps and volcanic cones, in which
hillslope segment must store the difference, re-
the initial topography can be relatively well con-
sulting in an increase in the average elevation.
strained at a known time in the past. Certain
Conversely, if more sediment leaves the segment
conditions must be met for the diffusion equa-
downslope than enters the segment upslope (as
tion to apply to hillslope evolution, however. The
in the section between x1 and x2 ), there is a
diffusion equation is applicable to gently slop-
net loss of sediment and the elevation must de-
ing hillslopes comprised of relatively uniform,
crease. The diffusion equation applies this logic
unconsolidated alluvium or regolith subject to
to in¬nitesimally-small hillslope segments in the
creep, rainsplash, and bioturbation. The accuracy
same way that the ¬rst derivative of calculus cal-
of the model breaks down as hillslope gradients
culates the average slope h/ x of a function
increase and mass-movement processes become
h(x) in the limit as x goes to zero. Mathemati-
cally, conservation of mass requires that the in-
The applicability of the diffusion equation has
crease or decrease in the elevation be equal to the
two requirements: slope-proportional transport
change in ¬‚ux per unit length, divided by the bulk
and conservation of mass. First, the ¬‚ux of sedi-
density ρ:
ment per unit length, q, must be proportional to
the hillslope gradient:
‚h 1 ‚q
q = ’ρκ (2.2)
‚t ρ ‚x

where t is time. Substituting Eq. (2.1) into Eq. (2.2) time, however, the rate of erosion and deposition
gives the classic diffusion equation: decreases and the widths of the top and bottom
of the scarps where erosion and deposition oc-
‚h ‚ 2h
=κ 2 curs increases (Figure 2.1f).
‚t ‚x
Many physical processes can be classi¬ed as
If, alternatively, q is expressed as a volumet- either diffusive or advective. Therefore, one rea-
ric ¬‚ux per unit length, the ρ terms are re- son for studying the simple diffusion and advec-
moved from both Eqs. (2.1) and (2.2). We will tion equations in detail is that they provide us
use both the mass and volumetric ¬‚uxes in this with tools for studying more complex processes
chapter. that can be broadly classi¬ed as one of these two
Equation (2.3) is applicable to hillslopes with a process types. Conceptually, diffusive processes
constant cross-sectional topographic pro¬le along result in the smoothing of some quantity (tem-
strike (i.e. a fault scarp). More generally, hillslopes perature, elevation, concentration of some chem-
evolve according to the 2D diffusion equation, ical species, etc.) over time. For instance, if one
given by injects a colored dye into a clear ¬‚uid being
vigorously mixed, the average dye concentration
‚h ‚ 2h ‚ 2h
=κ +2 = κ∇ 2 h (2.4) smooths out over time (rapidly at ¬rst when con-
‚t ‚ x2 ‚y
centration gradients are high, then more slowly)
Figure 2.1a illustrates a hypothetical fault according to the diffusion equation. Similarly, if
scarp 2 m high and 10 m wide after some ero- faulting offsets an alluvial fan, the scarp eleva-
sion has taken place. Figures 2.1b and 2.1c illus- tion smooths out over time as sediment is eroded
trate the gradient and curvature of the scarp, from the top of the scarp and deposited at the
respectively. The diffusion equation states that base.
the sediment ¬‚ux q at any point along a hill- Advection, in contrast, involves the lateral
slope is proportional to the hillslope gradient translation of some quantity without spreading.
(i.e. Figure 2.1b). The magnitude of the ¬‚ux is If a colored dye is injected into a laminar ¬‚uid
illustrated in Figure 2.1a at several points along moving down a pipe, for example, the dye will
the pro¬le using arrows of different lengths. At be transported down the pipe with the same
the top of the scarp, the ¬‚ux increases from left velocity as the ¬‚uid but with little spreading.
to right, indicating that more material is mov- Similarly, if a fault offsets a bedrock pediment
ing out of the section than is being transported instead of an alluvial fan, the fault will evolve
into it from upslope. This results in erosion along primarily by weathering-limited slope retreat in-
the top of the scarp where the change in gradient stead of smoothing. In such cases, slow bedrock
along the pro¬le (i.e. the curvature) is negative. weathering acts to maintain a steep fault scarp
Conversely, ¬‚ux decreases from left to right at over time because transport of the weathered
the bottom of the scarp, indicating that more material is rapid enough that sediment does
material is moving into that segment than out not accumulate at the scarp base and no depo-
of it. The result is an increase in surface eleva- sition occurs. This type of slope evolution can
tion (i.e. deposition) along the base of the scarp be modeled by the advection equation, where
where curvature is positive. The rate of erosion the erosion rate is proportional to the hillslope
or deposition varies with time (Figures 2.1d--2.1f) gradient:
according to the magnitude of curvature. Imme-
‚h ‚h
diately after faulting, for example, the change in (2.5)
‚t ‚x
¬‚ux per unit length (i.e. the curvature in Figure
2.1c) is very large and concentrated right near The variable c is the advection velocity in Eq.
the top and bottom of the scarp. The result is very (2.5), and it equals the rate of horizontal re-
rapid, concentrated erosion and deposition at the treat of the scarp. In general, transport-limited
top and bottom of the scarp, respectively. Over slopes and channels evolve diffusively while

(a) (d)
1 1

h h
(m) erosion
x1 x2


x3 x4
0 0

0.0 0.0
(b) (e)

kt = 25 m2

kt = 25 m2
= 0.8
= 2.5

(c) (f)




5 0 5 10 15
0 5 10 15
x (m)
x (m)
In many applications, diffusion is one com-
Fig 2.1 Evolution of a topographic scarp, illustrating (a)
ponent of a more complex model. The transport
elevation, (b) slope, and (c) curvature. In (a), arrows of
of contaminants in groundwater ¬‚ow, for ex-
varying length represent the sediment ¬‚ux at each point. In
the diffusion model, the ¬‚ux is proportional to the local ample, can be quanti¬ed using a combination
slope, and the resulting raising or lowering rate of the surface of advection (i.e. the translocation of contami-
is proportional to the change in ¬‚ux per unit length, which, in nants as they are carried along with the ¬‚uid)
turn, is proportional to the curvature. (d) and (e) Graphs of
and diffusion (i.e. the spreading of contaminant
elevation, slope, and curvature for ¬ve times following scarp
concentration from hydrodynamic dispersion as
offset (κt = 0.01, 0.03, 0.1, 0.3, 1.0.)
the contaminant travels along variable ¬‚ow paths
through a porous aquifer). Similarly, the concen-
tration of dust downwind from a playa is the re-
weathering-limited slopes and channels evolve
sult of both advection and turbulent diffusion.
advectively. Analytic and numerical techniques
Advection describes the motion of the dust as it
for solving the advection equation will be pre-
is carried downwind, while diffusion describes
sented in Chapter 4.

the spreading of the plume as it is mixed by at-
mospheric turbulence. In many of the examples
in this chapter, we will assume a ¬xed chan-
nel head at the base of a diffusing hillslope. Of
course, channels have their own complex evolu-
tion and are rarely ¬xed at the same elevation transport
over time. As such, a complete understanding of
hillslopes or channels requires a fully coupled
model of the ¬‚uvial system of which hillslope
diffusion is one component. Although the tech- linear (slope-proportional)
transport (i.e. diffusion eqn.)
niques used in this chapter are framed around
solving a single equation that has limited ap-
plication in geomorphology, the techniques we
will discuss form the basis for solving many real-
world problems in which diffusion is one part of Fig 2.2 Schematic graph of ¬‚ux versus slope. The diffusion
a more realistic model. equation assumes a linear relationship between ¬‚ux and
The applicability of the diffusion model to slope. More generally, sediment ¬‚ux increases nonlinearly
hillslope evolution depends on the processes act- with slope, diverging at a critical slope value Sc . However,
even the most general ¬‚ux relationship can be approximated
ing to move sediment on the hillslope. Creep
as linear for small slopes. As such, the predictions of linear
and rain splash are examples of hillslope pro-
and nonlinear transport models are the same for
cesses that are accurately represented by the dif-
gently-sloping topography.
fusion equation for gentle slope angles (Carson
and Kirkby, 1972). In these processes, particle
movement takes place by a series of small jumps
the divide (Carson and Kirkby, 1972):
triggered by freeze/thaw cycles, wetting/drying,
or momentum imparted from rain drops. In these ‚h
q = ’κ x (2.6)
jumps, particles may move both upslope and ‚x
downslope but the distance traveled by a parti-
cle moving downslope is slightly greater than the The evolution of hillslopes governed by this
distance of those traveling upslope. Moreover, the ¬‚ux relationship is considered in Section 2.3.2.
increase in distance traveled by a particle mov- Mass movements are not diffusive because the
ing downslope increases linearly with the hill- sediment ¬‚ux increases rapidly and nonlinearly
slope angle or gradient for relatively small slopes. with small increases in hillslope gradient near
the point of slope failure Sc , as illustrated in
Bioturbation is another important hillslope pro-
Figure 2.2.
cess that can often be approximated as diffusive,
The value of κ controls how fast diffusion
particularly for gravelly hillslopes of arid envi-
takes place. In the Basin and Range Province of
ronments where processes such as slope wash
the western US, an approximate value of κ =
are relatively ineffective due to the coarse tex-
1 m2 /kyr has often been used based on studies of
ture of the slope material. Slope wash, rilling,
pluvial shorelines and other landforms of known
and mass movements are examples of hillslope
age (Hanks, 2000). It is widely recognized, how-
processes that are not diffusive. Slope wash is
ever, that the 1 m2 /kyr value varies spatially ac-
not diffusive because the shear stress exerted by
cording to climate, vegetation, soil texture, and
overland ¬‚ow increases with distance from the di-
other variables, even if that variation is not yet
vide as ¬‚ow accumulates along the slope. While
well characterized. One drawback of the diffusion
hillslopes dominated by slope wash are not gov-
approach in geomorphology is that we often do
erned by the diffusion equation, they can be de-
not know the value of κ for a given location with
scribed by a more general, spatially variable dif-
any precision. This is a valid criticism, but it is
fusion equation in which the ¬‚ux is proportional
still important to study hillslope evolution even if
to the hillslope gradient and the distance from

the absolute rates of change cannot be precisely uplifting mountain range. In such cases, the
determined. diffusion equation simpli¬es to the time-
Three classic publications form the basis of independent equation:
many of the results in this chapter. Carslaw and ‚ 2h
κ +U = 0 (2.7)
Jaeger (1959) is the de¬nitive text on analytic so- ‚ x2
lutions to the diffusion equation. Although the
where U is the uplift rate. Solving for the hill-
book was written primarily for engineering ap-
slope curvature, Eq. 2.7 gives:
plications to heat conduction in solids, the solu-
‚ 2h U
tions presented in that book can be applied to =’ (2.8)
‚ x2 κ
diffusion problems that arise in any ¬eld of sci-
ence. Culling (1963) was among the ¬rst authors Integrating Eq. (2.8) gives an equation for the hill-
to apply the diffusion equation to hillslope evo- slope gradient:
lution. His paper is still the most complete ref- ‚h U
= ’ x + c1 (2.9)
erence on the subject. Culling™s paper includes ‚x κ
solutions appropriate for many different initial
where c 1 is an integration constant. In order to
conditions as well as a stochastic treatment of
constrain the value of c 1 , it is necessary to spec-
the relationship between the random motion of
ify the value of ‚h/‚ x at one end of the hillslope
sediment particles on the hillslope and the diffu-
pro¬le. The upslope end of the hillslope is cho-
sive evolution of the hillslope. Carson and Kirkby
sen to coincide with a divide in this case. Divides
(1972) is the most complete reference on the rela-
are defined as locations across which no water or
tionship between speci¬c hillslope processes and
sediment passes. As such, a no-flux boundary con-
their signature forms.
dition is appropriate for the upslope end of this
In the following sections we present solutions
pro¬le. The no-¬‚ux boundary condition is given
to the diffusion equation primarily within the
context of hillslope evolution. As such, the dif-
fusing variable h will represent elevation along =0 (2.10)
a topographic pro¬le. It should be emphasized, x=0

which implies c 1 = 0. Integrating Eq. (2.9) gives
however, that the solutions we describe are gen-
eral. Therefore, any solution that we obtain for h
h(x) = ’ x + c2 (2.11)
(for a particular initial condition and boundary

conditions) can be applied to any other diffusion
where c 2 is another integration constant. At the
problem, whether the physical process is heat
downslope end of the pro¬le, the elevation is as-
conduction, contaminant transport, etc. That is
sumed to be constant. This boundary condition is
part of the power of applied mathematics: once
appropriate for a hillslope-channel boundary in
a set of solutions or a solution method has been
which all of the sediment delivered to the chan-
learned it can often be rapidly applied to many
nel is transported away from the slope base with
different physical systems.
no erosion or deposition. If the hillslope has a
length L and an elevation of zero at the slope
2.2 Analytic methods and
h(L ) = 0 (2.12)
then c 2 = U L /2κ and the hillslope pro¬le is given
2.2.1 Steady-state hillslopes
Steady-state landscapes have received a great deal h(x) = (L ’ x 2 ) (2.13)

of attention in geomorphology in recent years.
The topography of a steady-state landform is time- Equation (2.13) shows that diffusive hillslope pro-
independent, such as when uplift and erosion ¬les are parabolas in cases of steady-state up-
are in perfect balance at every point in a rapidly lift. The hillslope relief in Eq. (2.13) is given by

U L 2/2κ, so the relief of diffusive hillslopes is pro- Fourier showed that nearly any function could be
represented as a series of sinusoidal functions:
portional to uplift rate and slope length, and
inversely proportional to diffusivity. Equations ∞
nπ x
f (x) =
(2.9)--(2.13) illustrate how boundary conditions an sin (2.16)
are used to solve differential equations. n=1

Equation (2.8) suggests that in regions of sim- where f (x) is any function that goes to zero
ilar climate and hillslope processes, hillslope cur- at x = 0 and x = L . Given a function f (x) that
vature can be a relative measure of tectonic up- obeys these boundary conditions, the values of
lift rates. Equation (2.8) may be applicable to the the Fourier coef¬cients an can be calculated using
early stage of uplift when hillslope gradients are
2 nπ x
still relatively low and weathering rates are suf- an = f (x ) sin dx (2.17)
¬ciently rapid for hillslopes to be regolith cov- 0

ered. In general, however, Eq. (2.8) is of limited The solution to the diffusion equation with ini-
use in many areas of rapid uplift because the re- tial condition given by Eq. (2.16) is the summa-
quirements for the diffusion equation generally tion of Eq. (2.15) from n = 1 to n = ∞:
do not apply. Later in this chapter we will solve

nπ x ’κn2 π 2 t/L 2
the solution to the nonlinear diffusion equation.
h(x, t) = an sin e (2.18)
That equation usually provides a more accurate n=1
representation of hillslope evolution in rapidly
Substituting Eq. (2.17) into Eq. (2.18) gives the gen-
uplifting terrain.
eral solution

2 nπ x
2.2.2 Fourier series method h(x, t) = f (x ) sin
Two general classes of solutions exist for the dif- 0 n=1
nπ x ’ κn2 π 2 t
fusion equation: series solutions and similarity
— sin e L2 dx (2.19)
solutions. First consider the one-dimensional (1D) L
diffusion equation in a region bounded by x = 0
This approach suggests a powerful means of
and x = L . By 1D, we mean a pro¬le that varies
solving the diffusion equation in bounded re-
with only one independent spatial variable, x.
gions: using the initial condition f (x), solve for
Let™s assume that the initial condition at t = 0
the coef¬cients an of the Fourier series repre-
is given by a function f (x) and that the value of
sentation of f (x); then, plug the an values into
h is set to zero at both ends of the region for
Eq. (2.18) to obtain the solution. Other boundary
all time t (i.e. h(0, t) = h(L , t) = 0). If the initial
conditions (e.g. values of h(0, t) and h(L , t) that
condition f (x) is given by
are nonzero, or boundary conditions speci¬ed by
nπ x hillslope gradients) can also be handled by using
f (x) = an sin (2.14)
L cosine, linear, and/or constant terms in the sum-
mation. For example, if we wish to place a divide
then the solution to the diffusion equation is
at x = 0 then the series solution is
nπ x ’κn2 π 2 t/L 2
h(x, t) = an sin ∞
e (2.15) L
1 nπ x
h(x, t) =
L f (x ) cos
L 2L
’L n=1
The correctness of this solution can be checked
nπ x ’ κn2 π22 t
— cos
by differentiation and substitution of the re- e 4L dx (2.20)
sults into the diffusion equation. This solution
is only valid for the sinusoidal function given by Equation (2.20) is the same as Eq. (2.19) ex-
Eq. (2.14), so it may not appear to be of much use. cept that cosine terms are used and the hill-
slope domain is assumed to be from x = ’L
Equation (2.15), however, is actually a powerful
to x = L . If f (x) is chosen to be symmetric (i.e.
building block that can be used to solve the dif-
f (x) = f (’x)) then Eq. (2.20) is consistent with
fusion equation for nearly any initial condition
by using a Fourier series representation for f (x). the divide boundary condition given by Eq. (2.13).

The Fourier series approach is only applicable integration. The approach is to introduce a new
to bounded regions. Of course, all geomorphic variable that is a combination of x and t:
applications are bounded, because the Earth is x
·= √ (2.22)
2 κt
not in¬nite in extent! In some cases, however, it
is useful to de¬ne mathematical models in in¬- and rewrite the diffusion equation for θ and its
nite or semi-in¬nite domains. For example, if one boundary conditions in terms of ·. This can be
is interested in modeling the ¬rst few years of done (following the approach in Turcotte and
radionuclide migration into an alluvial deposit Schubert (1992)) using the chain rule for differ-
that is hundreds of meters thick, a semi-in¬nite entiation:
model is more convenient to work with and prac-
‚θ ‚θ ‚· ‚θ 1·
1x1 dθ
= = ’√ = ’
tically speaking more accurate than a Fourier
‚t ‚· ‚t ‚· 4 κt t d· 2t
series approach.
The Fourier series approach works because
‚θ dθ ‚· dθ 1
= = √
the diffusion equation is linear (i.e. it does not (2.24)
‚x d· ‚ x d· 2 κt
involve squared or higher-order terms of h or
‚ 2θ 1 d2 θ ‚· 1 1 d2 θ
its derivatives). Only in linear equations can dif- =√ = (2.25)
‚ x2 2 κt d·2 ‚ x 4 κt d·2
ferent solutions be superposed to obtain other
The diffusion equation for θ becomes
1 d2 θ

’· = (2.26)
2.2.3 Similarity method 2 d·2

Series solutions with ¬xed-elevation boundary
and the boundary conditions become θ(∞) = 0
conditions are appropriate for hillslopes with a
and θ (0) = 1. Equation (2.26) is a second-order or-
well-de¬ned base-level control (such as a nearby
dinary differential equation (ODE) that can be
channel that carries sediment away from the
reduced to a ¬rst-order ODE by introducing the
slope base to maintain a constant elevation). In
some cases, however (such as when sediment is

allowed to deposit and prograde at the slope base) φ= (2.27)

it is more appropriate to assume that diffusion
occurs over an in¬nite or semi-in¬nite region. In Equation (2.26) then becomes
such cases, similarity solutions are a powerful 1 dφ
’·d· = (2.28)
approach. Let™s consider the semi-in¬nite region 2φ
given by x = 0 to x = ∞ with an initial condition
Equation (2.28) can be integrated to obtain
h(x, 0) = h 0 . Let™s also assume that at t = 0 the el-
’·2 = ln φ ’ ln c 1
evation h is instantaneously lowered to zero at (2.29)
x = 0. In the context of hillslope evolution, this
Taking the exponential of both sides of Eq. (2.29)
problem corresponds to instantaneous base level
drop of an initially planar hillslope or terrace.

In order to solve this problem, it is useful to 2
φ = c 1 e’· = (2.30)
introduce a dimensionless elevation given by d·
Integrating Eq. (2.30) gives
θ= ’1 (2.21) ·
h0 2
e’· d· + 1
θ = c1 (2.31)
The diffusion equation for θ is identical to the 0

where · is an integration variable and the con-
diffusion equation for h. The boundary condi-
tions on θ are now given by θ(x, 0) = 0, θ(0, t) = 1, dition θ (0) = 1 was used to constrain the integra-
and θ(∞, t) = 0. Similarity solutions make use of tion constant. Application of θ(∞) = 0 requires
a mathematical trick that reduces the diffusion that

equation (a partial differential equation) to an 2
e’· d· + 1
0 = c1 (2.32)
ordinary differential equation that is solved by 0

The integral in Eq. (2.36) is given by π/2, so the on h. The initial condition for w is given by

constant c 1 has the value ’2/ π. Equation (2.31)
w(x, 0) = ’ (L ’ x 2 )
becomes (2.39)

2 2
e’· d·
θ =1’ √ (2.33)
π The solution for w is found by substituting
Eq. (2.39) into Eq. (2.20) to obtain
The integral in Eq. (2.33) has a special name called
the error function: ∞
(2n + 1)π x
32 22
’ κn π t
w(x, t) = 3 cos e
· 2
2 π (2n + 1) 3
’· 2 2L
erf(·) = √ e d· (2.34) n=0
π 0
The solution for θ, therefore, is
Substituting Eq. (2.40) into Eq. (2.38) and re-
θ = 1 ’ erf(·) = erfc(·) (2.35)
arranging for h gives
where erfc(·) is called the complementary error

U L2 x2 (’1)n
function. Written in terms of the original vari- h(x, t) = 1’ 2 ’ 3 2
πL (2n + 1)3
2κ L
ables, Eq. (2.35) becomes n=0

(2n + 1)π x 22
’ κn π t
x — cos e (2.41)
h(x, t) = h 0 erf √
(2.36) 2L
2 κt
Most solutions of the diffusion equation in in¬- Examples of Eq. (2.41) are plotted in Fig-
ure 2.3a for values of κt ranging from 1000 m2
nite and semi-in¬nite domains can be written in
to 2500 m2 . These results provide an estimate for
terms of error functions.
Series solutions and similarity solutions are the time scale necessary to achieve steady state.
If we assume κ = 1 m2 /kyr, then steady state is
the two basic types of analytical approaches to
solving the diffusion equation. In the following achieved after approximately 2.5 Myr of uplift for
sections we will apply these equations to speci¬c a hillslope 50 m in length. Interestingly, the time
problems that arise in geomorphic applications. scale required to develop steady state does not
depend on the uplift rate.
2.2.4 Transient approach to steady state Figure 2.3b presents the same results as Fig-
In Section 2.2.1 we found the steady-state diffu- ure 2.3a except that the hillslope length has been
sive hillslope to be given by Eq. (2.13). Here we normalized to 1. In this approach the results de-
solve the time-dependent diffusion equation that pend only on a single nondimensional parameter
κt/L 2 . This nondimensional framework is more
describes the transient approach to steady state.
We consider a diffusing hillslope of length L with ¬‚exible because we do not need to plot differ-
constant uplift rate U : ent solutions for hillslopes of different lengths.
While Figure 2.3a is speci¬c to hillslopes 50 m
‚h ‚ 2h
’κ 2 =U (2.37) in length, Figure 2.3b allows us to estimate the
‚t ‚x
time required to achieve steady state for any hill-
The boundary conditions for this problem are
slope length. If κ = 1 m2 /kyr and L = 100 m, for
given by Eqs. (2.13) and (2.12). The initial condi-
example, steady state is achieved (i.e. the relief
tion is h(x, 0) = 0. The constant uplift term on
has achieved 90% of the steady-state relief) after
the right side of Eq. (2.37) can be eliminated by
κt/L 2 = 1 or t = 10 Myr.
introducing a new variable w that quanti¬es the
In this chapter we will usually present results
difference from steady state. The variable w is de-
in terms of actual length scales in order to de-
¬ned by
velop intuition about the rates of landscape evo-
U2 lution. In practice, however, it is often best to
w(x, t) = h(x, t) ’ (L ’ x 2 ) (2.38)
2κ solve the problem in a fully nondimensional way
and then apply speci¬c values of κ, L , and t to
The evolution of w is governed by the diffusion
equation with the same boundary conditions scale the results appropriately.

(a) 1.0 (a)
2500 m2
2250 m2

UL2 Q2a
0.4 Q2b
= 250 m2
0 10 20 30 40 50 (b)
x (m)
(b) 1.0
kt/L2 =
Q2a Q2b

kt/L2 = 0.1
0.0 0.2 0.4 0.6 0.8 1.0


. 2
( 14)