. 5
( 14)


classic dilution-mixing model scour-dilution-mixing model
0 100
contaminant concentration
estimated scour depth in each channel pixel. As
Fig 3.9 Schematic diagram of (a) classic dilution-mixing
model, and (b) scour-dilution-mixing model. (a) The basic a simple example, assume that the ¬‚ood-envelope
operation of the classic model is illustrated in the lower left. curve shows the peak unit discharge for extreme
Here, lines with arrows represent upstream contributing
¬‚oods to be linearly proportional to contributing
area, a proxy for sediment yield. Concentrations downstream
area. Then, if the scour depth is observed to be
of a con¬‚uence are equal to the average of the upstream
1 m (e.g. using a scour chain emplaced prior to
concentrations weighted by contributing area. In a drainage
a recent extreme ¬‚ood) at a channel pixel with
basin, concentrations decrease rapidly with distance
a contributing area of 100 km2 , then the scour
downstream from a localized contaminant source, with
depth at any channel pixel in the basin will be
abrupt decreases at major con¬‚uences (illustrated by
equal to the square root of the normalized con-
grayscale pattern). (b) In the scour-dilution-mixing model, the
vertical and cross-channel distributions of contaminant are tributing area at that pixel (i.e. the contributing
also considered. At each pixel (i , j ), the total volumes of area divided by 100 km2 ) multiplied by 1 m.
upstream contaminant and uncontaminated channel-bed
The model routes the contaminant (expressed
sediments are mixed with the local volumes of contaminant
as an equivalent thickness for each pixel) down-
and uncontaminated sediments. The contaminant distribution
stream using the same ¬‚ow-routing algorithm
at that pixel is de¬ned by that concentration value and the
used for routing contributing area. Finally, as the
local scour depth. The resulting mixture is transported
contaminant is routed downstream, the model
downstream where the same operation is repeated.
calculates the concentration at each pixel as the
ratio of total contaminant volume upstream di-
vided by the total volume of scour depth for all
values less than 1/ X are made zero in the grid,
active channels upstream. Mathematically, the
leaving ¬nite values in the remaining pixels.
model states that the concentration at a pixel (i,j)
Fourth, the contributing-area grid is converted
is given by
to a unit-discharge grid using empirical data (i.e.
a ¬‚ood-envelope curve), which, in turn, is con-
h l,m + h i, j · Hl,m + Hi, j
verted to a scour/mixing-depth grid using the min up up
C i, j =
empirical square-root relationship documented
Hl,m + Hi, j
by Leopold et al. (1966). Following this series up

of operations, the resulting grid represents the (3.4)

where h l,m is the effective contaminant thickness Zone (CFVZ) (Crowe and Perry, 1990). Its relative
at pixel (l,m), Hl,m is the scour depth, and the youth and location with respect to Yucca Moun-
sums are computed over all pixels (l,m) upstream tain has made it a focus of intense study as a
from pixel (i, j). possible analog for an eruption through the pro-
The fundamental tool for the scour-dilution- posed nuclear waste repository.
mixing model is the MFD algorithm. In this al- There are two principal drainages that trans-
gorithm, contributing area, contaminants, and port tephra from Lathrop Wells (Figure 3.10).
scour depths are routed from each pixel to all The western drainage system transports material
of its neighboring pixels weighted by topographic from the exposed tephra sheet on the northwest
slope. This algorithm is implemented by initializ- side of Lathrop Wells cone west and south into
ing the grid with the quantity to be routed. In the Amargosa Valley. The eastern drainage system
case of contributing area, the grid is initialized heads near the northern margin of the tephra
with the area of that pixel, x 2 (e.g. 900 m2 for a sheet and transports material around the east-
30-m resolution DEM). In the case of contaminant ern side of the Lathrop Wells cone and adja-
routing, each pixel is equal to the effective thick- cent lava ¬‚ows. Twenty ¬ve samples of channel-
ness of contaminant in channels of the source re- bed sediments were collected along these two
gion. The MFD or bifurcation-routing algorithm channel systems in order to evaluate the signif-
then ranks all grid points in the basin from high- icance of the dilution process and validate the
est to lowest in elevation. Routing is then calcu- scour-dilution-mixing model. In small channels
lated successively at grid points in rank order on the tephra sheet, the bottom of the scour
from highest to lowest, thereby ensuring that zone was clearly visible and ranged from 12 to
all incoming areas of contaminant have been 29 cm in depth in channel-bed-sediment sam-
accounted for before downstream routing is ple pits (Figure 3.11). In larger channels down-
performed. stream from the tephra sheet, the bottom of the
To test the model, we consider the concen- scour zone was not visible within our sample pits
tration of basaltic sediment in channels drain- (which were limited to 30--40 cm in depth). Sam-
ing the Lathrop Wells tephra sheet. This tephra ples were taken by uniformly sampling pit-wall
sheet provides a concentrated source of ash and material to a depth of 30 cm, or to visible scour
lapilli that acts as a contaminant source for depth if less than 30 cm. Samples were sieved
downstream channels. The tephra found in chan- and material larger than silt size was separated
nels downstream from the tephra sheet is com- into basaltic and non-basaltic fractions with a
parable in texture to the non-basaltic sand of magnet and visually checked for complete sep-
downstream channel beds. As such, the model as- aration using a stereomicroscope. The basaltic
sumption that contaminant material and uncon- material was strongly magnetic and total sepa-
taminated channel sediments are transported at ration error was estimated to be less than 5%
similar rates is likely to be a good approximation by mass. Non-basaltic material was composed of
in this case. In addition, the proximity of this Miocene welded tuff and eolian sand. The frac-
area to Yucca Mountain (i.e. 18 km south of the tion of basaltic sediment by mass was calculated
proposed repository location) makes it an ideal for each sample. This value was converted to a
analog (e.g. similar climate, hydrology, and geo- fraction by volume (for comparison to the volu-
morphology) for tephra redistribution following metric dilution-mixing model) using an average
measured basalt density of 1.5 g/cm3 (a relatively
a hypothetical volcanic eruption at Yucca Moun-
tain. The physical volcanology of the Lathrop low value for rock because of its high vesicular-
ity) and a tuff density of 2.7 g/cm3 .
Wells volcanic center was described in detail by
Mixing occurs along both the eastern and
Valentine et al. (2005). At 77 000-years old (Hei-
western drainages according to the trend ev-
zler et al., 1999), the Lathrop Wells cone is the
ident in the plots of basalt concentration as
youngest basaltic volcano in the Yucca Mountain
a function of distance from the channel head
region. It is the southern-most surface expression
in Figure 3.10b. These data indicate that the
of the Plio--Pleistocene-age Crater Flat Volcanic

116.51° W

measured data
36.70 N
1505.02 3D model
1505.01 1505.04 1505.05 1505.06
1405.04 1405.01
1405.02 tephra sheet
east channel 0.1
1405.06 1405.03
BG BG 1505.08
major tributary
1505.07 junctions
east channel
3 5
Lathrop Wells
x (km)
36.69 1405.10 lava flows
Volcanic Cone
west channel
measured data

1505.09 3D model

