equal to zero.

incised reaches and wide, shallow reaches. The

This analysis predicts the growth curve shown

most abrupt slope break is where the channel

in Figure 7.10b. The growth rate is positive

changes from distributary to incised. This break

(i.e., perturbations are unstable) for small wave-

is similar to the headcuts often observed in ar-

lengths, rises to a steep maximum at

royos and discontinuous ephemeral streams. The

basic two-step model evolution is a robust feature

w0h0

»=π (7.36) of the model and was observed for a wide range

S0

of model parameters.

and quickly becomes negative (i.e., perturbations As the instability develops from random ini-

decay to zero) for larger wavelengths. The growth tial conditions, the dominant wavelength in-

rate is a function of κ (i.e., channels with larger creases, and the channel geometry becomes

values of κ develop oscillations more rapidly), but more regularly periodic. The increase in oscilla-

the wavelength corresponding to the maximum tion wavelength indicates that nonlinear, ¬nite-

growth rate (Eq. (7.36)) is independent of κ. This amplitude effects are an important part of the

independence is important for testing Eq. (7.36) model behavior. This result is con¬rmed by Fig-

against ¬eld measurements because κ is not well ure 7.11b, which compares the oscillation wave-

constrained, but the average channel width, bank lengths predicted by the linear-stability analy-

height, and slope can be readily measured for any sis and the direct numerical solution for dif-

ferent values of (w 0 h 0 /S0 )1/2 . The lower solid

channel.

The two-step Lax--Wendroff scheme (e.g., Press line corresponds to the linear-stability predic-

et al., 1992) was used to integrate Eqs. (7.28) and tion (Eq. (7.36)), and the ¬lled squares and up-

(7.29) for the direct numerical solution. Code per solid line are the fully nonlinear numerical

that implements this solution is presented in Ap- results.

pendix 6. The width was not allowed to go be- To test the model predictions for oscillation

low 10 m or above 300 m in the model. Without wavelength, we constructed a GIS database of

bounds, the model develops in¬nitely small and DOQQs, recti¬ed false-color Landsat imagery, and

large channel widths that are unrealistic. The 30-m resolution DEMs for all of southern Arizona.

model results are not sensitive to the speci¬c val- This database was used to measure oscillation

ues of the lower and upper bounds as long as wavelengths and average channel widths, depths,

they are small and large, respectively, compared and slopes for 15 oscillating channels. Our data

to the equilibrium channel width. set includes examples of each of the three chan-

Figure 7.10c shows plots of bed elevation nel types, including channels over a broad range

h (top) and channel width w (bottom) for the of sizes. Figure 7.11b gives the average oscillation

early (left) and late (right) stages of the numer- wavelength measured for each channel as a func-

√

tion of w 0 h 0 /S0 . The data points corresponding

ical model. We assumed an initial width, bank

height, and slope of 100 m, 2 m, and 0.01, re- to 10 yr ¬‚ow depths are plotted with circles; 25 yr

spectively, with small (1%) random variations su- ¬‚ow depths are plotted with triangles. The agree-

perimposed on the initial channel width. The ment between the observed and predicted trend

solutions are plotted for a temporal sequence, in the data are quite good. This agreement pro-

with thicker lines representing later times. The vides validation for the quantitative prediction of

early stage of the model is characterized by the nonlinear model.

174 INSTABILITIES

drumlins often parallels the topographic surface,

(a)

suggesting that they formed by localized, upward

¬‚ow of sediment. If so, what were the driving

S V1

Bo Ce H NA

forces for this ¬‚ow, and what factors controlled

the morphology of the resulting drumlins? For

DM

V2

example, why are drumlins so variable in size

and shape, even over distances as short as a few

C

kilometers (Figure 7.12)? Hindmarsh (1998) and

WB CO

100 km

Fowler (2001) proposed that shearing of an in-

V

Q Bi

compressible viscous till layer could be responsi-

P

ble for producing drumlins even in the absence

(b) NA

of a bedrock core. Schoof (2002), however, has

P

DM V2

5 V

shown that bedforms produced by these models

V1 Ce

H

Q

are not consistent with many aspects of natural

C

Bo

3

drumlins.

numerical

l S

Most conceptual drumlin models invoke ei-

solution

(km)

ther subglacial hydraulic processes or deform-

1 linear stability

ing bed processes. This section seeks to unite

Bi CO

l = p (w 0h 0/S 0)1/2

these two hypotheses by modeling sediment de-

WB

0.5

formation caused by porewater migration. In

0.3 this model, porewater migration generates buoy-

0.3 0.5 1

0.03 0.05 0.1

(w 0h 0/S 0)1/2 (km) ancy forces that drive converging, upward ¬‚ow

of the matrix to form hummocky moraine (un-

Fig 7.11 Database of oscillating channel geometries in

der ¬‚at ice sheets) and drumlins (under slop-

southern Arizona. (a) Location map. Clockwise from upper

ing ice sheets) that have a characteristic scale

left: (Bo) Bouse Wash; (Ce) Centennial Wash; (H)

that depends on the sediment texture and thick-

Hassayampa River; (S) Sycamore Creek; (NA) North Airport

ness of the till layer. This hypothesis builds on

Wash; (V1) Vail Wash; (V2) Vail-Dead Mesquite Wash; (DM)

Menzies™ (1979) conceptual model of drumlin for-

Dead Mesquite Wash; (C) Cottonwood Wash; (WB) Wild

Burro Wash; (CO) Canada del Oro Wash; (P) Penitas Wash; mation by porewater expulsion. Menzies showed

(Bi) Bobaquivari Wash; (V) Vamori Wash; (Q) East La Quituni that sedimentary microstructures often indicate

√

Wash. (b) Plot of oscillation wavelength » vs. w 0 h 0 /S0 .

that drumlin sediments had a higher concentra-

Lower line is linear-stability prediction. Upper line and ¬lled

tion of porewater than interdrumlin sediments.

squares are numerical results. Observed data for channels in

Menzies used this observation to propose that

(a) are plotted with ¬lled circles (10 yr ¬‚ood depth) and open

drumlin formation was driven by stresses associ-

triangles (25 yr ¬‚ood depth). Modi¬ed from Pelletier and

ated with porewater migration and expulsion. To

Delong (2005).

date, however, this process has not been quanti-

¬ed, nor is it known what drumlin morphologies

result from this process.

7.6 How are drumlins formed? The modeling of sediment deformation with

porewater migration involves a complex set of

Drumlins are subglacial bedforms elongated par- diffusive and advective equations. The equations

allel to the ice-¬‚ow direction and composed pri- that describe deformable porous media such

marily of subglacial till or sediment. Some drum- as subglacial sediment were ¬rst developed in

lins have a bedrock core. These drumlins most the geophysics literature to describe the ¬‚ow of

likely formed by sediment accretion in the lee- magma through porous rock (Scott and Steven-

side low-pressure zone caused by streaming ice son, 1986; Spiegelman and McKenzie, 1987) and,

¬‚ow around the core (Boulton, 1987). More enig- later, to ¬‚uid ¬‚ow in sedimentary basins (e.g.

matic, however, are the drumlins composed en- Dewers and Ortoleva, 1994). Deformable porous

tirely of sediment. Subsurface bedding in these media are comprised of two ¬‚uids with different

7.6 HOW ARE DRUMLINS FORMED? 175

(a) (b)

Rogen Moraine

Rogen Moraine

20 km

20 km

43.4°N

Columbia 43.6° N

Dodge

Oswego

Wayne

43.2

43.4

43.0

Dane

43.2

Jefferson

Onondaga

42.8

Ontario Cayuga

Seneca

76.0

76.4

76.8 89.2° W

77.2° W 88.8

phases, Darcy™s Law for the migration of the liq-

Fig 7.12 Shaded-relief maps of drumlin ¬elds in (a)

north-central New York and (b) east of Madison, Wisconsin. uid in the matrix, Stokes™ Law for the viscous ¬‚ow

Inset regions illustrate the variety of drumlin shapes in each of the matrix, and an empirical power-law rela-

area. Rogen moraines (bedforms oriented perpendicular to

tionship between permeability and porosity:

the ice-¬‚ow direction) are present at the margins.

‚φ

+ ∇ · φv = 0 (7.37)

‚t

