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