1505.12 0.3
N west channel
x (km)
0.5 km
(a) (b)
are relatively high to the north of the cone where
Fig 3.10 (a) Sample location map for tephra concentration
most of the tephra was deposited (11.1% and 6.8%
in the east and west channels draining the Lathrop Wells
in samples SD091405.06 and SD091505.05, respec-
tephra sheet. (b) Comparison plots of measured and modeled
predictions for tephra concentration as a function of distance tively) and lower to the south (3.0% in sample
along-channel (note logarithmic scale on vertical axis, needed SD091505.08). These background samples were
to visualize data over 1“2 orders of magnitude). Also graphed collected in large drainages that provide repre-
are best-¬t exponential curves. Modi¬ed from Pelletier et al.
sentative average values over large areas. In or-
der to resolve the two background values in the
model, a ternary mask (Figure 3.12a) was used
with a source concentration of 73%, a distal
concentration of basaltic tephra is reduced by source value of 9% (north of the primary tephra
more than 50. sheet) and background value of 3% (everywhere
The Lathrop Wells tephra sheet was mapped else).
by Valentine et al. (2005), providing a prelimi- Additional inputs to the model are the thresh-
nary source map for input to the scour-dilution- old stream power and a 10-m resolution US Geo-
mixing model for this application. The source logical Survey DEM (Figure 3.12a). The value for
map is complicated in this case by the fact that the stream-power threshold was chosen to be
tephra from the 77 ka eruption dispersed tephra 0.015 km by forward modeling. In this forward-
in all directions. Samples collected in the study modeling procedure, a map of incised channels
reveal three separate groups of values. On the was made by comparing the stream power at
tephra sheet itself, contaminant concentrations each point in the grid to a threshold value. This
vary between 61% and 83%, with an average value threshold was varied in value to produce differ-
of 73%. Therefore, 73% was used as the value ent incised-channel maps. Maps were then over-
for the source concentration in the model. Off lain and visually compared with US Geological
of the tephra sheet, background concentrations Survey orthophotographs to identify which value

Fig 3.11 Field photos illustrating
(a) (b)
the three-dimensional pattern of
contaminant dilution near Lathrop
Wells tephra sheet. (a) and (b)
Channel pits on the tephra sheets
scour zone
expose a ¬‚uvially mixed scour zone
76% tephra
15 cm scour zone ranging from 12 to 29 cm in
29 cm SD091405.08
70% tephra
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
debris flow
tephra predominantly by Miocene volcanic
sheet tuff and eolian silt and sand
transported from the upper tephra
sheet. (c) The effects of dilution
middle sheet lower sheet visible as high-concentration
channels draining the tephra sheet
(c) (channel at right, 76% tephra) join
with low-concentration channels (at
“low”-concentration “high”-concentration
left, 26% tephra). Tephra
tephra sheet
channel channel
concentration in these channels
26% tephra (SD091405.07) 76% tephra (SD091405.08)
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. For color version, see
plate section. From Pelletier et al.
(2008). Reproduced with permission
of Elsevier Limited.

“mixed” channel
tributary junction
41% tephra (SD091405.09)

of the critical stream power produced the most of drainage basins in the Yucca Mountain region.
accurate map. Maximum discharge is used because the largest
In the ¬rst step of the model run, all pix- ¬‚oods control the long-term scour/mixing depth.
els of the contributing area grid were initial- The unit-discharge map was then used to create
ized to 10’4 km2 (i.e. the pixel area for a 10-m- a relative map of scour/mixing depth in channels
resolution DEM). The MFD algorithm was then (Figure 3.12c) using the square-root dependence
used to calculate the contributing-area grid. In of Leopold et al. (1966) and a stream-power
the second step, the contributing-area grid was threshold of 0.015 km. Only a relative grid of
converted to a unit discharge map using the re- scour/mixing depth was needed in this case be-
gional ¬‚ood-envelope curve of Squires and Young cause the scour depths in the source area were
(1984). These authors found the maximum dis- used to calculate the effective volume of contam-
charge to be proportional to the contributing inant routed downstream, and the scour depths
area to the 0.57 power using available stream downstream of the source region were used to
gage data and paleo¬‚ood estimates from a range dilute that contaminant. Since the contaminant

101 km
116.52 W 116.50 10
116.48 10
N1 km
36.74° N


tephra sheet
36.70 slope*area1/2


LathropWells cone
relativete phra
0.01 0.1 1.0

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

Fig 3.12 Model prediction for tephra concentration and
tration (as opposed to an effective thickness), the
scour/mixing depth downstream from the tephra sheet of the
absolute value of the scour depth is not needed to
Lathrop Wells volcanic center. (a) Model inputs include a
predict concentration values downstream; only
10-m resolution US Geological Survey Digital Elevation Model
relative values are needed. The model predicted
(DEM) of the region and a grid of source and background
tephra concentrations (expressed as a fraction) in
concentrations. Input concentration values are: 73% (Lathrop
Wells tephra sheet), 9% (distal source region), and 3% all channels draining the Lathrop Wells tephra
(background). (b) Map of stream power in the vicinity of the sheet. This map has a quadratic gray scale vary-
source region. Scale is logarithmic, ranging from 10’2 to ing from < 0.01 (black) to 1.0 (white).
101 km. (c) Map of scour/mixing depth (scale is quadratic,
The numerical model prediction was com-
ranging from 10 cm to 1 m). (d) Map of tephra concentration
pared to the measured tephra concentrations by
(scale is quadratic, ranging from 0.001 to 1). For color
extracting tephra concentrations along channel
version, see plate section. From Pelletier et al. (2008).
pro¬les of the western and eastern channels in
Reproduced with permission of Elsevier Limited.
Figure 3.10b. Along both channels, concentra-
tions are near their maximum values for ap-
proximately the ¬rst 1.0 km downstream from
concentration is a ratio of these quantities,
the channel head. Along the eastern channel
the concentration was unaffected by scaling the
(Figure 3.10b, top graph), abrupt concentration
scour depth up or down by a constant factor. In
declines are associated with major tributary
all cases where the source is speci¬ed as a concen-

junctions at 1.6 and 2.7 km from the channel servations at the San Francisco volcanic ¬eld sug-
head. The ¬rst major step down is associated gest that hillslopes above a threshold gradient
of approximately 0.3 m/m (17—¦ ) are susceptible to
with an increase in discharge and scour depth
associated with a relatively uncontaminated dis- mass movements and/or rilling within 1 kyr af-
tributary channel entering the channel from the ter an eruption. The threshold stream power for
northeast just upstream of the U-shaped chan- the Fortymile Wash drainage basin was estimated
nel bend. Another abrupt drop in concentration to be 0.05 km (corresponding to a drainage den-
sity of X = 20 km’1 ), based on forward modeling
occurs approximately 2.7 km from the channel
head (just downstream of sample SD091515.07) and comparison with USGS orthophotographs.
where another large tributary enters. As the For simplicity, contaminated tephra mobilized
channel curves around the southeast corner of from steep slopes and channels was assumed to
the Lathrop Wells lava ¬eld, it bifurcates into two be transported into the channel system instantly
distributary channels. Each pair of samples lo- following the eruption, although in nature this
cated at similar distances from the channel head process will take place over decades to centuries.
were averaged to comprise the ¬nal two points As an example, input to the tephra redistribu-
on the curve. The overall downstream-dilution tion model was constructed based on the fallout
pattern can be roughly approximated by an ex- from a hypothetical eruption with a northeast-
ponential function with a characteristic length directed plume computed using the ASHPLUME
of 0.9 km. For every 0.9 km downstream, there- numerical model (Figure 3.13a) (Jarzemba, 1997;
fore, the tephra concentration decreases by ap- Jarzemba et al., 1997; BSC 2004a,b). The hypothet-
proximately a factor of 1/e. Along the western ical fallout distribution used in Figure 3.13a cor-
responds to an eruption power of 9.55 — 109 W
channel, concentrations follow a similar pattern
in the bottom graph of Figure 3.10b. Note that (a moderate-to-large magnitude eruption for the
this pro¬le extends only about half as far from Crater Flat volcanic ¬eld), a duration of approx-
the channel head as the eastern channel, due to imately 3 days, a wind speed of approximately
the presence of the mine road which disturbs the 5 m/s, and an eruption velocity of 2 m/s. The
western channel about 2.5 km from the channel contributing-area grid is used to calculate the
head. scour-depth grid in several steps for this appli-
The ability of the model to represent ob- cation. First, contributing area is related to peak
served trends of downstream dilution at Lathrop discharge on a pixel-by-pixel basis using the re-
Wells provides a basis for using the model to gional ¬‚ood envelope curve of Squires and Young
predict dilution patterns following a hypothet- (1984), as in the Lathrop Wells example. Sec-
ical eruption at Yucca Mountain. To do this, a ond, discharge is related to scour depth on a
spatially distributed model was used to calcu- pixel-by-pixel basis using unit discharge with the
late the volume of tephra fallout mobilized from square-root dependence of Leopold et al. (1966).
steep slopes and active channels (Figure 3.13). Together, these relationships imply that scour
The model accepts input from the ASHPLUME depth is proportional to contributing area to the
numerical model (e.g. Figure 3.13a, assuming a 0.29 power. To determine the proportionality con-
southwesterly wind). The ASHPLUME model is stant needed to relate contributing area to a spe-
a numerical advection-dispersion-settling model ci¬c value of scour depth, the model uses mea-
that considers variable grain sizes and wind- surements of the US Geological Survey during the
speed pro¬les from the ground to 10 km altitude. 1995 ¬‚ood at the Narrows section of Fortymile
The model then calculates the DEM slopes and Wash (Beck and Hess, 1996) located close to the
contributing areas (Figures 3.13b and 3.13c). Fall- fan apex. The contributing area at this location is
the ratio of the basin area, 670 km2 , to the width,
out tephra is assumed to be mobilized if it lands
44 m. The maximum scour depth at this location
on pixels with hillslope gradients greater than
was measured to be at least 1.14 m, and had an
a threshold value equal to S or with a stream
average value of 0.72 m. The model runs of this
power greater than a threshold value equal to
paper used 1 m as a representative average scour
1/ X , where X is the drainage density. Field ob-