‚φ

densities and viscosities. The low-density, low-

’ ∇ · (1 ’ φ)V = 0 (7.38)

viscosity ˜˜liquid™™ phase migrates according to ‚t

Darcy™s Law within a ˜˜matrix™™ of high-density, kφ

φ(v ’ V ) = ’ ∇ P (7.39)

μ

high-viscosity ¬‚uid governed by Stokes™ Law. De-

∇ P = ’·∇ · ∇ · V

formable porous media exhibit two kinds of

behavior not observed in incompressible ¬‚uids. 4

+ ζ + · ∇(∇ · V ) ’ (1 ’ φ) ρgk

First, the liquid can become segregated from the 3

matrix by a positive feedback between porosity, (7.40)

permeability, and matrix expansion. In this feed- kφ = aφ n (7.41)

back, preferential migration of the liquid into

where φ is the porosity, v is the liquid velocity,

regions of the matrix with higher porosity (and

V is the matrix velocity, kφ is the permeability, μ

hence permeability) expand the matrix to further

increase the porosity. In a layer undergoing com- is the liquid kinematic viscosity, P is the excess

pressure (above hydrostatic), ζ and · are the bulk

paction, this feedback leads to alternating zones

and shear viscosities of the matrix, ρ is the dif-

of high and low porosity with a characteristic

spacing that depends on the layer thickness. Sec- ference in density between the liquid and matrix,

ond, the expansion and contraction of the matrix g is the gravitational acceleration, k is the unit

generates buoyancy forces that drive convective vector in the vertical direction, and a and n are

¬‚ow of the matrix. empirical parameters.

McKenzie (1984) ¬rst proposed a set of equa- Equations (7.37)--(7.41) have been applied to 1D

tions to describe partially molten rock as a de- (vertical) compaction problems, but they are not

formable porous medium. His equations include: easily solved in 2D or 3D. However, by introduc-

conservation of mass for the liquid and matrix ing stream functions and U for the incom-

176 INSTABILITIES

pressible and compressible components of the length is the length scale at which matrix

matrix velocity V (i.e. V = ∇ · + ∇ · U ), Spiegel- deformation and liquid migration occur at simi-

man (1993) developed a set of equations readily lar rates and neither is a limiting factor. The com-

solved in 2D and 3D: paction length is broadly analogous to the ¬‚exu-

ral parameter ±, which characterizes lithospheric

φ0

2

∇4 =’ ∇ · φk (7.42) strength (Turcotte and Schubert, 2002). The ¬‚ex-

·

ural parameter represents a similar balance, in

’kφ ∇ 2 C ’ ∇kφ · ∇C + C that case between ¬‚exural rigidity and buoyancy

· under crustal loading.

= ∇ · kφ (∇ · ∇ 2 ) ’ (1 ’ φ0 φ)k (7.43)

An order-of-magnitude estimate for δ in sub-

φ0

glacial sediments can be obtained using a ma-

∇ 2 U = φ0 C (7.44)

trix shear viscosity of 1010 Pa s (Murray, 1997), a

‚φ permeability of 10’10 m2 (appropriate for sand-

+ (∇ · + ∇U ) · ∇φ = (1 ’ φ0 φ)C (7.45)

‚t rich deposits; Freeze and Cherry, 1979), and the

kinematic viscosity of water: 10’2 Poise. The bulk

kφ = aφ n (7.46)

viscosity of subglacial sediments (i.e. their vis-

where φ is the initial porosity, ζ is the ra- cous resistance to expansion and contraction) is

tio of the shear to bulk viscosities, and C is not well constrained, but for simplicity we as-

the compaction rate (de¬ned as C = ∇ · v ). Al- sume it to be equal to the shear viscosity. These

though originally developed for partially-molten values yield an order-of-magnitude estimate δ ≈

rock, Eqs. (7.42)--(7.46) are applicable to any de- 100 m. Till thicknesses in New York and Wiscon-

formable porous medium governed by Darcian sin are between 0 and 50 m, or between l = 0 and

and Stokes™ ¬‚ow. Spiegelman™s equations may ap- l = 0.5 when scaled to this compaction length.

pear formidable, but their solution is not funda- The model results described below do not depend

mentally more dif¬cult than solving the diffu- on the precise value of l or δ, as long as l < δ. If

sion or wave equations. The strategy is to solve this inequality holds, l is the length scale con-

Eqs. (7.42)--(7.46) sequentially for the porosity at trolling the instability.

time t1 given the porosity at an earlier time t0 , Spiegelman et al. (2001) solved Eqs. (7.42)--

subject to the boundary conditions of the prob- (7.46) in 2D and showed that partially molten

lem. In this study, Eqs. (7.42)--(7.44) were solved us- rock undergoes an instability in which vertically-

ing the Alternating-Direction-Implicit (ADI) tech- migrating melt focuses into regularly-spaced,

nique (Press et al., 1992), and Eq. (7.45) was solved high-porosity channels alternating with low-

using upwind differencing. porosity channels. Channels develop in Eqs.

Equations (7.42)--(7.46) are scaled to the fun- (7.42)--(7.46) by a positive feedback in which as-

damental length scale of the problem: the com- cending melt migrates preferentially into regions

paction length, δ, given by of higher initial porosity (and hence permeabil-

ity), causing the matrix to expand locally to make

kφ · + 4 ν

δ= 3

room for the melt. Matrix expansion increases

(7.47)

μ

the local porosity, further enhancing the focus-

The compaction length is the spatial scale at ing of melt in a positive feedback. Although par-

which compaction, requiring both viscous ¬‚ow tially molten rock is governed by very different

of the matrix and Darcian ¬‚ow of the liquid, physics than subglacial sediments, both are com-

takes place most rapidly. At small spatial scales, posed of liquid embedded in a deformable ma-

liquid migrates through the matrix readily, but trix, and hence both may be expected to ex-

compaction occurs slowly because of viscous re- hibit the same fundamental instability mech-

sistance of the matrix. At large spatial scales, anism during compaction. Can this instability

the matrix can readily deform, but compaction mechanism be responsible for drumlins?

is limited by the time required for the liq- Here we explore the results of three nu-

uid to migrate long distances. The compaction merical experiments designed to illustrate how

7.6 HOW ARE DRUMLINS FORMED? 177

(a) (b) (c)

deformable, flat ice sheet sloping ice sheet

ice sheet permeable

2D till

3D

3D

bedrock rigid, impermeable

Fig 7.13 Model geometry and boundary conditions for the

This boundary condition assumes the existence

three model experiments. The upper boundary is deformable

of a meltwater channel at the sediment--ice in-

and permeable, and the lower boundary is rigid and

terface that conducts porewater out of the sys-

impermeable.

tem as it leaves the matrix. The standard free-

slip boundary condition for incompressible ¬‚ows

gives ∇ · n = 0. The sediment--ice interface is a

subglacial sediment deforms under different

moving boundary in the model: during each time

sediment-thickness and ice-surface-slope condi-

step, displacement occurs by an amount V t,

tions. In the ¬rst experiment (Figure 7.13a), 2D

which raises or lowers the interface locally, de-

compaction was considered under uniform pres-

pending on the sign of V .

sure. In the second experiment (Figure 7.13b), 3D

The model parameters for all the numerical

compaction was considered under uniform pres-

experiments in this section are φ0 = 0.5, n = 3,

sure. In the third experiment (Figure 7.13c), a

and ζ = 1. The sediment thickness was varied

sloping ice sheet was considered in 3D.

between l = 0.2 (˜˜thin™™ sediment) and l = 0.6

In order to solve Eqs. (7.42)--(7.46), the values or

derivatives of , U , C , and φ must be speci¬ed at (˜˜thick™™ sediment) in the experiments. Porosities

in subglacial sediments typically vary from 0.2 to

the upper and lower boundaries of the grid. In all

0.5. A relatively high porosity value was chosen to

calculations, we assume an impermeable, rigid,

clearly illustrate the model drumlins in cross sec-

no-slip lower boundary (the bedrock-sediment

tion (i.e. the initial porosity controls the relief of

interface) and a permeable, deformable, free-

the bedforms relative to the sediment-layer thick-

slip upper boundary (the sediment--ice interface).

ness). The value n = 3 is a representative value for

