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