width. Over time, large depressions deepen and

to compare spatial domains (e.g. lakes, landslides,

expand more quickly than smaller depressions

vegetation clusters) modeled on a computer with

nearby. Large depressions expand faster because

154 NON-NEWTONIAN FLOW EQUATIONS

(b)

(a)

104

(c)

103

N = cA

N (>A)

102

(d) (e) 101

100

101

100 102

10 10

A (km 2)

Fig 6.25 Comparison of modeled and observed topography

resulting from glacial erosion beneath ice-sheet interiors. The

Great Lakes. In this section we refer to these

pixel size in the model is assumed to be 100 m. (a) Shaded

lakes collectively as the ˜˜great™™ lakes of Canada.

relief of bedrock topography with ¬lled depressions. (b) A

binary image of ¬lled depressions (i.e. lakes) from (a). (c) The These lakes are roughly uniformly-spaced in a

cumulative frequency-size distribution of lakes in (b), well ¬t ring around Hudson Bay, with distances between

by the same relationship observed for Canadian lakes less lake centers ranging from 370 km to 1000 km.

than 104 km2 in area: N(> A ) ∝ A ’1 . (d) Shaded relief of

The map of Figure 6.26a was created by extract-

topography in central Canada for qualitative comparison with

ing all the lakes from a 1-km resolution DEM

(a). (e) Binary image of lakes in (d).

of Canada and the northern US. The shading

of each lake is proportional to elevation, with

black indicating lakes near sea level and light

those of populations observed in nature. The cu- gray shades representing lakes as high as several

mulative frequency-size distribution of lakes is hundred meters above sea level. The shading val-

the number of lakes larger than a given area. ues indicate, for example, that the great lakes of

Figure 6.25c illustrates the cumulative frequency- western Canada are separated from Hudson Bay

size distribution for the lakes of Figure 6.25e. The by a broad ridge of topography several hundred

data ¬t a power-law relationship N (> A) ∝ A ’1 . meters above sea level and parallel to the Hudson

This compares well with the frequency-size dis- Bay shoreline.

tribution of small and medium-sized lakes in To investigate the role that ¬‚exure may have

Canada shown in Figure 1.29. in modifying the basic glacial-erosion instability,

we performed a numerical experiment of glacial

6.6.2 Length scales > ¬‚exural wavelength erosion beneath an initially semi-circular ice

The largest lakes in Canada and the US are il- sheet. The 2D ¬‚exure of the lithosphere beneath

lustrated in Figure 6.26a, including (from north- the ice sheet was included in this model ex-

periment. The grid for this model is 3000 km —

west to southeast) Great Bear Lake, Great Slave

Lake, Lake Athabaska, Lake Winnipeg, and the 3000 km with a pixel resolution of 10 km. The

6.6 GLACIAL EROSION BENEATH ICE SHEETS 155

elastic thickness is 50 km. As in Section 6.6.1,

(a) Great Bear

we assumed a ¬‚at initial surface with white

noise variations as the initial bedrock surface.

600 km

Great Slave

The ice-surface topography, ¬‚ow, and bedrock

de¬‚ection in the ¬rst model iteration are given

Athabasca

370

in Figures 6.27a--6.27c. The accumulation pattern

Reindeer

was assumed to be uniform for elevations greater

than 1 km.

1000

The ice-surface topography, ¬‚ow, bedrock de-

Winnipeg

¬‚ection, and bedrock topography of the model

after ten iterations are given in Figures 6.27d--

Superior

950

6.27g. The bedrock topography does not include

500

400 the de¬‚ection under the ice load. Therefore, it

represents what the landscape would look like

400

if the ice load was removed and only the effect

of glacial erosion remained. In this experiment,

(b)

the largest lakes are uniformly spaced in a ring

around the center of the ice sheet, directly un-

derneath the equilibrium line. In the center of

100

the ice sheet the original rough topography re-

80

mains because erosion is minimal there. Near the

60

120 km

equilibrium line the glacial instability is ampli-

¬ed by the de¬‚ection of the lithosphere at wave-

lengths equal to ≈ 500 km. Bedrock depressions

that form at wavelengths of 500 km are ampli¬ed

120 km

faster than smaller wavelengths because their ice

100

80

in¬ll de¬‚ects the lithosphere, further depressing

60

the ice surface and focusing ice drainage into the

depressions.

A grayscale image and contour map of the

elastic thickness of central Canada is given in Fig-

ure 6.26b. The elastic-thickness estimates were ob-

Fig 6.26 The “great” lakes of Canada and their relationship

tained by Tony Watts (personal communication,

to the ¬‚exure of the lithosphere. (a) Topographic depressions

2002) and are based on the coherence method be-

of Canada (all non-white areas) with grayscale shading