The ice sheet plays only a passive role in this

granular media. The value ζ = 1 is constrained

model, providing hydrostatic pressure but no

by the assumption of equal bulk and shear vis-

shear. This contrasts with the models of Hind-

cosities. Small (5%) random variations in porosity

marsh (1998) and Fowler (2001), in which drum-

were superimposed on the uniform initial poros-

lins were formed by the shearing of an incom-

ity to represent small-scale heterogeneity of the

pressible viscous till layer by the overriding ice

sheet. At the lower boundary, the condition = sediment matrix. These random variations play

an important role in the model as ˜˜seeds™™ for

0 speci¬es that both the horizontal and vertical

the instability.

components of the matrix velocity are zero (this

Figure 7.14 illustrates the evolution of the

combines both rigid and no-slip boundary con-

model using grayscale images of the compaction

ditions). The rigidity of the lower boundary also

rate, C . Areas of white and yellow in these images

requires setting (i.e. the component of the ma-

are undergoing rapid expansion while red and

trix velocity perpendicular to a rigid boundary

is zero). C = 0 and φ = 0 de¬ne the imperme- black areas are undergoing contraction. Porosity

and compaction rate are closely correlated in the

ability of the lower boundary (Spiegelman, 1993).

model, so white and yellow regions also indicate

A deformable, permeable upper boundary re-

quires that U = 0 and that the gradients of com- high porosity. The 2D model evolution can be di-

vided into three basic stages. In the early stage

paction and porosity normal to the upper bound-

ary are zero: ∇C · n = 0 and ∇φ · n = 0, where (Figure 7.14a), regions of slightly higher poros-

ity (and hence permeability) near the top of the

n is the unit normal vector (Spiegelman, 1993).

178 INSTABILITIES

sit higher for the same reason. In the ¬nal stages

early

(a) of the model, high-porosity regions conduct the

40 — vert. remaining porewater out of the upper boundary

l = 0.2

exag.

(Figure 7.14c). Deformation stops when all of the

liquid has been squeezed from the matrix.

The time required for this instability to oc-

expansion

middle contraction

cur in natural drumlins can be estimated as

(b)

T ≈ l/K , where K is the sediment hydraulic

conductivity and l is the layer thickness. Using

a conductivity of K = 10’9 m/s (in the middle

range of glacial till values, ranging from 10’12 to

late

10’6 m/s; Freeze and Cherry, 1979) and l = 50 m,

porewater

(c)

T ≈ 1000 yr. This estimate depends sensitively on

the assumed value for conductivity; drumlins in

coarse-grained sediments could develop in less

sediment

than a year, while drumlins in ¬ne-grained sed-

Fig 7.14 2D model evolution, with 40— vertical

iments could require tens of thousands of years

exaggeration. Color plots of compaction rate C shown at (a)

to form. The overlying ice may also play a role

early, (b) middle, and (c) late times in the model. In the early

in controlling the time scale for drumliniza-

phase, alternating zones of matrix expansion and contraction

tion. The Rayleigh--Taylor instability model of salt

develop near the upper boundary as zones of initially higher

domes, considered in Section 7.2, provides a use-

porosity expand, capture upwardly migrating porewater, and

ful analogy here. It is clear that halite is less

further expand in a positive feedback. In the second stage,

viscous and has a lower yield stress than typi-

zones of expansion ascend and drive converging ¬‚ow of the

matrix into the high-porosity drumlin core. In the third stage, cal clastic sedimentary rocks. Nevertheless, this

porewater is squeezed from the matrix until the sediment is does not prevent the salt from intruding into the

fully compacted. For color version, see plate section.

overlying rock, nor does it prevent the instability

wavelength from being controlled by the thick-

layer expand as upwardly migrating porewater ness of the salt. The thickness and viscosity of

is focused into those areas. Viscous stresses sup- the salt control the wavelength because it is that

press expansion and contraction near the rigid layer that undergoes boundary-layer ¬‚ow along

lower boundary of the model. The spatial distri- the rigid lower boundary of the model. The over-

bution of compaction is largely random in this lying rock, in contrast, plays only a passive role,

early stage, re¬‚ecting random variations in ini- just as the ice plays a passive role in drumlin for-

tial porosity. As the instability develops, how- mation. Nevertheless, the rock must still ¬‚ow in

ever, regions of matrix expansion and contrac- order for the instability to occur. Therefore, the

tion grow to a maximum wavelength equal to viscosity of the rock (or ice, in the case of drum-

about ten times the layer thickness. As expand- lins) can control the time for the instability to

ing regions increase in porosity, their buoyancy occur.

increases with respect to the surrounding sedi- In the second model experiment, a 3D ge-

ment. This buoyancy imparts an upward velocity ometry was considered under uniform pressure

to both the liquid and matrix in these regions, de- (Figures 7.15a and 7.15b, ¬nal model topography

forming the ice--sediment interface above them. shown in color shaded relief). This experiment

Conversely, compacting regions become heavier results in ˜˜egg-carton™™ topography without a pre-

than surrounding regions and sink. This response ferred orientation. This topography is not simi-

can be thought of in terms of an isostatic bal- lar to drumlins, but it does closely resemble the

ance. Mid-ocean ridges sit higher, topographi- hummocky moraine found in large parts of the

cally, than the surrounding ocean ¬‚oor because Canadian prairie (Boone and Eyles, 2001). The for-

of the lower density associated with focused melt mation of hummocky-moraine topography is not

beneath spreading centers. Subglacial bedforms well understood, but this experiment suggests

7.6 HOW ARE DRUMLINS FORMED? 179

Fig 7.15 3D model evolution, map

(a) (c)

30 d view. Shaded-relief and color maps

of the model topography following

compaction. (a) Flat ice sheet, thin

sediment; (b) ¬‚at ice sheet, thick

sediment; (c) sloping ice sheet, thin

sediment; (d) sloping ice sheet, thick

sediment. For a ¬‚at ice sheet,

hummocky moraine is formed, with

the size of the hummocks controlled

by the sediment thickness. For

sloping ice sheets, drumlins are

sloping ice sheet, thin till: l = 0.2

flat ice sheet, thin till: l = 0.2

formed, with the width of drumlins

controlled by the initial

(d)

(b)

sediment-layer thickness. A sloping

ice sheet with thin sediment

reproduces bedforms most similar

to drumlins. With thicker sediment,

wide bedforms are created, some of

which join laterally to form

barchanoid drumlins and Rogen

moraine.

flat ice sheet, thick till: l = 0.6 sloping ice sheet, thick till: l = 0.6

that it may be formed by porewater expulsion un- drumlins can become wide enough to join lat-

der the nearly uniform pressure conditions found erally to form barchanoid drumlins and Ro-

beneath the interiors of great ice sheets. gen moraine (bedforms oriented perpendicular to

The 3D model runs were varied by changing the ice-¬‚ow direction). This suggests that Rogen

the initial till thickness to determine the effect moraine can be formed, in part, by the same in-

on bedform morphology (Figures 7.15a and 7.15b). stability mechanism as drumlins and hummocky

Varying the thickness results in a proportionate moraine.

change in the width of hummocky moraine. In The instability of this model occurs when the

the third model experiment (Figures 7.15c and effective pressure exceeds the yield stress, initiat-

7.15d), an ice sheet with a uniform slope of 0.01 ing viscous ¬‚ow of the matrix. The effective pres-

was assumed. To model a sloping ice sheet, the sure near the base of the sediment layer is given

vector k was replaced with (k + S j) in Eqs. (7.42) by

and (7.43), where S is the ice-surface slope and j is

pe = (1 ’ φ) ρgl (7.48)

the unit vector in the down-ice direction. Sloping

or ≈ 105 Pa for a till layer with a porosity be-

ice results in elongated, tapered bedforms very

similar to drumlins (Figure 7.15b). In this model, tween 0.2 and 0.5 that is at least 20 m thick. Typ-

the blunt up-ice end of the drumlin is associated ical yield stress values for subglacial sediments

are on the order of 104 --105 Pa, with higher values

with the early phase of sediment upwelling. As

the upwelling plume migrates down ice and pore- for well-drained sediments. These values indicate

water is expelled from the matrix, the plume that saturated tills should commonly exceed the

