. 9
( 14)


as a function of their size. This curve can be used
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


N = cA

N (>A)

(d) (e) 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

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
in Figures 6.27a--6.27c. The accumulation pattern
was assumed to be uniform for elevations greater
than 1 km.
The ice-surface topography, ¬‚ow, bedrock de-

¬‚ection, and bedrock topography of the model
after ten iterations are given in Figures 6.27d--
6.27g. The bedrock topography does not include
400 the de¬‚ection under the ice load. Therefore, it
represents what the landscape would look like
if the ice load was removed and only the effect
of glacial erosion remained. In this experiment,
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
the ice sheet the original rough topography re-
mains because erosion is minimal there. Near the
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
in¬ll de¬‚ects the lithosphere, further depressing
the ice surface and focusing ice drainage into the
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


ice-surface topography
iteration 1 ice-surface topography

(b) (e)

basal sliding velocity
basal sliding velocity iteration 10
iteration 1

ice margin

lighting direction
lithospheric deflection
lithospheric deflection iteration 10
iteration 1


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
velocities, (f) bedrock de¬‚ection, and (g) bedrock topography
(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
50 km.
where D is the ¬‚exural rigidity given by
E T e3
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, ρ=
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.

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
150 km
50 km 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).

(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
high in zones of diverging ice.

iteration 10 ice-surface topography


iteration 10 basal sliding velocity


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

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

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
covers topography only above the ELA.
6.5 The Boussinesq equation
‚h kρg ‚ ‚h
= h (E6.1)
‚t μφ ‚ x ‚x
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


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

y = ’b

r1, m

w x

r2, m

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
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
in Figure 7.2. A ¬‚uid with a thickness b and den- (7.2)
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

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

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

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
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 » » » »
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
4μ » »
„= 2πb
place the other. The pressure difference between ’1
(ρ1 ’ρ2 )gb 2
tanh 2πb ’ sinh 2πb
cosh 2π b
» » »
the two sides of the interface resulting from this
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

ratio of meander wavelength and channel width
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.
(r1 ’ r2)gbt

In this section we present a linear stability
analysis which predicts channel migration rates

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
100 101
10 the cross-sectional circulation or secondary ¬‚ow
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
„= latent heat of a curved bank. A linear stability
(ρ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

t Fig 7.4 Schematic diagram of the
w2 meander instability model
y calculation.


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
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
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.
In other words, the bank will erode in locations
„ (y) = „b,s 1 ’ y + „c where the variation in shear stress downriver is
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 »
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
„b,c ≈ „b,s 1 + = „b,s (1 + w K ) (7.17) which can be simpli¬ed to
™ 1 2π w
where K is the planform curvature of the cen-
∝ 1’ (7.22)
» »
terline (equal to the inverse of R ). Combining

is dif¬cult, however, to compare migration rates
from many different bends without ¬rst quantify-
lmax = 10.88w
ing the effects of bank-material texture. Lumping
migration rate (m/yr)

data from meander bends with different erodi-
bilities would increase the scatter compared to
the Nanson and Hickin data, even if individual
bends closely follow a universal growth curve. Mi-
gration rate data from within individual bends
would help to resolve this question.
The analysis that led to Eq. (7.22) also ap-
plies to channels carved by thermal erosion.
0.0 Supraglacial meltwater streams (Parker, 1975)
10 20
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
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

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;

wind Fig 7.6 A schematic diagram of
l Werner™s model for eolian dune
formation and evolution. After
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


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

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

(a) Wild Burro DEM
(c) Vamori
300 m
>100 cm
1988 flood depths
w/h long./width
120/ w profile B
1500 x (m)
2500 0
1000 500
(b) Dead Mesquite headcuts headcuts
channel fans
300 m

5 km B

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

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
qs = ’B (7.25) is less than that value. Expressing this relation-
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
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
‚h 1 ‚ 2h 1 ‚w ‚h
= κw 0 ’2 S0 ’ Equations (7.28) and (7.29) may be solved by us-
‚t w ‚x w ‚x ‚x
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

(a) 1/2
l = p(w0 h0/S0)
w 0.8

d 0.4

unstable stable

(c) 3 4 1/2
1 2 6
l ((w0 h0/S0) )
h backfilling
instability growth phase propagation phase
downstream downstream

w 2

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
maximum value at » = π (w 0 h 0 /S0 )1/2 . (c) Numerical
solution. Initial phase (left) is characterized by growth of
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 »
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
w 0 and slope S0 . The variables h and w refer to 2
κ κ S0 1

δ= 1’ (7.33)
the variations of h and w from their average val- » h0w0 δ
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)
» »

2π κ π π
where δ is the growth rate and » is the wave- S0
δ= ’ ’ (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

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


. 9
( 14)