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