. 11
( 14)


Dust accumulation on desert alluvial fan ter-
head entrenchment is the surface ˜˜age™™ and it
races can be affected by changes in either depo-
corresponds to the time interval of eolian dust
sition or erosion. Previous studies on dust accu-
accumulation. Figure 8.7 plots silt accumulation
mulation of alluvial fans terraces have generally
rates versus time interval on logarithmic scales
argued that climatically controlled dust deposi-
for these locations. Error bars represent time
tion plays the predominant role in controlling

2 1 10
(a) (b) (c)


1 0.5

bounded random walk
0 0 10
101 105
0 5 15 103 104
4 102
2 T (k.y.)
T (k.y.)
T (yr)
(d) (e) (f)


log R
random walk

0 10
10 log T
0 101 105
103 104
100 102
T (k.y.) T (yr)

intervals of 103 --104 yr, consistent with an early
Fig 8.8 (a) Plot of dust accumulation versus time for two
Holocene dust pulse.
runs of the bounded random walk model. (b) Inset of (a)
illustrating the scale-invariance of the bounded random walk In order to interpret the temporal scaling of
model by rescaling time by a factor of 4 and accumulation by dust accumulation rates with time interval, con-
a factor of 2. (c) Plot of accumulation rate versus time sider two end-member models. The ¬rst model
interval illustrating power-law scaling behavior. (d) Plot of
assumes that the magnitude of dust deposi-
dust accumulation versus time in the cyclical-climate model.
tion and erosion is constant through time (with
(e) Plot of accumulation rate versus time interval in the
a rate equal to 0.01 g/cm2 /yr in this example)
cyclical-climate model, illustrating an inverted S-shaped curve.
and that erosion or deposition takes place with
(f) Schematic diagram summarizing the results of the two
equal probability during each time step (1 yr).
models for accumulation rate versus time interval. Modi¬ed
Figure 8.8a illustrates two runs of this ˜˜bounded
from Pelletier (2007a).
random walk™™ model. The walk is bounded be-
cause the net accumulation is always positive
interval (i.e., surface age) uncertainties. The ¬rst
(in an unbounded random walk, negative values
four plots in Figure 8.7 exhibit a power-law
are also allowed). At the longest time interval
trend (i.e., a straight line on log--log scales) given
(104 yr), both model runs have accumulated ap-
proximately 1.5 g/cm2 of dust (equal to a thick-
ness of 1.5 cm if a density of 1 g/cm3 is assumed).
R ∝ T ’± (8.9)
On shorter time intervals (103 and 102 yr), net
accumulations are smaller than the 1.5 value,
where R is the accumulation rate, T is the time
interval of accumulation, and ± is an exponent but they are not linearly proportional to the
close to 0.5. Values of ± for each data set were time interval. For example, net accumulation is
0.3--0.5 g/cm2 over 103 yr for the two model runs
obtained by a linear ¬t of the logarithm of accu-
shown, a much larger change per unit time than
mulation rate to the logarithm of time interval,
the 1.5 g/cm2 observed over 104 yr. Figure 8.8b il-
weighted by the age uncertainty for each point.
lustrates the temporal scaling of accumulation
Accumulation rates plotted in Figures 8.7e and
rates versus time interval in this model. The re-
8.7f (lower and upper fans near Silver Lake Playa)
sults in Figure 8.8b illustrate that accumulation
show only a weak dependence on time interval,
rates decrease, on average, at a rate proportional
however. Figure 8.7e, for example, has a nearly
to one over the square root of the time inter-
constant rate across time intervals, while Figure
val, similar to that in Figures 8.7a--8.7d. This
8.7f shows a peak in accumulation rates at time

power-law scaling behavior can also be derived The scaling behavior of accumulation rates
theoretically. The average distance from zero in a shown in Figures 8.7a--8.7d has important impli-
bounded random walk, < h >, increases with the cations for the completeness of eolian deposits.
square root of the time interval, T (van Kampen, Sedimentary deposits are said to be ˜˜incomplete™™
2001): if unconformities or hiatuses with a broad dis-
tribution of time intervals are present (Sadler,
< h >= aT 1/2 (8.10)
1981). The accumulation rate curve can be used
where the brackets denote the value that would to quantify the completeness of a given deposit
be obtained by averaging the result of many dif- using the ratio of the overall accumulation rate
ferent model runs and a is the magnitude of ero- to the average rate at a particular time inter-
sion or deposition during each time step. As ap- val T (Sadler and Strauss, 1990). In the bounded
plied to eolian dust accumulation, the average random walk model, for example, a complete
distance < h > represents the total accumulation stratigraphic section (i.e., one with no hiatuses)
of eolian dust (in g/cm2 or cm) on the surface, of 10 000 yr duration would have a total accu-
mulation of 100 g/cm2 based on an annual accu-
including the effects of episodic erosion. The av-
mulation rate of 0.01 g/cm2 . The bounded ran-
erage accumulation rate is obtained by dividing
Eq. (8.10) by T to give dom walk model, in contrast, has an average
accumulation of 1 g/cm2 over 10 000 yr based
< R >= aT ’1/2 (8.11)
on Eq. (8.11). Therefore, a 10 000 yr sequence is
In the second end-member model, dust depo- only 1% complete (i.e., 99% of the mass that was
sition is assumed to have a high value during once deposited has been eroded). A power-law
interglacial periods (0.01 g/cm2 /yr), a low value accumulation-rate curve (e.g., Eq. (8.10)) speci¬-
during glacial periods (0.002 g/cm2 /yr) and no ero- cally implies a power-law or fractal distribution
sion. Figure 8.8d illustrates the accumulation of hiatuses, with few hiatuses of very long dura-
rate versus time interval for this model. The tion and many of short duration.
plot follows an inversed S-shaped curve shown
in Figure 8.8e. This relationship is characteris-
8.6 Stochastic trees and
tic of accumulation rate curves for all periodic
diffusion-limited aggregation
functions (Sadler and Strauss, 1990) and it is
markedly different from the observed scaling in
Figures 8.7a--8.7d.
Stochastic models have featured prominently in
The results of Figure 8.7 suggest that eolian
our understanding of drainage networks. Hor-
dust accumulation on alluvial fans located close
ton (1945) was among the ¬rst to quantitatively
to playa sources (Figures 8.7e--8.7f) is fundamen-
study the geometry of drainage networks using
tally different from that on fans located far from
Strahler™s classi¬cation scheme. In this scheme,
playa sources (Figures 8.7a--8.7d). The fan study
channels with no upstream tributaries are classi-
sites close to Silver Lake Playa exhibit a linear
¬ed as order 1. When two streams of like order
accumulation trend on the lower fan site and
join together they create a stream of the next
a cyclical-climate trend on the upper fan site
highest order. Horton showed, remarkably, that
while the remaining fan sites exhibit accumu-
the ratios of the number, length, and area of
lation trends consistent with the bounded ran-
streams within a drainage basin nearly always
dom walk model. These results suggest that the
follow common values independent of Strahler
bounded random walk model is representative of
order, i.e.
fans located far (i.e., greater than a few km) from
Ni A i+1 L i+1
playas, and that the classic linear-accumulation ≈ 4, ≈ 4, ≈2 (8.12)
N i+1 Ai Li
or cyclical-climate models are representative of
accumulation on fans proximal to playa sources where i is the Strahler order. In other words the
where dust deposition is rapid enough to over- ratio of the number of order 1 streams to order
whelm the effects of episodic disturbances. 2 streams is nearly the same as the number of

