lithosphere to loading and unloading can gener- ify the boundary conditions. The value of w must

ally be neglected over time scales of interest in go to zero as x goes to ∞ and dw/dx must go to

mountain building. zero at x = 0. This implies that c 1 and c 2 are zero

The equation for the 1D ¬‚exural-isostatic dis- and c 3 = c 4 to give

placement of a uniformly rigid lithosphere to ver- x x

w(x) = c 3 e’x/± cos + sin (5.10)

tical loading and unloading is: ± ±

d4w The value of c 3 is proportional to the applied load

+ (ρm ’ ρc )gw = q(x)

D (5.6)

dx 4 V . From Eq. (5.3) and Eq. (5.10) we have

This equation combines ¬‚exural bending stresses d3 w

1 4D c 3

V=D = (5.11)

(¬rst term on left side), the upward buoyancy ±3

dx 3

2 x=0

force exerted on the bottom of the lithosphere

We used 1 V for the load in Eq. (5.11) because we

(second term on left side), and the applied load 2

are solving for only half of the pro¬le (i.e. the

q(x) (right side).

other half of the pro¬le supports the other half

of the load). Substituting Eq. (5.11) into Eq. (5.10)

5.2 Methods for 1D problems gives

V ± 3 ’x/± x x

w(x) = cos + sin

e (5.12)

± ±

5.2.1 Displacement under line loading 8D

for x ≥ 0 (Figure 5.2b).

In some cases the topographic load is nar-

row enough compared to the ¬‚exural wave- In the case of de¬‚ection of a diving board,

length of the lithosphere that the load can be we found that the solution is characterized by

112 FLEXURAL ISOSTASY

thickness of 40 km, for example, the ¬‚exural pa-

200

rameter under crustal loading is ≈ 130 km, while

for ice loading it is approximately 90 km.

crust

150

(km) 5.2.2 Series solutions and Fourier ¬ltering

The ¬‚exure equation is a linear equation. As such,

100

ice we can use Fourier series and superposition to de-

velop general solutions for the de¬‚ection of the

50 lithosphere under any topographic load, given

the solution for de¬‚ection under a periodic load.

Let™s assume a periodic topography with ampli-

0

60

20 40 80 tude h 0 and wavelength »:

0

T e (km)

2π x

h(x) = h 0 sin (5.16)

Fig 5.3 Plot of ¬‚exural parameter ± as a function of elastic »

thickness Te for crustal and ice loads, assuming E = 70 GPa

The vertical load corresponding to this topogra-

and ν = 0.25.

phy is given by

a cubic polynomial with no characteristic length 2π x

qa (x) = ρc gh 0 sin (5.17)

»

scale. In the case of ¬‚exural-isostasy, however, the

de¬‚ection is characterized by the length scale ±.

and hence the ¬‚exural-isostatic equation we must

Also, the zone of depression is accompanied by

solve is

a smaller zone of uplift ¬‚anking the margin of

d4 w

the load. This uplift zone is the forebulge. The + (ρm ’ ρc )gw = qa (x)

D (5.18)

dx 4

location of the forebulge is determined by dif-

ferentiating w with respect to x and setting the The displacement of the lithosphere will follow

result equal to zero: the same periodic form as the loading -- the mag-

nitude of the de¬‚ection will depend sensitively on

dw 2w 0 ’x/± x

=’ sin = 0

e (5.13) the wavelength of the load and ¬‚exural param-

± ±

dx

eter of the lithosphere, but there will be no off-

Therefore,

set or shift between the load and the response.

xb = ± sin’1 (0) = π± (5.14) Therefore, we look for solutions of the form

The amplitude of the forebulge can be found by 2π x

w = w 0 sin (5.19)

»

substituting xb into Eq. (5.12) to obtain

w b = ’w 0 e’π = ’0.0432w 0 (5.15) where w 0 is yet to be determined. Substituting

Eq. (5.19) into Eq. (5.18) gives

Therefore, the height of the forebulge is equal

h0

to about 4% of the maximum de¬‚ection of the

w0 = (5.20)

ρm 2π 4

’1+

lithosphere beneath the load. This calculation as- D

ρc ρc g »

sumes that the lithosphere is entirely rigid. In

For small wavelengths, i.e. where

some cases, the loading causes the lithosphere

to break under the stresses of the load. In these 1

D 4

»

cases, the sin term in Eq. (5.12) becomes zero and 2π (5.21)

ρc g

the width of the zone of subsidence becomes nar-

rower (Turcotte and Schubert, 2002). the denominator in Eq. (5.20) is very large

The relationship between ± and T e is plotted and hence w 0 h 0 . This means that small-

in Figure 5.3 for loading by crustal rocks (i.e. wavelength topographic variations cause very lit-

ρ ≈ 600 kg/m3 ) and ice (i.e. ρ ≈ 2300 kg/m3 ), as- tle displacement of the crust and we say that the

suming E = 70 GPa and ν = 0.25. For an elastic crust is rigid for loads of this scale. At large scales,

5.3 METHODS FOR 2D PROBLEMS 113

where (so called because it ¬lters out low wavenumbers)

before the inverse Fourier transform is applied to

1

D 4

bring the function back to real space.

» 2π (5.22)

ρc g Appendix 4 provides code that implements

