<< . .

. 18
( : 51)

. . >>

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

<< . .

. 18
( : 51)

. . >>