(a) (b) (c)


Fig 8.10 TDCNs for Shreve magnitude (a) 1, (b) 2, (c) 3,
and (d) 4.

to transform this structure into that of Figure
8.9a without lifting one of the strings off of the
table. What is nice about the TDCN framework
Fig 8.9 Illustration of the concept of Topologically Distinct
is that the entire set of all possible TDCNs can
Channel Networks (TDCNs). (a) A hypothetical channel
network of magnitude 4. (b) Two channel networks that are be identi¬ed for a given Strahler order or Shreve
topologically identical to that in (a) because they can be magnitude. In the Shreve classi¬cation scheme,
transformed into (a) by moving stream segments around
channels without tributaries are de¬ned to be
within the plane. (c) A network that is topologically distinct
magnitude 1 and when two streams join they
from (a).
form a stream with magnitude equal to the sum
of the two tributaries. Figure 8.10, for example,
shows all of the TDCNs up to Shreve magnitude
order 2 streams to order 3 streams (and so on).
4. The number of TCDNs of Shreve magnitude M
The fact that drainage basins have common ra-
is given by
tios of number, length, and area across Strahler
orders (i.e. spatial scales) is the basis for fractal (2N ’ 2)!
N= (8.13)
models of drainage network geometry.
M !(M ’ 1)!
Following the establishment of Horton™s Laws,
where M ! is equal to M — (M ’ 1)— (M ’ 2) . . . 2 — 1.
stochastic models of drainage network geometry
were proposed. One such framework is Topologi- An analysis of tree structures randomly se-
cally Distinct Channel Networks (TDCNs). TDCNs lected from all possible TDCNs indicates that
are tree structures that cannot be topologically these structures obey Horton™s Laws. Since TD-
transformed into one another within their plane. CNs are a general mathematical model for tree
This concept is best illustrated with an example. structures and do not involve any geomorphic
In Figure 8.9a, a possible tree structure of order processes, the natural conclusion is that Horton™s
2 is shown. Now, imagine that this network is Laws are not really a consequence of drainage
made of a series of strings knotted together at network evolution, but instead are a general
junctions, and that the strings are laying on a feature of tree structures (Kirchner, 1993). Fur-
table. A structure is topologically indistinct from ther studies of drainage network geometry have
another tree structure if the strings can be rear- discovered other kinds of statistical relation-
ranged without lifting any of the strings off of ships commonly found in real drainage networks.
the table. Figure 8.9b illustrates two possible tree Tokunaga (1984) and Peckham (1996) expanded
structures that are topologically identical to the the Horton analysis to quantify not just the num-
structure in Figure 8.9a. Figure 8.9c, in contrast, ber of streams of a given order, but the num-
is topologically distinct because there is no way ber of ˜˜side branches™™ of a given order i that

feed into streams of order j. This approach recog-
nizes the fact that small streams that join with Seed cells

small streams in the headwaters of the Missis-
Accreted cells
sippi River basin in Minnesota are fundamentally
different from small streams that join with the
Newly added cell
mainstem Mississippi River in Louisiana.
Diffusion-Limited Aggregation (DLA) provides Random walk
one simpli¬ed stochastic model of drainage net-
work evolution that matches the observed Toku- Prohibited sites

naga statistics of real drainage networks. In the
Other allowed sites
DLA model for drainage networks (Masek and for accretion
Turcotte, 1993), incipient channels are de¬ned on
Fig 8.11 Illustration of the mechanism for network growth
a rectangular grid using ˜˜seed™™ pixels along one
in the DLA model. A particle is randomly introduced to an
boundary (Figure 8.11). Then, a particle is intro-
unoccupied cell. The particle undergoes a random walk until
duced into an unoccupied cell chosen randomly
a cell is encountered with one (and only one) of the four
from all unoccupied cells on the grid. The par-
nearest neighbors occupied (hatched cell). The new cell is
ticle then undergoes a random walk until it en- accreted to the drainage network. From Masek and Turcotte
counters a cell with one of the four nearest neigh- (1993). Reproduced with permission of Elsevier Limited.
bors occupied with a channel pixel. The new cell
is then accreted to the network at this position
(Figure 8.11). In Dunne™s model of drainage network evolution,
Clearly, the DLA model is a very abstract headward-growing channels cause groundwater
model of the processes involved in the evolution to be de¬‚ected towards the channel head, thereby
of real drainage networks. So, how are we to inter- increasing the headward growth rate of ˜˜master™™
pret the fact that the model gives rise to drainage streams at the expense of smaller streams nearby.
networks that look realistic and also match some The DLA model works in a broadly similar fash-
of the statistical features of real drainage net- ion. Random walkers dropped randomly on the
works (Figure 8.12)? One explanation is that surface will tend to accrete new channels close
the random walk process mimics the particle to the tips of existing channels because these
pathways taken by groundwater ¬‚ows that drive channels extend out farther into the grid. The
headward growth of channels by spring sapping. behavior of random walkers on the grid acts as a

(b) Fig 8.12 (a) Drainage network of
the Volfe and Bell Canyons, San
Gabriel Mountains, near Glendora,
California, obtained from ¬eld
mapping. (b) Illustration of one
realization of the DLA model. From
Masek and Turcotte (1993).
Reproduced with permission of
Elsevier Limited.

0.5 km