loses buoyancy and narrows. The streamlined threshold stress required for viscous ¬‚ow, as as-

drumlin ˜˜tail™™ represents the last vestige of sumed in the model.

porewater expulsion. The model result obtained The fundamental model prediction is a

with thick till (Figure 7.15d) also suggests that linear relationship between drumlin width and

180 INSTABILITIES

(a) (b)

l

1 w

2

threshold-filtered

shaded relief

(c)

1 w

1 w

2

2

Fig 7.16 Drumlin mapping algorithm. (a) Step 1: Construct shadow that just barely covered the shaded side

a grayscale shaded-relief image from a DEM with the

of drumlins completely. Several illumination di-

illumination direction chosen to be perpendicular to the

rections were used in different portions of each

predominant drumlin axis. (b) Step 2: Apply a threshold ¬lter

drumlin ¬eld to properly illuminate each drum-

to the shaded-relief image to make all shadows black and all

lin from a perpendicular direction. Second, the

other areas white. Step 3: Analyze the binary image for all

shaded-relief image was converted to a black and

connected clusters of black pixels.

white image using a threshold ¬lter. In this step,

shadows remain black while shades of gray are

initial sediment-layer thickness. This prediction converted to white. Third, the black-and-white

is testable using morphological analyses of drum- image was input into a custom software pro-

lins and groundwater well data that include gram that identi¬ed all connected clusters of

depths to bedrock. To do this, drumlin widths black pixels and computed their centroid posi-

were ¬rst mapped using a semi-automated al- tions, areas, lengths, and widths using standard

gorithm (Figure 7.16). This algorithm uses the de¬nitions (Chorley, 1959; Smalley and Unwin,

perpendicular shadow cast by each drumlin in 1968). Drumlin widths were taken to be twice

a shaded-relief image as a proxy for its shape. the cluster width because each cluster repre-

First, a shaded-relief image was constructed using sents only one side of the drumlin in shadow.

a high-resolution Digital Elevation Model (DEM). This algorithm leaves many small clusters re-

For the north-central New York area, 10-m reso- maining (e.g. Figure 7.16b), so the dataset was

lution US Geological Survey DEMs were used. In ¬ltered to remove clusters that are too small

(i.e. less than 0.01 km2 ) to be drumlins. This

Wisconsin, 30-m resolution DEMs were the best

available. The illumination direction was chosen procedure yielded a dataset of approximately

to be perpendicular to the predominant drum- 5800 drumlins in north-central New York, and

lin axis, and the azimuth was chosen to cast a 2900 drumlins in the Wisconsin study area.

7.6 HOW ARE DRUMLINS FORMED? 181

(a) (c)

43.4° N

43.6° N

43.2

43.4

43.0

43.2

42.8

drumlin width

89.2° W

76.0 88.8

76.4

76.8

77.2° W

0 500 m

0 1000 m

(b) (d)

Genessee

depth-to-bedrock

River Valley

60 m

0

0 60 m

areas of thick till and

= wide drumlins

100 m to 500 m in the Wisconsin study area

Fig 7.17 Grayscale maps for average drumlin width and

(Figure 7.17c).

bedrock depth in (a) and (b) north-central New York and (c)

The shadow-mapping algorithm does not pro-

and (d) Wisconsin, east of Madison. (a) and (c) Maps of

duce an ideal representation of the drumlin

drumlin width constructed by averaging widths within a 2-km

form, but it is a robust method. The principal

square moving window. (b) and (d) Maps of bedrock depth

also constructed with a 2-km square moving window and limitation of the shadow-based algorithm is that

using USGS groundwater well data. Curves are drawn to shadows created by other landforms may be in-

highlight areas in which thick sediment and wide drumlins advertently included in the analysis. This was not

(including Rogen moraine) coincide. For color version, see

a signi¬cant problem for most of the study area,

plate section.

however. Two problematic areas were found in

the southwest and southeast corners of the New

York study area. In these areas near the Finger

Finally, drumlin widths were spatially-averaged

using a 2 km — 2 km moving window (Figures Lakes, bedrock slopes are comparable to drum-

lin slopes, resulting in some bedrock ridges and

7.17a and 7.17c). The resulting values for aver-

valleys being included as drumlins. These areas

age drumlin width (including both drumlins and

(shown as shaded in Figure 7.17a) were excluded

Rogen moraine) vary from 100 m to 1000 m in

from the analysis.

the New York study area (Figure 7.17a) and from

182 INSTABILITIES

Till thicknesses were mapped using depth-to- topographic surface

(a)

bedrock (DTB) data from all publicly available

groundwater well records. Data were obtained

from the USGS New York of¬ce and the Wiscon-

sin State Geological and Natural History Survey.

In the New York study area, 2786 wells were avail-

conformable

able for six counties. In Wisconsin, 1349 wells

were available for six counties. DTB data were av- bedding surfaces

eraged using a 2 km — 2 km moving window (Fig- (b)

ures 7.17b and 7.17d). This averaging was done to

minimize the small-scale spatial variability asso-

ciated with well placement and to make the DTB

maps directly analogous to the maps of drumlin

width. DTB data in both study areas vary from 0

to 60 m. DTB data are a good proxy for till thick- unconformable

ness except in areas with signi¬cant ¬‚uvial depo-

sition. In the Genesee River Valley, for example,

(c)

¬‚uvial sediment thicknesses are comparable to

or greater than glacial sediment thicknesses in

the area, making any correlation with drumlin

widths unreliable.

Maps of average drumlin width and sediment unconformable

thickness show a strong correlation (except the

Fig 7.18 Drumlin stratigraphic end members (shown in

Genesee River Valley where ¬‚uvial sediments are

cross section, perpendicular to the drumlin axis). Three types

present). White curves in Figure 7.17 are used

are possible: (a) subsurface bedding broadly parallels the

to highlight several regions in each study area

topography (i.e. “concentric” drumlins), (b) subsurface

where wide drumlins and thick till coincide. In

bedding is undeformed, and (c) subsurface is deformed, but

the New York study area (Figures 7.17a and 7.17b) the deformation is poorly correlated with the topography.

wide drumlins (and, in places, Rogen moraine) oc- Only (a) is consistent with the proposed model.

cur in the northeast corner of the study area and

along a swath trending southwest-to-northeast

near the Lake Ontario shoreline. In the Wiscon- (1971) proposed that till-fabric data are consis-

sin study area (Figures 7.17c and 7.17d), the widest tent with a low-pressure zone along the drum-

drumlins and thickest till occur in the northwest, lin axis. Without a bedrock core to initiate the

south, and east portions of the study area. These low-pressure zone, however, it is unclear how

results prove that the geometry of the till layer low pressure develops in Evensen™s model. The

plays a dominant role in controlling the drum- buoyancy generated by concentrating porewater,

lin geometry and they lend string support to the however, provides a natural explanation for the

proposed model. alternating pressure zones and converging ma-

Sedimentological and stratigraphic observa- trix ¬‚ow inferred by Evensen and others using

tions provide additional constraints on drumlin- till-fabric analyses.

forming mechanisms. Analyses of till fabric (the Figure 7.18 summarizes three possible end

predominant orientation of elongated pebbles in members for drumlin stratigraphy (shown in

the till) provide evidence for convergent ¬‚ow of cross section, perpendicular to the drumlin axis).

sediment upward and inward toward the drum- Figure 7.18a illustrates the case in which the sub-

lin axis (Evensen, 1971; Stanford and Mickelson, surface bedding broadly parallels the topogra-

1985). Based on observations in the Wisconsin phy. Drumlins with this stratigraphy are referred

study area of Figure 7.12b, for example, Evensen to as ˜˜concentric™™ (Hart, 1997). In Figure 7.18b,

7.7 SPIRAL TROUGHS ON THE MARTIAN POLAR ICE CAPS 183

the subsurface bedding bears no indication of to form a self-sustaining poleward-migrating to-

deformation. In Figure 7.18c, sediment deforma- pographic wave. The relationship between this

tion is pervasive, but the deformation is uncorre- model and the spiral morphology of troughs has

lated with topography. Of these three end mem- not been fully established, but Fisher (1993) in-

bers, only Figure 7.18a is consistent with the pro- troduced an asymmetric ice-velocity ¬eld into

posed model, although reworking of drumlins Howard™s model and obtained spiral forms. Re-