U slope*area1/2
tephra thickness slope

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 (17o) > 0.05 km
0.05 km
Fig 3.13 Digital grids used in the modeling of tephra depth for this reach. Combining these measure-
redistribution following a volcanic eruption at Yucca ments yields
Mountain. (a) Shaded relief image of DEM and Fortymile
Wash drainage basin (darker area). Tephra from only this A i, j 44
Hi, j = Hapex (3.5)
drainage basin is redistributed to the RMEI location.
670 x
(b) Grayscale image of DEM slopes within the Fortymile
Wash drainage basin. (c) Black-and-white grid of areas in the where Hi, j is the scour depth at any point in the
drainage basin with slopes greater than 17—¦ .
basin upstream from the apex, Hapex = 1 m is the
(d) Black-and-white grid of active channels in the drainage
representative scour depth at the fan apex based
basin (de¬ned as pixels with contributing areas greater than
on measurements by the US Geological Survey,
0.05 km2 ). All tephra deposited within the black areas of
and A i, j is the contributing area at point (i, j).
(c) and (d) are assumed to be mobilized by mass movement,
The scour-dilution-mixing model predicts rel-
intense rilling, or channel ¬‚ow. For color version, see plate
atively high tephra concentrations along the
section. Modi¬ed from Pelletier et al. (2008). Reproduced
northeast-directed plume track (Figure 3.14b). As
with permission of Elsevier Limited.
tephra from these channels moves downstream,
tephra in channels that are integrated into the

mately nearly 2% for the southwesterly wind case
scour/mixing depth tephra concentration
illustrated in Figure 3.14. The primary control on
1.0 m 10%
1% 100%
0.1 0.2
tephra concentration at the outlet is the topog-
(a) raphy in the primary fallout region (which con-
trols the volume of tephra mobilized on steep
slopes and in active channels). The highest out-
let tephra concentrations occur for southerly and
westerly winds (Figure 3.15), because these wind
conditions transport tephra north and east from

the repository location to areas with a relatively
high density of tephra-transporting steep slopes
and active channels within the Fortymile Wash
drainage basin.

3.1 Download a DEM of any area dominated by ¬‚uvial
processes. Implement the fillinpitsandflats
concentration routine (Appendix 2) to ¬lter out pits and ¬‚ats
southwesterly = 1.97% in the DEM. Implement the MFD ¬‚ow-routing
method to calculate the contributing area A for
each pixel and plot the results as a grayscale or
N 3 km color map.
3.2 After completing the above exercise, develop a
program module that distinguishes hillslope and
Fig 3.14 Digital grids output by the modeling of tephra
channel pixels in the DEM based on a threshold
redistribution following a potential volcanic eruption at Yucca
value of contributing area (i.e. assume channel
Mountain. (a) Map of scour-depth grid, calculated from
pixels are those with A > A c ). Create maps of the
contributing area grid. (b) Grayscale map of tephra
channel network corresponding to several differ-
concentration in channels, calculated using the
ent values of A c and compare the map to observed
scour-dilution-mixing model. For color version, see plate
channel-head positions in the DEM.
section. Modi¬ed from Pelletier et al. (2008). Reproduced
3.3 After completing the above exercise, use the FTCS
with permission of Elsevier Limited.
or ADI technique to solve the 2D diffusion equa-
tion for the evolution of the hillslope pixels in
the DEM. Assume that the channel network has
main Fortymile Wash becomes progressively di- a ¬xed elevation through time. This exercise will
luted with each successive in¬‚ow of relatively require using a mask grid to identify the channel
uncontaminated tributary sediments. Note that pixels and to ensure that only hillslope pixels are
some of the tephras draining from channels in varied through time.
the vicinity of Yucca Mountain are locally de- 3.4 Modify the fillinpitsandflats routine to map
posited in alluvial ¬‚ats and do not reach the the area of all closed depressions in a DEM by
including a counter that is incremented during
main Fortymile Wash. At the basin outlet the ef-
each loop through the routine. Download a USGS
fects of dilution mixing and localized deposition
DEM of an area that includes depressions (e.g.
of tephra limit the tephra concentration to 1.97%
karst, glacial lakes, etc.) and use the modi¬ed rou-
in this example. The effects of variable wind di-
tine to map the depressions.
rection on contaminant concentrations can be
3.5 Download a USGS Seamless DEM of an alluvial
determined (Figure 3.15) by running the scour-
channel or alluvial fan. Use the successive ¬‚ood-
dilution-mixing model on a sequence of rotated ing algorithm with MFD routing to model the in-
plumes of the same size and shape as the ex- undated area associated with some hypothetical
ample in Figure 3.14. Outlet tephra concentra- ¬‚ood depth. This will require that you identify a
tions for these rotated plumes vary from a mini- pixel or set of pixels and prescribe an input ¬‚ow
mum of less than 0.1% to a maximum of approxi- depth.

southerly westerly northerly
winds winds winds
(a) (c)

concentration concentration concentration
at outlet = 0.73% at outlet = 0.48% at outlet = 0.076%

N 3 km
tephra concentration
Fig 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
the flow grid with a uniform depth of runoff (e.g.
winds. Tephra concentration at the basin outlet varies from a
1 cm) for a DEM of your choice. As runoff is routed
maximum of 0.73% (for southerly winds) to a minimum value
from each pixel to its neighbors (proceeding from
of 0.076% (for northerly winds) for this eruption scenario.
high elevations to low elevations in rank order),
Southerly winds result in higher concentrations because the
estimate the transit time between successive grid
high relief of the topography north of the repository is
points using a ¬‚ow velocity calculated by Man-
capable of mobilizing more tephra than other wind-direction
ning™s equation or assuming critical ¬‚ow. Assume
scenarios. For color version, see plate section. Modi¬ed from
no in¬ltration. At each pixel, multiple up-slope
Pelletier et al. (2008). Reproduced with permission of Elsevier
neighbors will contribute runoff, and each neigh-
bor will have a different transit time. Adopt the
transit time associated with the greatest incom-
3.6 Develop a numerical model based on the MFD ing discharge. Integrate the transit times for each
routing algorithm that estimates the time to peak step to estimate the time to peak discharge at the
discharge for a ¬‚ash ¬‚ood. Begin by initializing basin outlet.
Chapter 4

The advection/wave equation

