. 6
( 14)


max h
0 50 x (km) 100
3 0.1
mean h
0.3 calibrated with
Stock et al. (2004)
mean E
mean U
0 0.0
10 30
20 0
10 30
t (Myr) t (Myr)
Fig 4.6 Plots of erosion rate versus time for the best-¬t rates before and after knickpoint passage. The in-
sediment-¬‚ux-driven model (a) and (b) and stream-power
set graph in Figure 4.5d plots the erosion rate of
model (c) and (d). (a) Plot of erosion rate versus time for
the upland plateau along a linear transect shown
three locations along the Kings and South Fork Kings River
in Figure 4.5b. The hillslope erosion rate E h was
(locations in Figure 4.5), illustrating the passage of two
prescribed to be 0.01 mm/yr based on measured
knickpoints corresponding to the two uplift pulses. The
cosmogenic erosion rates (Small et al., 1997; Stock
late-Cenozoic pulse of incision at location 2 provides a
et al., 2004, 2005). The plot in Figure 4.6a shows
forward-model calibration based upon the cosmogenically
a range of erosion rates on the upland surface
derived incision rates (Stock et al., 2004, 2005) at this
location. The maximum value of the erosion rate at this from a minimum of 0.01 mm/yr (on hillslopes) to
location is matched to the observed rate of 0.3 mm/yr by a maximum value of 0.03 (in channels). Erosion
varying the value of K s . Inset into (a) is the upland erosion
rates in upland plateau channels are controlled
rate along the transect located in Figure 4.5b, illustrating
by both E h and K s . The fact that the model pre-
minimum erosion rates of 0.01 mm/yr and maximum values of
dicts an average erosion rate of ≈ 0.02 mm/yr,
0.03 mm/yr, consistent with basin-averaged cosmogenic
consistent with basin-averaged rates measured
erosion rates of ≈ 0.02 mm/yr. (b) Maximum and mean
cosmogenically (Riebe et al., 2000, 2001), provides
elevations of the model, illustrating the importance of
additional con¬dence in this approach.
isostatically-driven uplift in driving peak uplift. (c) and (d)
Figure 4.6b plots maximum and mean eleva-
Plots of erosion rate and elevation versus time for the
stream-power model, analogous to (a) and (b). Modi¬ed from tion values as a function of time in the model.
Pelletier et al. (2007c). Reproduced with permission of Mean surface elevation increases during active
Elsevier Limited.
uplift but otherwise decreases. Maximum eleva-
tion continually increases through time to a ¬-
River was preceded and followed by much lower nal value of over 4 km as a result of isostatic up-
rates of incision (≈ 0.02--0.05 mm/yr) from 5 to lift of Boreal Plateau remnants driven by canyon
3 Ma and 1.5 Ma to present. The model predicts cutting downstream. Mean erosion and uplift
the same order-of-magnitude decrease in incision rates are also plotted in Figure 4.5e, showing


h (km)
(c) Chagoopa

Kern R.
x (km) 20
Late Cretaceous knickpoint (c)

h (km)
Late Cenozoic knickpoint

1.0 n = 1
x (km)
0 40 80
Late Cretaceous knickpoint
h (km)

Late Cenozoic knickpoint

h (km) 1.0 n = 2
2.0 3.0 4.0
0.5 1.0
x (km)
0 40 80
Fig 4.7 Best-¬t sediment-¬‚ux-driven model results for the
this model occurs as two pulses: a 1-km pulse
North Fork Kern River basin (i.e. 1 km of late Cretaceous
in the late Eocene (t = 35 Ma) and a 0.5-km pulse
uplift and 0.5 km of late-Cenozoic uplift), illustrating broadly
in the late Miocene (t = 7 Ma). Figure 4.6c illus-
similar features to those of the actual North Fork Kern River
trates the knickpoints as they pass points along
basin illustrated in Chapter 1. Modi¬ed from Pelletier et al.
the Kings and South Fork Kings River. In this
(2007c). Reproduced with permission of Elsevier Limited.
model, knickpoint migration occurs at a simi-
lar rate in the ¬rst and second pulses of up-
lift, as illustrated by the equal durations between
a gradual increase in mean erosion rate over a
knickpoint passage following the ¬rst and second
40 Myr period to a maximum value of approxi-
pulses. The rates of vertical incision are lower in
mately 0.05 mm/yr. Isostasy replaces 80% of that
the second phase because the knickpoint has a
erosion as rock uplift, as prescribed by Eq. (4.27).
gentler grade. Model results are broadly compa-
Figure 4.7 illustrates the ability of the model to
rable to those of the sediment-¬‚ux-driven model,
predict details of the Sierra Nevada topography
except that the propagation of the initial knick-
using the North Fork Kern River as an example.
point requires only about half as long as the
The model predicts the elevations and extents of
sediment-¬‚ux-driven model to reach its present
the Chagoopa and Boreal Plateaux as well as the
approximate locations and shapes of the two ma-
The results of this model application sup-
jor knickpoints along the North Fork Kern River.
port two possible uplift histories corresponding
The best-¬t results of the stream-power model
to the best-¬t results of the sediment-¬‚ux-driven
are shown in Figures 4.5f--4.5h and 4.6c--4.6d,
with n = 1 and K w = 8 — 10’5 kyr’1 . Uplift in and stream-power models. Only one of these