by subsequent ice-sheet advances often leads to cent observations from Mars Global Surveyor,

erosional truncations and a more complex however, suggest that ice ¬‚ow near the troughs is

stratigraphy than the simpli¬ed case shown in not signi¬cant (Howard, 2000; Kolb and Tanaka,

Figure 7.18a. Nevertheless, concentric drumlins 2001) and it remains unclear precisely how the

are commonly observed in many drumlin ¬elds spirals form and what controls their spacing, ori-

(e.g. Newman and Mickelson, 1994; Zelcs and entation, and curvature.

Driemanis, 1997). The New York State drumlin In this section we explore a simple mathemat-

¬eld, in particular, is the classic location for con- ical model designed to capture the essential feed-

centric drumlins (Hart, 1997). back processes that couple the topography, the

As a composite of porewater embedded within distribution of solar heating on the surface, and

a matrix of sediment, subglacial sediments be- the resulting accumulation and ablation of ice.

have as deformable porous media. Porosity and The model is based on the processes in Howard™s

buoyancy play a signi¬cant role in the evolu- migrating scarp model and includes the simplest

tion of these systems. First, spatial variations mathematical descriptions of these processes in

in porosity variations can be enhanced during order to determine the necessary conditions for

compaction by a feedback between porosity, realistic spiral troughs. Lateral heat conduction is

permeability, and matrix expansion. Bouyancy as- also included in the model, and this element is

sociated with the resulting regions of high poros- crucial for obtaining realistic troughs. The model

ity provide the driving force for converging ¬‚ow does not include wind erosion or ice ¬‚ow. The

of the sediment into regularly spaced bedforms. equations are similar to those from the ¬eld of

Near the margins of former ice sheets, where excitable media, where solitary and spiral waves

the ice surface is sloping, sediment and porewa- resulting from the dynamics of two interacting

ter are advected down the ice-surface-slope direc- variables have been studied in detail.

tion, creating bedforms very similar to the classic The model equations consist of two coupled,

drumlins of New York and Wisconsin. scaled equations for the deviation of local ice-

surface temperature, T , from its equilibrium

value (the temperature at which no accumula-

tion or ablation takes place) and the deviation of

7.7 Spiral troughs on the Martian

local ice-surface topography, h, below its equilib-

polar ice caps rium value (i.e. troughs have positive h):

‚T

= κ∇ 2 T + f (T , h) (7.49)

The spiral troughs on the Martian polar ice caps

‚t

are one of the most fascinating landforms in the ‚h 1

= T (7.50)

solar system. Spiral troughs form with an insta- ‚t „f

bility in which the ice-surface temperatures of 1

f (T , h) = (T (T ’ T 0 )(1 ’ T ) ’ h) if ∇h · r > 0

steep, equator-facing slopes exceed 0 —¦ C during „i

the summer, sublimating the ice locally to form (7.51)

steeper, lower-albedo scarps (through exposure of 1

= h if ∇h · r < 0

ˆ (7.52)

subsurface dust-rich layers) in a self-enhancing „i

feedback (Howard, 1978). In this model, some of

the water vapor released from the equator-facing The deviations from equilibrium temperature

scarp may accumulate on the pole-facing slope do not include diurnal or annual cycles; these

184 INSTABILITIES

trough deepening

initial scarp formation

(a) Fig 7.19 Diagram of trough

ti < t < tf

0 < t < ti evolution and solution to Eq. (7.52)

high albedo

Sun Sun in 2D. (a) Trough initiation involves a

rapid increase in ice-surface

ice surface ice surface

temperature to its maximum value

as the scarp is steepened and dust

low albedo

warmer

(dust exposure) layers are exposed. During trough

onsun-facing

deepening, the value of h steadily

cooler

ti

slope:

T: T 0 1, h: 0 tf increases to its maximum value over

on sun-facing (shielded)

a time scale „ f = 5 Myr, while the

ti

slope:

T: 1 h: 1

tf temperature decreases in the

deeper portions of the trough due

(b) propagation to shielding. (b) Solution to (1) with

1 t1

Sun

T Runga“Kutta integration and

initial condition: x = 0.2 km, κ = 175 km2 ,

T = small pulse, h = 0

T 0.5 1/2

„i = 0.05, „ f = 1, and T0 = 0.2 at

l = (ktf )

and t0 t1 = 8. The initial condition is a

small pulse with T > T0 at the left

accumulation

0

side of the domain. Modi¬ed from

Pelletier (2004b).

ablation

30 40

20

10

0

x (km)

are assumed to be included in the equilibrium temperatures, bounded by a maximum value of 1.

temperature. Deviations from the equilibrium In the second step, the temperature remains near

temperature initiate ablation and accumulation its maximum value of 1 as the trough deepens

over a longer time scale, „ f . Trough deepening

at any time, t, in the model, although accu-

mulation and ablation of water ice is under- does not continue inde¬nitely, however, because

stood to be restricted to the summer months on deep portions of equator-facing scarps cool due to

Mars. partial shielding from the pole-facing slope. The

heat-absorption function f (T , h) includes an h

Equation (7.52) includes three basic processes:

lateral heat conduction (the diffusion equa- term to represent this negative feedback. As h in-

tion), accumulation and ablation (‚h/‚t = T /„ f ), creases in value, the temperature eventually falls

and enhanced heat absorption on equator-facing below the sublimation temperature and accumu-

scarps (‚ T /‚t = f (T , h)). The initiation and deep- lation begins. The heat-absorption term only ap-

plies to equator-facing slopes (i.e. ∇h · r , where r

ening of troughs is assumed to take place in two ˆ

steps illustrated in Figure 7.19a. In the ¬rst step, is the radial unit vector from the pole).

Equation (7.52) includes four parameters: κ, „i ,

slope steepening and dust exposure are initiated

„ f , and T 0 . The thermal diffusivity of ice is κ =

from an undissected ice surface through a pos-

35 m2 /yr. To estimate the vertical ablation rate

itive feedback between ice-surface temperature,

slope gradient, and albedo. In this step, the ice- on sun-facing slopes, Howard (1978) used Viking

surface temperature is assumed to increase from observations of summertime atmospheric water-

vapor concentration to obtain 10’4 m/yr, imply-

its threshold value for sublimation, T 0 , to a max-

imum value of 1 over a time scale „i . This in- ing „ f = 5 Myr for an average trough of 500 m

stability is represented mathematically in f (T , h) depth. The time of scarp initiation is uncertain,

by a cubic polynomial that generates a nega- but may be small compared with trough deep-

ening. Here we assume „i = 0.05, „ f = 0.25 Myr,

tive feedback if T is below the sublimation tem-

which implies that ≈ 25 m of ablation must

perature, T 0 , and a positive feedback for larger

7.7 SPIRAL TROUGHS ON THE MARTIAN POLAR ICE CAPS 185

occur before an incipient scarp is suf¬ciently cation to Mars, here we focus on the excitable

steepened for subsurface dust to be exposed regime.

and the maximum ice-surface temperature to be The solution of Eq. (7.52) in 1D starting with a

narrow pulse of localized sublimation (T = 0.25

reached. The model behavior is not sensitive to

the precise value of „i as long as it is small com- in the pulse, T = 0 elsewhere) is given in Figure

pared with „ f . The threshold temperature for sub- 7.19b. The pulse initiates a self-sustaining train

limation, relative to the equilibrium and maxi- of waves that propagate with the trough trailing

mum temperatures, is between 0 and 1 and varies the temperature wave. The width of the troughs

√

as a function of latitude. At 75—¦ latitude, the max- is given by » ≈ κ„ f , or 13 km for the parameters

imum summer temperature of a nearly ¬‚at, high- characterizing the Martian polar ice caps. This

albedo slope is about ’10 —¦ C (Howard, 1978). The value is in good agreement with observed trough

temperature difference between equator-facing widths on Mars. This behavior contrasts with that

troughs and nearby ¬‚ats is about 40 —¦ C (Howard, of previous models, in which troughs widened

1978). Using these values, T 0 , or 0 —¦ C, is equal to inde¬nitely without reaching a steady state. An-

0.2. We assume this value applies to the entire other similarity with the observed topography

ice cap for simplicity. Equation (7.52) can be fur- of spiral troughs is trough asymmetry: both in