advection models. In the diffusion model, which
4.1 Introduction is applicable to certain transport-limited hill-
slopes and alluvial channels, the rate of erosion
or deposition is proportional to the change in gra-
In this chapter we will focus on the advection
dient, or curvature. In the diffusion model, ma-
equation, which in its simplest form is given by
terial removed from the upper portion of the
‚h ‚h
=c (4.1) slope break is deposited on the lower portion,
‚t ‚x
leading to a smoother topography over time. In
In general, advection equations involve a lat- the advection model, in contrast, the rate of ero-
eral translation of some quantity (e.g. tempera- sion is directly proportional to elevation gradi-
ture, topography, chemical concentration, etc.). ent, and deposition is not allowed. In the ad-
The coef¬cient c in Eq. (4.1) has units of length vection model, slope breaks undergo a lateral
over time and represents the speed at which the translation upstream with no change in shape.
quantity h is advected. In the context of land- The basic stream-power model of bedrock chan-
form evolution, the advection equation is used nel evolution with m = 1/2 and n = 1:
to model retreating landforms, including cliffs,
‚h 1 ‚h
banks, and bedrock channel knickpoints. The ad- = KwQ 2 (4.4)
‚t ‚x
vection equation is closely related to the classical
wave equation which is given by is simply an advection equation with spatially
variable advection coef¬cient. Conceptually, the
‚ 2h ‚ 2h
= c2 2 (4.2) stream-power model says that the action of
‚t 2 ‚x
bedrock channel incision can be quanti¬ed by
Equation (4.1) looks like the square root of the
˜˜advecting™™ the topography upstream with a lo-
classical wave equation. In the classical wave
cal rate proportional to the square root of dis-
equation, disturbances propagate as waves in
charge. It should be noted that ‚h/‚t in Eq. (4.4)
both directions (i.e. +x and ’x). In the basic ad-
is always negative (i.e. erosion). Therefore, if x is
vection equation, disturbances propagate in one
the distance from the divide, ‚h/‚ x is always neg-
direction only. In 2D, the simple advection equa-
ative and no negative sign is required in front of
tion has the form
K w.
‚h Advection equations have many applications
= c · ∇h (4.3)
‚t in Earth surface processes. The transport of heat,
In this case, the advection coef¬cient is a vector sediment, and chemical species in any geologi-
with components in the x and y directions. cal ¬‚uid (e.g. the atmosphere, ocean, lakes, and
Figure 1.7 contrasts the evolution of a hill- rivers) includes some element of advective trans-
slope or channel governed by the diffusion and port. Many physical processes combine effects of

diffusion and advection. Heat is advected in the temperature ¬eld, in turn, controls the ¬‚uid ¬‚ow
ocean, for example, as a result of well-established through the effects of temperature on density
currents (e.g. the Gulf Stream). Superimposed on and hence buoyancy. In such cases, numerical ap-
that net transport is a spreading of heat due to proaches are generally required.
the cumulative effects of many small eddies su- The methods we will describe in this chap-
perimposed on the current. This spreading is dif- ter can be used to solve a large set of equa-
fusive in nature. Not surprisingly, the advection tions known as ˜˜¬‚ux-conservative™™ equations.
equation features prominently in nearly all cli- Flux-conservative equations have the form
mate models, where the transport of heat and ‚h ‚F
=’ (4.6)
chemical species around the globe must be mod- ‚t ‚x
eled with high precision. In the last few decades,
where F is the ¬‚ux. Equation (4.6) reduces to
some of the best new techniques for solving the
the simple advection equation, for example, if
advection equation (e.g. semi-Lagrangian meth-
F = ’c h and c is a constant.
ods) have originated in the climate and weather
modeling communities.
4.2.1 Method of characteristics
The method of characteristics can be used to solve
4.2 Analytic methods any linear advection equation and even some
nonlinear advection equations. Here we present
the method as it is applied to the general linear
The behavior of the simple advection equation
advection equation of the form
(Eq. (4.1)) is to propagate disturbances in one di-
‚h ‚h
rection with velocity c. As such, the analytic so- + b(x, t) + c(x, t)h = 0
a(x, t) (4.7)
‚x ‚t
lution can be written as
subject to the initial condition h(x, 0) = f (x). The
h(x, t) = f (x + ct) (4.5)
goal of the method of characteristics is to change
where h(x, 0) = f (x) is the initial condition for h. the coordinates from (x, t) to a new coordinate
system (x0 , s) in which the advection equation (a
In more complex advection equations where the
advection coef¬cient is not constant but is a pre- PDE) is transformed into an ODE along certain
curves in the x ’ t plane. If we choose a new co-
scribed function of space and time, the method
of characteristics can be used to determine an- ordinate s such that
alytic solutions in many cases. The idea of the dx
= a(x, t) (4.8)
method of characteristics is to take advantage ds
of the lateral translation inherent in advection
problems. Rather than using the local shape of
= b(x, t)
the function to evolve the system forward in (4.9)
time (as in explicit approaches to the diffusion
then the advection equation becomes, in the new
equation), the method of characteristics seeks to
coordinate system,
track the trajectories of individual segments of
the solution laterally and through time. In many dx ‚h dt ‚h ‚h ‚h
= + = a(x, t) + b(x, t)
cases, these trajectories are described by simple ds ‚ x ds ‚t ‚x ‚t
functions that can be quanti¬ed analytically. (4.10)
In general, the advection coef¬cient c in the
Therefore, along the ˜˜characteristics curves™™
advection equation is a function of space and
given by Eqs. (4.8) and (4.9), the advection equa-
time that coevolves with h. In the stream-power
tion becomes an ODE of the form
model, for example, the advection coef¬cient
K w A m varies through time as individual grid + c(x, t)h = 0 (4.11)
points compete for drainage. Thermally-driven
convection is another example. In that case, the The strategy for applying this series of equa-
¬‚uid ¬‚ow advects the temperature ¬eld, and the tions for speci¬c linear advection equations is as

10 2
(a) (c)

t/c t/c

5 1

0 0
1.5 1.5
(b) (d)

1.0 1.0
h h
8 6 4 2 2 1.5 1 0.5
0.5 0.5

0.0 0.0
t/c = 0
t/c = 0
’0.5 ’0.5
0 2 4 6 8 10 0 2 4 6 8 10
x x
characteristic equation (Eq. (4.9)) becomes
Fig 4.1 Application of the method of characteristics to
dt/ds = 1, and integration gives the solution
solve Eqs. (4.1) and (4.12). (a) and (b) Plots of (a)
t = s . Next, the characteristic equation (Eq. (4.8))
characteristic curves and (b) model evolution for the simple
advection equation (Eq. (4.1)) for a hypothetical channel becomes x (t) = c with the initial condition
knickpoint. (c) and (d) Plots of (c) characteristics curves and
x(0) = x0 . The solution, by integration, gives
(d) model evolution for the advection equation with
the characteristic curves x = x0 + ct. Next, the
coef¬cient c x . In this case, the velocity of the knickpoint
solution to the ODE (Eq. (4.11)), h (t) = 0 with
decreases and the gradient steepens as the knickpoint is
initial condition h(0) = f (x0 ), is h(t) = f (x0 ). The
propagated into the drainage headwaters.
solution of the PDE is, transforming back to the
original coordinates, h(x, t) = f (x0 ) = f (x + ct).
Figure 4.1a plots the solution along with the
follows. First, we solve the two characteristic
characteristic curves. By writing the character-
equations, Eqs. (4.8) and (4.9), using integration.
istic curves as t = ’x/c ’ x0 /c, we can see that
The constants of integration are constrained by
the characteristic curves are parallel lines with
setting x(0) = x0 and t(0) = 0. Second, we solve
slope ’1/c in the x ’ t plane.
the ODE (Eq. (4.11)) with the initial condition
Next we consider a case in which the ad-
h(0) = f (x0 ), where x0 are the initial points of the
vection coef¬cient varies with distance x. As
characteristic curves along the t = 0 axis in the
a concrete example, consider the evolution of
x ’ t plane. Third, we solve for s and x0 in terms
a bedrock channel evolving according to the
of x and t and substitute these values into h(x0 , s)
stream-power model with K w Q 1/2 = c x, where x
to get the solution of the original PDE.
is distance from the divide:
First we consider the simple advection
equation with constant coef¬cient c. In this ‚h ‚h
= cx
case, b(x, t) = 1 and c(x, t) = 0. Therefore, the (4.12)
‚t ‚x

