<< . .

. 17
( : 51)

. . >>

ample in Figure 3.14. Outlet tephra concentra- ¬‚ood depth. This will require that you identify a
tions for these rotated plumes vary from a mini- pixel or set of pixels and prescribe an input ¬‚ow
mum of less than 0.1% to a maximum of approxi- depth.

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

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

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

The advection/wave equation

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

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

10 2
(a) (c)

t/c t/c

5 1

0 0
1.5 1.5
(b) (d)

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

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

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

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

<< . .

. 17
( : 51)

. . >>