the Fourier ¬ltering technique for any input to-

Eq. (5.20) gives

pographic load. Figure 5.4 illustrates example

ρc

w0 = h0 output of that program for an east--west pro¬le

(5.23)

ρm ’ ρ c

through the central Andes. The east--west topo-

graphic pro¬le of the central Andes is character-

which is equivalent to the case of pure isostatic

ized by a high-elevation plateau (the Altiplano-

balance given by Eq. (5.2). In the large-scale limit,

Puna basin), high peaks on the margins of the

therefore, ¬‚exure plays no role and the topo-

Altiplano-Puna, and an active fold thrust belt to

graphic load is said to be fully compensated. At

the east of the Altiplano-Puna. In this example,

smaller scales, where the de¬‚ection is smaller

the topographic load of the Andes is considered

than the Airy isostatic limit, the topography is

to be all topography above 1 km in elevation. Fig-

said to be only partially compensated.

ure 5.4b plots the modeled de¬‚ection of the crust

The solution for periodic topography is use-

for ρc = 2.7 g/cm3 , ρm = 3.3 g/cm3 , and three hy-

ful for determining the degree of compensation

pothetical values of ±. The depth of total de¬‚ec-

of a topographic load. The real power of the pe-

tion beneath the Andes (≈ 15 km) depends lin-

riodic solution, however, is that it can be used

early on the relative density contrast between the

in conjunction with the Fourier series approach

crust and mantle, but only weakly on the ¬‚ex-

to determine the response of the lithosphere to

ural parameter for wide topographic loads. Low

any load using superposition. If the topographic

values of the ¬‚exural parameter predict narrow

load within a domain of length L is expressed as

foreland sedimentary basins, while higher val-

a Fourier series:

ues predict broad basins (Figure 5.4b--5.4c). The

∞

πx

h(x) = height of the ¬‚exural forebulge (Figure 5.4c) also

an sin (5.24)

L increases with higher values of the ¬‚exural pa-

n=1

rameter. Gravity modeling has been used to con-

with coef¬cients

strain the width of the Andean foreland basin by

L

2 nπ x taking advantage of the density contrast between

an = h(x ) sin dx (5.25)

L L the basin in¬ll and the crust below. Those con-

0

straints, coupled with forward modeling of ¬‚ex-

then, according to Eq. (5.20), the de¬‚ection

ural pro¬les, suggest that ± = 100 km is a good

caused by the load is given by

estimate for the ¬‚exural parameter of the central

∞ nπ x

Andes.

sin

w(x) = L

an (5.26)

ρm nπ 4

’1+ D

n=1 ρc ρc g L

Equation (5.26) can be used to obtain analytic ex-

5.3 Methods for 2D problems

pressions for the ¬‚exural-isostatic displacement

underneath loads with simple (rectangular, tri-

5.3.1 Integral method

angular, etc.) loads.

Most real topography, however, does not lend Thus far we have considered solutions in 1D only.

itself to simple analytic forms. Equation (5.26) These solutions are very useful for long moun-

can be implemented numerically for any topo- tain chains or other landforms where the loading

graphic load, however, using the Fourier ¬lter- is relatively constant in one direction. In many

ing technique. In this technique, the Fast Fourier cases, however, (e.g. an isolated seamount or curv-

Transform (FFT) is used to quickly calculate the ing mountain belt) a full 2D solution is needed.

Fourier coef¬cients an for an input load function. In such cases, 2D methods analogous to those we

Then, each coef¬cient is multiplied by the ¬lter developed for 1D are available.

114 FLEXURAL ISOSTASY

(a) 5

4

h3

(km)

(c) 0.4

2

w

0.2

1 (km)

0

0.0

(b) 0

± = 50 km

75 km

0.2

100 km

w5

0.4

(km) ± = 50 km

1200

800 1000 1400

75 km

100 km x (km)

10

15

600 800 1000 1200

400

200

0

x (km)

Solutions corresponding to a spatially dis-

Fig 5.4 Flexural model of the response of the lithosphere

to the distributed load of the central Andes. (a) Topographic tributed load can be obtained by representing

pro¬le of the central Andes. (b) De¬‚ection of the lithosphere the load as a collection of many point sources

beneath the central Andes for three hypothetical values of

and integrating the response to each point. For

± = 50, 75, 100 km. (c) Close-up plot of the ¬‚exural pro¬le

a spatially distributed topographic load given by

of the foreland basin and forebulge, illustrating the

h(x, y), the de¬‚ection is given by:

dependence of basin width and forebulge amplitude on ±.

∞ ∞

ρc

First we consider the solution to the 2D w(r ) =

ρm ’ ρc

¬‚exural-isostatic equation: 0 0

√

2

— kei (x ’ x )2 + (y ’ y )2

D ∇ 4 w + ρgw = q ±

(5.27)

— h(x , y )dx dy (5.29)

corresponding to a point load V at r = 0. The

solution is the Kelvin function kei:

√

V ±2 This approach is broadly analogous to the series

2r

w(r ) = kei (5.28)

πD ± solution method of Eq. (5.26).

Few analytic solutions can be obtained using

The shape of the kei function is broadly similar Eq. (5.29), both because topographic loads rarely

