<< . .

. 34
( : 51)

. . >>

(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

<< . .

. 34
( : 51)

. . >>