simpli¬ed model for the groundwater ¬‚ow path- dust emissions. Our model aims to quantify the
ways in Dunne™s model. range of water-table depths over which this tran-
sition from a wet to dry playa takes place. An
important assumption in this analysis is that pre-
8.7 Estimating total ¬‚ux based on a cipitation events are rare enough that the surface
is near its steady-state moisture value most of the
statistical distribution of events:
dust emission from playas The vertical transport of moisture in an un-
saturated soil is governed by Richards™ equation:
The relationship between the hydrological state ‚ψ ‚ ‚ψ
= ’1
Kψ (8.14)
of a playa and its dust-emitting potential is a nice
‚t ‚z ‚z
example of the coupling of hydrology and geo-
where ψ is the suction, t is time, z is height above
morphology. In this section we consider a model
the water table, and K ψ is the hydraulic conduc-
for the long-term average dust emission from
tivity. For steady upward ¬‚ow, Eq. (8.14) can be
playas. Speci¬cally, the goal of the model is to cal-
written as
culate the long-term average dust emission rate
‚ ‚ψ
and the role of the water table depth in control-
E= ’1
Kψ (8.15)
‚z ‚z
ling dust emission rate. This example illustrates
the use of probability distributions to calculate where E is the steady-state evaporation rate.
long-term geomorphic rates. Gardner (1958) proposed the following relation-
In nature, dust storms often originate when ship between K ψ and ψ to solve Eq. (8.15):
coarse sand from the playa margin enters into
Kψ =
saltation. As sand blows across the playa sur- (8.16)
ψn +b
face, the crust is disturbed, releasing both sand
where a and b are empirical constants for each
and ¬ne-grained sediments (silt and clay). Sand
soil. Gardner (1958) provided analytical solutions
from the playa surface becomes part of a self-
for suction pro¬les at several values of n. Most
sustaining saltation cloud that drives dust pro-
soils have n values between 2 and 3 (Gardner
duction from the surface. The model we consider
and Fireman, 1958). Here we assume n = 2, the
in this section does not resolve different compo-
most appropriate value for ¬ne-grained playa sed-
nents of that playa. Rather, it assumes that the
iments according to Gardner (1958). The solution
¬‚ux of blowing sand along the playa margin can
to Eq. (8.15) for n = 2 and with boundary condi-
be used as a proxy for dust emission from the
tion ψ = 0 at z = 0 is (Gardner, 1958):
playa. This assumption is supported by data show-
ing a strong correlation between horizontal salta-
E /a
z= ψ
tion ¬‚ux and the vertical dust ¬‚ux (Gillette et al., (8.17)
(E /a)b + 1
2003). However, the model is not applicable to a a
playas with no upwind sand available for trans-
The value of E is determined by the maximum
soil moisture ¬‚ux within the pro¬le, given by
Soil moisture affects dust emission from
Gardner (1958) and Warrick (1988) as:
playas through its effects on soil cohesion. In or-
E 2
der to model soil moisture at the surface of a
= (8.18)
a 2d
playa, we consider the physics of steady-state cap-
illary rise from a water table. For a shallow wa- where d is the water table depth. Substituting Eq.
ter table, capillary rise leads to a moist surface (8.18) into Eq. (8.17) and solving for the suction at
z = d gives
that suppresses dust emissions under all but the
most extreme wind conditions. For a deep wa- 2
π π
2d 2 2
ter table, the dry surface is at or near its max- ψz=d = b+ b+1
π 2 2d
imum dust-producing potential and changes in
subsurface moisture have little or no effect on (8.19)

Van Genuchten (1980) predicted the surface soil velocities (Takle and Brown, 1978; Bowden et al.,
moisture from the suction to be 1983), was used for the probability distribution
θz=d ’ θr of friction velocities:
’ »+1
= 1 + (±ψz=d )1+» (8.20)
θs ’ θr γ
fw (u— ) = γβ ’γ uγ ’1 e’(u— /β) (8.24)

where θr is the residual soil moisture, θs is the
where γ and β are parameters ¬t to measured
saturated soil moisture, ± is the inverse of the
bubbling pressure, and » is the pore-size distri- wind data. The long-term average saltation ¬‚ux
is obtained by integrating Eq. (8.22):
bution parameter. Equations (8.19) and (8.20) pro-
vide an analytical solution for the surface soil ∞
< q >= du— fw (u— )u— (u2 ’ u2 ) (8.25)
moisture as a function of water-table depth and — —t
g u—t
soil-hydrologic parameters.
Substituting Eq. (8.24) into Eq. (8.25) yields a
The second model component relates the sur-
closed-form solution for the long-term average
face soil moisture to the threshold friction veloc-
saltation ¬‚ux:
ity. Chepil (1956) developed the ¬rst empirical re-
lationship between these variables. He obtained: 3 u—t
<q >= 3β 2 1 + ,
g γ2 γ β
2 2
0.6 γ
u—t = u2 + 1 u—t
’ u2 ,
ρ θw (8.26)
γ β
where u—td is the dry threshold friction velocity, where (x, y) is the incomplete gamma function,
⎛ ⎛ ⎛ ⎞ ⎞2 ⎞ 1
⎛ ⎛ ⎞1+» ⎞’ »+1
» 2

⎜ 0.6 1 ⎜ ⎜ ⎟ ⎟⎟
π π
2d 2 2
u—t = ⎜u2 ⎜(θs ’ θr ) ⎜⎝1 + ⎝’± b + ⎟ + θr ⎟ ⎟
⎠ ⎠
+ b+1
⎝ —td ρ θw ⎝ ⎝ ⎠ ⎠⎠
2 2 2d


In many cases of interest, including the Soda Lake
ρ is the density of air (1.1 kg/m3 ), and θw is the
example described below, the values of the in-
wilting-point moisture content (typically 0.2--0.3
complete gamma function in Eq. (8.26) are nearly
for ¬ne-grained soils (Rawls et al., 1992)).
equal to 1. In such cases, Eq. (8.26) can be approx-
The ¬nal model component is the saltation
imated as
equation for transport-limited conditions, given
by Shao and Raupach (1993) as:
< q >≈ (6β 2 ’ u2 ) (8.28)
g γ2
q = u— (u2 ’ u2 ) (8.22)
— —t
g Relative changes in saltation ¬‚ux, given by
where q is the ¬‚ux in kg/m/s, g is the gravita- Eq. (8.26), are also expected to apply to dust
tional acceleration (9.8 m/s2 ), and u— is the fric- emissions because of the proportionality between
tion velocity. The friction velocity can be obtained saltation and dust emissions observed in many ar-
from measured velocities at a height zm above the eas (Gillette et al., 2003). Dust emissions can also
ground using a modi¬ed law of the wall: be explicitly calculated using Eq. (8.26) if an esti-
κ mate is available for the K factor, or the ratio of
u— = u (8.23)
the vertical dust ¬‚ux to the horizontal saltation
ln z0
¬‚ux. The value of K depends primarily on the
where u is the measured velocity, κ is von Kar- surface texture and must be determined empiri-
man™s constant (0.4), and z0 is the aerodynamic cally.
surface roughness. Three sites at the margins of Soda (dry) Lake,
A two-parameter Weibull distribution, com- California (Figure 8.13), were selected as an illus-
monly used to quantify the distribution of wind tration of the model calibration procedure. Soda

116.1 116.0 115.9
116.2 W
36.3 N

5 km


North Soda Lake


20 m
Soda Mtns.

Soda Lake
Old Dad Mtn.



MojaveR. Crucero Devil™s Playground
6m Balch

water table depths
CLIM-MET stations springs