This equation assumes that the square root of dis- 1.5
charge increases linearly with distance from the
divide, which is appropriate for a humid drainage
system (i.e. Q ∝ A) where drainage area A in- t/c = 0
creases as the square of distance from the divide
(i.e. an approximately semi-circular basin). Equa-
tion (4.12) corresponds to the general, ¬rst-order
PDE with b(x, t) = 1, c(x, t) = 0, and a(x, t) = c x.
As in the simple advection equation, the char-
acteristic equation (Eq. (4.9)) yields t = s. The
second characteristic equation to solve is x (t) =
c x with initial condition x(0) = x0 . Solving it gives
the characteristic curves x = x0 exp(’c x). As be- 4 10
0 2 6
fore, h(t) = f (x0 ). Finally, the characteristic curves
are given by x = x0 exp(c x), so x0 = x exp(’c x). Fig 4.2 Instability of the FTCS method applied to the
Transforming back to the original variables gives simple advection equation for a Gaussian (bell curve) initial
h(x, t) = f (x0 ) = f (x exp(’c x)). Figure 4.1b illus- condition (results obtained from modi¬ed version of Matlab
code distributed by Marc Spiegelman). The solution becomes
trates solutions of this equation with the char-
wildly unstable after the model has advected the Gaussian a
acteristic curves also shown. In this model, the
distance of only a few times greater than its width.
knickpoint propagates more slowly and becomes
steeper as it moves upstream into areas with
lower stream power. some speci¬c applications, Lagrangian methods
In the ¬rst set of papers to present the method are hard to use when the characteristic curves
of characteristics in the concept of landscape evo- are complex. To understand how a Lagrangian
lution, Luke (1972, 1974) showed how the method method works, imagine the mixing of dye in a
could be used to solve the nonlinear advection stirred ¬‚uid. As the dye becomes stretched and
equation mixed into the ¬‚uid, the Lagrangian method
produces a tortuous grid that follows the ¬‚uid.
‚h ‚h
=c (4.13) As the system evolves, determining neighboring
‚t ‚x
grid points (to model viscous interactions, for ex-
for a knickpoint initial condition. Luke (1972, ample) requires complex bookkeeping. Eulerian
1974) showed that in nonlinear equations like methods are the most versatile, provide the same
Eq. (4.13), the slope pro¬le develops ˜˜shocks™™ output format regardless of the speci¬c model
due to the fact that certain slope segments are runs, and are the method of choice for most
advected faster than others. In these nonlinear applications.
cases, the method of characteristics can be ex-
tended to incorporate these effects.
4.3.1 Failure of the FTCS method
In Chapter 2, we introduced the concept of ¬nite
differencing with the FTCS method. Applied to
4.3 Numerical methods
the basic advection equation, the FTCS method
Methods for solving the advection equation can
h i+1 ’ h i’1
h in+1 ’ h in n n
be classi¬ed as either Eulerian or Lagrangian. Eu- (4.14)
t 2x
lerian methods use a ¬xed grid while Lagrangian
methods use an adaptive grid that follows the Unfortunately, the FTCS method is inherently
characteristic curves of the system. In this book unstable when it is applied to advection equa-
we will consider Eulerian methods exclusively. tions. Figure 4.2 illustrates the evolution of a
Although Lagrangian methods are useful for Gaussian initial condition being advected to the

right. In diffusion problems, the smoothing na- it has the same form as Eq. (4.14):
ture of the equations results in a stable method h i+1 ’ h i’1
h in+1 ’ h in n n
(for suf¬ciently small time steps). The advection t 2x
equation lacks any inherent smoothness, there- 1 h i+1 ’ 2h in + h i’1
n n
+ (4.17)
fore small numerical errors introduced by the 2 t
FTCS method continually grow over time until
Equation (4.17) is the FTCS representation of the
they swamp the actual solution. Figure 4.2 also
illustrates what numerical instabilities often look
‚h ‚h ( x)2 ‚ 2 h
like: zig-zag functions that oscillate wildly from
=c + (4.18)
‚t ‚x 2 t ‚ x2
one grid point to the next.
Equation (4.18) illustrates that the Lax method
4.3.2 Lax method introduces a small diffusion term into the
The instability of the FTCS method can be cured
by a simple change: replace the term h in in the The problem of numerical diffusion turns out
to be less than fatal, however. Smolarkiewicz
time derivative by its average value:
(1983) noted that since we can quantify the mag-
h in ’ (h + h i’1 ) nitude of the numerical diffusion, we can make
2 i+1 a correction to reverse its effect. We describe this
˜˜Smolarkiewicz correction™™ procedure below in
With this change, the FTCS discretization be-
the context of upwind-differencing. In addition,
the numerical diffusion effect can be minimized
1n ctn
h in+1 = (h i+1 + h i’1 ) + (h i+1 ’ h i’1 )
n n
(4.16) by taking a half-step in time. This procedure is
2 2x
called the two-step Lax--Wendroff method.
This scheme turns out to be stable for suf¬ciently
small time steps (i.e. t ¤ x/c). The stability of 4.3.3 Two-step Lax“Wendroff method
the Lax method and other discretization schemes Thus far we have focused on the effects of dif-
can be analyzed using the von Neumann stabil- ferent methods of spatial discretization. In all of
ity analysis (Press et al., 1992). The stability cri- our methods, the updating of h i between time
terion t ¤ x/c (also called the Courant con- step n and n + 1 happens in a single ˜˜jump.™™ In
dition) is broadly applicable to many advection this section we explore an implementation of the
equations. In cases where c varies in space and Lax method that uses two steps to perform the
time, the Lax method will be stable everywhere updating (i.e. the temporal discretization). In this
if the value of t is chosen to ensure stability intermediate step, one de¬nes intermediate val-
where c has its largest value. Conceptually, the ues h i+1/2 at the half time steps tn+1/2 and half
Courant condition ensures that changes in the grid points xi+1/2 . These intermediate values are
value of h at one point on the grid travel along calculated using the Lax method:
the grid more slowly than the advection velocity 1n t
h i+1/2 = (h i+1 + h in ) ’ (F i+1 ’ F in )
c. If the time step is too large, changes propagate
2 2x
along the grid faster than the advection veloc-
Using Eq. (4.19), the user then calculates the
ity, and this unphysical condition will cause the
¬‚uxes at the half grid points and half time steps,
solution to become unstable. n+1/2
F i+1/2 . Then, the values of h in+1 are calculated by
It turns out that there is a price to be paid for
the space-centered expression
the stability introduced by the Lax method. The
t n+1/2
Lax method introduces a small but nonzero dif- n+1/2
h in+1 = h in ’ (F i+1/2 ’ F i’1/2 ) (4.20)
fusive component to the solution. In a sense, the
Lax method is not solving the advection equation The two-step Lax--Wendroff method eliminates
at all, but rather an advection--diffusion equa- most of the numerical diffusion of the Lax
tion. To see this, let™s rewrite Eq. (4.16) so that method.

4.3.4 Upwind-differencing method 1.5
Another simple but powerful numerical method
t/c = 0
for the advection equation is called upwind dif-
ferencing. Upwind differencing involves calculat-
ing the slope along the direction of transport:
§ hn ’ hn 8 6 4 2
⎪ i+1
⎨ if c in > 0
h in+1 ’ h in x
= ci
⎪ h in ’ h i’1
t © if c in < 0
x 0.0
In the FTCS technique, the ˜˜centered-space™™ gra-
dient is calculated by taking the difference be- upwind differencing
tween the value of the grid point to the left (i.e.
at i ’ 1) and the value to the right (i.e. at i + 1) (b)
of the grid point being updated. This approach
is prone to instability because it does not make
use of the value at the grid point i itself. As such,
the difference h i’1 ’ h i+1 can be small even if the
n n
8 6 4 2
value of h in is wildly different from the values on
either side of it. In effect, the FTCS method cre-
ates two largely decoupled grids (one with even
i, the other with odd i) that drift apart from
t/c = 0
each other over time. In the upwind method, upwind differencing
that problem is corrected by calculating gradi- with correction
ents using only one adjacent point (i.e. h i’1 ’ h in
4 10
0 2 6
or h in ’ h i+1 ). Which direction to choose? If the