mapped to elevation (lakes in black are at sea level, lighter tween gravity and topography. The greatest elas-

tones have water surfaces at higher elevations. The distances tic thickness in Canada exists just southwest of

between the largest lakes are indicated. (b) Grayscale and

Hudson Bay with a maximum value of nearly

contour map of the elastic thickness of the lithosphere from

140 km. To the northwest of this maximum value,

Tony Watts (personal communication, 2002), recti¬ed to (a).

following the lines adjoining the great lakes

The elastic thicknesses are greatest in the areas between

from Winnipeg to Athabasca to Great Slave and

Athabasca Lake, Lake Winnipeg, and Lake Superior, where

Great Bear, elastic thicknesses decrease gradu-

the distances between lakes are also the greatest. The

ally to about 100 km at Athabasca, then more

distances separating lakes are lower in the northwest and

rapidly to about 60 km near Great Slave. Moving

southeast portion of Canada where the elastic thickness

values fall to approximately one half their values in central in the southeast direction from Winnipeg, elas-

Canada. If Reindeer Lake is considered to be a great lake, tic thicknesses also decrease, to 100 km at Supe-

however, the relationship between lake spacing and elastic

rior and 60 km at Ontario. This pattern suggests

thickness is less robust.

that the distances between the great lakes may

correlate with the local elastic thickness, with a

156 NON-NEWTONIAN FLOW EQUATIONS

(d)

(a)

ice-surface topography

iteration10

iteration 1 ice-surface topography

(b) (e)

basal sliding velocity

basal sliding velocity iteration 10

iteration 1

(f)

(c)

ice margin

forebulge

lighting direction

lithospheric deflection

lithospheric deflection iteration 10

iteration 1

(g)

500 km

iteration10 bed topography

Fig 6.27 Numerical experiment illustrating the role of

¬‚exure in lake formation. (a) Shaded-relief image of

Is the relationship between great-lake spac-

ice-surface topography, (b) grayscale image of basal-sliding

ing and elastic thickness consistent with ¬‚ex-

velocities, and (c) shaded-relief image of lithospheric

de¬‚ection in the ¬rst iteration of the model. Basal sliding ure? The elastic thickness is related to the three-

velocities and bedrock erosion are diffuse in this early stage dimensional ¬‚exural parameter by the relation-

of the model. (d) Ice-surface topography, (e) basal-sliding

ship:

velocities, (f) bedrock de¬‚ection, and (g) bedrock topography

1

(no de¬‚ection component) after 10 iterations of the model

D 4

β= (6.53)

run. Domain is 3000 km in length and elastic thickness is

ρg

50 km.

where D is the ¬‚exural rigidity given by

E T e3

D=

proportionality factor of 10. For example, in ar- (6.54)

12(1 ’ ν 2 )

eas where the elastic thickness is 100 km, the dis-

tance between lakes is ≈ 1000 km. It should be E = 70 GPa, ν = 0.25, ρ=

Assuming

2300 kg/m , and g = 10 m/s , and an elastic

3 2

noted that if Reindeer Lake is considered as a

thickness T e = 50 km, the ¬‚exural parameter β

great lake in this analysis, the relationship be-

tween great-lake spacing and elastic thickness is is 77 km. The wavelength of de¬‚ection is related

to the ¬‚exural parameter by » = 2πβ = 480 km,

not as robust.

6.6 GLACIAL EROSION BENEATH ICE SHEETS 157

or about ten times the elastic thickness. This (a) profile 1

analysis indicates that the observed relationship

between great-lake spacing and elastic thickness

profile 2

is consistent with ¬‚exural control.

Although the relationship between lake spac-

ing and elastic thickness suggests that ¬‚exure

is the key process controlling great-lake forma-

tion, the ice margin may be home to many

feedback processes that lead to instabilities. In (b)

addition to the focusing of ¬‚ow by the bed topog-

raphy (through its in¬‚uence on the ice-surface

topography), feedbacks between ice temperature,

¬‚ow, meltwater content, and other variables may

10 km

also be at work. For example, the EISMINT model

comparison project noted that an ice-marginal

basal velocity

0.7 bars

instability driven by a feedback between ice tem-

perature, ¬‚ow, and ice-sheet geometry was re- (c)

produced in all of the models tested (Payne et

al., 2000). Since ice temperature is not a vari-

able in our model, the instability found in the

EISMINT experiments is independent from the

¬‚exural mechanism of our model. The wave-

5“7 km

length of the thermomechanical instability ob-

basal velocity

1.5 bars

served in the EISMINT experiments varied be-

1.5 km 960 m

tween model experiments but was generally profile 2

profile 1

several hundred kilometers. As such, a ther-

momechanical instability could also generate

˜˜great™™ lakes.

950 m

0

150 km

50 km 0

0

6.6.3 Near ice margins

Finger lakes are elongate glacially carved lakes Fig 6.28 Ice-surface topography and ¬‚ow above a rough

surface with „ = 0.7 bars. (a) Shaded-relief image of an ice

formed on the margins of ice sheets. The type ex-

sheet 150 km wide and 50 km from divide (top of image) to

ample is the Finger Lakes Region of New York

margin. Two pro¬les of the ice-surface topography are given:

State. In this chapter we use the term more

pro¬le 1 is along-dip from divide to margin; pro¬le 2 is

broadly to include elongate glacial troughs near

along-strike from left to right. (b) Grayscale image of basal

the ice margin whether or not they impound velocity corresponding to (a) with an ELA of 1 km. Ice-¬‚ow

water. For example, the Finger Lakes Region in- channels with a characteristic spacing of 10 km are present in

cludes ¬ve major lakes, several smaller lakes, and the ¬‚ow even though there is no characteristic scale to the

bed topography. (c) A thicker ice sheet („ = 1.5 bars) results

dozens of glacial valleys of similar shape. By ˜˜¬n-

in narrower ¬‚ow channels (assuming the ELA is the same as

ger lakes™™ we will be referring to these features

in (b)).

collectively. In this section we present model re-

sults for the ice-sheet geometry, basal velocities,

and bedrock erosion near linear and curved ice

margins. We show that uniformly-spaced ice-¬‚ow 6.28a and 6.28b. The basal shear stress in this

channels result from ¬‚ow over a rough surface case is assumed to be 0.7 bars and the ELA is

near the ice margin. 1 km. Also given in Figure 6.28a are topographic

The ice-surface topography and basal-sliding pro¬les parallel and perpendicular to the mar-

velocities for an ice sheet sliding over a rough gin. Pro¬le 1 shows the parabolic pro¬le charac-

surface near a linear margin are given in Figures teristic of ice sheets on ¬‚at beds (e.g. Nye, 1951).

158 NON-NEWTONIAN FLOW EQUATIONS

(a) Fig 6.29 Numerical experiment

illustrating the model behavior of

focused erosion near a curved ice

margin. (a) Shaded-relief image of

ice-surface topography, (b) grayscale

image of basal-sliding velocities at

iteration 1 ice-surface topography

the beginning of the model run. (c)

(b) Ice-surface topography, (d)

basal-sliding velocities, and (e) bed

topography of the model after ten

iterations. Finger lakes have formed

with orientations parallel to the

iteration 1 basal sliding velocity ice-¬‚ow directions and lake spacing

is low in zones of converging ice and

(c)

high in zones of diverging ice.

iteration 10 ice-surface topography

(d)

iteration 10 basal sliding velocity

(e)

iteration 10 bed topography

is to create strongly localized ¬‚ow near the ice

Pro¬le 2 shows the microtopography of the ice

margin even though the ice-surface topography

surface along a direction parallel to the margin.

is nearly uniform along the margin. The width

One of the most important ¬ndings of the chap-

of ice-¬‚ow channels depends on the slope of the

ter is illustrated in Figure 6.28: even though the

ice sheet at the ELA. A thicker ice sheet with the

bed topography has no characteristic scale (i.e. it

same ELA (1 km), for example, results in more

is white noise), the resulting ice-surface topogra-

closely-spaced channels, To illustrate this, Figure

phy has a characteristic scale of about 10 km. This

6.28c is a basal-velocity map for an ice-sheet re-

characteristic scale is even more apparent in the

construction with „ = 1.5 bars. Ice-¬‚ow channels

map of basal-sliding velocities (Figure 6.28b). Fo-

in this case are spaced at ≈ 5--7 km instead of

cusing of the ice ¬‚ow occurs as a result of the ad-

10 km.

ditive de¬‚ection of ice ¬‚ow along ¬‚ow paths from

How does the focusing effect differ for a

the divide to the margin. This focusing involves

curved ice margin? Figures 6.29a and 6.29b illus-

several steps. First, the rough bed is re¬‚ected

trate the ice-surface topography and basal-sliding

in the ice surface as microtopographic channels

velocities for a curved margin. The shape of this

(Figure 6.28a). Second, these channels act to focus

margin is identical to the 14 ka ice margin of the

the ¬‚ow a little bit as the ice is routed from pixel

Finger Lakes Region. The bed topography in this

to pixel. The cumulative effect of this focusing

6.6 GLACIAL EROSION BENEATH ICE SHEETS 159

(a) Fig 6.30 Spacing of ¬nger lakes

and troughs in upstate New York.

(a) Shaded-relief DEM image of

Finger Lakes area. White curve

indicates approximate position of

the Valley Heads Moraine (based on

the Finger Lakes and Niagara sheets

of the sur¬cial geologic map of New

York (New York Geological Survey,

1999)). Black lines indicate the

location of the topographic transect

(b)

5 km 10 5 5

20 given in (b). The transect shows

600 m

alternating zones of narrow (≈ 5 km)

and wide (up to 20 km) valley

elevation spacing. Zones of narrow spacing

correspond to regions of converging

¬‚ow in ice embayments. Conversely,

200 m

distance along profile wider spacings correspond with

270 km

0

areas of diverging ice ¬‚ow.

gin was coincident with this moraine, although

example was also chosen to be a ¬‚at surface with

several phases of glaciation may have contributed

white noise variations. The curved margin results

to their formation.

in ice-¬‚ow channels that are alternately closely

Figure 6.30b is the topographic pro¬le of the

spaced and widely spaced. Channels are more

region along the black line of Figure 6.30a. The

closely spaced in margin embayments where ice

pro¬le of Figure 6.30b shows the same pattern as

¬‚ow is convergent. Widely spaced channels occur

the model results of Figures 6.29c--6.29e: closely

beneath ice lobes where ¬‚ow is divergent. Fig-

spaced (≈ 5 km) troughs in margin embayments

ures 6.29c--6.29e illustrate the ice-surface topog-

and widely spaced (10--20 km), deeper troughs be-

raphy, basal velocities, and bedrock topography

neath ice-sheet lobes.

after 10 model iterations. The ice ¬‚ow channels

If ¬nger lake formation is an intrinsic feature

have carved distinct glacial troughs, depressing

of ice margins, why don™t ¬nger lakes occur more

the ice-surface topography and focusing ¬‚ow and

commonly in formerly glaciated terrain? One pos-

erosion into the troughs.

sible reason is that intense glacial erosion and a

Geomorphically, the Finger Lakes Region

stable ice margin are both required for ¬nger-

is dominated by subparallel, glacially scoured

lake formation. If the equilibrium line is rapidly

troughs with southernmost extents in Seneca

migrating during ice-sheet advance or retreat, for

and Cayuga Lakes, the two largest of the Fin-

example, there may not be suf¬cient time for the

ger Lakes. The ¬ve largest troughs of the region

instability between depressions, ice ¬‚ow, and ero-

comprise the Finger Lakes proper, but there are

sion to develop deep troughs. It may be that the

numerous other troughs cut into the Allegheny

Finger Lakes were formed because 14 ka was a

Plateau of smaller size that are not deep enough

period of suf¬cient margin stability for erosion

to impound water. Figure 6.30a is a shaded-

to carve the topography we see today. Alterna-

relief image of the topography of the region

tively, the margin may not have been unusually

with the location of the 14 ka ice margin (the

stable, but erosion may have been particularly in-

Valley Heads Moraine) shown in white. The

tense during this period. For example, meltwater

troughs vary in spacing from 5 to 20 km along

pulses may have driven rapid erosion just prior

strike with the greatest spacing between Seneca

to 14 ka followed by rapid ¬lling of the Finger

and Cayuga Lakes. The southern tips of the Fin-

Lakes between 14 and 13 ka. This hypothesis is

ger Lakes coincide with the Valley Heads Moraine

consistent with the lake stratigraphy of the re-

(14 ka). As such, the scouring of the Finger Lakes

gion (Mullins and Hinchey, 1989).

Region most likely took place when the ice mar-

160 NON-NEWTONIAN FLOW EQUATIONS

and compute pro¬les for each position by numeri-

cally integrating from the terminus up the valley.

6.3 Using the code in Appendix 5 as a guide,

4m model the evolution of a 1D gravity ¬‚ow with

temperature-dependent viscosity ¬‚owing down an

inclined plane of slope S.

6.4 Download a DEM of an area formerly covered by

alpine glaciers. Using the sandpile algorithm code

10 m given in Appendix 5 as a guide, map the thick-

ness of glacier ice assuming a basal shear stress

Fig 6.31 Schematic diagram of lava-¬‚ow margin in Exercise

of 1 bar. Choose an ELA value and assume that ice

6.1.

covers topography only above the ELA.

6.5 The Boussinesq equation

‚h kρg ‚ ‚h

= h (E6.1)

‚t μφ ‚ x ‚x

Exercises

6.1 Lava-¬‚ow margins in the vicinity of Lathrop Wells is used to quantify the evolution of a water table

cone, Amargosa Valley, Nevada, are approximately of height h in an uncon¬ned aquifer of perme-

ability k and porosity φ, where ρ is the density of

10 m wide and 4 m thick in areas of ¬‚at terrain

water, g is gravity, and μ is the dynamic viscos-

(Figure 6.31). Estimate the yield stress for these

¬‚ows based on a 1D perfectly plastic model. ity of water (Turcotte and Schubert, 2002). Plot

6.2 Use a topographic map or DEM to extract a lon- the shape of the steady-state water table in a 1D

gitudinal pro¬le of a steep ¬‚uvial valley. Import aquifer of length 10 km above an impervious stra-

tum with a dip of 1%. Assume k =10’12 m2 and

the pro¬le in two-column format x, h into Ex-

μ =10’3 Pa s. (Hint: map this problem onto the

cel. Using Excel, model the pro¬le of a slow-

moving debris ¬‚ow with yield stress 0.1 bars ¬‚ow- perfectly plastic model for alpine glaciers. See Sec-

ing down the valley bottom using the perfectly tion 9.5 of Turcotte and Schubert for more details

plastic model. Choose several terminus positions on the Boussinesq equation.)

Chapter 7

Instabilities

to spatial and temporal oscillations in arid al-

7.1 Introduction luvial channels. The study of instabilities does

not involve any fundamentally new types of

equations. Rather, most equations involving ge-

Many geomorphic processes involve positive feed-

omorphic instabilities combine diffusion- and

backs. The incision of a river valley into an

advection-type equations similar to those we

undissected landscape is an example of a posi-

have already considered. In most cases, how-

tive feedback process in which an incipient chan-

ever, two or more variables are coupled in such

nel erodes its bed and expands its drainage area

a way as to enhance each other in a positive

in a positive feedback. In some cases, positive-

feedback.

feedback processes give rise to landforms with

This chapter does include one fundamentally

no characteristic scale. Many drainage networks,

new tool: the linear-stability analysis. In this anal-

for example, are branching networks with no

ysis, the fundamental equations that describe

characteristic scale above the scale that de¬nes

the evolution of a particular system are ¬rst lin-

the hillslope-channel transition. In other cases,

earized (i.e. nonlinear terms are neglected) and

periodic landforms are created (Hallet, 1990).

then solved for the case of an initially small

Sand dunes, for example, form by a positive

perturbation of wavelength » (where » is a free

feedback between the height of the dune, the

parameter). Conceptually, a linear stability ap-

de¬‚ection of air ¬‚ow above it, and the result-

proach considers the evolution of the system

ing pattern of erosion and deposition (i.e. en-

with small, incipient landforms with a wide

hanced windward erosion and leeward deposi-

range of sizes (i.e. wavelengths) and determines

tion). Unstable geomorphic systems all include

whether and how fast each perturbation will

some type of positive feedback. In many arid al-

grow to become full-sized landforms. For those

luvial channels, for example, sediment becomes

wavelengths that do grow, the wavelength that

stored preferentially at certain points along the

grows the fastest provides an estimate of the char-

channel pro¬le, thereby causing the channel to

acteristic scale of the landform as a function of

aggrade and the channel bed to widen. Chan-

the system parameters. This fastest-growing wave-

nel widening further enhances sediment stor-

length often, but not always, provides a good es-

age and bed aggradation in a positive feedback.

timate of the spatial scale of full-sized landforms.

This instability continues until a suf¬ciently

As such, linear stability analyses often provide a

steep slope develops downstream of the aggrad-

useful starting point for understanding the emer-

ing zone, triggering channel entrenchment and

gence of characteristic scales in complex geomor-

narrowing. The positive feedback between chan-

phic systems.

nel width and localized erosion/deposition leads

162 INSTABILITIES

y = ’b

r1, m

w x

y=0

r2, m

y=b

y

Fig 7.2 Geometry of the Rayleigh“Taylor instability of a

dense ¬‚uid overlying a lighter ¬‚uid. After Turcotte and

Schubert (2002). Reproduced with permission of Cambridge

University Press.

Fig 7.1 Surface expression of salt domes in the Gulf of

Mexico, from an oblique perspective image of high-resolution

bathymetry. The topography has a hummocky shape that

layers is initially centered on y = 0 but has a si-

re¬‚ects the underlying characteristic scale and spacing of salt

domes. Image reproduced with permission of Lincoln Pratson. nusoidal variation given by

2π x

w(x) = A 1 cos (7.1)

»

7.2 An introductory example: the The goal of the linear stability analysis is to

Rayleigh“Taylor instability predict how w varies with time and to deter-

mine the wavelength » associated with the fastest

growth.

As an introduction to the study of instabilities

The ¬rst step of the Rayleigh--Taylor analysis

in geomorphology, we consider the Rayleigh--

is to solve for the ¬‚ow ¬eld in the two-layer sys-

Taylor instability of a dense ¬‚uid overlying a

tem given a growing interface given by Eq. (7.1).

lighter ¬‚uid. The Rayleigh--Taylor instability de-

The ¬‚ow ¬elds are characterized by four velocity

scribes the formation of salt domes in sedimen-

functions: u1 , v 1 , u2 , v 2 . These functions are the

tary basins (Figure 7.1) and it provides a simple

horizontal (u) and vertical (v ) components of the

model for the initiation of mantle plumes. When

velocity in the upper (1) and lower (2) layers. Tur-

a dense layer overlies a lighter layer, the system

cotte and Schubert (2002) used the stream func-

is gravitationally unstable, i.e. the lighter ¬‚uid

tion method to obtain expressions for these ¬‚ow

will rise and penetrate through the upper layer.

velocities. In this method, the velocity functions

Here we will quantify how this instability works

u and v are obtained by differentiating a single

following closely the approach of Turcotte and

˜˜stream function™™ ψ:

Schubert (2002).

The geometry of the two-layer system is shown

‚ψ

u=’

in Figure 7.2. A ¬‚uid with a thickness b and den- (7.2)

‚y

sity ρ1 overlies a lighter ¬‚uid of thickness b and

‚ψ

density ρ2 . Here we assume that both ¬‚uid layers

v= (7.3)

have the same viscosity μ. We also assume that ‚x

both the top and bottom of the system have rigid

surfaces. The interface between the two ¬‚uid Turcotte and Schubert™s (2002) analysis for this

7.2 AN INTRODUCTORY EXAMPLE: THE RAYLEIGH“TAYLOR INSTABILITY 163

problem gives

» »

2π x 2π y 2π x y 2π b 2π y y 2π y 2π y

ψ1 = A 1 sin + A 1 sin + ’

cosh tanh sinh cosh sinh

» » » » » » »

b 2π b b 2π b

’1

2

» »

1 1 2π b

— + ’ tanh (7.4)

»

sinh 2π b cosh 2π b sinh 2π b cosh 2πb

2πb 2π b

» » » »

The expression for ψ2 is obtained from Eq. (7.4)

Substituting Eq. (7.4) into Eqs. (7.3), (7.3), and (7.8)

by replacing y with ’y.

gives

Equation (7.4) can be used to solve for the evo-

2μA 1 2π x 2π

lution of the interface. The rate of change of the (P 1 ) y=0 = cos

» »

amplitude of the interface, ‚w/‚t must equal the b

»

vertical component of the ¬‚uid velocity at just 1

— +

above the interface. If not, a void would open up sinh » cosh 2π b

2π b

2π b »

in the system. This condition can be written as 2

»

1 2π b

— ’ tanh

‚w »

2πb 2π b

2π b

sinh » cosh »

= (v 1 ) y=0 (7.5)

‚t (7.9)

As an approximation, the velocity in Eq. (7.5) Substituting Eq. (7.9) into Eq. (7.7) gives

2

» »

4μA 1 2π 2π x 1 1 2π b

(ρ1 ’ ρ2 )gw = ’ + ’

cos tanh

» » »

sinh 2πb cosh 2πb sinh 2πb cosh 2π b

b 2π b 2π b

» » » »

(7.10)

is evaluated at the unperturbed position of the Solving Eq. (7.10) for A 1 and substituting the re-

interface (i.e. y = 0). The value of v 1 at y = 0 is sult into Eq. (7.6) gives

obtained by differentiating Eq. (7.4) with respect

‚w (ρ1 ’ ρ2 )gb

to x and evaluating the resulting expression at =w

‚t 4μ

y = 0. The result is

2π b ’1

»2

tanh 2π b ’ sinh 2π b cosh

‚w » » »

2π A 1 2π x — 2πb

= 2πb ’1

cos (7.6) »

+ sinh 2π b cosh

‚t » » » »

2πb

(7.11)

In order to constrain the value of the coef-

¬cient A 1 in Eq. (7.6), we must incorporate the Equation (7.11) describes an exponential growth

role of buoyancy on ¬‚uid motion. As the interface process with time scale „ given by

between the two ¬‚uids is disturbed downward

2π b ’1

»

(or upward) the less (more) dense ¬‚uid will dis- + sinh 2πb

cosh

4μ » »

„= 2πb

place the other. The pressure difference between ’1

(ρ1 ’ρ2 )gb 2

»

tanh 2πb ’ sinh 2πb

cosh 2π b

» » »

2πb

the two sides of the interface resulting from this

(7.12)

buoyancy force is given by

Equation (7.12) provides the necessary relation-

2(P 1 ) y=0 = ’(ρ1 ’ ρ2 )gw (7.7)

ship between the time scale of the instability and

The pressure on y = 0 in the upper layer can the wavelength ».

be found by using a force balance equation that Figure 7.3 plots the relationship between

equates the pressure gradient in the ¬‚ow with the dimensionless growth rate of a perturba-

the viscous stresses: tion (i.e. 1/„ ) and its dimensionless wavelength

‚P ‚ 2u ‚ 2u 2πb/». The fastest growing wavelength can be

’μ +2 =0 (7.8)

determined by differentiating Eq. (7.12) as a

‚x ‚ x2 ‚y

164 INSTABILITIES

ratio of meander wavelength and channel width

100

is a constant equal to about ten for these diverse

lmax= 2.568b

phenomena (Leopold et al., 1964). This similarity

suggests a common mechanism for meandering.

10

(r1 ’ r2)gbt

In this section we present a linear stability

analysis which predicts channel migration rates

4

and the proportionality between meander wave-

10 length and channel width using a very simple

geometric approach that applies to both river

channels and lava channels. Most models of allu-

vial channel meandering stress the importance of

10

102

100 101

10 the cross-sectional circulation or secondary ¬‚ow

2pb/l

in meander initiation (Callander, 1978; Rhoads

Fig 7.3 Plot of the inverse of the dimensionless time scale „ and Welford, 1991). This model does not replace

as a function of the dimensionless wavelength for the

those more sophisticated models, but rather es-

Rayleigh“Taylor instability (Eq. (7.12)).

tablishes suf¬cient conditions for meandering to

occur in a wide variety of channel ¬‚ows. Mean-

function of „ and setting the result equal to zero. ders of the correct scaling form when the bank

shear stress (in the case of sediment transport)

Turcotte and Schubert (2002) give the result of

or temperature (in the case of thermal erosion)

that equation as

is proportional to the curvature of the bank and

»max = 2.568b (7.13)

cross-sectional variations in shear stress or tem-

perature have a linear form. For alluvial chan-

where »max is the fastest-growing wavelength. Per-

nels, this curvature dependence arises from the

turbations of this wavelength have a growth time

centripetal force of ¬‚uid rounding a bend, and

scale given by

in thermal erosion it arises through the excess

13.04μ

„= latent heat of a curved bank. A linear stability

(7.14)

(ρ1 ’ ρ2 )gb analysis predicts that channels are unstable to

growth at all wavelengths larger than a critical

Equation (7.12) predicts that the fastest growing

wavelength given by » = 2π w, where w is chan-

wavelength is equal to about two and half times

nel width, and stable below it (provided that

the layer thickness. The Rayleigh--Taylor instabil-

some channel ¬‚ows exceed the critical condition

ity provides a classic example of a linear stability

required for erosion). The wavelength of fastest

analysis. We will use the linear stability analysis

growth predicted by this analysis is equal to

procedure again in Section 7.5, to explore the in- √

2 3πw or 10.88w, which is very similar to ob-

stability mechanism that creates oscillating arid

servations from alluvial and lava channels.

alluvial channels.

In alluvial rivers, the dependence of channel

migration on the local curvature of the bed has

been stressed in measurements of channel mi-

7.3 A simple model for river

gration by Nanson and Hickin (1983) and in the-

meandering oretical studies by Begin (1981) and Howard and

Knutson (1984). Nanson and Hickin (1983) found

Alluvial rivers, supraglacial meltwater streams, that channel migration rates were a maximum

within a narrow interval of R /w near 3, where

Gulf Stream meanders, and lava channels on

Earth, the Moon, and Venus (Komatsu and Baker, R is the radius of curvature of the meander,

1994) all exhibit meandering with similar propor- and w is the channel width. We will compare

tionality between meander wavelength and chan- our theoretical results to their data. Begin (1981)

nel width (Leopold et al., 1964; Stommel, 1965; showed, by computing the centripetal force nec-

Parker, 1975; Komatsu and Bakerm, 1994). The essary to accelerate ¬‚ow around a bend, that the

7.3 A SIMPLE MODEL FOR RIVER MEANDERING 165

v

t Fig 7.4 Schematic diagram of the

w2 meander instability model

y calculation.

w

0

l

shear stress at the bed is proportional to the cur- Eqs. (7.15) and (7.17) for the cross-sectional and

vature. Howard and Knutson (1984) developed a along-channel variation in shear stress gives

simulation model of meandering where the mi- 2

„ (x, y) = „b,s (1 + w K ) 1 ’ y + „c (7.18)

gration rates were proportional to local curvature w

up to some maximum value. They obtained real-

The shape of the curved-channel bank within

istic meandering channels with their model.

this linear-stability approach is given by yb =

In Figure 7.4 we present the geometry to be

sin(2π x/»). The curvature K is given by

considered in our model. We assume symmetry

d yb /dx 2 . Substituting the sinusoidal bank shape

2

within our 2D model, i.e. that the forces and ero-

into Eq. (7.18) gives

sion on one bank are equal in magnitude and

2

opposite in direction to those on the other bank. 2 2π

„ (x, y) = „b,s 1’ y 1+w

»

First we consider the distribution of boundary w

shear stress in a straight channel of width w. The 2π x

— sin + „c (7.19)

fundamental assumption of this analysis is that »

the cross-sectional shear stress pro¬le is linear,

We further assume that spatial variations in bed

decreasing from a maximum value at the bank

shear stress produce erosion and accretion of the

(y = w/2) to a minimum value at the channel

bank proportional to the variation in shear stress.

centerline:

In other words, the bank will erode in locations

2

„ (y) = „b,s 1 ’ y + „c where the variation in shear stress downriver is

(7.15)

w positive, and a point bar will prograde in places

where „b,s is the maximum shear stress at the where variations in shear stress are negative. The

bank for a straight channel and „c is the mini- rate of migration, therefore, is proportional to

the derivative of „ (x, y) with respect to x, evalu-

mum shear stress at the centerline. In a curving

ated at y = w/2:

channel, empirical data indicate that the bank

shear stress is equal to (Richardson, 2002) ‚„ d 2π x

= sin

‚x »

dt

w 2

y=w/2

„b,c = „b,s +1 (7.16)

2R 2

w 2π 2 2π x

= „b,s ’ sin

where „b,c refers to the bank shear stress for a » »

2 w

curved channel, and R is the radius of curvature (7.20)

of the channel centerline. For incipient mean-

The normalized growth rate of perturbations of

ders, the radius of curvature of the centerline

the bed is given by

is much larger than the channel width. As such,

we can use the approximation (1 + )2 ≈ 1 + 2 2

™ 2π„b,s 2 w 2π

= ’ (7.21)

to rewrite Eq. (7.16) as » »

w 2

w

„b,c ≈ „b,s 1 + = „b,s (1 + w K ) (7.17) which can be simpli¬ed to

R

2

™ 1 2π w

where K is the planform curvature of the cen-

∝ 1’ (7.22)

» »

terline (equal to the inverse of R ). Combining

166 INSTABILITIES

is dif¬cult, however, to compare migration rates

0.8

from many different bends without ¬rst quantify-

lmax = 10.88w

ing the effects of bank-material texture. Lumping

migration rate (m/yr)

0.6

data from meander bends with different erodi-

bilities would increase the scatter compared to

the Nanson and Hickin data, even if individual

0.4

bends closely follow a universal growth curve. Mi-

gration rate data from within individual bends

would help to resolve this question.

0.2

The analysis that led to Eq. (7.22) also ap-

plies to channels carved by thermal erosion.

0.0 Supraglacial meltwater streams (Parker, 1975)

40

30

10 20

0

l/w and lava channels on the Earth, the Moon, and

Venus (Hulme, 1973) form by thermal erosion. In

Fig 7.5 Plot of the normalized growth rate of perturbations the case of lava channels, melting of the channel

as a function of wavelength » (Eq. (7.22)). Also shown (dots) bed by hot lava (and solidi¬cation on the opposite

are the bank-migration data of Nanson and Hickin (1983)

bank) causes channel migration to occur. Curved

along the Beatton River for comparison.

banks have a melting point that depends on the

curvature of the bank. The melting temperature

of the bank is given by

Equation (7.22) is plotted in Figure 7.5 as the solid

curve. The maximum growth rate occurs where

T = Tm + c K (7.23)

the ¬rst derivative of ™/√ with respect to » is

zero. This occurs at » = 2 3πw = 10.88w, which

where T m is the melting point for a straight chan-

is close to the observed scaling relation (Leopold

nel, K is the planform curvature of the bank, and

et al., 1964).

c is a constant that depends on the latent heat of

For comparison, Figure 7.5 shows observed mi-

the bank material. The temperature required to

gration rate data for the Beatton River from Nan-

melt a curved bank is elevated because eroding

son and Hickin (1983). These authors employed

the bank by a given amount requires melting a

dendrochronological techniques to estimate the

larger volume of material per unit surface area

migration rates for sixteen bends of that river.

than for a straight channel. As such, the tempera-

They expressed their data as a function of the

ture boundary condition in lava channels has the

average radius of curvature of the bend. To com-

same curvature dependence as the shear stress

pare their results to our theoretical growth curve,

boundary condition in alluvial channels. In ad-

we used the relation » = π R which is applica-

dition, turbulent eddies will transport heat away

ble for an ideal shape of fully developed me-

from the hot center of the ¬‚ow towards the cool

anders: the sine-generated curve (Langbein and

sides of the ¬‚ow in the same way that turbu-

Leopold, 1966). Their results show the same de-

lent eddies transport shear stress through an al-

pendence with meander wavelength predicted

luvial channel, resulting in an approximately lin-

by our linear stability analysis. The migration

ear temperature pro¬le. These underlying simi-

rate rises rapidly for small »/w, reaches a maxi-

larities between alluvial and lava channels sug-

mum value, and then decreases. Migration rates

gests that the same geometrical instability mech-

decrease rapidly for tightly curved bends with

anism is at work in both cases.

»/w < 2π , consistent with our model™s prediction

that channel bends become stable below that

wavelength. It should be noted that some stud- 7.4 Werner™s model for eolian

ies following Nanson and Hickin (1983) showed

dunes

that migration rates are far less systematic that

the results obtained for the Beatton River. Most

of those studies have concluded that the Nanson Eolian ripples, dunes, and megadunes obey a

and Hickin results are not generally applicable. It striking periodicity. Many attempts have been

7.4 WERNER™S MODEL FOR EOLIAN DUNES 167

made to understanding the controls on eolian Folk, 1976). Wilson (1972) has proposed the same

