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