¬‚ux of material is moving from left to right, then
Fig 4.3 Solution to Eq. (4.1) for an initial condition of a
physically it makes sense that the value of h in+1 hypothetical knickpoint (a) with upwind-differencing and
n n
should depend on h i’1 , not on h i+1 . If the ¬‚ux of (b) including the Smolarkiewicz correction. Results obtained
material is in the opposite direction, h in+1 should using modi¬ed version of Matlab code distributed by Marc
depend on h i+1 . The upwind-differencing method Spiegelman. Without the correction, the knickpoint in
implements this approach. (a) gradually acquires a rounded top and bottom. With the
correction, the initial knickpoint shape is preserved almost
Upwind differencing also suffers from the nu-
exactly as it is advected upstream.
merical diffusion of the Lax method. In order to
reverse the effects of numerical diffusion, we fol-
low the procedure of Smolarkiewicz (1983). Given This formulation shows that the effects of nu-
a small numerical diffusion characterized by the merical diffusion can be written in terms of a
equation diffusion velocity. The Smolarkiewicz correction
involves taking a small advective step with anti-
‚h ‚ ‚h
= D num (4.22) diffusion velocity ’Vd with the upwind differenc-
‚t ‚x ‚x
ing method following the original upwind differ-
where D num is a numerical diffusivity, the effects encing step (Eq. (4.21)) with velocity c in . For the
of numerical diffusion can be written as upwind-differencing method, D num = (|c in | x ’
‚h ‚(Vd h) t(c in )2 )/2.
= (4.23)
‚t ‚x Figure 4.3 illustrates the solution to the ba-
sic advection equation with upwind differenc-
ing with and without Smolarkiewicz correction.
’ D num if h > 0 The advantage of upwind differencing with Smo-
Vd = ‚x
h (4.24)
0 if h = 0 larkiewicz correction is that it represents a good

compromise between accuracy and ease of imple- that river knickpoints have not yet reached the
mentation that makes it the method of choice upper plateau surface, the large-scale morphol-
for many advection problems. The primary limi- ogy of the Sierra Nevada suggests that 1.5 km of
tation of the method is the relatively small time late Cenozoic surface uplift occurred (Clark et al.,
step (given by the Courant condition), which ef- 2005). This conclusion has been challenged, how-
fectively limits the grid resolution that can be ever, by stable-isotope paleoaltimetry of Eocene
implemented in a reasonable amount of comput- gravels (Mulch et al., 2006) and the antiquity
ing time. of the Sierra Nevadan rain shadow (Poage and
In the context of landform evolution, solving Chamberlain, 2002), which both support the con-
the advection equation to model bedrock chan- clusion of a high (>2.2 km) Sierra Nevada in
nel evolution is signi¬cantly simpler than many early Eocene time. A high early Eocene Sierra
of the examples we have given here. In ¬‚uvial Nevada is not necessarily inconsistent with late
landscapes, slopes are ˜˜advected™™ upslope, so up- Cenozoic surface uplift if the range was sig-
wind differencing is simply a matter of calcu- ni¬cantly denuded in the intervening period.
lating the local channel gradient in the downs- Given the low erosion rates observed in the high
lope direction. The numerical diffusion of the Sierra today, however, these results are dif¬cult to
upwind-differencing method still applies, how- reconcile.
ever, so a Smolarkiewicz correction is needed for A large number of studies clearly indicate
precise results. When solving advection problems that signi¬cant late Cenozoic stream incision and
more generally, such as the transport of heat in rock uplift occurred throughout much of the
Sierra Nevada. Unruh (1991) documented 1.4—¦ of
a complex ¬‚ow ¬eld for long time periods, great
care must be taken in order to determine the ¬‚ux post-late-Miocene westward tilting of the Sierra
directions correctly in order to obtain stable, pre- Nevada based on the stratigraphy of the San
cise results. Joaquin Valley. Stock et al. (2004, 2005) utilized
cosmogenic isotope studies of cave sediments to
document a pulse of relatively high channel inci-
4.4 Modeling the sion rates (≈ 0.3 mm/yr) between 3 and 1.5 Ma in
¬‚uvial-geomorphic response of the South Fork Kings River and nearby drainages
of the southern Sierra Nevada. What is currently
the southern Sierra Nevada to
unknown is whether signi¬cant surface uplift ac-
uplift companied this late Cenozoic rock uplift. The re-
lationship between rock uplift and surface up-
In this section we consider a detailed case study lift remains uncertain for two principal reasons.
of the application of bedrock incision mod- First, surface uplift triggers knickpoint retreat
els to modeling the propagation of knickpoints along mainsteam rivers, but the time lag be-
and escarpments at mountain-belt scales. This tween incision at the range front and incision
application illustrates the use of the upwind- tens of kilometers upstream is not well con-
differencing method combined with the MFD strained. Thermochronologic data of Clark et al.
¬‚ow-routing method discussed in Chapter 3. This (2005), for example, provide a 32-Ma maximum
application also provides the motivation for an age for the onset of stream incision at their
in-depth comparison of the behavior of the sample localities, but range-wide surface uplift
stream-power and sediment-¬‚ux-driven bedrock could have occurred earlier. Second, rock uplift,
incision models. Codes for implementing the 2D stream-incision, and local surface uplift can all
bedrock drainage network model described in occur in the absence of range-wide surface up-
this section are given in Appendix 3. lift. As stream incision removes topographic loads
The uplands of the Sierra Nevada are a slowly from the crust, isostatic rebound raises slowly
eroding plateau perched 1.5 km above narrow eroding upland plateau remnants to elevations
river canyons. If we assume that a low-relief land- much higher than the original, regionally ex-
scape near sea level was uplifted recently and tensive surface. As such, no simple relationship