bedform spacing in terms of both microscale pro- mechanism for ripples and megadunes with the

cesses (e.g. the trajectories of individual grains bedforms on these scales created by different

in saltation and reptation) and macroscale pro- scales of atmospheric circulation. However, care-

cesses (e.g. the interaction between a growing ful ¬eld studies have not found evidence for the

bedform, the wind ¬‚ow above it, and the result- existence of these persistent atmospheric circula-

ing pattern of erosion and deposition). Despite tions (Livingstone, 1986).

the efforts of many scientists, the mechanisms In 1995, Werner made a signi¬cant advance in

responsible for the formation of eolian ripples, eolian geomorphology in a paper that described

dunes, and megadunes and the factors responsi- a very simple but powerful model of eolian dune

ble for their size and spacing are not fully un- formation (Werner, 1995). Werner™s model begins

derstood. Perhaps the most well studied type of with a set of discrete sand slabs on a rectangular

bedform is the eolian ripple. Bagnold (1941) per- grid. Slabs may represent individual sand grains

formed the pioneering work on the formation of or a collection of grains; the number of slabs on

eolian ripples. He proposed that ripples form by the grid is limited for computational reasons, but

a geometrical instability in which a perturbation the model is not sensitive to the size of each slab.

exposes the windward side of the perturbation During each time step of the model, a sand slab is

