<< . .

. 25
( : 51)



. . >>

600
0 200 400
200
100 400 2000 6000 m
0 x (km)
spatial distribution of glaciated regions in the
Fig 5.9 (a) Grayscale map of topography in the central
western United States. Porter et al. (1983) mapped
Andes indicating sediment-rich, strongly over¬lled basins in
the northern central Andes (near 20—¦ S) and areas of Pliocene--Quaternary glaciation in the
western United States, based on the elevations
sediment-starved, under¬lled and weakly over¬lled basins
(near 30—¦ S) corresponding to a humid-to-arid transition. of the lowest cirques. Their map identi¬es the
(b) Basin topographic pro¬les extracted from SRTM data. approximate local elevation of vigorous glacial
Also plotted using logarithmic scale (with base-level elevation
erosion. To construct a map of the extent of late
subtracted) to illustrate exponential trend and associated
Cenozoic glacial erosion, the contours of glacial-
best-¬t κ values. (c) Best-¬t model pro¬les with associated
limit elevation from Porter et al. (1983) were ¬rst
Q s values. Sediment supply decreases from north to south
digitized, projected to Lambert Equal-Area, and
with precipitation rates and as deformation becomes more
then interpolated to obtain a grid of glacial-
spatially distributed. Note that 100 m has been added to the
limit elevations. Second, a 1 km resolution DEM
northernmost pro¬le so that the plots can more easily be
(digital elevation model) of the western United
distinguished. Modi¬ed from Pelletier (2007). Reproduced
with permission of the Geological Society of America. States was recti¬ed to the glacial-limit-elevation
grid. Third, the glacial-limit-elevation and DEM
grids were compared pixel-by-pixel to determine
the width of the glaciated region, uplift will be whether local elevations were above or below
limited by ¬‚exural bending stresses. In contrast, the glacial-limit elevation. Pixels with elevations
if the ¬‚exural wavelength and the width of the above the glacial limit were considered areas of
eroded region are comparable, bending stresses glacial erosion and were given a value of 1. Con-
will be smaller and uplift will approach the Airy versely, pixels with elevations below the glacial
isostatic limit in which 80% of the eroded mate- limit were considered not to have been subjected
rial is replaced through rock uplift. to glacial erosion and were given a value of 0.
To map the compensation resulting from The resulting binary grid is a map of relative
glacial erosion, we ¬rst must determine the glacial erosion; an underlying assumption is that
122 FLEXURAL ISOSTASY


to gravity, and q(x, y) is the vertical load (Watts,
117° W 103° W
49° N
2001). In order to account for spatial variations
300 km in elastic thickness we retained D (x, y) inside the
parentheses in Eq. (5.41). This results in addi-
tional terms involving gradients of D that are
not present when D is assumed to be spatially
uniform.
42° N
The binary grids of glacial erosion were di-
YP-A
rectly input into the ¬‚exure equation as the
loading term, so that q(x, y) = 1 where erosion
occurred and q(x, y) = 0 where erosion did not
occur. The solution to the ¬‚exure equation with
SN
these loading terms is the ratio of the positive de-
¬‚ection of the lithosphere to the depth of glacial
SJ erosion. This is the compensation resulting from
glacial erosion. Here we use the alternating dif-
ference implicit (ADI) technique to solve Eq. (5.41),
following the work of van Wees and Cloetingh
109° W
Te (1994). Previous techniques we have discussed
0 100 km
(e.g. Fourier ¬ltering) cannot be used when D
Fig 5.10 Grayscale map of elastic thickness in the western varies spatially.
United States. Data from Lowry et al. (2000) and Watts For the analysis, we will use Lowry et al.™s
(2001). Boundaries of glaciated areas are shown in white.
(2000) data set for elastic thickness (T e ) in the
Labeled ranges are: SN (Sierra Nevada), SJ (San Juan), and
western United States. These data are based on
YP-A (Yellowstone Plateau-Absaroka). Modi¬ed from
the coherence method and provide a direct mea-
Pelletier (2004a).
sure of lithospheric strength that includes the
weakening effects of faults. These data have the
erosion is uniform within the glaciated regions. highest resolution available for the region, but do
The boundaries of the glaciated regions obtained not completely cover the western United States.
in this way are mapped in white in Figure 5.10 To cover the entire area, the Lowry et al. data was
(which also shows state boundaries for reference). augmented with the global compilation of Watts
Also shown in Figure 5.10 is a grayscale map of (2001), and the data sets joined smoothly with a
the elastic thickness in the western United States. low-pass ¬lter. The resulting T e map was recti¬ed
Ranges along the Paci¬c coast were not included to the glacial erosion maps and is shown as a
in Porter et al. (1983) and have also been ex- grayscale map in Figure 5.10. T e values were con-
cluded from this analysis. In addition, the north- verted to a grid of ¬‚exural rigidity, D (x, y), for
ern Cascades were excluded because they were use in Eq. (5.41), by means of the relationship
completely covered by the Cordilleran Ice Sheet
E T e3
and hence subject to a different style of glacia- D= (5.42)
12(1 ’ ν 2 )
tion than other ranges of the western United
States. assuming E = 70 GPa and ν = 0.25.
The de¬‚ection of a thin elastic plate with vari- A map of glacial-erosion compensation in the
able elastic thickness under a vertical load is western United States is given in Figure 5.11a.
given by Also shown in Figure 5.11b is the compensa-
tion determined by assuming that glacial ero-
∇ 2 (D (x, y)∇ 2 w) + ρgw = q(x, y) (5.41)
sion was concentrated near the equilibrium line,
where w is the de¬‚ection, D (x, y) is the ¬‚exural not uniformly within the glacial region as in Fig-
rigidity, ρ is the density contrast between the ure 5.11a. Although compensation is strictly de-
crust and the mantle, g is the acceleration due ¬ned to be between 0 and 1, we have mapped
EXERCISES 123