exists between the timing of local surface uplift channels are suf¬ciently overwhelmed with sed-
and range-wide surface uplift. iment to cause a reduction in erosion rates.
Cosmogenic nuclide measurements provide Eq. (4.25) also does not explicitly include the ef-
important constraints on rates of landscape evo- fects of grain size on abrasional ef¬ciency. There-
lution in the Sierra Nevada. Small et al. (1997) and fore, Eq. (4.25) is only an approximation to the
Stock et al. (2005) measured erosion rates of ap- considerable complexity of the saltation-abrasion
proximately 0.01 mm/yr on bedrock-exposed hill- process.
slopes of the Boreal Plateau. Riebe et al. (2000, Figure 4.4 compares the behavior of two nu-
2001) calculated basin-averaged erosion rates of merical landform evolution models incorporat-
≈ 0.02 mm/yr (largely independent of hillslope ing stream-power erosion (Figures 4.4a--4.4f) and
gradient) on the Boreal Plateau, increasing to sediment-¬‚ux-driven erosion (Figures 4.4g--4.4l)
≈ 0.2 mm/yr in steep watersheds draining directly for a vertically uplifted, low-relief plateau 200 km
in width. Uplift occurs at a constant rate U =
into mainstem river canyons. Stock et al. (2004,
2005) also measured mainstem incision rates 1 m/kyr for the ¬rst 1 Myr of each simulation.
The values of K w = 3 — 10’4 kyr’1 and K s =
along the South Fork Kings and nearby rivers
3 — 10’4 (m kyr)’1/2 used in these runs would
as they incised into the Chagoopa Plateau. These
authors found incision rates to be <0.07 mm/yr predict identical steady state conditions if up-
between 5 and 3 Ma, accelerating to 0.3 mm/yr lift was assumed to continue inde¬nitely. In the
between 3 and 1.5 Ma, and decreasing to model, hillslope erosion rates are prescribed to
≈ 0.02 mm/yr thereafter. These studies provide di- be 0.01 mm/yr and ¬‚exural-isostatic rebound is
rect constraints for parameters of hillslope and included. Isostatic rebound was estimated by as-
channel erosion in the numerical model. suming regional compensation over a prescribed
It is not currently known whether water or ¬‚exural wavelength and averaging the erosion
sediment acts as the primary erosional agent in rate over that wavelength. The ¬‚exural wave-
bedrock channels. The answer likely depends on length was assumed to be 200 km (the width of
the lithology of the eroding rock. Field studies, the model domain). The average erosion rate E av
for example, indicate that sediment abrasion is was computed for each time step and isostatic
likely to be most important in massive bedrock rock uplift U i rate was calculated using
lithologies like granite, while water can play a ρc
Ui = E av (4.26)
large role in plucking blocks out of jointed sedi-
mentary rocks (Whipple et al., 2000). The classic
where ρc and ρm are the densities of the crust
method for quantifying bedrock channel erosion,
and mantle, respectively.
the stream-power model, assumes that water is
First, consider the behavior of the stream-
the primary erosional agent in bedrock channels.
power model illustrated in Figures 4.4g--4.4l. Up-
Sklar and Dietrich (2001, 2004) proposed an
lift initiates a wave of bedrock incision in which
alternative ˜˜saltation--abrasion™™ model in which
channel knickpoints propagate rapidly through
bedrock channel erosion is controlled by the sed-
the drainage basin. In this model, knickpoints
iment ¬‚ux delivered from upstream. A simpli¬ed
reach the drainage headwaters after 25 Myr and
version of their model is obtained by replacing
the maximum elevation at that time is nearly
drainage area with sediment ¬‚ux in the stream-
3 km (Figure 4.4l). Following 25 Myr, the range
power model to obtain a sediment-¬‚ux-driven ero-
slowly erodes to its base level. In the sediment-
sion model:
¬‚ux-driven model (Figures 4.4a--4.4f), knickpoint
‚h 1 ‚h
= U ’ K s Q s2 (4.25) migration occurs at far slower rates due to the
‚t ‚x
lack of cutting tools supplied from the slowly-
and introducing a new coef¬cient of erodibil- eroding upland plateau. In this model, knick-
ity K s . This version of the sediment-¬‚ux-driven points require 60 Myr to reach the drainage head-
model incorporates the ˜˜tools™™ effect of saltat- waters and the maximum elevation at that time
ing bedload but not the ˜˜cover™™ that occurs when is nearly 4 km. Incision of lowland channels

sediment-flux-driven model
(a) 4
50 km
h 60 Myr
40 Myr
20 Myr

80 Myr
t = 20 Myr t = 40 Myr 40 x (km)
0 80
(c) 4 (f)
h max

t = 60 Myr t = 80 Myr 40 t (Myr)
0 80
stream-power model
(g) 4 (k)
10 Myr
20 Myr
30 Myr
40 Myr
t = 10 Myr t = 20 Myr 40 x (km)
0 80
(i) 4 (l)
h max

t = 30 Myr t = 40 Myr 20 t (Myr)
0 40
h (km)
0.0 1.0 2.0 3.0 4.0

Fig 4.4 Model results comparing the (a)“(f) sediment-¬‚ux-driven and (g)“(l) stream-power models following 1-km uniform block
uplift of an idealized mountain range. In the stream-power model, knickpoints rapidly propagate into the upland surface, limiting
the peak elevation to ≈ 3 km at 25 Myr following uplift. In the sediment-¬‚ux-driven model, knickpoints propagate slowly due to
the limited cutting tools supplied by the slowly eroding upland plateau, which raises plateau remnants to nearly 4 km 60 Myr
following uplift. Modi¬ed from Pelletier et al. (2007c). Reproduced with permission of Elsevier Limited.

drives regional rock uplift in both models. In has a mean elevation of approximately 0.5 km, a
the sediment-¬‚ux-driven model, however, the low maximum elevation of 1 km, and local relief sim-
rates of channel incision on the upland plateau ilar to that of the modern Boreal Plateau. All pix-
cause surface uplift of the relict surface rem- els draining to areas north, east, and south were
nants to higher elevations than in the stream- ˜˜masked™™ from the model domain so that only
power model. The magnitude of this ¬‚exural- the evolution of west-draining rivers was consid-
isostatic rebound effect is equal to the inverse of ered. The drainage divide at the crest was treated
the relative density contrast between the crust as a free boundary.
and the mantle, ρc /(ρm ’ ρc ), which is between Hillslopes and channels in the model are dis-
4 and 5 for typical continental crust (Turcotte tinguished on the basis of a stream-power thresh-
and Schubert, 2002). Typical ¬‚exural wavelengths old (Tucker and Bras, 1998; Montgomery and
in mountain belts are 100--200 km. Therefore, if Dietrich, 1999). If the product of slope and the
the extents of individual surface remnants are square root of contributing area is greater than
equal to or smaller than this wavelength, canyon a prescribed threshold given by 1/ X , where X
cutting nearby will cause substantial local sur- is the drainage density, the pixel is de¬ned to
face uplift of these remnants up to four times be a channel, otherwise it is a hillslope. The
value of X was chosen to be 0.005 km’1 on the
higher than that of the original, regionally ex-
tensive plateau, even as the mean elevation of basis of observed drainage densities on the up-
the range is decreasing. In the limit where the land Boreal Plateau observed on a 30-m resolu-
upland erosion rate goes to zero, a vertical es- tion DEM. Hillslope erosion in the model occurs
carpment would form in the sediment-¬‚ux-driven at a prescribed rate E h independent of hillslope
model that would retreat laterally at a rate equal gradient. Diffusion or nonlinear-diffusion mod-
to the rate of bedrock weathering (i.e. typically < els are commonly used to model hillslope evolu-
m/kyr). At such rates, knickpoint retreat of tion in numerical models, but cosmogenic ero-
100 km requires more than 100 Myr. Appendix 3 sion rates measured on the Boreal Plateau (Riebe
provides codes for implementing the sediment- et al., 2000) support the slope-independent ero-
¬‚ux-driven model with ¬‚exural-isostatic rebound. sion model used here. The value of E h is con-
Given the predominantly granitic lithology of strained by cosmogenic erosion rates on bare
the Sierra Nevada, the results of the sediment- bedrock surfaces to be approximately 0.01 mm/yr
¬‚ux-driven model are most relevant to this ap- (Small et al., 1997; Stock et al., 2004, 2005). This
plication. Here we assume an initially low-relief, value applies only to the upland low-relief sur-
low-elevation landscape similar to the Boreal face, but in practice hillslope pixels in the model
Plateau of the modern high sierra. However, no only occur on the upland landscape because of
explicit assumption is made about the timing of the large gradients and/or drainage areas that
initial uplift of the range. The initial topography develop in the incised river canyons below the
input to the model respects the modern drainage upland landscape. Hillslopes play an important
architecture of the southern Sierra Nevada but role in the sediment-¬‚ux-driven model by supply-
replaces the modern stepped topography with ing sediment to propagate knickpoints, but the
a uniformly low-relief, low-elevation landscape model results are not sensitive to the particular
(Figure 4.5a). To construct the initial topography, value of X as long as it is within a reasonable
range of 0.002 km’1 to 0.01 km’1 .
downstream ¬‚ow directions were ¬rst computed
Bedrock channel erosion occurs in all chan-
using a 500-m resolution DEM of the southern
nel pixels during each time step of the model ac-
Sierra Nevada. Second, a constant base-level ele-
cording to Eqs. (4.4) or (4.25) implemented with
vation of 200 m a.s.l. was identi¬ed based on the
upwind differencing. Active uplift in the model
modern range-front elevation. Starting from pix-
was assumed to occur as uniform, vertical, block
els in the DEM with this base-level elevation, DEM
uplift of U = 1 m/kyr during each uplift pulse.
pixels with higher elevations were then ˜˜reset,™™
Slope--area relationships in steady-state bedrock
enforcing a uniform hillslope gradient of 5%. The
rivers indicate that the ratio m/n = 0.5 is typical.
initial model topography produced in this way