to more impacts and faster surface creep than picked up from a randomly chosen pixel. If there

the leeward side, resulting in more grains being is no sand at that pixel, nothing happens. If a

ejected into saltation on the windward slope. In slab is present, it is picked up and moved down-

his theory, the spacing of ripples is related to a wind a distance of l grid point and deposited

characteristic saltation distance that is constant back on the surface with a probability that de-

in time. Sharp (1963), however, observed that in- pends on the presence or absence of sand at the

cipient ripples increase in spacing over time until new location. If sand is present at the new loca-

a steady-state spacing is achieved. In addition, ex- tion, the slab is deposited with probability ps . If

perimental studies of grain impacts suggest that no sand is present, it is deposited with a lower

a wide distribution of energies are imparted to probability pns . In nature, sand is more likely to

grains on the bed during each impact, resulting be deposited on a patch of sand than on a bare

in a wide distribution of saltation path lengths desert surface because the boundary layer above

(Mitha et al., 1986). As such, Bagnold™s hypothesis a sandy surface is usually rougher than above a

is not well supported by available data, and the bare, smooth surface, and because energy is ab-

controls on ripple spacing are still not well un- sorbed from the saltating grain by the granular

derstood. Anderson and Bunas (1992) developed bed. If the slab is not deposited, it is transported