histories is consistent with the conclusion of high of bedrock channel and coupled bedrock-alluvial
(> 2.2 km) Sierra Nevada in early Eocene time channel modeling using the stream power and
(Poage and Chamberlain, 2002; Mulch et al., 2006), sediment-¬‚ux-driven erosion models.
however. That constraint, taken together with the Analyses of sediment-load data have shown
fact that sediment abrasion is the dominant ero- that erosion rates are approximately proportional
sional process in massive granitic rocks (Whipple to the mean local relief or mean elevation of a
et al., 2000) support the greater applicability of drainage basin. In tectonically inactive areas, for
the sediment-¬‚ux-driven model in this case. example, Pinet and Souriau (1988) obtained
The results of this section are consistent with
E av = 0.61 — 10’7 H (4.28)
the basic geomorphic interpretation of Clark et al.
(2005) that the upland plateaux and associated
where E av is the mean erosion rate in m/yr
river knickpoints of the southern Sierra Nevada
and H is the mean drainage basin elevation
likely record two episodes of range-wide surface
in m. Similar studies have documented approx-
uplift totaling approximately 1.5 km. The results
imately linear correlations between mean ero-
presented here, however, suggest that the ini-
sion rates and mean local relief, relief ratio, and
tial 1-km surface uplift phase occurred in late
basin slope both within and between mountain
Cretaceous time, not late Cenozoic time. In the
belts (e.g. Ruxton and MacDougall, 1967; Ahnert,
sediment-¬‚ux-driven model, the slow geomorphic
1970; Summer¬eld and Hulton, 1994; Ludwig and
response to the initial uplift phase is caused by
Probst, 1998). Equation (4.28) can be expressed as
the lack of cutting tools supplied by the slowly-
a differential equation for mountain-belt topog-
eroding upland Boreal Plateau. If the behavior
raphy following the cessation of tectonic uplift:
of the sediment-¬‚ux-driven model is correct, the
model results suggest that 32 Myr (i.e. the time ‚H 1
=’ H (4.29)
since onset of Sierra Nevadan uplift, according to ‚t „d
Clark et al. (2005)) does not afford enough time
with „d = 16 Myr using the correlation coef¬cient
to propagate the upland knickpoint to its present
location. The self-consistency of the model re- in Eq. (4.28). Equation (4.29) has the solution
H = H0 e ’t/„d , where H0 is the mean basin ele-
sults provide con¬dence in this interpretation.
The model correctly reproduces details of the vation immediately following uplift. Isostatic re-
bound will increase „d by a factor of ρc /(ρm ’ ρc ),
modern topography of the range, including the
where ρm and ρc are the densities of the man-
elevations and extents of the Chagoopa and Bo-
tle and crust, thereby giving „d ≈ 50--70 Myr. This
real Plateaux and the elevations and shapes of
value is about ¬ve times smaller than the age
the major river knickpoints.
of Paleozoic orogens, several of which (e.g. Ap-
palachians, Urals) still stand to well over 1 km in
peak elevation. This is the paradox of persistent
4.5 The erosional decay of
mountain belts.
ancient orogens One possible reason for persistent mountain
belts is the role that piedmonts play in raising
One of the great paradoxes in geomorphology the effective base level for erosion. As a moun-
concerns the ˜˜persistence™™ of ancient orogens. tain belt ages, sediment can be deposited at the
Analyses of modern sediment-load data imply footslope of the mountain, causing the base level
that mountain belts should erode to nearly sea of the bedrock portion of the mountain belt to
level over time scales of tens of millions of years. increase (Baldwin et al., 2003; Pelletier, 2004c).
Several mountain belts that last experienced ac- To explore this effect we need to develop mod-
tive tectonic uplift in the Paleozoic (e.g. Ap- els for the coupled evolution of bedrock chan-
palachians, Urals) still have mean elevations at nels and alluvial piedmonts. Here we consider
or near 1 km in elevation. In this section we ex- a model that couples the stream power model
plore this question using 1D and 2D modeling (Eq. (4.4) to a simple model for alluvial-channel

large-k case: k = 4 — 10’3 km2/yr
K = 10’7/yr

h/h0 h/h0
t = 50 Myr
t = 100 Myr
t = 100
t = 150
h(L) = 0
hb ha

small-k case: k = 2 — 10’3 km2/yr
K = 10’7/yr
(c) K = 10’7/yr
K = 3 — 10’8/yr
k = 10’3 km2/yr
h/h0 t = 50 Myr
k = 3 — 10’3 km2/yr h/h0
k = 10’2 km2/yr t = 100
t = 150


condition is applied a small distance L h from the
Fig 4.8 (a) Model geometry and key variables. (b) Plots of
divide (i.e. the hillslope length) to give h b (L h , 0) =
elevation vs. downstream distance for two values of κ. Each
plot shows a temporal sequence from t = 50 Myr“250 Myr in h 0 . Sea level is assumed to be constant, giving
50 Myr intervals. (c) Plots of mean basin elevation vs. time for h a (L , t) = 0. The two remaining boundary condi-
L m = 100 km, L = 500 km, and a range of values of K and κ.
tions are continuity of elevation and sediment
From Pelletier (2004c).
¬‚ux at the mountain front, x = L m :

‚h b ‚h a
evolution based on the diffusion equation (Paola h b (L m ) = h a (L m ), dx = κ
‚t ‚x
et al., 1992): Lh x=L m
‚h ‚ 2h
=κ 2 (4.30)
‚t ‚x Equations (4.12) and (4.30) were solved using
where κ is the piedmont diffusivity. The diffu- upwind differencing and FTCS techniques for the
sion equation is a highly simpli¬ed model of pied- bedrock and alluvial portions of the basin, re-
mont evolution, but enables us to make prelimi- spectively. Codes for modeling coupled bedrock-
nary conclusions regarding the role of piedmont alluvial channel evolution are given in Appendix
aggradation on the time scale of mountain-belt 3. The effect of the piedmont diffusivity value on
denudation. mountain-belt denudation is illustrated in Figure
The model geometry and boundary condi- 4.8b. In these two experiments, the model param-
eters are identical (K = 10’7 yr’1 , L m = 200 km,
tions are illustrated in Figure 4.8a. h a and h b refer
L = 500 km) except for the values of κ, which dif-
to the alluvial and bedrock portions of the pro-
fer by a factor of 2. A smaller value of κ results in
¬le. Tectonic uplift is assumed to occur as a rigid
block between x = 0 and x = L m prior to t = 0. a steeper piedmont, lower bedrock relief, and a
h 0 is the maximum elevation immediately follow- smaller denudation rate. After 250 Myr, the small-
κ case has a mean elevation about twice as large
ing the cessation of active uplift. This boundary

as the large-κ case. Varying the piedmont length bedrock and alluvial components, and a transcen-
L for a ¬xed value of κ has a similar effect. dental equation is obtained for „d :
Figure 4.8c illustrates the decline in mean
basin elevation with time on a semi-log plot for Ar /(K „d )
h b (x, t) = h 0
a piedmont of ¬xed length (L = 500 km) and a (4.32)
range of values of K and κ. Following an ini-
L ’x
Ar /(K „d ) sin √
tial rise to a maximum value, the decay in mean κ„d
h a (x, t) = h 0 (4.33)
basin elevation is exponential in each case, but Lm L ’L m
sin √

varies as a function of both the bedrock and pied-
κ„b L ’ Lm
1’ = tan
mont parameters, illustrating their dual control √ (4.34)
K „d κ„b
on the tempo of denudation.
In nature, values of K , κ, L m , and L vary from
one mountain belt to another. Persistent moun- Equation (4.34) cannot be further reduced
tain belts can be expected to have extreme values and must be solved numerically for speci¬c val-
of one or more of these parameters, resulting in ues of the model parameters. As an example,
large values of „d compared to the global average. for the model parameters of the Appalachian-
type case of Figure 4.9a, Eq. (4.34) gives „b = 186
L and L m can be determined from topographic
and geologic maps. Values of K and κ are not well Myr. In the limit of no piedmont (i.e. L m ’ L ),
constrained, but qualitative data on climate, rock Eq. (4.34) has the appropriate limiting behavior
„b ’ Ar /K .
type, and sediment texture can be used to vary
those parameters around representative values. Shepard (1985) found actual bedrock pro¬les
The effects of unusually resistant bedrock, for in nature to be best ¬t by a power law when slope
example, are illustrated in the Appalachian-type was plotted versus downstream distance, consis-
example of Figure 4.9a. Hack (1957, 1979) argued tent with Eq. (4.32). The power-law exponent in
that the resistant quartzite of the Appalachian the model is a function of the bedrock erodibil-
highlands was a key controlling factor in its evo- ity (more resistant bedrock leads to more con-
lution. In this example, the values of L m and L cave pro¬les) but the exponent is also a func-
were set to 100 km and 500 km. The values K = tion of the basin diffusivity, κ (through „d ), il-
5 — 10’8 yr’1 and κ = 6 — 10’4 km2 /yr led to pro- lustrating the explicit coupling of the bedrock
¬les most similar to those observed (Figure 4.9a). and piedmont pro¬le shapes. This result suggests
In contrast, Figure 4.8b illustrates the Ural-type that piedmont deposition must be included to
example (i.e. an unusually broad piedmont). Here, correctly model bedrock drainage basin evolu-
L m and L were set to 600 km and 2000 km. The tion. The sine function is a poor approximation
value of K most consistent with observed pro¬les to many observed piedmont pro¬les, but this dis-
in the Urals (K = 10’7 yr’1 , κ = 2 — 10’2 ; Figure crepancy is largely the result of neglecting down-
4.9b) is a factor of two greater than the value for stream ¬ning. A more precise piedmont form
could be achieved by making κ a function of
the Appalachians. The Ural example shows that
a broad piedmont can diminish bedrock relief, downstream distance.
reducing the average basin denudation rate com- Montgomery and Brandon (2002) recently
pared to a narrower piedmont of the same slope. questioned the linear relationship between mean
The eastern Appalachian piedmont is relatively erosion rate and mean elevation (or mean lo-
modest in size, so resistant bedrock is the most cal relief) that lies at the heart of the persis-
likely explanation for its persistence and strongly tent mountain belt paradox. Their study is best
concave headwater pro¬les. The values of K and known for documenting a rapid increase in ero-
κ in Figure 4.9 are not unique; a range of values sion rates in terrain with a mean local relief
is consistent with the observed pro¬les. greater than 1 km. However, these authors also
The model equations can also be solved an- expanded upon earlier data sets that established
alytically for the decay phase of mountain-belt a linear relationship between mean erosion rate
evolution using separation of variables. Power- and mean local relief. Montgomery and Brandon
law and sinusoidal solutions are obtained for the (2002) found that the expanded data set was best

(b) 1.0
82 80
83° W
h model
36° N
(a) 250 Myr

K = 5 — 10’8 yr’1
k = 6 — 10’4 km2/yr
0.5 250 Myr

50 km 400
0 200
x (km)
48° E 60
54 0.7
60° N
250 Myr


K = 10’7 yr’1
k = 2 — 10’2 km2/yr
52 model
0.5 250 Myr

200 km
0 1000
x (km)
represented by a power law:
Fig 4.9 (a) and (b) Appalachian-type example. (b)
Longitudinal pro¬les of channels in the Santee and Savannah
E av = 1.4 — 10’6 R z
drainage basins (location map in (a)), plotted with a model
pro¬le (K = 5 — 10’8 yr’1 , κ = 6 — 10’4 km2 /yr) at where E av is in units of mm/yr and R z is
t = 250 Myr for comparison. Evolution of the model pro¬le the mean local relief (de¬ned as the average
in intervals of 50 Myr is also shown. (c) and (d) Urals-type difference between the maximum and minimum
example. (d) Pro¬les in the Volga drainage basin (location map
elevation over a 10 km radius) in m. This re-
in (c)), plotted with a model pro¬le (K = 10’7 yr’1 ,
sult necessitates a reassessment of the persistent
κ = 2 — 10’2 km2 /yr) for comparison. From Pelletier
mountain-belt paradox.

(a) (b) Eav ∝ H

Eav ∝ H 2
10’2 1.3
Eav∝ Rz
(mm/yr) SP

10’3 SFD
Eav∝ Rz SFD

10’1 10’3 10’2 10’1 100
Rz (km) H (km)
uplift of 0.05 mm/yr for the ¬rst 100 Myr of the
Fig 4.10 (a) Plots of mean erosion rate vs. mean local relief
simulation (suf¬cient to develop steady state), fol-
following uplift for the stream-power and
sediment-¬‚ux-driven models. The stream-power model lowed by a long period of erosional decay and
isostatic response to erosion. The values of K =
closely approximates a power-law relationship with an
2 — 10’3 kyr’1 and K s = 4 — 10’2 (m kyr)’1/2 used
exponent of 1.3 (plots closely overlap), while the
sediment-¬‚ux-driven model also follows a power law but with
in these runs predict identical steady state condi-
an exponent of 1.8. (b) Plots of mean erosion rate versus
tions at the cessation of tectonic uplift at 100 Myr.
mean elevation following uplift for the stream-power and
Figure 4.10a plots the relationship between mean
sediment-¬‚ux-driven models. Following a transient phase,
erosion rate E av and mean local relief for the two
erosion rates in the stream-power model are proportional to
models following tectonic uplift. The sediment-
mean elevation H while in the sediment-¬‚ux-driven model
¬‚ux-driven model closely follows the power-law
they are proportional to the square of H .
trend documented by Montgomery and Brandon
(2002) with an exponent of 1.8 (model results and
power-law trend strongly overlap). The stream-
In order to determine the relationships be-
power model also follows a power-law trend but
tween erosion rate, mean local relief, mean el-
with a smaller exponent of 1.3. Figure 4.10b plots
evation, and speci¬c models of bedrock channel
the relationships between mean erosion rate and
erosion, we must work with a model designed
mean elevation for the stream-power (SP) and
to use the stream-power model and a sediment-
sediment-¬‚ux-driven (SFD) models. Following an
¬‚ux-driven model interchangeably. The model we
initial transient phase, mean erosion rates in the
will consider is similar to that of Figure 4.4 but
stream-power model are proportional to mean el-
it assumes that each pixel in the model is above
evation while in the sediment-¬‚ux-driven model
the threshold for channelization. As such, it does
they are proportional to the square of mean ele-
not include hillslope processes. The stream-power
vation. The fact that the trends between mean
version of the model uses Eq. (4.4) while the
erosion rate and mean local relief differ from
sediment-¬‚ux-driven version uses Eq. (4.25). The
those between mean erosion rate and mean eleva-
values of m and n used in the model are 0.5 and
tion indicates that declines in mountain-belt re-
1.0, respectively (Kirby and Whipple, 2001). Iso-
lief and mean elevation are not identical. For the
static rebound in the model is estimated using
purposes of modeling erosional decay it is mean
Eq. (4.26). Sediment ¬‚ux Q s in Eq. (4.25) is com-
elevation, not relief, that is most signi¬cant
puted in the model by downslope routing of the
because mean erosion rate is the derivative of
local erosion rate computed during the previous
mean elevation with respect to time, while ero-
time step, multiplied by the pixel area.
sion rates and relief are not as simply related.
Model runs reported here use a square do-
Nevertheless, observed relationships between
main of width 128 km subject to uniform vertical

100 max, sediment-flux-
h driven (SFD)
(km) 50 km

mean (H), SFD

max, stream-power
mean (H), SP
0 600
t (Myr)
(b) (c)
10 mean relief (Rz )
mean relief (Rz )
h max
(km) max
mean (H)
mean (H)
H ∝ 1/t
H ∝ ed’t/t
td = 30 Myr
10 stream-power model (SP) sediment-flux-driven model (SFD)
200 300
0 800
100 100 t (Myr)
t (Myr)
Fig 4.11 (a) Semi-log plots of maximum and mean elevation where c is a constant. Equation (4.36) has the so-
versus time following the cessation of tectonic uplift for the lution
stream-power and sediment-¬‚ux-driven model of a uniformly
uplifted square domain of width 128 km. The bedrock
H = ct + (4.37)
erodibilities of each model have been scaled to predict the H0
identical steady-state landscapes. The sediment-¬‚ux-driven
where H0 is the initial mean elevation following
model exhibits persistent mountain belts relative to the
tectonic uplift.
exponential decay of the stream-power model. (b) Semi-log
The inverse power-law dependence in Eq.
plots of maximum and mean elevations and mean local relief
versus time following uplift for the stream-power model (4.37) predicts a fundamentally different his-
(same model output as in (a)). Also plotted is an exponential tory of mountain-belt decay than the exponen-
decay model with a time constant „d = 30 Myr. (c) Log“log tial model that follows from Eq. (4.29). The in-
plot of maximum and mean elevations and mean local relief
verse power law has a broad tail in which to-
versus time for the sediment-¬‚ux-driven model, illustrating
pography decays much more slowly than pre-
the asymptotic approach to 1/t behavior at large time scales.
dicted by the exponential model at long time
scales. Figure 4.11 compares the topographic de-
mean erosion rates and mean relief are still cay of the stream-power and sediment-¬‚ux-driven
very useful for distinguishing between different models. Figure 4.11a plots the mean and max-
models of topographic evolution. imum elevations and mean local relief of each
Data plotted in Figure 4.10b suggest that topo- model on a semi-log scale as a function of time
graphic decay in the sediment-¬‚ux-driven model following tectonic uplift. For large times, the
is described at large times by mean elevation in the stream-power model de-
cays as e’t/„d with „d ≈ 30 Myr. The sediment-¬‚ux-
= ’c H 2 (4.36) driven model, however, decays much more slowly

of 1 cm and a diffusivity of κ = 0.1 cm2 /yr, model
with time. One can think of the inverse power-
the relative concentration of an abrupt pulse
law model as an exponential model in which the
effective time scale (i.e. „d in Eq. (4.29)) is con- of radionuclides into a semi-in¬nite soil pro¬le
with time. Plot the concentration pro¬les for t =
tinually increasing as the topography decays. As
10, 100, and 1000 yr following the fallout event.
such, the inverse-power-law model suggests that
4.2 Assume that the precipitation in a mountain range
there is no single time scale for mountain belt
increases linearly with elevation. In such a case,
decay, but rather a spectrum of time scales with
the stream-power model for a semi-circular basin
increasing values for older mountain belts. The in a humid environment (Eq. (4.12)) becomes
maximum and mean elevations in the stream-
‚h ‚h
power and sediment-¬‚ux-driven models are plot- h
= cx 1 + (E4.1)
‚t ‚x
ted in Figures 4.11b and 4.11c on semi-log and h0
log--log plots, respectively, with the exponential
where h 0 is a characteristic length scale. Using the
and power-law asymptotic trends of each model
method of characteristics, plot the evolution of
also shown.
knickpoints in this system. How do the knickpoints
Erosion in natural bedrock channels likely
change in shape and speed as they propagate up
falls between the two end-member models of
into the headwaters?
stream-power and sediment-¬‚ux-driven erosion,
4.3 Extract a ¬‚uvial channel pro¬le from a topographic
but, as previously noted, ¬eld studies suggest that
map or a DEM starting at the channel head. Using
the sediment-¬‚ux-driven model is more appro- Excel, evolve the channel forward in time accord-
priate for massive bedrock lithologies common ing to Eq. (4.12). Neglect uplift. Choose a value of
in ancient mountain belts. The results described c that produces realistic erosion rates (e.g. 1 m/kyr
in this section illustrate that the stream-power for very steep terrain).
model in its basic form and a sediment-¬‚ux- 4.4 Construct a simple model of bank retreat assum-
ing that the rate of bank retreat is inversely pro-
driven model predict power-law relationships be-
portional to the channel width w:
tween mean erosion rate and mean local re-
lief with exponents of 1.3 and 1.8, respectively.
‚h c ‚h
The results of the sediment-¬‚ux-driven model are (E4.2)
‚t w ‚x
more consistent with the analysis of Montgomery
and Brandon (2002). The sediment ¬‚ux-driven with c = 0.001 yr’1 . Assume a channel with an ini-
model predicts that mountain-belt topography tial width of 100 m. As the banks retreat and the
decays proportionally to 1/t, resulting in persis- channel widens, update the value of w in Eq. (E4.2).
tent mountain belts relative to the exponential- How long is required before the channel doubles
in width?
decay model. Conceptually, mountain belts per-
4.5 Cliff retreat can be modeled with the advec-
sist in the sediment-¬‚ux-driven model because
tion equation. Consider a cliff and talus slope
the cutting tools responsible for abrading chan-
as illustrated in Figure 4.12. Initially, the cliff
nel beds decrease over time, thereby decreasing
channel incision rates in a positive feedback.
This mechanism provides a new hypothesis for
mountain belt persistence in massive bedrock
lithologies dominated by the saltation-abrasion

4.1 The migration of radionuclides into a soil by diffu-
sive processes was considered in Chapter 2. Modify
the model to include advective transport. This ra-
tio of the diffusivity to the advection velocity is
Fig 4.12 Schematic diagram of Exercise 4.5.
often called the dispersivity. Given a dispersivity

and talus slope both have height H /2. Assume Assume that the rock and talus slope have equal
that the cliff retreat can be modeled advec- densities.
tively with retreat rate c = 1 m/kyr. Model the 4.6 Repeat the above for a case in which the channel
evolution of the cliff and talus slope through crosses a major lithologic boundary (e.g. the Kaibab
time. Impose conservation of mass at the base limestone at the rim of the Grand Canyon). Choose
of the cliff (i.e. the volume of rock removed a smaller value of c to represent the more resistant
from the slope must be deposited on the slope). unit.
119° W 118° W
along-dip profile
(b) 1°
profile in (b)
0 50
x (km)
San Joaquin R.
profile in (b)
37° N

inset in (c)
Kings R.

Kaweah R.

along-strike profile
36° N
Kern R.

100 150
0 50
x (km)
3.0 Boreal
profile in (f)
2.0 Chagoopa
Kern R.
x (km)
0 30
3.0 knickpoints
linear 2.0
Kern R.
profile in (e)
x (km)
0 100
S (m/m)
h (km)
0.0 0.5 >1.0
3.0 4.0
1.0 2.0

Plate 1.8 Major geomorphic features of the southern Sierra Nevada. (a) Shaded relief map of topography indicating major rivers
and locations of transects plotted in b. (b) Maximum extents of the Chagoopa and Boreal Plateaux based on elevation ranges of
1750“2250 m and 2250“3500 m a.s.l. Also shown are along-strike and along-dip topographic transects illustrating the three levels
of the range in the along-strike pro¬le (i.e. incised gorges, Chagoopa and Boreal Plateaux) and the westward tilt of the Boreal
Plateau in the along-strike transect. (c) and (d) Grayscale map of topography (c) and slope (d) of the North Fork Kern River,
illustrating the plateau surfaces (e) and their associated river knickpoints (f). Modi¬ed from Pelletier (2007c). Reproduced with
permission of Elsevier Limited.
116° W
117° W Plate 2.18 (a) Location map and
Yucca LANDSAT image of Eagle Mountain
Eagle piedmont and adjacent Franklin
Nevada Lake Playa, southern Amargosa
Valley, California. Predominant wind
N direction is SSE, as shown by the
e Amargosa
Ca vad Valley
wind-rose diagram (adapted from
Funeral li a
36.5° N Mtns. fo
January 2003“January 2005 data
Eagle from Western Regional Climate
Mtn. Center, 2005). Calm winds are
Death de¬ned to be those less than 3 m/s.
Valley Qa2
(b) Soil-geomorphic map and
oblique aerial perspective of Eagle
Mtns. Qa5“Qa7 Mountain piedmont, looking
36.0° N N 10 km Soil-geomorphic map and predominant playa
southeast. Terrace map units are
N pits
sample-pit locations wind direction
based on the regional classi¬cation
by Whitney et al. (2004).
outline of active
wind (c)
Approximate ages: Qa2 “ middle
modern playa
Pleistocene, Qa3 “ middle to late
3“6 m/s
Pleistocene, Qa4 “ late Pleistocene,
6“9 m/s

Qa5“Qa7 “ latest Pleistocene to
Franklin Lake
active. (c) Map of eolian silt
1% Playa
thickness on Qa3 (middle to late
Pleistocene) surface, showing
3% maximum thicknesses of 80 cm
close to the playa source,
Amargosa Eagle
River decreasing by approximately a
factor of 2 for each 1 km
secondary hotspot
downwind. Far from the playa,
background values of approximately
1 km 20 cm were observed. Modi¬ed
0 10 20 40 80 cm
Silt thickness on Qa3
from Pelletier and Cook (2005).
(a) 80 (b)
u = 5 m/s
K = 5 m2/s
p = 0.05 m/s
(x', y')
u = 5 m/s
K = 10 m2/s source
p = 0.075 m/s
wind x
2 3 4
0 direction
distance downwind from playa (km)
2 km
0.1 0.2 0.4 1.0

wind direction

numerical model results
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s
wind direction
(d) (e)
numerical model results
numerical model results
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s
1 km
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s

q = 10
q = ’10 q = 30
10 20 40
0 80 cm

Plate 2.19 (a) Plot of eolian silt thickness versus downwind distance, with analytic solutions for the two-dimensional model for
representative values of the model parameters. (b) Schematic diagram of model geometry. Depositional topography shown in this
example is an inclined plane located downwind of source (but model can accept any downwind topography). (c) Color maps of
three-dimensional model results, illustrating the role of variable downwind topography. In each case, model parameters are
u = 5 m/s, K = 5 m2 /s, and p = 0.05 m/s. Downwind topography is, from left to right, a ¬‚at plane, an inclined plane, and a
triangular ridge. Width and depth of model domain are both 6 km. (d) Color maps of three-dimensional model results for
deposition downwind of Franklin Lake Playa, illustrating the role of variable wind direction. From left to right, wind direction is
θ = ’10—¦ (0—¦ is due south), 10—¦ , and 30—¦ . (e) Map of three-dimensional model results obtained by integrating model results over a
range of wind conditions, weighted by the wind-rose data in Figure 2.18a. Modi¬ed from Pelletier and Cook (2005).
(a) Plate 3.1 Comparison of steepest-descent and MFD
¬‚ow-routing algorithms. (a) Steepest descent: a unit of
precipitation that falls on a grid point (shown at top) is
successively routed to the lowest of the eight nearest
neighbors (including diagonals) until the outlet is reached.
Precipitation can be dropped and routed in the landscape in
steepest descent MFD method
any order. MFD: all incoming ¬‚ow to a grid point is
117.05° W 117.00 116.95 116.90
(b) distributed between the down-slope pixels, weighted by bed
2 km slope. To implement this algorithm ef¬ciently, grid points
36.225° N
should be ranked from highest to lowest, and routing should
be done in that order to ensure that all incoming ¬‚ow from
up-slope is accumulated before downstream routing is
calculated. (b) Map of contributing area calculated with
steepest descent for Hanaupah Canyon, Death Valley,
California, using USGS 30 m DEMs. Grayscale is logarithmic
and follows the legend at ¬gure bottom. (c) Map of
steepest descent
contributing area computed with bifurcation routing, for the
(c) same area and grayscale as (b). Multiple-¬‚ow-direction
routing results in substantially different and more realistic
¬‚ow distribution, particularly for hillslopes and areas of
distributary ¬‚ow. Modi¬ed from Pelletier (2004d).


106 107 m2
100 104

scour/mixing depth tephra concentration

1% 100%
1.0 m
0.1 0.2
(a) (b)


Plate 3.14 Digital grids output by the modeling of tephra outlet
redistribution following a potential volcanic eruption at Yucca concentration
southwesterly = 1.97%
Mountain. (a) Map of scour-depth grid, calculated from
contributing area grid. (b) Grayscale map of tephra
concentration in channels, calculated using the
scour-dilution-mixing model. Modi¬ed from Pelletier et al.
N 3 km
(2008). Reproduced with permission of Elsevier Limited.
(a) (b)

(c) (d)

Plate 3.7 Illustration of a successive ¬‚ow-routing algorithm applied to Fortymile Wash alluvial fan in Amargosa Valley,
Nevada. (a) Shaded-relief image of a high resolution (1 m/pixel) photogrammetric DEM of Fortymile Wash, (b)“(d) ¬‚ow depths
predicted by the model for different values of the prescribed upstream discharge: (b) 300 m3 /s, (c) 1000 m3 /s, and (d) 3000 m3 /s.
(a) (b)


Plate 3.11 Field photos illustrating the three-dimensional pattern of contaminant dilution near Lathrop Wells tephra sheet.
(a) and (b) Channel pits on the tephra sheets expose a ¬‚uvially mixed scour zone ranging from 12 to 29 cm in thickness. Two types
of deposits occur beneath the scour zone: the tephra sheet itself (exposed on the upper and lower sheet) and debris-¬‚ow
deposits comprised predominantly by Miocene volcanic tuff and eolian silt and sand transported from the upper tephra sheet.
(c) The effects of dilution visible as high-concentration channels draining the tephra sheet (channel at right, 76% tephra) join with
low-concentration channels (at left, 26% tephra). Tephra concentration in these channels correlates with the darkness of the
sediments, with the dark-colored channel at right joining with the (larger) light-colored channel at left to form a light-colored
channel downstream. From Pelletier et al. (2008). Reproduced with permission of Elsevier Limited.
116.52° W 116.50° 10 ’2 10’1
116.48° 101 km
N1 km


tephra sheet
36.70° slope*area1/2


Lathrop Wells cone
relative tephra
0.01 0.1 1.0

topography and source mask
scour/mixing depth 0.1 0.3 1.0
900 1000 m

Plate 3.12 Model prediction for tephra concentration and scour/mixing depth downstream from the tephra sheet of the
Lathrop Wells volcanic center. (a) Model inputs include a 10-m resolution US Geological Survey Digital Elevation Model (DEM) of
the region and a grid of source and background concentrations. Input concentration values are: 73% (Lathrop Wells tephra sheet),
9% (distal source region), and 3% (background). (b) Map of stream power in the vicinity of the source region. Scale is logarithmic,
ranging from 10’2 to 101 km. (c) Map of scour/mixing depth (scale is quadratic, ranging from 10 cm to 1 m). (d) Map of tephra
concentration (scale is quadratic, ranging from 0.001 to 1). From Pelletier et al. (2008). Reproduced with permission of Elsevier
tephra thickness slope

10 ’3 10 ’2 10 ’1 100 km
0.01 0.1 1m 0.01 0.03 0.1 0.3 1.0

(a) (c)



N 3 km
slope > 0.3 (17°) > 0.05 km
Plate 3.13 Digital grids used in the modeling of tephra redistribution following a volcanic eruption at Yucca Mountain.
(a) Shaded relief image of DEM and Fortymile Wash drainage basin (darker area). Tephra from only this drainage basin is
redistributed to the RMEI location. (b) Grayscale image of DEM slopes within the Fortymile Wash drainage basin.
(c) Black-and-white grid of areas in the drainage basin with slopes greater than 17—¦ . (d) Black-and-white grid of active channels in
the drainage basin (de¬ned as pixels with contributing areas greater than 0.05 km2 ). All tephra deposited within the black areas of
(c) and (d) are assumed to be mobilized by mass movement, intense rilling, or channel ¬‚ow. Modi¬ed from Pelletier et al. (2008).
Reproduced with permission of Elsevier Limited.
(b) (c)

Plate 3.15 Maps of tephra concentration in channel sediments following a hypothetical volcanic eruption at Yucca Mountain
with (a) southerly, (b) westerly, and (c) northerly winds. Tephra concentration at the basin outlet varies from a maximum of 0.73%
(for southerly winds) to a minimum value of 0.076% (for northerly winds) for this eruption scenario. Southerly winds result in
higher concentrations because the high relief of the topography north of the repository is capable of mobilizing more tephra than
other wind-direction scenarios. Modi¬ed from Pelletier et al. (2008). Reproduced with permission of Elsevier Limited.
Plate 4.5 Maps of best-¬t
model results for the
sediment-¬‚ux-driven model
(b)“(d) and stream-power
model (f)“(h) starting from
the low-relief, low-elevation
surface illustrated in (a). The
actual modern topography
is shown in (e). The best-¬t
uplift history for the
sediment-¬‚ux-driven model
occurs for a 1-km pulse of
uplift starting at 60 Ma
(t = 0) and a 0.5 pulse
starting at 10 Ma (t = 50 yr
out of a total 60 Myr). The
best-¬t uplift history for the
stream-power model
occurs for a 1-km pulse of
uplift at 35 Ma and a 0.5 km
pulse at 7 Ma. Also shown
are transect locations and
stream location identi¬ers
corresponding to the
erosion-rate data in
Figure 4.6. Modi¬ed from
Pelletier et al. (2007c).
Reproduced with
permission of Elsevier
103° W
103° W 117° W
117° W 49° N
49° N
(a) (b)
300 km

42° N

uniform erosion above glacial-limit erosion concentrated near glacial-limit
109° W
109° W
0 compensation 1
Plate 5.11 Map of compensation of glacial erosion assuming (a) uniform erosion within glaciated areas and (b) erosion
concentrated within 500 m altitude of glacial limit. Modi¬ed from Pelletier (2004a).
(a) Wild Burro
(c) Vamori
300 m
>100 cm
1988 flood depths
w/h long./width
120/ w profile B
2500 0
1000 500
0.3 2000 1500 x (m)

(b) Dead Mesquite headcuts
channel fans
300 m

5 km B

Plate 7.9 Examples of oscillating channels in southern Arizona. Alternating reaches: B“braided and I“incised. (a) Wild Burro
Wash data, including (top to bottom) a high-resolution DEM (digital elevation model) in shaded relief (PAG, 2000), a color Digital
Orthophotoquadrangle (DOQQ), a 1:200-scale ¬‚ow map corresponding to an extreme ¬‚ood on July 27, 1988 (House et al., 1991),
and a plot of channel width w (thick line) and bed elevation h (thin line, extracted from DEM and with average slope removed).
(b) Dead Mesquite Wash shown in a color DOQQ and a detailed map of channel planform geometry (after Packard, 1974). (c) Vamori
Wash shown in an oblique perspective of a false-color Landsat image (vegetation [band 4] in red) draped over a DEM. Channel
locations in Figure 7.11a. Modi¬ed from Pelletier and Delong (2005).
40 — vert.
l = 0.2

middle expansion

late porewater


Plate 7.14 2D model evolution, with 40— vertical exaggeration. Color plots of compaction rate C shown at (a) early,
(b) middle, and (c) late times in the model. In the early phase, alternating zones of matrix expansion and contraction develop near
the upper boundary as zones of initially higher porosity expand, capture upwardly migrating porewater, and further expand in a
positive feedback. In the second stage, zones of expansion ascend and drive converging ¬‚ow of the matrix into the high-porosity
drumlin core. In the third stage, porewater is squeezed from the matrix until the sediment is fully compacted.
(a) (c) Plate 7.17 Grayscale
maps for average drumlin
43.4° N
width and bedrock depth in
43.6° N
(a) and (b) north-central
New York and (c) and (d)
Wisconsin, east of Madison.
43.4 (a) and (c) Maps of drumlin
width constructed by
averaging widths within a
2-km square moving
window. (b) and (d) Maps of
bedrock depth also
constructed with a 2-km
square moving window and
drumlin width using USGS groundwater
well data. Curves are drawn
76.8 89.2° W
77.2° W 88.8
to highlight areas in which
0 1000 m 0 500 m thick sediment and wide
drumlins (including Rogen
(b) (d)
moraine) coincide.

River Valley
60 m
0 60 m
areas of thick till and
= wide drumlins
(b) (c)

(d) (e)

Plate 7.20 Numerical model results and north-polar topography. (a) Shaded-relief image of Martian north-polar ice cap DEM
constructed using MOLA topography. The large-scale closeup indicates examples of gullwing-shaped troughs, bifurcations, and
terminations. Highest elevations are red and lowest elevations are green. (b)“(e) Shaded-relief images of the model topography,
’h, for (b) t = 10, (c) t = 100, (d) t = 1000, and (e) t = 2000 starting from random initial conditions. The model parameters are
L = 250 (number of pixels in each direction), x = 0.4, T0 = 0.3, „i = 0.05, „ f = 1. In the Barkley approximations to Eq. (7.52),
T0 is combined into two parameters, a = 0.75 and b = 0.01. The model evolution is characterized by spiral merging and
alignment in the equator-facing direction. A steady state is eventually reached with uniformly rotating spirals oriented clockwise or
counter-clockwise depending on the initial conditions. Modi¬ed from Pelletier (2004b).
Chapter 5

Flexural isostasy

Isostasy refers to the buoyant force created
5.1 Introduction when crustal rock displaces mantle rock. Isostatic
balance requires that, over geologic time scales,
the weight of the overlying rock must be uni-
Flexural isostasy is the de¬‚ection of Earth™s litho-
form at any given depth in the mantle. Because
sphere in response to topographic loading and
the mantle acts as a ¬‚uid over long time scales,
unloading. When a topographic load is generated
any lateral pressure gradient will initiate man-
by motion along a thrust fault, for example, the
tle ¬‚ows to restore equilibrium. Isostasy simply
lithosphere subsides beneath the load. The width
says that the hydrostatic pressure gradient in
of this zone of subsidence varies from place to
the mantle must be zero, otherwise the mantle
place depending on the thickness of the litho-
would correct the imbalance by ¬‚owing. Isostatic
sphere, but it is generally within the range of 100
balance requires that the hydrostatic force pro-
to 300 km. Conversely, a reduction in topographic
duced by the topographic load, ρc gh (where ρc is
load causes the lithosphere to rebound, driving
the density of the crust, g is the acceleration due
rock uplift. Flexural-isostatic uplift in response to
to gravity, and h is the elevation of the mountain
erosion replaces approximately 80% of the eroded
belt) be equal to the buoyancy force (ρm ’ ρc )gw
rock mass, thereby lengthening the time scale
(where ρm is the density of the mantle and w is
of mountain-belt denudation by a factor of ap-
the depth of the crustal root) produced by the
proximately ¬ve because erosion must remove
displacement of low-density crustal rocks with
all of the rock that makes up the topographic
higher-density mantle rocks (Figure 5.1):
load and the crustal root beneath it in order
to erode the mountain down to sea level. Given
ρc gh = (ρm ’ ρc )gw (5.1)
the ubiquity of erosion in mountain belts, it is
reasonable to assume that ¬‚exural-isostasy plays Solving for w, the thickness of the crustal root,
a key role in nearly all examples of large-scale gives:
landform evolution. Flexural isostasy also plays
an important role in the evolution of ice sheets (5.2)
ρm ’ ρ c
because the topographic load of the ice sheet
For typical crust and mantle densities, e.g. ρc =
causes lithospheric subsidence, thereby in¬‚uenc-
2.7 g/cm3 and ρm = 3.3 g/cm3 , the depth of the
ing rates of accumulation and ablation on the
crustal root is approximately ¬ve times the
ice sheet. In this chapter, we will discuss three
height of the topographic load.
broadly-applicable methods (series and integral
Flexure refers to the forces and displacements
solutions, Fourier ¬ltering, and the Alternating-
involved in bending the elastic lithosphere. Flex-
Direction Implicit (ADI) method) for solving the
ure affects isostasy because small-scale variations
¬‚exural-isostatic equation in geomorphic applica-
in the topographic load (e.g. peaks and valleys)




0.6 1.0
0.2 0.4 0.8


Fig 5.1 Isostatic balance requires that a topographic load,
given by ρc g h, be balanced by a much larger crustal root,
which extends to a depth given by Eq. (5.2).

can be supported by the rigidity of the litho- 2 4 6
sphere. It is only at length scales larger than the
¬‚exural wavelength (i.e. ≈ 100--300 km) that iso- Fig 5.2 (a) Nondimensional de¬‚ection of an elastic beam of
static balance is achieved. The equations used to rigidity D and length L subject to an applied force V at the
describe ¬‚exure were originally developed in the end of the beam. (b) Nondimensional de¬‚ection of an elastic
mechanical engineering literature to describe the beam overlying an inviscid ¬‚uid of different density, subject to
a line load at x = 0. The de¬‚ection of the elastic beam in this
response of elastic beams and plates to applied
case has a characteristic length scale given by ±.
forces. The de¬‚ection of a diving board under
the weight of a diver is a simple example of 2D
¬‚exure. If a diver stands still at the end of the
of a diver at the end of a diving board (x = 0)
board, all of the forces and torques on the diving
of length L can be obtained by integrating four
board must be in balance, otherwise the board
times to obtain the cubic function (Figure 5.2a):
would accelerate. A force-balance analysis of the
diving board indicates that the fourth derivative
V (L ’ x)3
of the displacement w is proportional to the ap- w(x) = (5.4)
plied load q(x) (Turcotte and Schubert, 2002):
The ¬‚exural rigidity of the lithosphere is, in turn,
= q(x)
D (5.3) controlled by the elastic thickness T e as well as
dx 4
the elastic properties of rock by the relationship
where D is a coef¬cient of proportionality that
de¬nes the ¬‚exural rigidity of the board. The so- E T e3
D= (5.5)
lution to Eq. (5.3) corresponding to the weight V 12(1 ’ ν 2 )

where E is the elastic modulus and ν is the Pois- approximated as acting at a single point along
son ratio of the rocks in the crust. Equation (5.5) a 1D pro¬le. Within a 1D model, such a load is
indicates that the ¬‚exural rigidity is very sensi- referred to as a ˜˜line™™ load. A narrow mountain
tive to the elastic thickness T e . One common tech- range that extends for thousands of kilometers
nique for mapping T e in a region involves compar- along-strike is one example of a model well ap-
ing the topography and gravity ¬elds at different proximated by line loading. In such cases, the
length scales. At small length scales, topography term q(x) on the right side of Eq. (5.6) is equal
to zero except at x = 0 where it is equal to some
and gravity are poorly correlated, while at large
scales they have a high degree of coherence. Lo- prescribed value V . In order to solve Eq. (5.6) un-
cal values of T e can be inferred by mapping the der a line load it is easiest to set q(x) equal to
smallest length scales at which topography and zero for all x and then introduce the load V as
a boundary condition for x = 0. With q(x) = 0,
gravity are strongly correlated.
In this chapter our primary focus is on ver- Eq. (5.6) becomes
tical (i.e. topographic) loads. In some cases, such
d4 w
+ (ρm ’ ρc )gw = 0
D (5.7)
as when a subducting slab exerts a horizontal
dx 4
frictional force on the overriding plate, horizon-
The general solution to Eq. (5.7) is obtained by
tal forces play an important role in lithospheric
¬‚exure. Such cases are rare, however, in geo-
x x
w(x) = ex/± c 1 cos + c 2 sin
morphic applications. We will also ignore the
± ±
time-dependent response of the mantle to load- x x
+ e’x/± c 3 cos + c 4 sin (5.8)
ing, which can be important for some problems ± ±
(e.g. ¬‚exural isostatic response to Quaternary ice
where c 1 through c 4 are integration constants,
sheets). Once a load is applied, there is necessar-
ily some time delay before the crust can respond
by uplift or subsidence. This time scale is gener- ±= (5.9)
(ρm ’ ρc )g
ally on the order of 104 yr, and generally longer


. 6
( 14)