103° W
103 W 117° W
117° W 49 N
49 N
(a) (b)
300 km




42 N




uniform erosion above glacial-limit erosion concentrated near glacial-limit
109° W
109° W
0 compensation 1
Fig 5.11 Map of compensation of glacial erosion assuming Plateau--Absaroka Range. The high values in this
(a) uniform erosion within glaciated areas and (b) erosion region re¬‚ect the extensiveness of the glaciated
concentrated within 500 m altitude of glacial limit. For color
area, which is nearly 200 km in width. In ad-
version, see plate section. Modi¬ed from Pelletier (2004a).
dition, however, nearby ranges are also being
eroded, reducing the bending stresses that would
otherwise be concentrated at the edges of the
the full solution to the ¬‚exure equation in Fig-
Yellowstone area, limiting its potential uplift. In-
ure 5.11, including areas of minor subsidence (in
stead, nearby ranges, including the Beartooth,
black). All of the major ranges in the western
Gallatin, Wyoming, and Wind River Ranges, ac-
United States have compensation values greater
commodate the bending stresses to produce a
than or equal to 0.4, indicating that nearly half
˜˜bulls-eye™™ over the central part of the area -- the
of all glacially eroded mass is returned to the
Yellowstone Plateau--Absaroka Range.
region™s ranges as erosionally driven rock uplift.
The pattern of compensation varies little between
Figures 5.11a and 5.11b because ¬‚exure is insen-
sitive to small-scale variations in loading.
Exercises
The greatest difference between Figures 5.11a
and 5.11b is the mean values: values of compen- 5.1 Calculate the degree of isostatic compensation
for continental topography with a wavelength of
sation obtained by assuming concentrated ero-
20 km and an elastic thickness of 5 km, assum-
sion near the glacial limit are 20% lower than
ing E = 70 GPa, nu = 0.25, ρm = 3300 kg/m3 , and
for uniform erosion. These lower values re¬‚ect
ρc = 2700 kg/m3 .
the additional bending stresses present in the
5.2 In addition to driving rock uplift, the isostatic re-
case of localized erosion; i.e., for the same eroded
sponse to erosion also causes subsidence in some
mass, a circular area (i.e., uniform erosion within
locations. Consider a 1D model of periodic line un-
the range) results in less intense bending stresses loading (Figure 5.12). Each line unloading might
than a ˜˜ring-shaped™™ area (i.e., localized erosion represent glacial erosion concentrated in the high
near the glacial limit). peaks of a Basin and Range-style topography. Calcu-
The highest compensation values in the west- late the distance between unloading, L , that pro-
ern United States are found in the Yellowstone duces the largest value of subsidence in spaces
124 FLEXURAL ISOSTASY



w(x)
0
x
L LIS (9 ka)

Fig 5.12 Schematic diagram of the 1D periodic line loading
calculation used in Exercise 5.2.