a model that produced realistic ripples based on another distance l and again deposited with a

probability given by ps or pns . One additional rule

the microscale processes of grain movement and

controls the probability of deposition: if a sand

bed impact, but the spacing of ripples in their

slab lands in a ˜˜shadow zone,™™ it is deposited with

model was controlled by an ad hoc ˜˜ceiling™™ in

100% probability. A shadow zone is de¬ned as the

the model domain that does not occur in nature.

domain downwind of any topographic high that

The mechanisms controlling the formation

lies below a plane with an angle of 15—¦ to the

and spacing of eolian dunes are also incom-

horizontal. To map the shadow zone, a path is

pletely understood. Eolian dunes represent a dis-

traced from every local ridge down a 15—¦ slope

tinct bedform type from eolian ripples since

in the downwind direction until the plane in-

they occur at a larger scale with no transitional

tersects the surface. The 15—¦ angle of the shadow

bedforms. It has been proposed that dunes re-

zone can be varied, but its value should be signif-

sult from the presence of pre-existing wave-like

icantly less than the angle of repose. Physically,

motions or secondary circulations in the atmo-

the shadow zone represents the zone of recircu-

spheric boundary layer. Periodic motions intrin-

lation on the lee side of a growing dune. In this

sic to the ¬‚ow or resulting from the response