to Eq. (5.12), i.e. it has a deep zone of subsidence conform to analytic expressions and because the

beneath the load and a smaller forebulge ¬‚ank- kei function is dif¬cult to integrate. Appendix 4

ing the load. provides code for series approximations to the kei

5.3 METHODS FOR 2D PROBLEMS 115

ral parameter, the de¬‚ection follows the disk-like

shape of the load except at the edges where ¬‚ex-

0 0.2

ural bending stresses limit the curvature of the

w/h de¬‚ection pro¬le.

1 R/a =

1.0

4.0

5.3.2 Fourier ¬ltering

2 In the previous section we showed how the de¬‚ec-

tion beneath spatially distributed loads could be

calculated by dividing the load into many small

3

point sources and summing the contributions of

each point. This is a powerful approach, but it

4

can be computationally challenging when calcu-

’1 0 1

lating the response at high spatial resolution over

x/a

a large region. Such calculations can require eval-

uating the kei function 104 --106 times for each of

Fig 5.5 De¬‚ection of disk-shaped loads for various values

104 --106 points. The Fourier-¬ltering method pro-

of R /±.

vides a faster way by taking advantage of the

computational speed of the Fast Fourier Trans-

function that can be used to integrate Eq. (5.29)

form algorithm.

numerically for any topographic load.

In order to implement the Fourier ¬ltering al-

One analytic expression is available, how-

gorithm in 2D, the FFT of the 2D load h(x, y) over

ever, that nicely illustrates the response of the

a domain of width L x and length L y is computed

lithosphere to loads of different size. The re-

to obtain a set of Fourier coef¬cients an,m . Each

sponse of the lithosphere to a disk-shaped load

Fourier coef¬cient is then multiplied by the fac-

depends on the ratio of the disk radius R to the

tor

¬‚exural parameter ±. The solution is given by

⎛ ⎛ ⎞4 ⎞’1

(Hetenyi, 1979): 2 2

⎝ ρm ’ 1 + D ⎝π n m ⎠⎠

+

hρc R R r

ρc ρc g

w(r ) = Lx Ly

ker ber

ρm ’ ρ c ± ± ±

(5.32)

R R r

’ kei +1

bei (5.30)

± ± ± Then, the inverse Fourier transform is applied to

the ¬ltered data set.

for r < R and

Appendix 4 provides a code for implementing

hρc R R r

w(r ) = the Fourier ¬ltering method in 2D. Figure 5.6 il-

ber ker

ρm ’ ρ c ± ± ±

lustrates example output from this program for

R R r

the central Andes for ± = 50, 75, and 100 km. Fig-

’ bei kei (5.31)

± ± ±

ures 5.6b--5.6d are grayscale plots of the absolute

for r > R . value of the de¬‚ection on a logarithmic scale. As

the value of ± increases, the pattern of de¬‚ection

Figure 5.5 plots Eq. (5.31) for loads of differ-

ent radii, expressed as a ratio of the ¬‚exural pa- becomes smoother. In each grayscale image the

rameter ±. For cases in which the load radius forebulge is represented by the ¬rst ring around

is much less than the ¬‚exural parameter, the the central zone of subsidence. Subsequent rings

load is only partially compensated. For cases in are the backbulge and other smaller zones of al-

which the radius is equal to the ¬‚exural parame- ternating uplift and subsidence with increasing

ter, the maximum de¬‚ection is equal to the Airy distance from the load. Figure 5.6e illustrates the

isostatic limit (i.e. Eq. (5.2)), but the maximum de- variations in forebulge amplitude in the vicinity

¬‚ection is attained only beneath the center of the of the Bolivian Orocline using a contour map. The

load. For load radii much larger than the ¬‚exu- forebulge amplitude varies from 160 m to 400 m

116 FLEXURAL ISOSTASY

basin environments are also attractive environ-

Fig 5.6 2D ¬‚exural-isostatic response of the lithosphere to

ments for studying the linkages of climate, tec-

the topographic load of the central Andes (a) for ± = 50 (b),

tonics, and erosion, because the tectonic param-

75 (c), and 100 km (d). In (b)“(d), the grayscale images

eters controlling foreland basin geometry can be

represent the absolute value of the de¬‚ection on a

unusually well constrained. Speci¬cally, thrust-

logarithmic scale. This color scale clearly shows that the

width of the forebulge increases with increasing ±. (e) A belt migration rates can be determined using

contour map of the forebulge in the vicinity of the Bolivian a mass-balance framework that relates crustal

Orocline. This example shows that the height of the shortening to thrust-belt migration (DeCelles and

forebulge changes dramatically (from a minimum of 160 m to

DeCelles, 2001) and basin subsidence is often

a maximum of 400 m) in the vicinity of a curved load.

well characterized using ¬‚exural modeling, grav-

ity measurements, and seismic stratigraphy (e.g.

within a distance of ≈ 100 km, illustrating the

Watts et al., 1995). Sediment transport in ¬‚uvial

importance of 2D effects in the vicinity of a bend

channels has a linear relationship with channel

in the topographic load.

slope, resulting in a diffusion model for the basin

topographic pro¬le (Begin et al., 1981; Paola et al.,

1992).

5.4 Modeling of foreland basin

Foreland basin evolution can be modeled as