sites has a residual value of θr = 0.05, and that
Fig 8.13 LANDSAT image of Soda (dry) Lake study site,
with locations of CLIM-MET stations and water table depths the soil moisture is within a few percent of
indicated. Modi¬ed from Pelletier (2006). Reproduced with
that value most of the time. For this reason, the
permission of Elsevier Limited.
steady-state approximation is an appropriate esti-
mate for long-term-average moisture conditions,
although it does not represent transient effects
Lake is located at the western edge of the Mo-
of moisture in the days to weeks following rare
jave National Preserve just south of Baker, Cal-
precipitation events.
ifornia. CLIM-MET data collected from 1999 to
Wind speeds measured 2 m above the ground
present by a team from the US Geological Sur-
were converted to friction velocities using Eq.
vey are available for these areas. The water table
(8.23) with z0 = 0.005 m. This value was obtained
beneath Soda Lake is at or near the ground sur-
from Marticorena et al.™s (1997) relationship z0 =
face along the western edge of the playa, deepen-
0.005h c , together with a canopy height of h c =
ing to the east to 20 m and lower beneath the
1 m. Friction velocities inferred from the mea-
alluvial fans on the eastern side of the basin.
sured wind speed data were ¬t to the two-
Water-table depths indicated in Figure 8.13 were
parameter Weibull distribution. The best-¬t pa-
obtained from the US Geological Survey National
rameters obtained for the Balch and Crucero sta-
Groundwater Database.
tions are β = 0.398 and 0.405, and γ = 2.17 and
Figure 8.14 illustrates the CLIM-MET station
2.15, respectively. Based on these numbers, we
data for hourly rainfall, soil moisture, peak wind
chose β = 0.40 and γ = 2.15 as representative val-
speed, and the number of particles in saltation
ues for the meteorological conditions at Soda
measured at the Balch and Crucero stations near
Lake. The dry threshold friction velocities var-
the southern playa margin. Figures 8.15b and
ied from 0.5 to 0.8 m/s at the three CLIM-MET
8.15e illustrate that the soil moisture at both

(a) (e)





(— 104)
(— 104)

tively, for a range of water table depths up to
Fig 8.14 Hourly CLIM-MET data from (a)“(d) Crucero and
10 m. Figure 8.15c indicates that saltation (and
(e)“(h) Balch stations. (a) and (e) Rainfall, (b) and (f)
volumetric soil moisture, (c) and (g) peak wind speed at 2 m hence dust emission) is essentially absent for wa-
above ground, (d) and (h) number of particles in saltation, as ter tables within 4 m of the surface for this set of
measured by SENSIT instruments. Modi¬ed from Pelletier
parameters. As the water table depth increases,
(2006). Reproduced with permission of Elsevier Limited.
saltation ¬‚ux also increases, asymptotically ap-
proaching a maximum value of nearly 6 g/m2 /s.
stations. We chose a value of u—td = 0.5 for our For water tables deeper that 10 m, the surface is
dry enough that soil moisture plays a relatively
calculations, but clearly a spatially distributed
minor role. The shape of the saltation curve de-
model that resolves variations in threshold fric-
pends primarily on the value of the wilting point
tion velocity would be the most accurate ap-
moisture θw . Figure 8.15d, for example, illustrates
proach for characterizing the playa basin as a
¬‚ux curves for representative end-member values
of 0.2 and 0.3. A larger wilting-point moisture
Soil-hydrologic parameters were estimated
means a smaller role for soil moisture in sup-
from values given in Gardner and Fireman (1958)
pressing emissions. The θw = 0.3 case, therefore,
and Rawls et al. (1992). Gardner and Fireman
has saltation initiated at a shallower water ta-
(1958) studied soils with ¬ne sandy loam and
ble depth (2 m instead of 4 m), and rises more
clay textures, and obtained b values of 0.04 and
0.055 m2 , respectively. Rawls et al. (1992) provided rapidly to its asymptotic value as a function of
estimates of ± and » appropriate for ¬ne-grained water table depth. In contrast, saltation curves
are relatively insensitive to the values of b, ± and
soils. Representative values of ± = 0.3 m’1 , » =
». Varying each of these parameters by 10%, for
0.2, and θw = 0.2 were chosen.
example, results in ¬‚ux changes of only a few
Figures 8.15a, 8.15b, and 8.15c illustrate the
percent for any particular water table depth.
solution to Eqs. (8.19), (8.26), and (8.25), respec-

104 10
(a) <q>
yz=d 5
u*td = 0.45
u*td = 0.5
2.0 10
(b) <q>
u*t 1.5 (g/m2/s)
b = 0.45
(m/s) 5
b = 0.4
10 (c)
(c) <q>
g = 1.95
g = 2.15
0 4
0 6 8 10
d (m)
Fig 8.16 Sensitivity of saltation ¬‚uxes to changes in (a) dry
threshold friction velocity, (b) Weibull scale factor, and (c)
qw = 0.3
Weibull shape factor, for a range of water table depths. From
qw = 0.2 Pelletier (2006). Reproduced with permission of Elsevier
0 Limited.
2 6
d (m)

Fig 8.15 Solutions to (a) Eq. (8.19), (b) Eq. (8.26), and (c)
dependence on saltation ¬‚ux on friction velocity.
Eq. (8.25) for the model parameters. (d) Sensitivity to the
A 10% change in the Weibull shape factor γ leads
value of the wilting-point moisture. Modi¬ed from Pelletier
to a small but signi¬cant ¬‚ux response (≈ 20%)
(2006). Reproduced with permission of Elsevier Limited.
for deep water tables.
Field measurements of saltation activity in-
Figure 8.16 illustrates the sensitivity of the dicate that threshold friction velocities are vari-
saltation curve to changes in the dry threshold able in time (e.g. Gillette et al., 1980, 1997). This
friction velocity and the mean and variability of variability can occur for many reasons, including
wind speeds. Although likely future changes in limited availability of transportable material,
regional wind speeds cannot be easily quanti¬ed, short-term bursts in near-surface turbulence,
many global climate models suggest that wind temporal changes in microsur¬cial characteris-
speeds will become more variable in the future. tics (i.e. formation and disturbance of protec-
This type of scenario should be considered as a tive crusts), and the selective entrainment of ¬ne
possibility within future air-quality management particles and subsequent surface armoring with
plans. Figure 8.16a illustrates the saltation ¬‚ux coarse lag deposits. In the absence of detailed
as a function of water-table depth corresponding quantitative observations of controlling sur¬cial
to a 10% decrease in the dry threshold friction characteristics, a stochastic model component
velocity. This 10% decrease causes a proportion- is necessary to represent this variability. Here
ate increase in saltation ¬‚ux for all water table we generalize the model equations to include a
depths. A 10% increase in the Weibull scale factor range of threshold wind velocities.
β (closely related to mean wind speed), shown in Figures 8.17a--8.17c illustrate the variability
Figure 8.16b, results in a near doubling of salta- in threshold velocities in the CLIM-MET station
tion ¬‚ux for most water-table depths. This result data. In this ¬gure, soil moisture is plotted as a
is not surprising given the nonlinear, threshold function of wind speed, with saltating conditions

North Soda Lake

(a) (d)

30 0.6
no saltation
20 0.4

10 0.2
North Soda Lake

Chepil™s model
(b) (e)
30 0.6
20 0.4