zone, air¬‚ow velocities are greatly reduced and

of the ¬‚ow to a topographic disturbance upwind

any sand in saltation will likely be deposited.

may produce periodic bedforms (Wilson, 1972;

168 INSTABILITIES

wind Fig 7.6 A schematic diagram of

l Werner™s model for eolian dune

shadow

formation and evolution. After

zone

Werner (1995).

angle of repose

Finally, when the sand slab is deposited it may the example of Figure 7.7) and the direction of

create an oversteepened condition in which the transport (e.g. star dunes form if the sand bed

local slope is higher than the angle of repose. If is initially thin and the wind direction is pe-

this occurs, the model moves the slab down the riodically reversed). Werner™s model provides a

direction of steepest descent until a site is found nice illustration of the power of simpli¬ed mod-

that does not lead to an oversteepened condition. els for understanding the behavior of complex

Figure 7.6 illustrates the model schematically. A systems with positive feedback. Werner™s model

code that implements Werner™s model is provided can easily be modi¬ed to account for variable

in Appendix 6. bed topography. Figure 7.8, for example, illus-

Initially, a rough surface develops in Werner™s trates a dune ¬eld migrating through a crater. On

model from an initially ¬‚at surface due to the Mars, the crestlines of dunes are deformed down-

random entrainment of grains from the surface wind of craters, similar to the effect shown in

(i.e. sand from some areas will be entrained more Figure 7.8.

than in others because of the stochastic nature of The height and spacing of dunes in Werner™s

entrainment in the model). Incipient ridges then model increases continuously over time with-

develop small shadow zones that encourage de- out ever reaching a steady-state value. Numeri-

position on their lee sides. This deposition leads cal models that incorporate the interaction be-

