ńņš. 11 |

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

8.5 UNSTEADY EROSION AND DEPOSITION IN EOLIANENVIRONMENTS 195

2 1 10

(a) (b) (c)

(g/cm2)

h

RT

10

(g/cm2/yr)

R

1 0.5

10

bounded random walk

model

0 0 10

101 105

10

0 5 15 103 104

100

4 102

2 T (k.y.)

0

T (k.y.)

T (yr)

2

(d) (e) (f)

10

(g/cm2)

(g/cm2/yr)

h

log R

glacials

R

interglacials

1

bounded

random walk

cyclical-climate

cyclical-climate

model

0 10

15

10 log T

5

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-

by

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

196 STOCHASTIC PROCESSES

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

8.6 STOCHASTIC TREES AND DIFFUSION-LIMITED AGGREGATION 197

(a) (b) (c)

(a)

(b)

(d)

(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

198 STOCHASTIC PROCESSES

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

(a)

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

8.7 ESTIMATING TOTAL FLUX BASED ON A STATISTICAL DISTRIBUTION OF EVENTS 199

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:

time.

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

a

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

1

tanā’1

z= Ļ

tion ļ¬‚ux and the vertical dust ļ¬‚ux (Gillette et al., (8.17)

(E /a)b + 1

b+1

E E

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

port.

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

tan

Ļ 2 2d

imum dust-producing potential and changes in

subsurface moisture have little or no effect on (8.19)

200 STOCHASTIC PROCESSES

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 Ī³ Ī²

1

2 2

Īøz=d

0.6 Ī³

uā—t = u2 + 1 uā—t

(8.21)

ā—td

ā’ u2 ,

Ļ Īøw (8.26)

ā—t

Ī³ Ī²

where uā—td is the dry threshold friction velocity, where (x, y) is the incomplete gamma function,

and

ā ā ā ā ā2 ā 1

ā ā ā1+Ī» āā’ Ī»+1

Ī» 2

ā 0.6 1 ā ā ā āā

2

Ļ Ļ

2d 2 2

uā—t = āu2 ā(Īøs ā’ Īør ) āā1 + āā’Ī± b + ā + Īør ā ā

ā ā

+ b+1

tan

ā ā—td Ļ Īøw ā ā ā ā ā

Ļ

2 2 2d

(8.27)

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)

ā—t

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

zm

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

8.7 ESTIMATING TOTAL FLUX BASED ON A STATISTICAL DISTRIBUTION OF EVENTS 201

116.1 116.0 115.9

116.2 W

36.3 N

5 km

8m

North Soda Lake

36.2

20 m

Soda Mtns.

Soda Lake

Old Dad Mtn.

1m

36.1

7m

6m

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

202 STOCHASTIC PROCESSES

(a) (e)

(f)

q

q

(b)

(g)

(c)

(h)

(d)

(Ć— 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

whole.

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-

8.7 ESTIMATING TOTAL FLUX BASED ON A STATISTICAL DISTRIBUTION OF EVENTS 203

104 10

(a)

(a) <q>

(g/m2/s)

102

yz=d 5

u*td = 0.45

(m)

u*td = 0.5

100

0

2.0 10

(b)

(b) <q>

u*t 1.5 (g/m2/s)

b = 0.45

(m/s) 5

1.0

b = 0.4

0

0.5

10

10 (c)

(c) <q>

<q>

(g/m2/s)

(g/m2/s)

5

5

g = 1.95

g = 2.15

0

0 4

2

0 6 8 10

d (m)

10

(d)

<q>

Fig 8.16 Sensitivity of saltation ļ¬‚uxes to changes in (a) dry

(g/m2/s)

threshold friction velocity, (b) Weibull scale factor, and (c)

qw = 0.3

5

Weibull shape factor, for a range of water table depths. From

qw = 0.2 Pelletier (2006). Reproduced with permission of Elsevier

0 Limited.

10

8

4

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

204 STOCHASTIC PROCESSES

50

1.0

North Soda Lake

40

0.8

(a) (d)

30 0.6

no saltation

q

p

saltation

(%)

20 0.4

10 0.2

North Soda Lake

0

0

50

1.0

Balch

Chepilā™s model

40

0.8

(b) (e)

30 0.6

q

p

(%)

20 0.4

10 0.2

Balch

0

0

50

1.0

Crucero

40

0.8

(c) (f)

30 0.6

q

p

(%)

20 0.4

10 0.2

Crucero

0 0

15

10 20

16 18

20 14

0 8 10 12

5

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.

8.8 THE FREQUENCY-SIZE DISTRIBUTION OF LANDSLIDES 205

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)

approach.

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

uā—td,min

number of authors have examined the frequency-

ā

+ duā— fw (uā— )uā— (uā— ā’ uā—t ) (8.30) size distribution of landslides (Whitehouse and

2 2

uā—td,max

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-

(8.31)

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.

206 STOCHASTIC PROCESSES

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

(b)

(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

8.8 THE FREQUENCY-SIZE DISTRIBUTION OF LANDSLIDES 207

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)

1962):

Ļ„ f = Ļ„0 + (Ļ ā’ u) tan Ļ (8.35)

101

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

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,

(8.34)

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

208 STOCHASTIC PROCESSES

shows that there are no contours smaller than

(a)

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.

8.8 THE FREQUENCY-SIZE DISTRIBUTION OF LANDSLIDES 209

N ā Aā’1.6

102

N (>A)

101

(a) (b)

100

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-

(8.36)

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

210 STOCHASTIC PROCESSES

103

(a) 103 (b)

NāA

N (>A)

N (>A) NāA

102

NāA 102

101

101

100 100

100

10

10 10

10

10 10

A (km2) A (km2)

104

(c)

NāA

103

N (>A)

102

101

100

10 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

distribution.

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.

8.9 COHERENCE RESONANCE AND THE TIMING OF ICE AGES 211

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,

Tn

= T n ā’ T n3 + T nā’Ļ„ + D Ī·n

2002). The necessary components of a coherently- (8.37)

t

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

212 STOCHASTIC PROCESSES

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

T

models with a restoring or negative feedback. A

Lorentzian spectral signature was documented in

0.0

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),

t

North et al. (1981), Nicolis and Nicolis (1982), and

(b) Pelletier (1997), among others.

10

The climate model we consider here is a single

equation for global temperature analogous to Eq.

10

S( f ) (8.37). The Earth maintains a mean global tem-

S( f ) = C

perature of approximately 15 ā—¦ C through a bal-

10

ance between incoming solar radiation and long-

S( f ) ā f ā’2 wavelength outgoing radiation. Long-wavelength

10

outgoing radiation is characterized by the Stefan--

10 Boltzmann Law, in which outgoing radiated en-

ergy depends on the fourth power of absolute

10

temperature. For small temperature ļ¬‚uctuations

10 10

f

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

Tn

= ā’B T n

with an average period of 300 time steps. (b) Average power C (8.38)

t

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-

(2003).

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

8.9 COHERENCE RESONANCE AND THE TIMING OF ICE AGES 213

outgoing

(a) (b) incoming

Ī”T Ī”T

Ī”t Ī”T

Ī”T Ī”t

Ī”t

Ī”t

Th

interglacial

glacial glacial interglacial

=

+

T T

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

interglacial

glacial

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 |