geometry a combination of thrust-belt migration, ¬‚exural

subsidence, and diffusive ¬‚uvial sediment trans-

port (Flemings and Jordan, 1989). Modeling this

Flexure plays an important role in large-scale

system requires solving the diffusion equation

landform evolution because, as erosion and de-

with a moving boundary (the thrust front, where

position change the topographic load, ¬‚exural

sediment enters the basin) and moving sink

isostasy causes uplift and subsidence responses

term (the ¬‚exural depression). Typically, mov-

that vary as a function of spatial scale and litho-

ing boundaries and sinks in diffusion problems

spheric rigidity. In this section we explore the

require numerical approaches. However, in this

interaction between ¬‚exural isostasy and sedi-

case the problem can be solved analytically by

mentation in a foreland basin adjacent to a mi-

working in the frame of reference of the migrat-

grating thrust belt. Foreland basins provide a po-

ing thrust front and assuming steady-state con-

tentially useful record of tectonics and erosion

ditions. Conceptually, basins are in steady state

of the adjacent range, but our ability to inter-

when sediment progradation occurs at the same

pret that record uniquely is limited. Foreland

5.4 MODELING OF FORELAND BASIN GEOMETRY 117

rate as thrust-belt migration. Basins naturally where w(x) is the ¬‚exural pro¬le, given by Tur-

evolve toward this condition because when sed- cotte and Schubert (2002) as

iment progradation lags behind thrust-belt mi- x

w(x) = ’w 0 e’x/± cos (5.36)

±

gration at any point along its pro¬le, the basin

steepens locally until the two are in balance. where w 0 is the basin depth beneath the thrust

Conversely, if sediment progradation outstrips front and ± is the ¬‚exural parameter. Equa-

thrust-belt migration, slope decreases and progra- tion (5.36) corresponds to the ¬‚exural pro¬le cor-

dation slows. This approach is analogous to that responding to a line load. This is only an approx-

of the Kenyon and Turcotte model of deltaic sed- imation to the Andes, where loading is spatially

imentation considered in Chapter 2. In this case, distributed.

however, we will also consider the pattern of Combining Eqs. (5.33) and (5.35) and as-

¬‚exural subsidence that controls accommodation suming steady-state conditions gives the time-

space in the basin. independent equation:

Basin evolution in the model is governed by

d2 h dh dw

κ +F =F

the diffusion equation: (5.37)

2

dx dx dx

‚h ‚ 2h which can be written as a ¬rst-order differential

=κ 2 (5.33)

‚t ‚x equation for the slope S = dh/dx with solution

where h is the elevation, t is time, κ is the diffu-

’e’F x/κ x

dw

e’F x /κ

S(x) = Qs ’ F dx

sivity, and x is the distance from the thrust front.

κ dx

0

Boundary conditions are needed at both ends of

(5.38)

the basin to solve Eq. (5.33). At the thrust front,

Substituting Eq. (5.36) into Eq. (5.38) gives

a constant sediment ¬‚ux Q s (in m2 /yr) enters the

’e’F x/κ ±F ±F

F x x x

e’(F /κ+1/±)x

S(x) = Qs ’ ’ sin ’ 2sin ’

cos (5.39)

κ (± F /κ)2 + 2± F /κ + 2 κ ± ± ± κ

The basin topographic pro¬le is obtained by inte-

basin, providing a ¬‚ux boundary condition

grating Eq. (5.39) analytically or numerically and

‚h Qs enforcing the basin-outlet boundary condition to

= (5.34)

‚x κ constrain the integration constant.

x=0

Model solutions are plotted in Figure 5.7, rep-

In this 2D framework, the sediment ¬‚ux is given

resenting the effects of varying Q s on basin pro-

by ELd , where E is the basin-averaged erosion rate

¬les around a reference case with parameter val-

and L d is the upstream drainage basin length.

ues F = 10 mm/yr, ± = 150 km, w 0 = 5 km, κ =

At the basin outlet, a constant base-level eleva-

5000 m2 /yr, Q s = 10 m2 /yr, and assuming an in-

tion boundary condition is prescribed. Two cases

¬nite basin. Figure 5.7a illustrates three impor-

can be considered. For an in¬nite basin, h = 0 as

tant points. First, the basin pro¬le is highly sen-

x = ∞. For a basin of length L b and base-level ele-

sitive to incoming sediment supply: the topo-

vation h b (determined by sea level or a valley-¬‚oor

graphic relief of the basin doubles when Q s is

channel), the boundary condition is h(L b ) = h b . In

increased by 50% from 8 to 12 m2 /yr. This suggests

the moving reference frame, the foreland basin

that basin pro¬les may provide useful constraints

moves towards the thrust front and enters the

on upstream sediment supply through forward

¬‚exural depression with velocity F , the thrust-

modeling and comparison with observed pro¬les.

migration velocity. This motion is represented by

Second, as sediment supply decreases, the basin

an advection equation:

undergoes a transition from over¬lled to under-

‚h ‚w ‚h ¬lled (closed) basins. As such, the model provides

=F ’ (5.35)

quantitative criteria for topographic closure in

‚t ‚x ‚x

118 FLEXURAL ISOSTASY

0.8

(b)

reference case:

overfilled

e’Fx/k F = 10 mm/yr