to the development of taller dunes with longer tween dunes and the air¬‚ow above them, in

shadow zones in a positive feedback. Using dif- contrast, predict that dunes eventually achieve

ferent parameter values, Werner (1995) was able an optimum dune height and spacing. In na-

to reproduce a variety of dune types by vary- ture, it is unclear which model is correct. In lo-

ing the initial sand thickness (a relatively thin cations where dunes and megadunes are both

initial sand bed leads to barchan dunes, as in present, these two bedform types appear to

wind

100 — L2 timesteps 200 — L2 timesteps 300 — L2 timesteps

Fig 7.7 Results of Werner™s model for a grid of size

L = 500 and model parameters l = 5, pns = 0.4, ps = 0.6, occur at distinctly different spatial scales with

angle of repose = 30—¦ , and angle of shadow zone = 15—¦ for no intermediate bedforms (Wilson, 1972). In

three model durations. The initial condition is uniform sand other words, even in sand seas mature enough

with a thickness of 10 sand parcels on each pixel.

to have formed megadunes, dunes appear to

7.5 OSCILLATIONS IN ARID ALLUVIAL CHANNELS 169

the American West. In their model, channels ag-

wind grade and widen until a threshold slope and/or

width is reached, initiating incision in the over-

steepened reach and aggradation in the down-

stream reach. In this section, we explore a math-

ematical model for Schumm and Hadley™s arid cy-

cle of erosion and compare the model predictions

to the spatial pattern of oscillations commonly

observed in alluvial channels of the southwest-

ern US.

Oscillating alluvial channels in southern Ari-

zona are classi¬ed as one of three basic types

depending on their geomorphic position. Pied-

mont channels fed by montane drainage basins

are continuous channels that alternate between

braided and narrowly entrenched reaches with

wavelengths on the order of 1 km. Wild Burro

Wash on the Tortolita Piedmont is a classic ex-

ample of this type of oscillating channel (Figure

Fig 7.8 Results of Werner™s model modi¬ed to include

variable bed topography. In this example, a dune ¬eld interacts 7.9a). Channel width (measured from 1 m reso-

with a crater. The bed topography causes deformation of lution US Geological Survey Digital Orthopho-

dune crest lines within and downwind of the crater.

toquadrangles (DOQQs)) and bed elevation (mea-

sured from a high-resolution DEM with average

channel slope subtracted) in Wild Burro Wash are

reach a steady-state spacing of approximately 10--

plotted in Figure 7.9a. These data illustrate that

200 m, with stronger winds and larger grain sizes

oscillations in both bed elevation and channel

associated with greater spacing. This observation

width occur in Wild Burro and that the two os-

suggests that the aerodynamic interaction be-

cillations are in phase.

tween dunes and the air¬‚ow above them is essen-

The second type of oscillating channel is

tial for understanding the height and spacing of

the piedmont discontinuous ephemeral stream.

mature eolian dunes.

These channels are fed only by local piedmont

runoff, and they alternate between incised chan-

nels and sheet-¬‚ow-dominated channel fans with-

7.5 Oscillations in arid alluvial

out well-de¬ned banks (Bull, 1997). The transi-

channels tions between channel fans and incised reaches

occur abruptly at one or more steep ˜˜headcuts™™

Most ¬‚uvial-geomorphic studies that have mod- (Figure 7.9b). Dead Mesquite Wash in southeast

eled the evolution of channel longitudinal pro- Tucson is the classic example given by Bull (1997)

¬les have assumed channels that are uniform (Figure 7.9b). Packard (1974) documented in-phase

in width or a prescribed function of drainage oscillations in bed elevation and width in Dead

area (e.g. Begin et al., 1981; Slingerland and Mesquite, similar to the behavior of Wild Burro.

Snow, 1987; Sinha and Parker, 1996; Dade and The third type of oscillating channel is the

Friend, 1998). Actual channels often vary substan- classic southwestern valley-¬‚oor arroyo. These

tially in width, however, and these variations in systems alternate between narrow, deeply incised

width can strongly in¬‚uence channel evolution. channels and broad, shallow channel fans, but

Schumm and Hadley (1957), for example, devel- they typically drain one or more broad basins

oped their conceptual ˜˜arid cycle of erosion™™ on and adjacent ranges and have wavelengths on the

the basis of the episodic cut-and-¬ll histories and order of 10 km. Abrupt headcuts also character-

oscillating geometries observed in the arroyos of ize distributary-to-tributary transitions in these

170 INSTABILITIES

(a) Wild Burro DEM

(c) Vamori

downstream

300 m

B

downstream

photo

I

N B I

B I B I

>100 cm

1988 flood depths

50“10

30“50

10“30

0“10

w/h long./width

h

120/ w profile B

0.3

60/

0

I

0/

1500 x (m)

2500 0

1000 500

2000

(b) Dead Mesquite headcuts headcuts

incised

B

channel fans

300 m

downstream

I

N

5 km B

N

B I B I B I

Fig 7.9 Examples of oscillating channels in southern

two ways. First, a numerical model for the

Arizona. Alternating reaches: B“braided and I“incised.

coupled evolution of the channel longitudinal

(a) Wild Burro Wash data, including (top to bottom) a

high-resolution DEM (digital elevation model) in shaded relief pro¬le and cross section is presented. This model

(PAG, 2000), a color Digital Orthophotoquadrangle (DOQQ), predicts that arid alluvial channels are un-

a 1:200-scale ¬‚ow map corresponding to an extreme ¬‚ood on

stable over a range of wavelengths controlled

July 27, 1988 (House et al., 1991), and a plot of channel width

by channel width, depth, and slope. Second,

w (thick line) and bed elevation h (thin line, extracted from

a database of channel geometries in southern

DEM and with average slope removed). (b) Dead Mesquite

Arizona is constructed that con¬rms the model

Wash shown in a color DOQQ and a detailed map of channel

predictions and places the three channel types

planform geometry (after Packard, 1974). (c) Vamori Wash

within a single continuum from steep, small-

shown in an oblique perspective of a false-color Landsat

image (vegetation [band 4] in red) draped over a DEM. wavelength-oscillating channels to gently slop-

Channel locations in Figure 7.11a. Modi¬ed from Pelletier ing, long-wavelength-oscillating channels.

and Delong (2005). For color version, see plate section. The equation describing conservation of mass

for sediment states that erosion and deposition

systems. In contrast to the other oscillating- in an alluvial-channel bed are proportional to the

channel types, not all arroyos are periodically gradient of total sediment discharge, or

alternating; some are tectonically controlled by

‚h 1 ‚(wqs )

=’

the elevations of basin divides. In southern Ari- (7.24)

‚t C 0 ‚x

zona, Cooke and Reeves (1976) identi¬ed Vamori

Wash as perhaps the best example of an oscil- where h is the elevation of the channel bed, t

lating arroyo system (Figure 7.9c). Entrenched is time, C 0 is the volumetric concentration of

reaches are associated with low vegetation den- bed sediment, w is the channel width, qs is the

sity, whereas channel fans support a broad, dense speci¬c sediment discharge (i.e. the sediment dis-

riparian zone in Vamori Wash. charge per unit channel width), and x is the dis-

The common oscillatory pattern in these tance downstream. In most bedload transport re-

channel types suggests a common underlying lations, speci¬c sediment ¬‚ux is proportional to

mechanism that acts across multiple spatial the 3/2-power of the bed shear stress. This im-

scales. This section explores this hypothesis in plies that speci¬c sediment discharge is a linear

7.5 OSCILLATIONS IN ARID ALLUVIAL CHANNELS 171

function of channel gradient and a nonlinear are assumed to widen as they aggrade and nar-

function of the discharge per unit channel width, row as they incise. Bull™s model states that chan-

or nels incise if the stream power is greater than a

threshold value and aggrade if the stream power

b

‚h

Q

qs = ’B (7.25) is less than that value. Expressing this relation-

‚x

w

ship in terms of channel widening (aggradation)

where B is a mobility parameter related to grain and narrowing (incision) gives

size and Q is water discharge. The value of b is

‚w wqs ’ w 0 qs0

constrained by sediment rating curves and is be- =’ (7.29)

‚t h0 ’ h

tween 2 and 3 for both suspended-load and bed-

load transport. Here we consider the case b = 2. where qs0 is the equilibrium sediment ¬‚ux

Equation (7.25) assumes that the bed shear stress or stream power, and h 0 is the equilibrium

is much larger than the threshold for particle bank height. Equation (7.29) represents a cross-

entrainment. sectional mass balance: sediment removed from

If the channel width is assumed to be uniform the bank contributes to the local sediment-¬‚ux

along the longitudinal pro¬le, the combination de¬cit wqs ’ w 0 qs0 (i.e. the amount of sediment

of Eqs. (7.24) and (7.25) gives the classic diffusion that cannot be transported out of the reach), pro-

equation: moting further aggradation in the reach.

‚h ‚ 2h The nonlinear term in Eq. (7.28) alters the dy-

=κ 2 (7.26)

namics of channels markedly compared to the

‚t ‚x

diffusive behavior expressed in Eq. (7.26). Along a

where the diffusivity is given by κ = (B Q 2 )/

reach with uniform discharge and grain size, the

(C 0 w 0 ), and w 0 is the uniform channel width.

diffusion equation smoothes out curvatures in

If the channel width is not uniform along

the pro¬le over time. The nonlinear term in Eq.

the longitudinal pro¬le, the chain rule must be

(7.28), however, has the opposite effect through

used when differentiating qs in Eq. (7.24). This ap-

a positive feedback between channel width and

proach introduces an additional nonlinear term

slope. In this feedback, spatial variations in chan-

into Eq. (7.26), to give

nel slope generate variations in width via Eq.

‚h 1 ‚ 2h 1 ‚w ‚h (7.29). Large gradients in channel width, in turn,

= κw 0 ’2 (7.27)

‚t w ‚x w ‚x ‚x

2

increase the nonlinear behavior in Eq. (7.28), fur-

ther localizing erosion and deposition to com-

It is mathematically easier to de¬ne h as the

plete the feedback cycle. This cycle is balanced

difference between the local bed elevation and

by the diffusive term, and the balance between

that of a straight, equilibrium channel with uni-

these two terms controls the oscillation wave-

form discharge. Equation (7.27) then becomes

length.

‚h 1 ‚ 2h 1 ‚w ‚h

= κw 0 ’2 S0 ’ Equations (7.28) and (7.29) may be solved by us-

(7.28)

‚t w ‚x w ‚x ‚x

2

ing linear stability analysis and direct numerical

where S0 is the equilibrium channel slope. Equa- solution. Linear stability analysis works by solv-

tion (7.28) can be used to study the evolution of ing the linear approximation to Eqs. (7.28) and

perturbations from the equilibrium geometry. (7.29) for the growth rate of a small-amplitude

Channel widening and narrowing occurs by a oscillation superimposed on the initial channel

complex set of processes including bank retreat geometry (a channel with speci¬ed equilibrium

and bed scouring. Scouring often leads to chan- width w 0 , bank height h 0 , and slope S0 ).

nel narrowing as ¬‚ow is focused into the scour Linear-stability analysis involves solving for

zone and parts of the former channel bed effec- the growth rate of spatially periodic, small-

tively become part of the bank. Bull™s ˜˜threshold- amplitude perturbations from an equilibrium ge-

of-critical-power™™ concept (Bull, 1979) can be used ometry by retaining only linear terms in the

to relate the rate of channel widening and nar- model equations. In this linear stability analysis,

rowing with the excess stream power if channels we assume that the channel has an average width

172 INSTABILITIES

1

(b)

(a) 1/2

l = p(w0 h0/S0)

h

w 0.8

0.6

d 0.4

x

0.2

unstable stable

0

5

0

(c) 3 4 1/2

1 2 6

l ((w0 h0/S0) )

1

incision

0

h backfilling

random

l

h0

instability growth phase propagation phase

downstream downstream

3

w 2

w0

1

0

2 x (km) 3 4 5 2 x (km) 3 4 5

0 1 1

retaining only linear terms gives

Fig 7.10 Model behavior. (a) Schematic diagram of model

geometry. (b) Growth curve for linear-stability analysis,

) = 2π κ 2π ) + S0 w ei ( 2π x )