ther simpli¬ed by scaling time to the value of „ f . the model and on Mars, equator-facing scarps are

With this rescaling, the model is reduced to three about twice as steep as pole-facing scarps. Code

independent parameters: κ = 175 km2 , „i = 0.05, for implementing the 1D spiral trough model is

and T 0 = 0.2. given in Appendix 6.

Equation (7.52), without the directional de- Spiral formation in Eq. (7.52) is illustrated

pendence in f (T , h), is one example of the in Figures 7.20b--7.20e for a circular region

FitzHugh--Nagumo (FHN) equations, a well- with a permanently-frozen region near the pole

studied class of equations that describe excitable and starting from random initial conditions. A

media. The FHN equations have been applied to shaded-relief image of a polar-stereographic DEM

predator-prey systems (Murray, 1990), electrical of MOLA topography is given in Figure 7.20a for

conduction in nerve ¬bers (Winfree, 1987), chem- comparison. In these calculations we have used

ical oscillators (Fife, 1976), and electrical trans- the Barkley (1991) approximations to Eq. (7.52),

mission lines (Nagumo et al., 1965). The FHN equa- which use a semi-implicit algorithm for the cu-

bic function in f (T , h) in order to gain computa-

tions come in several forms that all reproduce

solitary waves in 2D (i.e. h(x)) and spiral waves tional ef¬ciency and investigate spiral evolution

in 3D (h(x, y)) for a range of parameter values. over long time scales. The early-time evolution

Excitable media are composed of two variables: of the model is characterized by the initiation,

an ˜˜activator™™ (T in Eq. (7.52)) and an ˜˜inhibitor™™ growth and merging of spirals. It should be noted

(h in Eq. (7.52)). Triggering or excitation of the that heat conduction in the neighborhood of an

medium occurs when a positive feedback is ini- equator-facing scarp may facilitate the sublima-

tiated above a threshold value of the activator. tion of nearby regions, even if those regions face

Local triggering may initiate the triggering of ad- toward the pole. The late-time evolution of the

jacent zones through diffusive spreading of the model is characterized by continued spiral merg-

activator. The value of the inhibitor increases dur- ing, alignment of troughs to face the equator,

ing excitation and eventually returns the system and a gradual increase in trough spacing as ˜˜de-

to the unexcited state through the negative cou- fects™™ are removed from the system. Mature spi-

pling between h and T in Eq. (7.52). The FHN equa- rals face the equator at low latitudes but curve

tions exhibit both excitable behavior, character- towards the pole at higher latitudes, just as those

ized by wide swings in h and T when „i „f, on Mars do. In the model, this poleward orienta-

and oscillatory behavior, characterized by smaller tion is required for the maintenance of steady-

¬‚uctuations when „i ≈ „ f (Gong and Christini, state, rigidly rotating spiral forms. By developing

2003). Although both regimes may have appli- poleward orientations, high-latitude scarps face

186 INSTABILITIES

(c)

(b)

(a)

t = 100

t = 10

(d) (e)

200 km

termination

bifurcation

gullwing

gullwing

wave-like undulations

t = 1000 t = 2000

Fig 7.20 Numerical model results and north-polar

troughs, bifurcations, and terminations (Figures

topography. (a) Shaded-relief image of Martian north-polar ice

7.20a, 7.20c, and 7.20d). The steady state of the

cap DEM constructed using MOLA topography. The

model after t = 2000 is a set of seven rigidly-

large-scale closeup indicates examples of gullwing-shaped

troughs, bifurcations, and terminations. Highest elevations are rotating spirals. The sense of spiral orientation

red and lowest elevations are green. (b)“(e) Shaded-relief has an equal probability of being clockwise or

images of the model topography, ’h, for (b) t = 10, (c)

counter-clockwise in the model depending on

t = 100, (d) t = 1000, and (e) t = 2000 starting from

the random number seed used to generate the

random initial conditions. The model parameters are

initial conditions. The spiral orientation of Fig-

L = 250 (number of pixels in each direction), x = 0.4,

ure 7.20e is opposite to those on the north pole

T0 = 0.3, „i = 0.05, „ f = 1. In the Barkley approximations to

of Mars, but 50% of random initial conditions

Eq. (7.52), T0 is combined into two parameters, a = 0.75 and

b = 0.01. The model evolution is characterized by spiral

merging and alignment in the equator-facing direction. A

steady state is eventually reached with uniformly rotating

spirals oriented clockwise or counter-clockwise depending

on the initial conditions. For color version, see plate section.

Modi¬ed from Pelletier (2004b).

perpendicular to the equator, thereby reducing

their ablation rates and decreasing their migra-

tion speeds. The decreased migration speed off-

sets the smaller distance required to complete

a revolution at high latitudes, enabling spirals

to rotate at a uniform angular velocity and pre-

serve their steady-state form. The model also re-

Fig 7.21 Vegetation bands of Ralston Playa, Nevada. Dirt

produces a number of secondary features of po-

road provides scale. Individual bands are 1“3 m wide.

lar troughs on Mars, including gullwing-shaped

EXERCISE 187

lead to spirals that rotate in the same sense of the western US and are comprised of brushy, salt-

tolerant vegetation. These vegetation bands most

as those on Mars (i.e. clockwise away from the

likely form by a positive feedback in which vege-

pole). This model behavior suggests that a uni-

tation growth promotes mounding of the surface

form sense of spiraling need not represent an

beneath the vegetation. Mounds trap occasional

asymmetrical process of trough evolution. In-

runoff into the playa, which promotes more vegeta-

stead, it could be the result of a minor asymme-

tion growth. Bands evolve because a line of shrubs

try in the initial pattern of sublimation, as in the

oriented along a contour pond the most water for

model. the least biomass. Develop a model for the forma-

tion of these bands based on a simpli¬ed model for

the positive feedback between vegetation, mound-

ing, and ponding. Actual vegetation bands follow

Exercise contours closely and are more closely spaced on

more-steeply dipping portions of the playa. Is your

7.1 Figure 7.21 illustrates the curious phenomenon of

model consistent with these observations?

vegetation bands. These bands form on some playas

Chapter 8

Stochastic processes

observed statistical distributions of rainfall

8.1 Introduction events, for example, can be used to generate syn-

thetic storms that match the statistical behav-

ior of actual storms. In cases where there is sig-

Thus far we have considered models in which

ni¬cant uncertainty in model input parameters,

the future behavior of the system can be de-

Monte Carlo methods can be used to determine a

termined using equations and boundary condi-

range of model outputs corresponding to a range

tions known at some initial time. In many Earth

of input parameters. Monte Carlo methods are

surface systems, however, the future behavior of

based on deterministic models; they work by run-

the system cannot be predicted with certainty,

ning many successive examples of a deterministic

either because the system behavior is sensitive

model using different input parameters sampled

to small-scale processes that cannot be fully re-

from prescribed probability distributions.

solved and/or because there is signi¬cant uncer-

tainty in model input parameters. The climate

system is a good example of a system that de-

8.2 Time series analysis and

pends on small-scale processes (i.e. turbulence)

that cannot be fully resolved in any numerical

fractional Gaussian noises

model. As a result, climate and weather models

are inherently limited in their ability to predict

Climatological and hydrological time series pro-

the details of future climate and weather pat-

vide many good examples of phenomena that re-

terns. Soil permeability is a good example of a

quire a stochastic modeling approach. Time se-

model input parameter with signi¬cant uncer-

ries data for rainfall and discharge at a point in

tainty. Soil permeability depends on the detailed

space are the result of complex processes, includ-

structure and composition of the soil at a range

ing convective instabilities in the atmosphere. Al-

of spatial scales, making an exact determination

though some elements of the climate system can

of permeability very dif¬cult over large areas. As

be quanti¬ed deterministically, convective insta-

such, numerical models that require soil perme-

bilities are very dif¬cult to model more than a

ability as an input (e.g. models for runoff, in¬l-

few days into the future. The occurrence of con-

tration, aquifer recharge, etc.) are limited in their

vective instabilities is one of the fundamental

precision no matter how ¬nely they resolve the

limitations on accuracy of weather forecasts.

underlying physics of the problem.

Figure 8.1 plots the daily mean discharge