Lake Agassiz
between the locations where unloading is concen-
trated.
5.3 Calculate the lithospheric de¬‚ection associated
with a narrow seamount of radius 10 km and
height 1 km using the Kelvin function solution.
5.4 Use the Fourier series approach to construct an
expression for the ¬‚exural response of the litho- Fig 5.13 Outline of the Laurentide Ice Sheet at 9 ka,
sphere to a periodic topographic load of trian- showing the extents of proglacial lakes at the ice sheet
gular mountain ranges with width L and height margins.
H.
5.5 Proglacial lakes formed at the margins of the
this model, assuming that the lake ¬lls up to sea
Laurentide Ice Sheet. The largest of these, Lake
level?
Agassiz, had a width of several hundred kilome-
5.6 Following the Andean example, download topo-
ters (e.g. Figure 5.13). Model the 1D bathymetric
graphic data for a portion of the Main Central
pro¬le of Lake Agassiz perpendicular to the ice
Thrust of the Himalaya. Using published estimates
margin assuming that the lake formed as a ¬‚ex-
for the elastic thickness of the lithosphere, model
ural depression. Assume that the proglacial to-
the 1D ¬‚exural pro¬le and determine the ap-
pography was initially ¬‚at and at sea level, and
proximate location of the Himalayan forebulge in
that the ice sheet acts as a rectangular load with
India.
height 2 km. How wide is the lake according to
Chapter 6




Non-Newtonian ¬‚ow equations


extend in¬nitely far into and out of the cross sec-
6.1 Introduction tion. The horizontal force F exerted on the left
side of a thin vertical slice of the ¬‚uid of width
x is equal to the total hydrostatic pressure in
Chapter 1 introduced the concept of non-
the ¬‚uid column. The hydrostatic pressure at a
Newtonian ¬‚uids and highlighted their impor-
depth y below the surface is given by ρgy. The
tance for our understanding of the mechanics
total hydrostatic pressure is calculated by inte-
of glaciers. In this chapter we explore methods
grating ρgy over the ¬‚uid column from y = 0 to
for modeling non-Newtonian ¬‚uids with applica-
y = h:
tions to ice sheets, glaciers, lava ¬‚ows, and thrust
sheets. In most cases we will make the simplify-
ρgh 2
h
F (x) = ρgydy =
ing assumption that the rheology of the ¬‚ow is (6.1)
2
0
constant in time. In nature, however, many types
of non-Newtonian ¬‚ows are complicated by the On the right side of the vertical slice (i.e. at po-
sition x + x), the pressure is lower because the
fact that the rheology of the ¬‚ow depends on
parameters that evolve with the ¬‚ow itself. The ¬‚ow is thinner:
viscosity of lava, for example, depends on its tem-
ρg(h ’ ρg 2
h)2
F (x + x) = = (h ’ 2h h
perature, which, in turn, depends on the ¬‚ow
2 2
behavior (e.g. slower moving lava will cool more ρg 2
+ ( h)2 ) ≈ (h ’ 2h h) (6.2)
quickly, increasing its viscosity and further slow- 2
ing ¬‚ow in a positive feedback). Here we focus
The ( h)2 term in Eq. (6.2) is much smaller than
primarily on modeling ¬‚ows with a constant, pre-
the 2h h term for ¬‚ows with small (i.e. < 0.1) sur-
scribed rheology. In cases where the rheology de-
face slopes, and hence the second-order term can
pends strongly on the state of the ¬‚ow the meth-
be neglected in such cases. This ˜˜small-angle™™ ap-
ods introduced in this chapter will serve as the
proximation is commonly used when modeling
foundation for more complex models which fully
ice sheets, glaciers, and lava ¬‚ows on Earth. The
couple the rheology with the ¬‚ow.
difference between the horizontal force exerted
on the opposite sides of the ¬‚uid slice must be
equal to the shear stress „0 exerted over the width
6.2 Modeling non-Newtonian and of the slice, x:
perfectly plastic ¬‚ows
’„0 = F (x + x) ’ F (x) = ρgh h (6.3)

First we consider the forces acting on a ¬‚uid Rearranging Eq. (6.3) gives
spreading over a ¬‚at surface (Figure 6.1). We be-
„0 = ’ρghS
gin with a 1D example, so the ¬‚uid is assumed to (6.4)
126 NON-NEWTONIAN FLOW EQUATIONS