10 0.2

(c) (f)
30 0.6
20 0.4

10 0.2
0 0
10 20
16 18
20 14
0 8 10 12
peak wind speed peak wind speed
(m/s) (m/s)

indicated by black dots and non-saltating condi-
Fig 8.17 (a)“(c) Control on saltation activity by wind speed
tions indicated by gray dots. These data illustrate
and soil moisture at (a) North Soda Lake, (b) Balch, and (c)
Crucero CLIM-MET stations. Saltation activity is indicated that high wind speeds occasionally fail to pro-
with black dots; gray dots indicate no saltation. Chepil™s duce saltating conditions. Conversely, low wind
quadratic relationship between soil moisture and threshold
speeds can sometimes trigger saltation. Subset-
friction velocity is also shown. (d)“(f) The frequency of
ting the CLIM-MET data at time intervals of a few
saltation activity at a given wind speed shows a linear increase
months does not signi¬cantly reduce this over-
from a minimum threshold velocity to a maximum threshold
lap. The range of threshold friction velocities at
velocity. Modi¬ed from Pelletier (2006). Reproduced with
these locations, therefore, appears to be related
permission of Elsevier Limited.
to microsur¬cial changes that occur over intra-
annual time scales.

The CLIM-MET station data support a linear
8.8 The frequency-size distribution
relationship between excess friction velocity and
of landslides
the probability of saltation, as shown in Fig-
ures 8.17d--8.17f. Mathematically, this can be ex-
pressed as a linear increase in the probability Landslides occur in regions of steep slopes and
of saltation p from a minimum value of 0 at are often triggered by intense rainfall, rapid
u = u—td,min to a maximum value of 1 at u—td,max : snowmelt, and ground motion from earthquakes
u— ’ u—td,min and volcanic eruptions. Given the spatial and
p= (8.29)
u—td,max ’ u—td,min temporal complexity of these triggering mecha-
nisms and the underlying topography, it is not
where u—td,min and u—td,max de¬ne a range
surprising that stochastic models play a promi-
of threshold friction velocities. For a surface
nent role in our understanding of landslide pop-
that experiences minor microsur¬cial variabil-
ulations. Although our ability to forecast land-
ity through time, the difference between u—td,min
slides is limited, we can gain understanding of
and u—td,max will be small. Conversely, large sur-
landslide population statistics by combining de-
¬cial changes will result in a correspondingly
terministic models of landslide occurrence with
large range of values. Among the three CLIM-MET
stochastic models of the underlying spatial and
stations, North Soda Lake and Crucero (Figures
temporal quantities such as topographic slope
8.17d and 8.17f) exhibit a relatively small range
and soil moisture. This section illustrates that
of threshold velocities, while Balch (Figure 8.17e)
exhibits a relatively large range.
In this section we investigate the frequency-
The range of threshold friction velocities can
size distribution of landslides in areas where dif-
be included in the model by integrating Eq. (8.29)
ferent triggering mechanisms dominate. We ar-
in a piecewise manner:
gue that fractional Gaussian noises are useful
ρ u—td,max
du— p(u— ) fw (u— )u— (u2 ’ u2 ) for modeling the variations in soil moisture and
<q >= — —t
g topography that trigger landslide instability. A

number of authors have examined the frequency-

+ du— fw (u— )u— (u— ’ u—t ) (8.30) size distribution of landslides (Whitehouse and
2 2
Grif¬ths, 1983; Ohmori and Hirano, 1988;
Hovius et al., 1997). These authors have consis-
to obtain
tently noted that the distribution is power-law
<q >≈ 6β 2 ’ u2 —td,max above some threshold size. This is precisely anal-
g γ2
ogous to the Gutenberg--Richter law for earth-
+ 2(u—td,max ’ u—td,min )2
quakes, which states that the cumulative fre-
quency of events with a seismic moment greater
The incomplete gamma functions have been than M o is a power-law function of M o .
Studies suggest that a threshold dependence
approximated as unity in Eq. (8.31) to yield an
expression analogous to Eq. (8.28). The additional on soil moisture is appropriate for slope insta-
terms in Eq. (8.28) were found to have a minor blity since landslides often occur when the prod-
effect on the dependence of saltation on water- uct of rainfall duration and intensity exceeds a
table depth (e.g. Figures 8.15 and 8.16) because threshold value (Caine, 1980; Wieczoreck, 1987).
the additional terms in Eq. (8.31) are related di- In situ monitoring of pore pressure has shown
rectly to wind speed rather than soil moisture. that increases in pore pressure resulting from
Equation (8.31) does, however, provide a more ac- heavy precipitation are coincident with landsli-
curate estimation for the absolute value of the des (Johnson and Sitar, 1990). Several authors
saltation ¬‚ux than Eq. (8.28), because it explicitly have shown that variations in the frequency of oc-
represents the range of threshold wind velocities currence of landslides respond to climatic chan-
ges, with higher rates of landslides associated
observed in CLIM-MET station data.

with wetter climates (Grove, 1972; Innes, 1983; (a)
Pitts, 1983; Brooks and Richards, 1994). Further
evidence for correlations between landslides and
soil moisture is the observation that soil drainage
is the most successful method of landslide pre-
vention in the United States. Seismic triggering of
landslides has also been studied by many authors
(Keefer, 1984; Jibson, 1996). Correlations between
landslides and earthquakes can be inferred in a
number of ways. Many earthquakes have been
known to trigger large numbers of landslides and
landslide densities have a strong correlation with
active seismic belts.
The Washita experimental station has pro-
vided very useful information on the evolution Fig 8.18 (a) Grayscale map of microwave remotely sensed
of soil moisture in space and time. By character- estimates of soil moisture on June 17, 1992 at the Washita
experimental watershed. (b) Grayscale map of a synthetic soil
izing the spatial variability of soil moisture us-
moisture ¬eld with the power spectrum S(k) ∝ k ’1.8 with
ing observed data sets, we can construct statis-
the same mean and variance as the remotely sensed data.
tical models that match the observed statistics.
Modi¬ed from Pelletier et al. (1997). Reproduced with
For this purpose, microwave remotely-sensed soil
permission of Elsevier Limited.
moisture data collected from June 10 to June 18,
1992 at the Washita experimental watershed can
be used (Jackson, 1993). The watershed received
In this equation, the evapotranspiration rate ·
heavy rainfall preceding the experiment but did
is assumed to be constant in space and time.
not receive rainfall during the experiment. The
Soil moisture disperses diffusively in the soil,
data are gridded estimates of soil moisture in the
and rainfall input ξ (x, t) is modeled as a random
top ¬ve centimeters of soil calculated using the
function in space and time. Variations in rain-
algorithm of Jackson and Le Vine (1996). Each
pixel represents an area of 200 m — 200 m and fall from place to place cause spatial variations
in soil moisture that are damped by the effects
the total area considered is 45.6 km by 18.6 km.
of diffusion and evapotranspiration. Without spa-
The soil moisture values do not correlate with re-
tial variations in rainfall input, there are no vari-
lief in the watershed, which is very small. Figure
ations in soil moisture according to this model.
8.18a shows a grayscale image of the soil mois-
As pointed out by Rodriguez-Iturbe et al. (1995),
ture at Washita on June 17, 1992. White spaces
the variability at the small scale of the Washita
indicate areas where the watershed is interrupted
watershed is not likely to be the result of spatial
by roads or lakes. Figure 8.18b is a synthetic
variations in rainfall.
two-dimensional fractional Gaussian noise with
β = 1.8 constructed using the Fourier-¬ltering An alternative approach to this problem
which will generate spatial variations in soil
method. The mean and variance of the synthetic
moisture at these small scales assumes that the
moisture ¬eld have been chosen to match those
evapotranspiration rate is a random function in
of the observed image. The synthetic image re-
space and time. This models spatial and temporal
produces the correlated structure of the real soil
variations in evapotranspiration resulting from
moisture image.
variable atmospheric conditions and heterogene-
Entekhabi and Rodriguez-Iturbe (1994) pro-
ity in soil, topography, and vegetation character-
posed a partial differential equation for the dy-
istics. The resulting equation is
namics of the soil moisture ¬eld s (x, t):