In some applications where deterministic

in the Gila River near the town of Winkel-

models cannot predict the future behavior of

man, Arizona from 1942 to 1980. Discharge in

a system precisely, stochastic models can be

this case is plotted on a logarithmic scale be-

useful for understanding the range of possible

cause of its highly skewed distribution, with peak

system behaviors. Stochastic models based on

8.2 TIME SERIES ANALYSIS AND FRACTIONAL GAUSSIAN NOISES 189

annual peak

1000

Q

10

(m3/s)

S(f )

100

S(f ) ∝ f

1

10

1 0.1

t 8000 12000

4000

0 100

10

10 10 10

(days, 1942“1980)

f (1/day)

Fig 8.1 Stream¬‚ow of the Gila River near Winkelman,

Fig 8.2 Power spectrum of the Gila River discharge plotted

Arizona, from 1942“1980.

in Figure 8.1.

values many times larger than average values. because they are controlled by a single meteoro-

Many kinds of datasets that arise in climatolog- logical event. At longer time scales, a seasonal

ical and Earth surface studies involve spatial or cycle can be expected and is apparent in Fig-

temporal series data. Figure 8.1 is a time series ure 8.1. At still longer time scales, interannual

because it is a quantity that varies as a function and decadal variations in ocean temperatures (re-

of time. Spatial series data, such as topographic sulting from known periodic processes such as

pro¬les, can be analyzed using many of the same ENSO and also aperiodic variations) lead to clus-

tools used to analyze time series data. Time se- ters of wet years and clusters of dry years. Hurst

ries data can be characterized by their probabil- et al. (1965) quanti¬ed this multi-scale autocor-

ity distributions and autocorrelation functions. relation pattern in hydrological data sets, now

The probability distribution of a dataset quanti- called the Hurst phenomenon. In this section we

¬es its mean value and its variability from the explore a particular type of stochastic model, i.e.

mean. The Gila River near Winkelman is charac- the fractional Gaussian noise, that reproduces

terized by a mean stream¬‚ow value of approxi- the autocorrelation structure common to many

mately 100 m3 /s. The mean value is not very rep- hydrologic and climatic time series.

resentative of the typical stream¬‚ow in this case, The autocorrelation structure of a time series

however, because the river also has extended dry can be quanti¬ed in a number of ways, but the

periods with discharge values less than 10 m3 /s power spectrum is one of the simplest. The power

and occasional large storm events that trigger spectrum is de¬ned as the square of the coef¬-

discharges greater than 1000 m3 /s. This river, like cients in a Fourier series representation of the

most in the southwestern US, has a large variance time series. The power spectrum quanti¬es the

in stream¬‚ow values due to the episodic nature variance of a time series at different frequencies.

of precipitation in semi-arid environments. The Figure 8.2 plots the power spectrum of the Gila

autocorrelation function of a time series quan- River stream¬‚ow dataset on logarithmic scales.

ti¬es how strongly points in the series are cor- This spectrum shows a prominent seasonal cycle

related over different time scales. In the case of (i.e. a large peak at a period of 1 year) superim-

river discharge, signi¬cant correlations can be ex- posed on a broad spectrum that decreases with

pected over many time scales due to the nature larger frequencies. Aside from the seasonal cy-

of rainfall and hydrologic response in a large cle, the background spectrum follows a power-

law relationship with frequency given by S( f ) ∝

river basin. During an individual ¬‚ood event

f ’1/2 . This relationship turns out to be nearly

that lasts several days or more stream¬‚ow values

can be expected to have signi¬cant correlation universal for stream¬‚ow datasets worldwide. In a

190 STOCHASTIC PROCESSES

this observed underlying structure. Figure 8.4

10’1

plots examples of three synthetic time series with

power spectra given by S( f ) ∝ f ’β , with β = 0,

1/2, and 1. Fractional Gaussian noises are simply

S(f ) time series that have Gaussian probability distri-

(yr) ’1/2

S(f ) ∝ f butions and power spectra have a power-law de-

pendence on frequency f . As the magnitude of

the power-spectral exponent increases, the cor-

relations between adjacent points in the series

become stronger. Fractional Gaussian noises can

10’2

be generated with the Fourier-¬ltering technique.

This technique proceeds in several steps. First, a

10’1

10’3 10’2 100 101

Gaussian white noise sequence is generated. Sec-

f (1/yr)

ond, a discrete Fourier transform of the Gaus-

Fig 8.3 Average of power spectra of monthly average sian white noise sequence is computed. The re-

sulting Fourier spectrum will be ¬‚at, i.e. β = 0

discharge data from 636 gage stations with the annual

periodicity ¬rst removed.

in the expression S( f ) ∝ f ’β . Except for the sta-

tistical scatter, the amplitudes of the Fourier co-

ef¬cients, |Ym | will be equal. Third, the Fourier

study of 636 gage stations, Pelletier and Turcotte coef¬cients are multiplied by a power-law func-

(1997) found that, when the annual peak was tion of frequency:

removed, the average power spectrum followed

’β

m 2

Ym =

a nearly straight line on log--log scales with a Ym (8.1)

N

slope of ’0.50 (Figure 8.3). The data for this anal-

The exponent β/2 is used because the power

ysis were chosen from all complete records with

durations greater than or equal to 512 months, spectrum is proportional to the square of the

including 636 records in the analysis. Pelletier Fourier coef¬cients. Fourth, an inverse discrete

and Turcotte (1997) took advantage of the large Fourier transform is computed using the new

number of available stations to investigate the coef¬cients Ym . This technique is called Fourier

possible regional variability of the power spec- ¬ltering because Eq. (8.1) has the effect of re-

tra. All of the regions exhibited the same power- moving, or ¬ltering out, some of the variability

at high frequencies (if β > 0) or low frequencies

law power-spectral form with an average expo-

nent of ’0.52 and a standard deviation of 0.03, (if β < 0). Two-dimensional fractional noises can

indicating little variation. Proxy data for hydro- also be constructed (Figure 8.5). These noises have

logic datasets (e.g. tree ring data in semi-arid en- applications to stochastic models of topography,

vironments) also exhibit that same spectrum out soil moisture, and the porosity structure of sed-

to time scales greater than 1000 years. Pelletier imentary basins. The details of this technique,

and Turcotte (1997) also showed similar power- together with codes for its implementation, are

spectral behavior in precipitation time series at given in Appendix 7.

time scales greater than 10 years, indicating that Fractional noises have many applications in

the correlation structure in stream¬‚ow data are geosciences (Pelletier and Turcote, 1999). For ex-

primarily a result of correlations in the climate ample, generating synthetic fractional noises can

system on decadal time scales and greater, while help to assess the likelihood of future hydrolog-

shorter-term correlations are at least partially the ical events that depend on the cumulative dis-

result of surface and subsurface ¬‚ow pathways. charge over a series of years. A drought of a given

Thus far we have shown that hydrological duration, for example, is the result of several con-

data sets have a common underlying autocorre- secutive years of below-average ¬‚ow. The correla-

lation structure. It is useful for many applica- tion structure of hydrological time series plays a

tions to generate synthetic data sets that match very signi¬cant role in the frequency of drought

8.3 LANGEVIN EQUATIONS 191

4.0

b = 1/2

b = 0 (white noise) 4.0

2.0

2.0

0.0

0.0

’2.0 ’2.0

(a) (b)

1000 2000

0

1000 2000

0

b=2

4.0 b = 1

2.0

2.0

1.0

0.0 0.0

’1.0

’2.0

(d)

(c)

1000 2000

0

1000 2000

0

on a stochastic diffusion model of convective

Fig 8.4 Examples of fractional Gaussian noises with β = 0,

transport. Fluctuations in precipitation, river dis-

1/2, 1, and 2.

charge, and tree ring growth are related in com-

plex but direct ways to variations in temperature

occurrence. If hydrological time series were com- and/or humidity.

pletely uncorrelated, then the probability of 10 The simplest approach to the turbulent trans-

years of below-average stream¬‚ow would be 1/210 port of heat and water vapor vertically in the at-

or 1/1024. However, the fact that hydrologic data mosphere is to assume that the ¬‚ux of heat or

sets have positive correlation over a range of time water vapor concentration is proportional to its