2π x 2π x

δh 0 ei ( h 0 ei (

’

where δ is the nondimensional growth rate. This analysis » » »

» » w0 0

predicts unstable behavior for small wavelengths, with a

(7.31)

maximum value at » = π (w 0 h 0 /S0 )1/2 . (c) Numerical

solution. Initial phase (left) is characterized by growth of

and

small, random perturbations and an increase in oscillation

)=’ κ

wavelength with time (line thickness increases with time). 2π

2π x +φ+ π 2π x + π

δw 0 ei ( h 0 ei ( ) (7.32)

» »

2 2

Steady-state phase (right) is characterized by solitary-wave »

w0h0

propagation of oscillations in upstream direction. Incision and

Equation (7.31) implies φ = ’π/2 (i.e. h and w

channel narrowing on leading edge of each wave are balanced

are 90—¦ out of phase). Substituting Eq. (7.32) into

by back¬lling and widening on trailing edge. Modi¬ed from

Pelletier and Delong (2005). Eq. (7.31) gives the following equation for growth

rate:

w 0 and slope S0 . The variables h and w refer to 2

κ κ S0 1

2π

δ= 1’ (7.33)

the variations of h and w from their average val- » h0w0 δ

2

w0

ues, and have solutions of the form

The solution to Eq. (7.33) is

2π x

) eδt , w = w ei ( 2π x +φ ) eδt

h = h 0 ei ( (7.30)

» »

0

2π κ π π

where δ is the growth rate and » is the wave- S0

2

δ= ’ ’ (7.34)

w0» » » w0h0

length. These solutions assume a spatially peri-

odic form, and an exponential increase in am-

plitude with time as the instability grows. Sub- Equation (7.34) is plotted in Figure 7.10b.

stituting Eqs. (7.30) into Eqs. (7.28) and (7.29) and This ¬gure shows that the growth rate of

7.5 OSCILLATIONS IN ARID ALLUVIAL CHANNELS 173

perturbations has a maximum at the ampli¬cation of spatial variations in width

within a range of wavelengths close to 1 km. In

w0h0 the latter stage, oscillations achieve a steady-state

»=π (7.35)

S0 amplitude and propagate upstream as a train of

solitary (i.e., nondispersive) topographic waves. At

This result can also be obtained by differentiating

any given instant, the channel geometry is char-

Eq. (7.34) with respect to » and setting the result