. 10
( 14)


acterized by alternating zones of narrow, deeply
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
»=π (7.36) of the model and was observed for a wide range
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
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.

drumlins often parallels the topographic surface,
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
example, why are drumlins so variable in size
and shape, even over distances as short as a few
kilometers (Figure 7.12)? Hindmarsh (1998) and
100 km
Fowler (2001) proposed that shearing of an in-
Q Bi
compressible viscous till layer could be responsi-
ble for producing drumlins even in the absence
(b) NA
of a bedrock core. Schoof (2002), however, has
5 V
shown that bedforms produced by these models
V1 Ce
are not consistent with many aspects of natural
l S
Most conceptual drumlin models invoke ei-
ther subglacial hydraulic processes or deform-
1 linear stability
ing bed processes. This section seeks to unite
l = p (w 0h 0/S 0)1/2
these two hypotheses by modeling sediment de-
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

(a) (b)

Rogen Moraine
Rogen Moraine

20 km
20 km
Columbia 43.6° N


Ontario Cayuga
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)
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-

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
∇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-
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
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

(a) (b) (c)
deformable, flat ice sheet sloping ice sheet
ice sheet permeable

2D till

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

sit higher for the same reason. In the ¬nal stages
(a) of the model, high-porosity regions conduct the
40 — vert. remaining porewater out of the upper boundary
l = 0.2
(Figure 7.14c). Deformation stops when all of the
liquid has been squeezed from the matrix.
The time required for this instability to oc-
middle contraction
cur in natural drumlins can be estimated as
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
10’6 m/s; Freeze and Cherry, 1979) and l = 50 m,
T ≈ 1000 yr. This estimate depends sensitively on
the assumed value for conductivity; drumlins in
coarse-grained sediments could develop in less
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

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

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

(a) (b)

1 w

shaded relief


1 w
1 w

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.

(a) (c)
43.4° N
43.6° N



drumlin width
89.2° W
76.0 88.8
77.2° W
0 500 m
0 1000 m
(b) (d)

River Valley
60 m
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

Till thicknesses were mapped using depth-to- topographic surface
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-
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,
¬‚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,

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):
= κ∇ 2 T + f (T , h) (7.49)
The spiral troughs on the Martian polar ice caps
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

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
(dust exposure) layers are exposed. During trough
deepening, the value of h steadily
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
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
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
side of the domain. Modi¬ed from
Pelletier (2004b).
30 40
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

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


t = 100
t = 10

(d) (e)
200 km


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

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

annual peak

S(f )
S(f ) ∝ f


1 0.1
t 8000 12000
0 100
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

this observed underlying structure. Figure 8.4
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
be generated with the Fourier-¬ltering technique.
This technique proceeds in several steps. First, a
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)
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

b = 1/2
b = 0 (white noise) 4.0



’2.0 ’2.0
(a) (b)
1000 2000
1000 2000
4.0 b = 1


0.0 0.0

1000 2000
1000 2000
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

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

point is proportional to the number of particles
at the grid point. The average rate at which par-
ticles enter a grid point is proportional to the
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
0 100
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

sin2 (kL ) 1
dk ∝ ω’ 2
S E (ω) ∝ (8.7)
‚T ‚J 2 k 4 + ω2
ρc =’ ’∞
‚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
J = ’σ + ·(x, t) (8.4)
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
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)
By expected, we mean that this is the average be-

havior obtained by many realizations of the ran-
E (ω) = sin(kL )J (k, ω) dk (8.6)
dom walk process.
(2π) 2 ω ’∞

10 (c)
(a) (b)




Fortymile Wash fan Whipple Mtns. fan (clay)
Coyote Mtns. fan
10 (d) (e)
R T 0.64




Silver Lake fan (lower)
Yucca Wash fan Silver Lake fan (upper)
106 102 103
106 102 103 106
104 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
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-


. 10
( 14)