k=

a = 150 km

Qs = 12 m2/yr

0.6 w0 = 5 km 10 000 m2/yr

k=

Qs = 10 m2/yr k = 5000 m2/yr

5000 m2/yr

Qs = 8 m2/yr Qs = 10 m2/yr

h (km)

k=

1500 —

0.4 2500 m2/yr

underfilled

vert. exag.

(closed basin)

0.2 forebulge k=

1000 m2/yr

sensitivity to Qs sensitivity to k

0.0 infinite basin

infinite basin

(a) no forebulge erosion

no forebulge erosion

0.8 Ld

sensitivity to Qs

h

finite basin

(c) (d)

no forebulge erosion

0.6 Qs = 12 m2/yr

Qs = 10 m2/yr flux boundary

Q

h (km)

Qs = 8 m2/yr

dh/dx = Qs /k sediment flux s

F k

0.4 diffusivity

migration

basin-fill profile

rate

(diffusion)

decollement

x=0

0.2 slope

x

flexural profile

w0

0.0

w

800 1000

400 600

200

0

x (km)

cause sediment back¬lling at the thrust front to

Fig 5.7 (a) and (b): Model solutions for the topographic

produce steep, short basins. Surprisingly, higher

pro¬le of an in¬nite basin assuming no forebulge erosion,

diffusivities (e.g. more humid conditions) pro-

showing sensitivity to variations in Q s and κ, respectively,

around a reference solution with F = 10 mm/yr, ± = 150 km, mote basin closure if sediment supply is held

w 0 = 5 km, κ = 5000 m2 /yr, and Q s = 10 m2 /yr. (c) Solutions constant. Naively, we might expect more humid

for a ¬nite basin of length L b = 1000 km for the same model conditions to result in steeper, more over¬lled

parameters as in (a). (d) Schematic diagram of the model.

basins. However, to understand the in¬‚uence of

Modi¬ed from Pelletier (2007b).

climate on basin geometry we must consider

the effects of climate on sediment supply (i.e.

foreland basins. Third, the distal basin pro¬le is weathering) and basin transport rates indepen-

approximated by a simple exponential function dently. Greater sediment supply puts more total

with length scale κ/F . This result suggests that sediment into the basin and therefore promotes

distal basin pro¬les can be used to uniquely and steeper and more over¬lled basins. Diffusivity val-

easily constrain κ values if thrust-belt migration ues, on the other hand, control how uniformly

rates are known. that sediment is spread across the basin. Higher

Figure 5.7b illustrates the sensitivity of model diffusivity values promote closed basins because

solutions to variations in κ around the reference a larger fraction of the total sediment is trans-

value of κ = 5000 m2 /yr. Diffusivity values are ported past the foredeep, leaving behind a depres-

sion if the ratio Q s /κ is suf¬ciently small. These

controlled by upstream basin length, precipita-

tion, and sediment texture (with longer basins, results emphasize that the key to understanding

higher precipitation rates, and ¬ner textures pro- basin geometry is not simply the sediment sup-

moting greater diffusivity values) (Paola et al., ply to the basin, but the ratio of sediment supply

1992). Low diffusivity values (e.g. κ = 1000 m2 /yr) to transport capacity (e.g. as emphasized in the

5.4 MODELING OF FORELAND BASIN GEOMETRY 119

geomorphic context by Bull (1979)). Figure 5.7c on crustal shortening rates. Flexural modeling

illustrates model solutions for a ¬nite basin us- coupled with geophysical and geomorphic obser-

ing the same model parameters as in Figure 5.7a. vations constrain the ¬‚exural parameter to be

Shorter basins have lower basin relief, all else be- approximately 150 km (resulting in a foredeep ap-

ing equal. proximately 250 km wide and a forebulge at a

Equation (5.40) assumes no forebulge erosion distance 400 km from the thrust front) (Coud-

(i.e. that the forebulge is buried or, if exposed, ert et al., 1995; Watts et al., 1995; Horton and

that the erosion rate is low). In some cases, it may DeCelles, 1997; Ussami et al., 1999). These data

be more appropriate to assume that the forebulge provide relatively ¬rm constraints on three of

the ¬ve model parameters: F = 10 mm/yr, ± =

erodes completely to sea level. In this alternative

150 km, w 0 = 4 km.

end-member scenario, Eq. (5.36) is amended to

w(x) = 0 for x > 3π±/4 (i.e. the point at which Figure 5.9a illustrates a map of the topog-

the forebulge rises above sea level). In this case, raphy of the central Andes and the adjacent

foreland basin from 17—¦ S to 31—¦ S. This region

the slope is given by Eq. (5.36) for x < 3π±/4,

and by

’e’F x/κ ±F ±F

F 1 3π ±

F

1.414e’( κ + ± ) 4

S(x) = Qs ’ 1+ + (5.40)

κ /κ)2 + 2± F /κ + 2 κ κ