scales makes the probability of droughts, espe- gradient. The mixing of heat energy and humid-

cially long droughts, far more frequent (Pelletier ity will then be governed by the diffusion equa-

and Turcotte, 1997). tion. Experimental and theoretical observations

support the hypothesis that diffusive transfer ac-

curately models vertical transport in the atmo-

8.3 Langevin equations sphere. Kida (1983) has traced the dispersion of

air parcels in a hemispheric GCM and found the

dispersion of those air parcels to obey the diffu-

If a broad range of climatological and hydrolog-

sion equation in the troposphere. Also, Hofmann

ical time series data have power spectra of the

form S( f ) ∝ f ’1/2 from time scales of decades and Rosen (1987) argued that turbulent diffusion

of aerosols from the El Chichon volcanic erup-

to centuries, we must look for some general un-

tion was consistent with the exponential decay

derlying climatological mechanism. In this sec-

of aerosol concentrations with time (i.e. turbu-

tion we present one possible model of ¬‚uctua-

lent diffusion and gravitational settling combine

tions in local temperature and water vapor based

192 STOCHASTIC PROCESSES

(a) (b)

b=0 b=1

(c) (d)

b=2 b=3

Fig 8.5 Examples of 2D fractional Gaussian noises with β =

(a) 0, (b) 1, (c) 2, and (d) 3.

is given by

c in+1 ’ c in ∝ c i+1 ’ 2c in + c i’1

n n

(8.2)

to predict exponential concentrations in time).

Vertical atmospheric turbulent diffusion is super- where i is an index for space and n is an in-

imposed, as it is in oceanic diffusion, on large- dex for time. In our model we establish a one-

scale Hadley and Walker circulations (Peixoto and dimensional lattice of 32 sites with periodic

Oort, 1992). boundary conditions at the ends of the lattice.

To see how time series with power-law power At the beginning of the simulation, we place

spectra arise from a stochastic diffusion pro- ten particles on each site of the lattice. At each

cess, we explore the behavior of a discrete, one- timestep, a particle is chosen at random and

dimensional stochastic diffusion process. A dis- moved to the left with probability 1/2 and to the

crete version of the diffusion equation for the right if it does not move to the left. In this way,

density of particles, c, on a one-dimensional grid the average rate at which particles leave a grid

8.4 RANDOM WALKS 193

point is proportional to the number of particles

20

at the grid point. The average rate at which par-

ticles enter a grid point is proportional to the

10

number of particles on each side multiplied by

one half (since the particles to the left and right

of site i move into site i only half of the time). 0

This is a stochastic model satisfying Eq. (8.2). The

probabilistic nature of this model causes ¬‚uctu- 10

ations to occur in the local density of particles

that do not occur in a deterministic model of dif- 20

fusion. This simple model produces ¬‚uctuations

40 60

20

0 100

80

in the number of particles in the central site of

the 32-site lattice that have a power-law power

Fig 8.6 Five examples of a simple random walk process.

spectra with β = 1/2 (Pelletier and Turcotte, Plotted in Gnuplot by mdoege@compuserve.com.

1997). Reproduced from http://en.Wikipedia.org/wiki/Image:

A stochastic diffusion process can be studied Random Walk example.png.

analytically by adding a noise term to the ¬‚ux of

The power spectrum of variations in E (t),

a deterministic diffusion equation (van Kampen,

S E (ω) =< |E (ω)|2 > is

1981):

∞

sin2 (kL ) 1

dk ∝ ω’ 2

S E (ω) ∝ (8.7)

‚T ‚J 2 k 4 + ω2

D

ρc =’ ’∞

(8.3)

‚t ‚x

for low frequencies. Since T ∝ E , S T (ω) ∝

ω’1/2 also. This is the same result as the discrete

where the ¬‚ux, J , is given by

model.

‚T

J = ’σ + ·(x, t) (8.4)

‚x

8.4 Random walks

The variable T represents the ¬‚uctuations in

temperature from equilibrium and ·(x, t) is a The simplest type of random walk is a stochastic

process consisting of a sequence of discrete steps

Gaussian, white noise. Equations (8.3)--(8.4) are

of ¬xed length. Figure 8.6 illustrates ¬ve exam-

one example of a Langevin equation (i.e. a partial

ples of a random walk which begins at zero and

differential equation with stochastic forcing).

involves unit jumps along the y axis up or down

We will use Eqs. (8.3)--(8.4) to calculate the

with equal probability for each step along the x

power spectrum of temperature ¬‚uctuations in

axis. Each realization of the random walk is dif-

a layer of width 2L exchanging heat with an in¬-

ferent, but the fact that the distance from the

nite, one-dimensional, homogeneous space (Voss

origin tends to increase with each step in all of

and Clarke, 1976). The Fourier transform of the

the examples gives the random walk a certain ele-

heat ¬‚ux of the stochastic diffusion equation

ment of predictability. What is remarkable about

is

random walks and other stochastic processes is

iω·(k, ω) that the average behavior of a large collection of

J (k, ω) = (8.5)

D k 2 + iω random walks is entirely predictable. In a simple

random walk, the expected distance from the ori-

The rate of change of heat energy in the layer

gin, | y| is proportional to the square root of the

will be given by the difference in heat ¬‚ux

number of steps along the x axis, x:

out of the boundaries, located at ±L: dE (t)/dt = √

J (L , t) ’ J (’L , t). The Fourier transform of E (t) is < | y| >∝ x (8.8)

then

By expected, we mean that this is the average be-

∞

havior obtained by many realizations of the ran-

1

E (ω) = sin(kL )J (k, ω) dk (8.6)

1

dom walk process.

(2π) 2 ω ’∞

194 STOCHASTIC PROCESSES

10 (c)

(a) (b)

RT

(g/cm2/yr)

10

R T

RT

R

10

10

Fortymile Wash fan Whipple Mtns. fan (clay)

Coyote Mtns. fan

10

(f)

10 (d) (e)

RT

RT

R T 0.64

(g/cm2/yr)

10

R

10

10

Silver Lake fan (lower)

Yucca Wash fan Silver Lake fan (upper)

10

106 102 103

106 102 103 106

104

104 105

105

102 103 104 105

T (yr)

T (yr)

T (yr)

dust accumulation on alluvial fan terraces. Low

Fig 8.7 Plots of silt accumulation rates versus time interval

rates of long-term (105 --106 yr) dust accumulation

on alluvial fan terrace surfaces from data in Reheis et al.

(1995): (a) Fortymile Wash fan, (b) Coyote Mtns. fan, (c) in the southwestern US, for example, have been

Whipple Mtns. fan (clay fraction shown), (d) Yucca Wash fan, interpreted as a direct result of low dust de-

(e) lower fan near Silver Lake Playa, and (f) upper fan near position rates during the predominantly cool,

Silver Lake Playa. Best-¬t lines are also shown indicating

wet Pleistocene, whereas high rates of shorter-

power-law scaling with exponents close to ’0.5 for (a)“(d).

term (103 --104 yr) dust accumulation have been

Trends in (e) and (f) are more consistent with linear

interpreted as an early Holocene pulse of dust

accumulation or cyclical-climate models. Modi¬ed from

deposition related to pluvial-lake dessication

Pelletier (2007a).

(McFadden et al., 1986; Chadwick and Davis, 1990;

Reheis et al., 1995). In this section we compare the

observed eolian dust accumulation rates at six al-

8.5 Unsteady erosion and luvial fan study sites with the predictions of a

climatically controlled deterministic model and

deposition in eolian

a simple random-walk model of eolian dust ac-

environments cumulation. The results suggest that eolian ero-

sion and deposition on alluvial fans is a highly

episodic process that closely resembles a random

The accumulation of dust on alluvial fan terraces

walk.

in desert environments provides a good example

Reheis et al. (1995) synthesized data on eo-

of the application of a simple random walk in

lian dust accumulation from alluvial fan terrace

geomorphology. In this example, the y axis in

study sites in the southwestern US with opti-

Figure 8.7 represents the accumulation of wind-

mal age control. These sites are located on gen-

blown dust (as a thickness or mass per unit area)

tly sloping, planar alluvial fan terraces that have

underneath the desert pavement, and the x axis

not been subject to ¬‚ooding since abandonment

represents time.

by fan-head entrenchment. The time since fan-