5
”x (a)
h
h4 t0 = 1.5 bar
1.0 bar
(km)
S
F 3
F+ ”F
2
x
1
0.5 bar
Fig 6.1 Schematic diagram of the forces acting on a gravity
¬‚ow atop a ¬‚at surface. 0
4
(b)
1.5 bar
where S is equal to the local ice-surface gradient, 3
h
h/ x. Field studies of modern ice sheets and 1.0 bar
(km)
glaciers have found the basal shear stress calcu- 2
lated with Eq. (6.4) to be typically between 0.5
and 1.5 bars (Nye, 1952), a surprisingly narrow 1
0.5 bar
range considering the spatial variability of ob-
served basal conditions, including gradients in 0
temperature, meltwater content, basal debris, till
rheology, and other variables.
200 300
100 400 500
0
6.2.1 Analytic solutions x (km)
Equation (6.4) is a general result for any ˜˜slowly™™
moving gravity ¬‚ow (i.e. slow enough that acceler- Fig 6.2 Plots of ice sheet pro¬les predicted by the perfectly
plastic model with (b) and without (a) isostatic adjustment
ations can be neglected). In cases where the ¬‚uid
included, for a range of threshold basal shear stresses.
moves rapidly when a threshold shear stress is
exceeded, Eq. (6.4) can be used to model the ¬‚ow
geometry by treating „0 as a constant. Taking the the right side of Eq. (6.7) has units of length
limit as x goes to zero and expressing Eq. (6.4) and is often written as h 0 (not to be confused
as a differential equation gives with an initial, or maximum thickness). For an
„0 ice sheet with „0 = 1 bar (105 Pa), ρ = 920 kg/m3 ,
dh
=
h (6.5)
and g = 9.8 m2 /s, h 0 ≈ 12 m. It should be noted
ρg
dx
that this approach breaks down near the divide,
Multiplying both sides of Eq. (6.5) by dx and in-
where additional stress components become im-
tegrating yields
portant.
„0
h2 = x+C (6.6) The parabolic solution of Eq. (6.7) is usually
ρg
associated with perfectly-plastic deformation of
The integration constant C is equal to the termi- near-basal ice. Perfectly plastic deformation, how-
nus position. If the ¬‚ow has a half-width of L , for ever, is an idealization of just one of the com-
example, then C = L , and Eq. (6.6) becomes plex basal ¬‚ow mechanisms we recognize today,
including debris-controlled frictional sliding and
2„0
h= (L ’ x) (6.7) till deformation. So, how can Eq. (6.7) provide
ρg
even an approximate model for real ice sheets
Figure 6.2a illustrates several examples of this and glaciers? The reason is that many deforma-
simple parabolic ice-sheet model solution. Equa- tion mechanisms mimic the threshold-controlled
tion (6.7) is usually referred to as a ˜˜perfectly ¬‚ow of the perfectly plastic model. Many ¬‚ow
plastic™™ model for ice sheets. The ratio „0 /ρg on mechanisms are threshold controlled and result
6.2 MODELING NON-NEWTONIAN AND PERFECTLY PLASTIC FLOWS 127


in low effective friction when a threshold shear Plots of Eqs. (6.12) and (6.13) are shown in Figure
6.2b for the same cases as Figure 6.2a with δ =
stress is exceeded. Threshold behavior, for exam-
ple, characterizes the plastic deformation of sub- 0.4. Equation (6.13) indicates that the volume of
glacial till (e.g. Tulaczyk et al., 2000) and the shear an ice sheet increases when isostatic effects are
stress necessary to overcome kinetic friction dur- taken into account. In essence, the ˜˜hole™™ created
ing basal sliding (e.g. Lliboutry, 1979). The per- by the weight of the ice provides an additional
fectly plastic model, therefore, can be generally con¬ning force that enables the ice to be thicker
applied to model the geometry of ice sheets and than it would be on a ¬‚at bed.
glaciers that move at low effective friction when Equation (6.4) can be generalized to cases with
a threshold shear stress is exceeded, whether rel- variable bed topography. In such cases, the form
ative motion occurs internally or by sliding over of Eq. (6.4) does not change, but the surface slope
deforming till or directly over the bed. To em- S is no longer simply equal to dh/dx. Instead, S
is equal to dz/dx = dh/dx + db/dx where b(x) is
phasize this broader applicability, we will also
use the term ˜˜threshold-sliding™™ to refer to the the bed elevation. Equation (6.4) then becomes
perfectly plastic model.
„0
dh db
+ =
The results of Figure 6.2a are unrealistic in h (6.14)
ρg
dx dx
at least one major respect: the subsidence of the
crust beneath the ice sheet was neglected. For
Equation (6.14) can be used to derive an ana-
ice sheets larger than a few hundred kilometers,
lytical expression for the approximate shapes of
the weight of the overlying ice will cause signif-
alpine glaciers with a uniform bed slope. In such
icant bedrock subsidence beneath the ice sheet
cases, db/dx = Sb and Eq. (6.14) becomes
in order to achieve isostatic equilibrium. This ef-
fect can be quanti¬ed with a simple adjustment dh
+ Sb = h0
h (6.15)
to the threshold-sliding model. If the bed eleva- dx
tion is given by b(x), then isostatic equilibrium
where h 0 = „0 /ρg. Dividing both sides of Eq. (6.15)
requires that
by h gives
ρi z = ’ρm b (6.8)
dh h0
= Sb ’
where ρi is the density of ice, z = b + h is the (6.16)
dx h
ice-surface elevation, and ρm is the density of the

<< . .

. 25
( : 51)



. . >>