ńņš. 5 |

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)

3.5 CONTAMINANT TRANSPORT IN CHANNEL BED SEDIMENTS 79

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

80 FLOW ROUTING

116.50

116.51Ā° W

0.8

ex

concentration

tephra

measured data

1505.03

36.70 N

BG

1505.02 3D model

1505.01 1505.04 1505.05 1505.06

prediction

1405.04 1405.01

1405.02 tephra sheet

1405.05

east channel 0.1

1405.06 1405.03

BG BG 1505.08

major tributary

1405.07

0.05

1505.07 junctions

1405.08

east channel

1405.09

3 5

4

2

1

0

Lathrop Wells

x (km)

36.69 1405.10 lava flows

Volcanic Cone

0.8

west channel

1405.11

measured data

concentration

1505.09 3D model

1405.12

prediction

tephra

1505.11

1505.10

36.68

1505.12 0.3

N west channel

1505.13

2

1

0

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-

(2008).

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

3.5 CONTAMINANT TRANSPORT IN CHANNEL BED SEDIMENTS 81

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

SD091405.03

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

deposits

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

82 FLOW ROUTING

(b)

101 km

100

116.52 W 116.50 10

116.48 10

(a)

N1 km

36.74Ā° N

3%

9%

(d)

36.72

tephra sheet

36.70 slope*area1/2

(c)

36.68

73%

LathropWells cone

relativete phra

concentration

36.66

0.01 0.1 1.0

0.001

topography and source mask

relative

scour/mixing depth 0.1 0.3 1.0

900 1000 m

800

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-

3.5 CONTAMINANT TRANSPORT IN CHANNEL BED SEDIMENTS 83

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-

84 FLOW ROUTING

U slope*area1/2

tephra thickness slope

100 km

10

10

10

0.01 0.1 1m 0.01 0.03 0.1 0.3 1.0

(b)

(a) (c)

vent

location

southwesterly

winds

N 3 km

slope*area1/2

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

0.29

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

EXERCISES 85

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

0.1 0.2

tephra concentration at the outlet is the topog-

(b)

(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

vent

conditions transport tephra north and east from

location

the repository location to areas with a relatively

high density of tephra-transporting steep slopes

and active channels within the Fortymile Wash

drainage basin.

Exercises

3.1 Download a DEM of any area dominated by ļ¬‚uvial

processes. Implement the fillinpitsandflats

outlet

concentration routine (Appendix 2) to ļ¬lter out pits and ļ¬‚ats

southwesterly = 1.97% in the DEM. Implement the MFD ļ¬‚ow-routing

winds

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.

86 FLOW ROUTING

southerly westerly northerly

winds winds winds

(b)

(a) (c)

concentration concentration concentration

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

N 3 km

tephra concentration

100%

10%

1%

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-

Limited.

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

88 THE ADVECTION/WAVE EQUATION

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

and

problems. Rather than using the local shape of

dt

= b(x, t)

the function to evolve the system forward in (4.9)

ds

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

dh

= + = a(x, t) + b(x, t)

cases, these trajectories are described by simple ds ā‚ x ds ā‚t ā‚x ā‚t

ds

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

dh

K w A m varies through time as individual grid + c(x, t)h = 0 (4.11)

ds

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

4.2 ANALYTIC METHODS 89

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

90 THE ADVECTION/WAVE EQUATION

This equation assumes that the square root of dis- 1.5

charge increases linearly with distance from the

4

divide, which is appropriate for a humid drainage

3

1.0

system (i.e. Q ā A) where drainage area A in- t/c = 0

2

1

h

creases as the square of distance from the divide

(i.e. an approximately semi-circular basin). Equa-

0.5

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-

0.0

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

8

0 2 6

x

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.

2

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

becomes

Methods for solving the advection equation can

h i+1 ā’ h iā’1

h in+1 ā’ h in n n

=c

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

4.3 NUMERICAL METHODS 91

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

=c

(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

equation

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

solution.

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-

1n

h in ā’ (h + h iā’1 ) nitude of the numerical diffusion, we can make

n

(4.15)

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,

comes

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

n+1/2

h i+1/2 = (h i+1 + h in ) ā’ (F i+1 ā’ F in )

n

(4.19)

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)

x

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.

92 THE ADVECTION/WAVE EQUATION

4.3.4 Upwind-differencing method 1.5

(a)

Another simple but powerful numerical method

t/c = 0

for the advection equation is called upwind dif-

1.0

ferencing. Upwind differencing involves calculat-

ing the slope along the direction of transport:

ā§ hn ā’ hn 8 6 4 2

0.5

āŖ i+1

āØ if c in > 0

i

h in+1 ā’ h in x

= ci

n

(4.21)

āŖ h in ā’ h iā’1

n

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.

1.5

at i ā’ 1) and the value to the right (i.e. at i + 1) (b)

of the grid point being updated. This approach

1.0

is prone to instability because it does not make

use of the value at the grid point i itself. As such,

h

the difference h iā’1 ā’ h i+1 can be small even if the

n n

8 6 4 2

0.5

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

0.0

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

n

4 10

8

0 2 6

x

or h in ā’ h i+1 ). Which direction to choose? If the

n

ļ¬‚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

n

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-

where

ing with and without Smolarkiewicz correction.

ā‚h

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

4.4 FLUVIAL-GEOMORPHIC RESPONSE OF THE SOUTHERN SIERRA NEVADA TO UPLIFT 93

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

94 THE ADVECTION/WAVE EQUATION

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-

Ļm

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

4.4 FLUVIAL-GEOMORPHIC RESPONSE OF THE SOUTHERN SIERRA NEVADA TO UPLIFT 95

sediment-flux-driven model

(b)

(a) 4

50 km

(e)

h 60 Myr

(km)

40 Myr

2

20 Myr

80 Myr

0

t = 20 Myr t = 40 Myr 40 x (km)

0 80

(d)

(c) 4 (f)

h max

(km)

2

mean

0

t = 60 Myr t = 80 Myr 40 t (Myr)

0 80

stream-power model

(h)

(g) 4 (k)

h

(km)

2

10 Myr

20 Myr

30 Myr

40 Myr

0

t = 10 Myr t = 20 Myr 40 x (km)

0 80

(j)

(i) 4 (l)

h max

(km)

2

mean

0

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.

96 THE ADVECTION/WAVE EQUATION

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

4.4 FLUVIAL-GEOMORPHIC RESPONSE OF THE SOUTHERN SIERRA NEVADA TO UPLIFT 97

sediment-flux (b) (d)

(c)

50 km

driven

1

crest 2 3

(e)

(a)

t = 60 Myr

t = 20 Myr t = 40 Myr

(f) (h)

(g)

1

observed

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

crest:

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

L

98 THE ADVECTION/WAVE EQUATION

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

4.4 FLUVIAL-GEOMORPHIC RESPONSE OF THE SOUTHERN SIERRA NEVADA TO UPLIFT 99

0.6 0.03 (a) (b)

E

E (mm/yr)

4 0.2

(mm/yr) 0.02

E, U

h

(mm/yr)

(km)

0.01

max h

0.3

0.3 0 50 100

calibrated with x (km)

1 2 0.1

mean h

Stock et al., Riebe et al.

2

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

E (mm/yr) 4 0.2

(mm/yr) 1 0.02

E, U

h

0.6

ńņš. 5 |