sediment-flux (b) (d)
50 km

crest 2 3


t = 60 Myr
t = 20 Myr t = 40 Myr

(f) (h)

2 3
t = 0 Myr

t = 35 Myr
t = 10 Myr t = 20 Myr
stream power
h (km)
0.0 1.0 2.0 3.0 4.0

Fig 4.5 Maps of best-¬t model results for the
channels using a bifurcation-routing algorithm
sediment-¬‚ux-driven model (b)“(d) and stream-power model
(f)“(h) starting from the low-relief, low-elevation surface that partitions incoming ¬‚uxes of sediment or
illustrated in (a). The actual modern topography is shown in drainage area between each of the downstream
(e). The best-¬t uplift history for the sediment-¬‚ux-driven
neighboring pixels, weighted by slope. Time steps
model occurs for a 1-km pulse of uplift starting at 60 Ma
in the model are variable (i.e. small time steps are
(t = 0) and a 0.5 pulse starting at 10 Ma (t = 50 yr out of a
needed when stream gradients are high following
total 60 Myr). The best-¬t uplift history for the stream-power
uplift in order to ensure numerical stability) but
model occurs for a 1-km pulse of uplift at 35 Ma and a 0.5 km
range from about 1 to 10 kyr.
pulse at 7 Ma. Also shown are transect locations and stream
Geomorphic and stratigraphic observations
location identi¬ers corresponding to the erosion-rate data in
suggest that isostatic rock uplift takes place by
Figure 4.6. For color version, see plate section. Modi¬ed from
Pelletier et al. (2007c). Reproduced with permission of westward tilting along a hinge line located west
Elsevier Limited. of the range front. Sediment loading of the San
Joaquin Valley and the presence of a major high-
angle fault along the eastern escarpment suggest
that isostatic rebound is a maximum at the range
This m/n ratio implicitly includes the effect of in-
crest and decreases uniformly towards the hinge
creasing channel width with increasing drainage
line (Chase and Wallace, 1988). Therefore, we as-
area. In this model experiment, n values of 1
sume here that isostatic uplift has a mean value
and 2 were used, scaling m accordingly to main-
tain a constant m/n = 0.5 ratio. Sediment ¬‚ux equal to Eq. (4.26) but is spatially distributed us-
ing a linear function of distance from the range
is computed in the sediment-¬‚ux-driven model
using erosion rates computed during the previ-
ous time step. Sediment ¬‚ux and drainage area ρc
L xc
Ui = 1’ E av (4.27)
are both routed downstream on hillslopes and in < xc > ρm ’ ρc

uplift was imposed from time t = 0 to 1 Myr (out
where xc is the distance from the range crest
along an east--west transect, L is the distance be- of a total duration of 60 Myr) and a second 0.5-km
tween the crest and the hinge line, and < xc > pulse of late Cenozoic surface uplift was imposed
from t = 50 to 50.5 Myr. Model experiments with
is the average distance from the range crest. For
the results of this section I assumed L = 100 km n = 2 were also performed and yielded a similar
(Unruh, 1991) and ρc (ρm ’ ρc ) = 5. best-¬t uplift history, but the knickpoint shapes
Calibration of K w and K s was performed for in that model were a poor match to the observed
each uplift history by forward modeling the re- knickpoint shapes along the North Fork Kern
sponse to a late Cenozoic pulse of uplift and River. Figures 4.5b--4.5d present color maps of
the model topography at t = 20, 40, and 60 Myr.
adjusting the values of K w and K s to match
the maximum incision rate of 0.3 mm/yr mea- Figure 4.5e illustrates the modern topography
sured by Stock et al. (2004, 2005) along the South of the southern Sierra Nevada with the same
Fork Kings River. For each uplift history consid- color scale for comparison. Figure 4.6a plots the
ered, the values of K s and K w were varied by channel incision rate at three points along the
trial and error until the maximum value of late Kings and South Fork Kings Rivers as a func-
Cenozoic incision matched the measured value tion of time. Location 2 is of particular signi¬-
of 0.3 mm/yr on the South Fork Kings River to cance since that point coincides with the South
within 5%. Late Cenozoic incision rates in the Fork Kings River where Stock et al. (2004, 2005)
model increase monotonically with K s and K w , measured maximum incision rates of 0.3 mm/yr.
so trial values of K s and K w were raised (or low- Each point along the Kings River experiences two
ered) if predicted incision rates were too low (or pulses of incision as knickpoints propagate up-
high). All uplift histories considered in the model stream as topographic waves. Figure 4.6a indi-
experiments assumed two pulses of range-wide cates that incision rates decrease with increas-
surface uplift with the late Cenozoic phase con- ing distance upstream and that knickpoint prop-
strained to occur on or after 10 Ma and on or be- agation is slower following the ¬rst phase of
fore 3 Ma to be consistent with ¬rm constraints uplift compared to the second phase of up-
imposed by previous studies. The magnitude of lift. Both of these behaviors can be associated
uplift during each pulse was tentatively con- with the ˜˜tools™™ effect of the sediment-¬‚ux-driven
strained using the offset of the two major river model. First, larger drainages have more cutting
knickpoints of the North Fork Kern River, which tools to wear down their beds, thereby increas-
suggests an initial pulse of ≈ 1 km surface uplift ing incision rates and knickpoint migration rates.
corresponding to the upstream knickpoint, fol- Second, knickpoint propagation occurs slowly fol-
lowed by a smaller (≈ 0.5 km) pulse correspond- lowing the ¬rst pulse of uplift because of the
ing to the downstream knickpoint. Forward mod- lack of cutting tools supplied by the low-relief
eling indicates that knickpoint relief during the upland plateau. Figure 4.5b, in particular, illus-
model remains nearly equal to the magnitude trates that canyon cutting has affected only a very
of the range-wide surface uplift that triggered small percentage of the landscape even as late as
knickpoint creation except for late stages of the 20 Myr following the initial uplift. In contrast,
model when the uplift of isolated plateau rem- knickpoint and escarpment retreat following the
nants can greatly outpace stream incision. Best- second phase of uplift is enhanced by the up-
¬t uplift histories were determined by visually stream sediment supply associated with the ¬rst
comparing the model-predicted topography with phase of uplift. As a result, the second knickpoint
those of the modern Sierra Nevada, including the travels almost as far upstream as the ¬rst knick-
extents of the Chagoopa and Boreal Plateaux and point in less than 20% of the time.
the major river knickpoints. The results of Figures 4.5b--4.5d and Figures
Figures 4.5b--4.5d illustrate results of the best- 4.6a--4.6b also provide several consistency checks
¬t sediment-¬‚ux-driven model experiment, with that aid in con¬dence building. First, the maxi-
n = 1 and K s = 2.4 — 10’4 (m kyr)’1/2 . In this ex- mum incision rate observed by Stock et al. (2004,
periment, a 1-km pulse of late Cretaceous surface 2005) from 3 to 1.5 Ma along the South Fork Kings

0.6 0.03 (a) (b)
E (mm/yr)
4 0.2
(mm/yr) 0.02
E, U
max h
0.3 0 50 100
calibrated with x (km)
1 2 0.1
mean h
Stock et al., Riebe et al.
3 calibrated with
mean E
Stock et al. (2004)

mean U
0.0 0 0.0
20 t (Myr) 40 60 20 t (Myr) 40 60
0 0
0.9 0.03 (c) (d)
E (mm/yr) 4 0.2
(mm/yr) 1 0.02
E, U


. 5
( 14)