(± F

for x > 3π ±/4. Model solutions assuming fore- corresponds to a humid-to-arid transition in the

Andes based upon mean winter precipitation val-

bulge erosion are given in Figure 5.8.

ues above 1500 m elevation. Associated with this

This model framework also provides a simple

climatic transition is a geomorphic transition

means to estimate the basin sediment delivery

from strongly over¬lled (700 m in relief) to weakly

ratio (the ratio of the sediment yield to erosion

over¬lled (100 m in relief) basins. Closed basins

rate). The sediment ¬‚ux at the entrance to the

(e.g. L. Ambargasta and L. Mar Chiquita) also oc-

foreland basin is Q s . For a basin of length L b ,

cur near the southernmost transect, but they are

the ¬‚ux leaving the basin is computed by evaluat-

ing Eq. (5.40) at x = L b and multiplying by κ. For located in wedgetop depozones of a broken fore-

land and hence are not well described by the

basins that have a base-level control located be-

model.

yond the forebulge distance, Eq. (5.40) reduces to

Three topographic pro¬les were extracted

a simple exponential function with length scale

κ/F . In the reference case, for example, a 500-km from the Shuttle Radar Topographic Mission

(SRTM) data and are plotted in Figure 5.9b. To

wide basin has a sediment delivery ratio of 1/e

= 0.37, and a 1000-km wide basin has a value model these pro¬les, it is most accurate to use

of 1/e2 = 0.14. This latter value is broadly consis- ¬nite-basin solutions constrained by values of h b

and L b for each pro¬le. The base level eleva-

tent with measured sediment delivery ratios for

tion h b of the Andean foreland is controlled by

large drainage basins in areas of active tectonism

Rio Paraguay, which ranges in elevation from 80

(Schumm, 1977).

to 10 m (decreasing from north to south) at dis-

As a test of the model, we consider the mod-

tances L b varying from 350 to 720 km from the

ern foreland basin of the central Andes. The

thrust front. Model results in Figure 5.7a sug-

latest phase of Andean deformation responsible

gest that the distal basin pro¬le is an exponential

for Eastern Cordilleran and Subandean uplift is

function with length scale κ/F . The inset plot in

generally considered to be late Miocene in age

Figure 5.9b showing hh b (note vertical logarith-

(Gubbels et al., 1993). The foreland basin is ap-

mic scale) as a function of x con¬rms this pre-

proximately 3.5--4.5 km thick beneath the thrust

diction for the central Andean foreland. Best-¬t

front based on seismic re¬‚ection data (Horton

exponential pro¬les (straight lines on this loga-

and DeCelles, 1997). DeCelles and DeCelles (2001)

rithmic scale) provide estimates of κ = 1500 m2 /yr

estimated the thrust-migration rate to be approx-

for the northern pro¬le and κ = 4000 m2 /yr for

imately 10 mm/yr for the central Andes based

120 FLEXURAL ISOSTASY

result of assuming uniform κ values. Down-

(a) 0.8 reference case:

stream ¬ning occurs in all basins, and therefore

F = 10 mm/yr

overfilled

± = 150 km

κ values are likely to increase somewhat down-

0.6 w0 = 5 km

h stream, resulting in gentler distal-basin slopes

Qs = 12 m2/yr κ = 5000 m2/yr

Qs = 10 m2/yr

Qs = 10 m2/yr

(km) compared to the model solutions of Figure 5.9b.

Qs = 8 m 2

/yr 1500x

0.4 vert. exag.

underfilled

(closed basin)

5.5 Flexural-isostatic response to

0.2 forebulge

glacial erosion in the western

sensitivity to Qs

0.0 US

infinite basin

forebulge erosion

(b) 0.8

The pace of alpine glacial erosion increased sig-

ni¬cantly starting 2--4 Myr ago when global cool-

0.6

h ing initiated the Plio-Quaternary era of large con-

κ=

10,000 m2/yr

(km) tinental ice sheet advances and retreats (Zhang

0.4 et al., 2001). This erosion must have triggered sig-

κ=

5000 m2/yr

ni¬cant ¬‚exural-isostatic unloading of glacially

κ=

2500 m2/yr carved mountain belts. Since glacial erosion is

0.2 κ=

strongly focused in the landscape, scouring out

1000 m2/yr

sensitivity to κ

valley ¬‚oors where ice ¬‚ow is focused but erod-

0.0 infinite basin

forebulge erosion

ing the highest points in the landscapes rela-

800 1000

400 600

200

0 tively slowly, it is likely that glaciated mountain

x (km) belts have increased in relief during the Plio-

Quaternary era due to isostatic rock uplift, even

Fig 5.8 (a) and (b) Model solutions for the topographic

in areas that are tectonically inactive. Determin-

pro¬le of an in¬nite basin assuming complete forebulge

ing precisely how much glacial erosion has taken

erosion, showing sensitivity to variations in Q s and κ,

place is dif¬cult to do in most cases. Nevertheless,

respectively, for the same model parameters as Figures 5.7a

it is still useful to ask how much rock uplift oc-

and 5.7b. These results show that model results are broadly

similar to those obtained by assuming no forebulge erosion; curred for a given amount of glacial erosion. By

the only signi¬cant difference is a minor reduction in basin using the geographic extent of glacial cover as a

relief when forebulge erosion is assumed.

proxy for the spatial distribution of erosion, we

can use the 2D ¬‚exural-isostatic model to com-

pute the ratio of late Cenozoic erosionally driven

the two southern pro¬les using F = 10 mm/yr. rock uplift to glacial erosion. In this section we

Given observed or inferred values for h b , L b , and compute that ratio for the western United States,

κ, as well as estimated values for F , ±, and w 0 following the approach of Pelletier (2004a). Ar-

from the published literature, a family of solu- eas where this ratio is large have likely under-

tions corresponding to a range of Q s values can gone the greatest relief production and should

be generated for each pro¬le and compared to be the focus of future efforts to identify signa-

the observed data. Figure 5.9c illustrates the ob- tures of erosionally driven rock uplift and re-

served pro¬les plotted with their corresponding lief production. In this application we will use

best-¬t pro¬les (κ values for the two southern- the Alternating-Direction-Implicit technique, tak-

most pro¬les were also varied to ¬nd the optimal ing into account spatial variations in ¬‚exural

¬t for both κ and Q s ). Model solutions match rigidity.

the observed data very closely for the strongly The ratio of rock uplift to erosion depends on

over¬lled pro¬le, and they match the ¬rst-order the area and shape of the eroded region and on

trends of the weakly over¬lled pro¬le. Discrep- the local ¬‚exural wavelength of the lithosphere.

ancies in the southernmost pro¬les may be the If the ¬‚exural wavelength is large compared to

5.5 FLEXURAL-ISOSTATIC RESPONSE TO GLACIAL EROSION IN THE WESTERNUS 121

68° W 57° W 0.8 1

(b)

k = 1500 m2/yr

17° S h ’ hb

(km) k = 4000 m2/yr

0.6

0.1

h

sediment-rich (km)

0.4

0.01

400

0 200 600

x (km)

0.2

0.0

600

0 200 400

x (km)

Rio Paraguay

(c)

k = 1500 m2/yr

Qs = 9 m2/yr

0.6

h ’ hb

(km)

sediment-starved 0.4 k = 2000 m2/yr

Qs = 6.5 m2/yr

k = 2000 m2/yr

Qs = 5 m2/yr

0.2

100 km (a)

31° S L. Mar Chiquita

L. Ambargasta

0.0

elevation

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

which integrates to

mantle. Substituting z = b + h, Eq. (6.8) becomes

ρi h h0

b=’ z = ’δz x’C = + 2 ln(h 0 ’ Sb h)

(6.9) (6.17)

ρm ’ ρi Sb Sb

where where C is the integration constant. Figure 6.3a

ρi illustrates solutions to Eq. (6.17) for C = 0, Sb =

δ= (6.10)

ρm ’ ρi 0.05 and 0.1, and h 0 = 12 m (i.e. „0 = 1 bar). In or-

der to better visualize the results, it is helpful to

The ice thickness can be written as

plot z = ’Sb x + h instead of simply h. Figure 6.3b

h = (1 + δ)z (6.11)

presents these results. The value of C can be var-

ied to shift the glacier terminus to the right or

Substituting Eq. (6.11) into Eq. (6.7) yields the

left, but otherwise the value of C does not affect

following equations for the ice-sheet elevation

the solution.

and thickness:

Another solution relevant to alpine glaciers

„0

z= (L ’ x)

2 (6.12) involves the variations in thickness caused by mi-

ρg(1 + δ)

nor variations in valley slope. If the bed topog-

raphy is slowly-varying and the region close to

„0 (1 + δ) the glacier terminus is not considered, the dh/dx

h= (L ’ x)

2 (6.13)

ρg term can be neglected relative to the db/dx term

128 NON-NEWTONIAN FLOW EQUATIONS

250 300

3.0

(a)

h

z

(m)

h 200 (km) h

Sb = 0.05

(m) 2.0 200

150

z

100 100

1.0

Sb = 0.1

b

50

0

0.0

20

10 30

0

0

x (km)

800

(b)

Fig 6.4 Plots of ice thickness variations in an alpine glacier

with slowly varying bed topography, according to Eq. (6.18),

h

z 400 with h 0 = 12m.

(m)

0 scribed, however, the geometry of the ice sheet

or glacier does not explicitly depend on mass bal-

ance. When considering ice bodies obeying Glen™s

Sb = 0.1

Flow Law, however, it is necessary to prescribe a

Sb = 0.05 mass balance rate M . Here we will assume M to

be uniform and constant, but in general it will

10 15

0 5

x (km) vary in both space and time.

The starting point for modeling ice-sheet or

Fig 6.3 Plots of (a) thickness and (b) ice-surface elevation

glacier geometries with a power-law ¬‚uid rheol-

for Sb = 0.05 and 0.1 for h 0 = 12 m (i.e. „0 = 1 bar).

ogy is the relationship between shear stress and

strain rate, given in Chapter 1 for the case of sim-

to give ple shear:

’1

db du

h = h0 = A„ n

(6.18) (6.19)

dx dy

Figure 6.4 plots solutions to Eq. (6.18) for h 0 = where u is the ¬‚uid velocity, y is the depth, „ is

12 m, a mean valley slope of 0.05, and an ideal- the shear stress, and A and n are rheological pa-

ized sinusoidal bed topography. The ice thickness rameters. Substituting Eq. (6.4) into Eq. (6.19) and

is inversely proportional to the bed slope, thick- integrating from 0 to h gives the velocity pro¬le

ening in areas of gentle slopes and thinning in

A

areas of steep slopes. u(y) = ’ (ρgS)n (h ’ y)n+1 ’ h n+1 (6.20)

2(n + 1)

When considering the geometry of ice sheets

and glaciers within the framework of the per-

Vialov (1958) was the ¬rst to obtain an analytic

fectly plastic model, it is not necessary to pre-

solution for the pro¬le of an ice sheet on a ¬‚at

scribe the accumulation or ablation rates because

bed evolving according to a power-law ¬‚uid rhe-

the shape of the ice sheet or glacier is con-

ology. Vialov assumed a radially symmetric ice

trolled entirely by the threshold shear stress „0

sheet with radius R . Equation (6.20) must be in-

(given standard values for ice density and gravity).

tegrated from 0 to h to obtain the total horizontal

Changes in the spatial distribution of accumula-

¬‚ux of ice:

tion may cause the ice terminus to migrate. As

such, the position of the terminus is controlled dh

U = ’C h 4 (6.21)

by mass balance. If the terminus position is pre- dr

6.2 MODELING NON-NEWTONIAN AND PERFECTLY PLASTIC FLOWS 129

where C is a single constant that combines pa-

rameters for viscosity, density, and gravity. Con- perfectly plastic

3.0

servation of mass is given by h

(km)

‚h ‚(r hU ) Vialov

= ’∇ · (H U ) + M = ’ +M (6.22)

‚t r dr 2.0

For an ice sheet in steady state,

d(r hU ) = M r dr (6.23)

1.0

which integrates to

Mr

U= 0.0

(6.24)

2h 1000

400 600 800

200

0

r (km)

Combining Eqs. (6.21) and (6.24) yields

1/n

M Fig 6.5 Comparison of the perfectly plastic model solution

dh = ’

1+2/n

r 1/n dr

h (6.25)

(with „0 = 1 bar) to the Vialov solution with Glen™s Flow Law

2C

with a maximum thickness of 2.5 km.

which can be integrated to obtain

1/n

n n M

(h 2+2/n ’h 2+2/n ) = ’ r 1+1/n (6.26) Figure 6.5 prescribes h c = 2.5 km. The Vialov solu-

c

2n+2 n+1 2C

tion predicts gentler slopes in the ice sheet inte-

where h c is the thickness at the center where

rior and steeper slopes near the margin relative

r = 0. Equation (6.26) can be rearranged to give

to the perfectly plastic model.

1/n

M

’ = ’2

2+2/n

h 2+2/n r 1+1/n

h (6.27)

c

2C 6.2.2 Numerical solutions

Thus far we have considered only ˜˜one-sided™™

The value of h c can be constrained by setting a

ice sheets and glaciers, i.e. those with elevations

boundary condition of vanishing ice thickness at

the margin, i.e. h(R ) = 0, to give that decrease monotonically from a known di-

vide position or increase from a known terminus

1/n

M

=2

h 2+2/n R 1+1/n position. Equation (6.14) is limited to one-sided

(6.28)

c

2C

ice sheets because dh/dx + db/dx must be posi-

Substituting Eq. (6.28) into Eq. (6.27) gives the full tive throughout the domain, increasing from the

margin (at x = 0) to the divide (at x = L ). In the

solution

1/(2n+2) more general case that includes variable bed to-

M

h=2 n/(2n+2)

pography, only the positions of the two margins

2C

are known and the topography on both sides of

n/(2n+2)

R 1+1/n ’ r 1’1/n (6.29)

the ice sheet must be solved for. In such cases, the

For n = 3, Eq. (6.29) becomes divide position cannot be prescribed but must

be solved for along with the geometry. For such

1/8

M 3/8

h=2 R 4/3 ’ r 4/3

3/8

(6.30) cases, Eq. (6.14) must be modi¬ed to include both

2C

positive and negative ice-surface slopes to obtain

It is dif¬cult to compare the Vialov solution „

dh db

+ = (6.31)

precisely to the perfectly plastic solution because ρgh

dx dx

one was derived for a 1D ice-sheet pro¬le and

The boundary conditions for Eq. (6.31) are h = 0

the other for a radially symmetric ice sheet. Nev-

ertheless, Figure 6.5 compares solutions obtained at the two ice margins.

for the two models. The ratio M/C acts as a free Equation (6.31) can be solved with stan-

parameter in the Vialov solution that relates to dard methods for boundary-value problems. The

the maximum height h c . The solution plotted in equation must ¬rst be modi¬ed to remove the

130 NON-NEWTONIAN FLOW EQUATIONS

absolute-value sign by squaring both sides: 2.5

“two-sided”

solution

h 2.0

2 2

„0

dh db

+ = “one-sided”

(6.32)

(km)

ρgh solution

dx dx

1.5

and then by differentiating both sides of Eq.

1.0

(6.32) to obtain:

„2 0.5

d2 h d2 b 1

+ 2 = ’ 20 2 (6.33)

ρ g h3 +

dh db

dx 2 dx dx dx

0.0

200 800 1000

400 600

0

x (km)

Equation (6.33) is nonlinear and singular at

divides (i.e. where dh/dx + db/dx = 0). As such, Fig 6.6 Comparison of analytic and numerical perfectly

care must be taken to solve it. One approach is plastic solutions for an ice sheet over a ¬‚at bed with „0 = 1