‚s ‚s
= D ∇ 2 s ’ ·s + ξ (x, t) = D ∇ 2 s ’ ·(x, y, t)s
(8.32) (8.33)
‚t ‚t

same distribution that was observed in real soil
moisture ¬elds is obtained.
Landslides occur when the shear stress ex-
102 ceeds a threshold value given approximately by
N ∝ A’0.8 the Mohr--Coulomb failure criterion (Terzaghi,
N (>A)
„ f = „0 + (σ ’ u) tan φ (8.35)
where „0 is the cohesive strength of the soil, σ
is the normal stress on the slip plane, u is the
pore pressure, and φ is the angle of internal
100 102
101 103 friction. Landslides are initiated in places where
A (km2)
„ f is greater than a threshold value. The move-
ment of the soil at the point of instability in-
Fig 8.19 Cumulative frequency-size distribution of patches
creases the shear stress in adjacent points on
of soil moisture larger than a threshold value for the synthetic
soil moisture ¬eld with power spectrum S(k) ∝ k ’1.8 . the hillslope causing failure of a connected do-
Modi¬ed from Pelletier et al. (1997). Reproduced with main with shear stress larger than the thresh-
permission of Elsevier Limited.
old value. To model the landslide instability and,
in particular, the frequency-size distribution of
landslides, it is therefore necessary to model
This equation is a variant of the Kardar--Parisi--
the spatial variations of „0 , σ , u, and φ. The
Zhang equation (Kardar et al., 1986), which has
variables „0 and u are primarily dependent on
received a great deal of attention in the physics
soil moisture for a homogeneous lithology and
literature. This equation is nonlinear and can
a slip plane of constant depth. The dependence
only be solved numerically.
of each of these variables on soil moisture has
Amar and Family (1989) have solved this equa-
been approximated using power-law functions
tion and have found that the solutions have tran-
(Johnson, 1984). The shear stress and normal
sects with power spectra that have a power-law
stress are linearly proportional to soil moisture
dependence on wave number with exponent β =
through the weight of water in the soil. The
1.8. The similarity between the ¬eld generated by
shear stress and normal stress are also trigono-
Eq. (8.33) and the soil moisture observed in the
metric functions of the local slope. Based on the
Washita images suggests that Eq. (8.33) may be
results given earlier, the soil moisture will be
capturing much of the essential dynamics of soil
modeled as a two-dimensional fractional Brow-
moisture at these scales.
nian walk with β = 1.8. Power spectral analyses
Rodriguez-Iturbe et al. (1995) quanti¬ed the
of one-dimensional transects of topography have
scale-invariance of soil moisture variations by
shown topography to have a power-law power
computing the cumulative frequency-size distri-
spectrum with β = 2 (e.g. Huang and Turcotte,
bution of patches of soil with moisture levels
1989). This behavior is applicable only over a cer-
higher than a prescribed value. They found that
tain range of scales, however. At scales smaller
the cumulative frequency-size distribution had a
than the inverse of the drainage density, the to-
power-law function on area with an exponent of
pography is controlled by hillslope processes and
approximately ’0.8. The number of soil patches
does not exhibit fractal behavior. To represent the
N depended on the area according to
small-scale smoothness of hillslope topography
N (> A) ∝ A ’0.8 as well as the fractal behavior at larger scales,
we can ¬rst construct a two-dimensional surface
with one-dimensional power spectrum S(k) ∝ k ’2
This is a power-law or fractal relation. We have
determined the equivalent distribution for the with the Fourier-¬ltering method. At small scales,
synthetic soil moisture ¬eld illustrated in Fig- however, the synthetic topography is linearly in-
ure 8.18b. The result is plotted in Figure 8.19. The terpolated. A shaded-relief example of the model

shows that there are no contours smaller than
8 — 8 pixels.
The shear stress necessary for failure is a com-
plex function of soil moisture and slope. How-
ever, to show how landslide areas may be asso-
ciated with areas of simultaneously high levels
of soil moisture and steep slopes, we assume
a threshold shear stress criterion proportional
to the product of the soil moisture and the
slope. In addition, we assume that slope and
soil moisture are uncorrelated. A grid of syn-
thetic soil moisture and topography of 512 —
(b) 512 grid points was constructed according to the
models described above. The domains where the
product of the soil moisture and the topogra-
phy were above a threshold value are shown in
Figure 8.21a. The threshold value was chosen
such that only a small fraction of the region
was above the threshold. Figure 8.21b shows the
cumulative frequency-size distribution of the re-
gions above threshold, our model landslides. It
can be seen that at large areas a power-law dis-
tribution with an exponent of ’1.6 exists. Actual
Fig 8.20 (a) Synthetic model of topography used in the
landslide distributions, discussed below, show a
stochastic landslide model. The topography is a scale-invariant
similar trend. Distributions obtained with differ-
function with S(k) ∝ k ’2 above a scale of eight lattice sites
ent realizations of the model yielded power-law
and is planarly interpolated below that scale. (b) Slope model
exponents of 1.6 ± 0.1 when ¬t to the landslides
corresponding to the topographic model of (a). The slopes
above A = 10. This form was independent of the
are white noise above a scale of 8 lattice points and constant
threshold value chosen for landslide initiation, as
below that scale. Modi¬ed from Pelletier et al. (1997).
Reproduced with permission of Elsevier Limited. long as the threshold was chosen such that only
a small fraction of the lattice is above threshold.
The exponent of this distribution is more neg-
ative than that of the soil moisture patches of
topography is illustrated in Figure 8.20a along Figure 8.19. This results from the less correlated
with its contour plot. The plot has 128 — 128 slope ¬eld ˜˜breaking up™™ some of the large soil
grid points with interpolation below a scale of 8 moisture patches. The effect of smooth topogra-
pixels. The result of the interpolation is clearly phy at the hillslope scales is to create a rolloff
identi¬ed as piecewise linear segments in the in the frequency-area distribution for small land-
transects along the boundaries of the plot. The slide areas. This is a consequence of the fact that
slope corresponding to this model of topography fractal scaling breaks down at hillslope scales,
is illustrated in Figure 8.20b. Below a scale of 8 and hence slopes tend to fail as a unit rather
pixels, the slopes are constant. Above this scale, than triggering many small, isolated landslides
the slopes are a two-dimensional Gaussian white on the same slope. The effect of strong ground
noise. This follows from the fact that our model motion from earthquakes is to lower the shear
for topography at large scales, a Brownian walk stress necessary for failure. This does not alter
with S(k) ∝ k ’2 , can be de¬ned as the summation the frequency-size distribution of landslides ac-
of a Gaussian white noise time series. White noise cording to this model since the form of the
means that adjacent values are totally uncorre- distribution is independent of the value of the
lated. The contour map of the slope function threshold.

N ∝ A’1.6
N (>A)


(a) (b)

100 101 102
A (km2)
were triggered by heavy rainfall. Large landslides
Fig 8.21 (a) Contour map of the product of the synthetic
in the two areas match power-law distributions
soil moisture ¬eld with the synthetic slope function. Areas
from Eq. (8.36) with exponents ’1.6 and ’2. Fig-
inside the contour loops represent model landslides. (b)
Cumulative frequency-size distribution of model landslides. ure 8.22c presents the cumulative frequency-area
The distribution compares favorably to the distributions of distribution of more than 11 000 landslides trig-
real landslides. Modi¬ed from Pelletier et al. (1997). gered by the 1994 Northridge earthquake over an
Reproduced with permission of Elsevier Limited.
area of 10 000 km2 (Harp and Jibson, 1995). Work-
ing with high-resolution aerial photography ac-
quired the day after the earthquake, Harp and
Do actual landslide distributions produce a
Jibson estimate that their catalog of landslides
similar distribution? Whitehouse and Grif¬ths
is complete for landslides larger that 25 m2 . At
(1983), Ohmori and Hirano (1988), and Hovius
smaller sizes, a signi¬cant number of landslides
et al. (1997), among others, have all presented
may have been missed. Large landslides in this
evidence that landslide frequency-size distribu-
dataset have a power-law dependence on area
tions are power-law functions of area for large
with an exponent of approximately ’1.6. As in
landslide areas. Figure 8.22 presents cumula-
the Japan and Bolivia datasets there is a roll
tive frequency-size distributions of landslide area
off in the power-law distribution for small areas.
from three areas ¬rst analyzed in Pelletier et al.
These results suggest that cumulative frequency-
(1997). In Figure 8.22a, the number of landslides
size distributions are remarkably similar despite
with an area greater than A are plotted as a func-
different triggering mechanisms, and that they
tion of A for seven lithologic units from a dataset
match the predictions predicted by a stochastic
of 3 424 landslides with areas larger than 104 m2
model for landslides triggered by soil moisture
in the Akaishi Ranges, central Japan (Ohmori and
and topography.
Sugai, 1995). For large landslides, the distribution
Some studies have suggested that rolloffs at
is well characterized by a power law with an ex-
ponent of approximately ’2, that is the small end of landslide distributions are the
result of an incomplete catalog (Stark and Hov-
N (> A) ∝ A ’2 ius, 2001). While it is certainly true that power-
law relationships describe only the upper tail of
the landslide size distribution, it is not true that
In this region landslides occur as a result of both
breaks in scaling at the lower end of the dis-
heavy rainfall and strong seismicity. Figure 8.22b
tribution are primarily due to incomplete sam-
presents two sets of landslide areas mapped in
pling. Consider the Northridge dataset shown in
adjacent watersheds in the Yungas region of the
Figure 8.22c. If power-law scaling held down to
Eastern Cordillera, Bolivia. All of these landslides

(a) 103 (b)
N (>A)
N (>A) N∝A
N∝A 102


100 100
10 10
10 10
A (km2) A (km2)

N (>A)


10 10
10 10
A (km2)

landslides clearly shows that there is some phys-
Fig 8.22 (a) Cumulative frequency-size distribution of
ical mechanism for limiting these small land-
landslides in six lithologic zones in Japan. The distributions are
well characterized by N(> A ) ∝ A ’2 above approximately slides. In the model of this section, the break
0.1 km2 . (b) Cumulative frequency-size distribution of in topographic scaling at the hillslope scale is
landslides in two areas in Bolivia. (c) Cumulative frequency-
the cause of the rolloff in the landslide size dis-
size distribution of landslides triggered by the 1994
tribution. The planarity of slopes at small scales
Northridge, California earthquake. This dataset is
means that when slides are triggered, large sec-
characterized by N(> A ) ∝ A ’1.6 above approximately
tions of the slope tend to fail, with fewer isolated
0.01 km2 .
sections of slope failure relative to a power-law
landslide areas as small as 10’4 km2 , for exam-
ple, the power-law extrapolation in Figure 8.22c
would predict approximately 1 million landslides 8.9 Coherence resonance and the
larger than 10’4 km2 . Working with detailed
timing of ice ages
aerial photographs, Harp and Jibson found that
they could resolve landslides as small as 25 m2 .
Therefore, if the break in power-law scaling is The behavior of Earth™s climate system in the last
due to undersampling, this implies that Harp and 2 Myr is a consequence of both deterministic and
Jibson missed approximately 990 000 landslides stochastic forces. Signi¬cant periodicities exist in
at least four times larger than the minimum area time series of Late-Pleistocene global climate (i.e.
they could resolve in their images. The fact that proxies for temperature and ice volume) near
so few small landslides occur in the Northridge 29 and 41 kyr. These periodicities are controlled
case relative to the power-law trend for larger by the precession and obliquity of Earth™s orbit.

There is little agreement, however, on which pro- can make a randomly kicked ball move back and
cesses and feedbacks control other characteristics forth periodically between the two hills if a de-
of Late-Pleistocene climate, including the 100-kyr layed feedback is introduced into the system. In
cycle and the asymmetry of glacial-interglacial this climate system, the growth and decay of ice
transitions. In this section we explore a simple sheets and the motion of the lithosphere beneath
model of Earth™s climate system in the Pleis- them provide that delayed feedback.
tocene as an example of how stochastic processes There are two time scales in this model sys-
can amplify nonlinear system behavior. Although tem: the delay time and the average time that
not strictly limited to surface processes, the cli- the ball spends in each valley before switching
mate system over these time scales is strongly over to the other valley. This latter time scale
controlled by both ice sheets and the lithospheric is controlled by the height of the hill and the
response to ice-sheet loading. As such, the cli- magnitude of the random variability. If these two
mate system involves key aspects of surface pro- time scales are close in value, a resonance occurs
cesses we have discussed in earlier chapters. in which the random switching between valleys
Many models have been introduced over the is made more regular by the delayed feedback,
past few decades that reproduce one or more of which encourages a repetition in the history of
the key features of Earth™s climate in the Pleis- the system. For example, if a series of random
tocene. Several of the most successful models kicks move the ball over the hill, the delay feed-
have focused on nonlinear feedbacks between back will act to push the ball up over the hill in
continental ice-sheet growth, radiation balance, the same manner at a later time by providing sys-
and lithospheric de¬‚ection (e.g. Oerlemans, tematic kicks in the uphill direction. This behav-
1980b; Pollard, 1983). In this section we will ex- ior is strongly self-reinforcing: the more regularly
plore a simple model of Late-Pleistocene climate a ball makes its way up over the hill, the more
that reproduces the dominant 100-kyr oscillation regularly will the delay feedback work to system-
of ice ages. The mechanism for producing this atically push it up the hill at the later time.
oscillation process is ˜˜coherence resonance.™™ Co- The standard equation that displays coher-
herence resonance has been investigated by sta- ence resonance is given by
tistical physicists for several years (e.g. Masoller,
= T n ’ T n3 + T n’„ + D ·n
2002). The necessary components of a coherently- (8.37)
resonating system are (1) a system with two stable
states, (2) a delayed feedback, and (3) random vari- where T is the system variable (e.g. elevation in
ability (Tsimring and Pikovsky, 2001). A number the case of a rolling ball, or temperature in the
climate system), n is the time step, „ is the delay
of physical, chemical, and biological systems have
time, and · is a white noise with standard de-
these components and exhibit coherence reso-
nance. A simple example of a system with two viation D . The ¬rst two terms on the right side
stable states is a ball that rolls back and forth be- of Eq. (8.37) represent bistability in the system.
tween two valleys when pushed over a hill. Ran- The delay feedback and noise are represented by
dom noise may be introduced into this system to the third and fourth terms, respectively. Systems
produce random ˜˜kicks™™ that push the ball over with both positive and negative feedback may oc-
the hill and into the other valley if enough kicks cur depending on whether is positive or nega-
work in the same direction. If the kicks are ran- tive. Coherence resonance occurs in both cases,
dom, they generally tend to cancel each other out but with a different periodicity.
to keep the ball trapped in one valley. From time An example of model output from Eq. (8.37) is
to time, however, a series of kicks in the same di- plotted in Figure 8.23a, with the corresponding
rection will occur that move the ball over the hill. power spectrum plotted in Figure 8.23b. Model
parameters for this example are „ = 600, =
Delayed feedback occurs if an additional force is
’0.1, and D = 0.1. The power spectrum in Fig-
exerted on the ball that is not random, but de-
pends on the elevation of the ball at a previous ure 8.23b was computed by averaging the spec-
time. The basic idea is that coherence resonance tra of 1000 independent model time series. The

spectrum in Figure 8.23b is constant for low fre-
(a) quencies and proportional to f ’2 above a thresh-
1.0 old frequency. This is called a Lorentzian spec-
trum and is a common feature of stochastic
models with a restoring or negative feedback. A
Lorentzian spectral signature was documented in
Late-Pleistocene paleoclimatic time-series data by
Komintz and Pisias (1979), suggesting a stochastic
element to the climate system limited by a neg-
ative feedback mechanism at long time scales.
This observation was the basis for many stochas-
1000 2000
0 tic climate models including Hasselman (1971),
North et al. (1981), Nicolis and Nicolis (1982), and
(b) Pelletier (1997), among others.
The climate model we consider here is a single
equation for global temperature analogous to Eq.
S( f ) (8.37). The Earth maintains a mean global tem-
S( f ) = C
perature of approximately 15 —¦ C through a bal-
ance between incoming solar radiation and long-
S( f ) ∝ f ’2 wavelength outgoing radiation. Long-wavelength
outgoing radiation is characterized by the Stefan--
10 Boltzmann Law, in which outgoing radiated en-
ergy depends on the fourth power of absolute
temperature. For small temperature ¬‚uctuations
10 10
around an equilibrium temperature, the T 4 Boltz-
mann dependence may be linearized to yield the
Fig 8.23 Example of coherence resonance in Eq. (8.37). (a)
¬nite-difference equation (North et al., 1981)
Time series of Eq. (8.37) with „ = 600, = ’0.1, and
D = 0.1. The time series of model output exhibits a cycle
= ’B T n
with an average period of 300 time steps. (b) Average power C (8.38)
spectrum of 1000 samples of the model illustrated in (a). The
power spectrum has a dominant peak with a period of „/2
where n is an index of the time step, C is
and smaller-amplitude periodicities at odd harmonics of the
the heat capacity per unit surface area of the
dominant periodicity, superposed on a “background”
atmosphere-ocean-cryosphere system, B is the co-
spectrum which is constant at small frequencies and
ef¬cient of temperature dependence of outgo-
proportional to f ’2 at larger frequencies. From Pelletier
ing radiation, and T n is the temperature differ-
ence from equilibrium. This equation character-
izes a negative feedback: global warming (cool-
time series in Figure 8.23a has two stable states ing) resulting from short-term climatic ¬‚uctua-
at T = +1 and T = ’1. The system switches be- tions results in more (less) outgoing radiation. C
tween these states with a period equal to „/2, may be estimated from the speci¬c heat of wa-
ter, 4.2 — 103 J/kg/—¦ C, and the mass of the oceans,
despite the absence of any periodic terms in the model
1.4 — 1021 kg, and the surface area of the Earth,
equation. In Figure 8.23 we have illustrated the
5.1 — 1014 m2 , to obtain C = 1.2 — 1010 J/—¦ C, as-
negative-feedback case by choosing < 0. Posi-
tive feedback produces a similar output, but with suming that the heat capacity of the climate sys-
a dominant period equal to „ rather than „/2 tem is dominated by the oceans. The value of B
(Tsimring and Pikovsky, 2001). Less-dominant pe- is constrained from satellite measurements to be
2.1 J/s/m2 /—¦ C (North, 1991). The ratio of B to C
riodicities also occur at odd harmonics of the
dominant periodicity (i.e. periods of „ /6, „ /10, is de¬ned to be c 1 in this model and de¬nes
„ /14, etc.) in Figure 8.23b. The ˜˜background™™ a radiative-damping time scale for the climate

(a) (b) incoming
”T ”T
”t ”T
”T ”t

glacial glacial interglacial

broad expansion
in northern Canada
narrow expansion
in southern Canada

(c) net change (d) W/bedrock deflection
”T ”T
”T = c |v|
”t ”t
”t 3

T |v| T

wavelength outgoing radiation. This variability
Fig 8.24 Schematic diagram of processes included in our
may be modeled by introducing a random noise
radiation-balance model. Each graph is a sketch of the
in Eq. (8.38) to obtain (Hasselman, 1971; North


. 11
( 14)