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

18 INTRODUCTION

readily transported across the fan. These distal-

100

fan eolian corridors are controlled by a thresh-

P20

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

10’1

threshold value, sand is trapped locally in low 102

100

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

(a)

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

R’1

stress necessary to move the particle past down- VP

A

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

(c)

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

1.3 A TOUR OF THE GLACIAL SYSTEM 19

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

20 INTRODUCTION

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

1971).

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

rates.

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

1.3 A TOUR OF THE GLACIAL SYSTEM 21

divide

accumulation

equilibrium

line

ablation

flowlines

ice margin

basal velocity

moraines

bedrock

till plane

predominant erosion little or no erosion or deposition

predominant

(frozen bed or low velocity)

deposition

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

‚h

= ’aV b

increase toward the margin, transporting all of the bedrock (1.19)

‚t

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

22 INTRODUCTION

a Newtonian ¬‚uid, strain rate is proportional to

4

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:

dz

2 n

„

du

=a (1.20)

„0

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

0

parameters. The values of a and n can be consid-

1.0 1.5 2.0

0.5

0.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

1.3 A TOUR OF THE GLACIAL SYSTEM 23

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

24 INTRODUCTION

(a)

105

(b) resolution limit

104

’1

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

103

102 uncompensated,

"great" lakes

compensated lakes

101

100

104 105

103

102

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.

1.3 A TOUR OF THE GLACIAL SYSTEM 25

drumlins

moraines

eskers

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

26 INTRODUCTION

high peaks

plateau surface

U-shaped valley

hanging valley glacial

moraines

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

1.3 A TOUR OF THE GLACIAL SYSTEM 27

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-

cirques

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

accumulation

ablation

evolves to

ELA

ice

bedrock

V

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

28 INTRODUCTION

P-T Granite

Ordovician

Metasediments

Icecover

30 km

Huayna Potosi

P-T Granite

Silurian

Metasediments

(a) (b)

modern ELA

LGM E LA

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

1.4 CONCLUSIONS 29

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

narrow

peak

ELA

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

processes.

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-

signi¬cant.

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

‚h

=’

q = ’ρκ (2.2)

(2.1)

‚t ρ ‚x

‚x

2.1 INTRODUCTION 31

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).

(2.3)

‚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

=c

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

32 THE DIFFUSION EQUATION

(a) (d)

q

1 1

h h

(m) erosion

(m)

x1 x2

deposition

x3 x4

0 0

0.0 0.0

(b) (e)

kt = 25 m2

kt = 25 m2

= 0.8

= 2.5

0.003

0.0005

(c) (f)

2

2

2

2

0.0

0.0

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.

2.1 INTRODUCTION 33

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-

nonlinear

q

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-

Sc

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

34 THE DIFFUSION EQUATION

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

by

context of hillslope evolution. As such, the dif-

‚h

fusing variable h will represent elevation along =0 (2.10)

‚x

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

U2

h(x) = ’ x + c2 (2.11)

(for a particular initial condition and boundary

2κ

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

base:

2.2 Analytic methods and

h(L ) = 0 (2.12)

applications

then c 2 = U L /2κ and the hillslope pro¬le is given

by

2.2.1 Steady-state hillslopes

U2

Steady-state landscapes have received a great deal h(x) = (L ’ x 2 ) (2.13)

2κ

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

2.2 ANALYTIC METHODS AND APPLICATIONS 35

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)

L

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

L

2 nπ x

still relatively low and weathering rates are suf- an = f (x ) sin dx (2.17)

L L

¬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)

L

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

∞

L

2 nπ x

2.2.2 Fourier series method h(x, t) = f (x ) sin

L L

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)

2L

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).

36 THE DIFFUSION EQUATION

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.

(2.23)

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

solutions.

1 d2 θ

dθ

’· = (2.26)

2.2.3 Similarity method 2 d·2

d·

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

variable

some cases, however (such as when sediment is

dθ

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

d·

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

gives

drop of an initially planar hillslope or terrace.

dθ

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

h

θ= ’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

2.2 ANALYTIC METHODS AND APPLICATIONS 37

√

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)

U2

w(x, 0) = ’ (L ’ x 2 )

becomes (2.39)

2κ

∞

2 2

e’· d·

θ =1’ √ (2.33)

π The solution for w is found by substituting

0

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

(’1)n

32 22

’ κn π t

w(x, t) = 3 cos e

· 2

4L

2 π (2n + 1) 3

’· 2 2L

erf(·) = √ e d· (2.34) n=0

π 0

(2.40)

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

32κ

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)

2

h(x, t) = h 0 erf √

4L

(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.

38 THE DIFFUSION EQUATION

(a) 1.0 (a)

=

kt

=

kt

2500 m2

2250 m2

0.8

0.6

2hk

UL2 Q2a

0.4 Q2b

Q2c

Q3

0.2

= 250 m2

kt

0.0

0 10 20 30 40 50 (b)

x (m)

(b) 1.0

=

kt/L2 =

kt/L2

1.0

0.8

0.9

0.6

2hk

UL2

Q2a Q2b

0.4

0.2

kt/L2 = 0.1

0.0

0.0 0.2 0.4 0.6 0.8 1.0

x/L