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

bar and L = 100, 300, and 1000 km.

to use the relaxation method with a Newton it-

eration step (Press et al., 1992). In this approach,

Eq. (6.33) is discretized as

1.2

1.0

h i+1 ’ 2h i + h i’1 + bi+1 ’ 2bi + bi’1

2 „0

2

2x h 0.8

= ’( x) 2 2 3

ρ g h i (h i+1 ’ h i’1 + bi+1 ’ bi’1 ) (km)

0.6

(6.34)

0.4

where x is the resolution cell size of the grid.

This equation is solved iteratively with 0.2

0.0

L i (h iold ) 20 40 100

60

0 80

h inew = h iold ’ (6.35)

‚ L i (h iold )/‚h i x (km)

Fig 6.7 Solutions for ice ¬‚ows on a tilted bed with

where

„0 = 1 bar, L = 100 km, and bed slopes of 0.01, 0.003, and

0.001. As part of this model solution, the divide position is

L (h i ) = h i+1 ’ 2h i + h i’1 + bi+1 ’ 2bi + bi’1 solved for simultaneously with the ice thickness.

„2 2

+ 222 3

ρ g h i (h i+1 ’ h i’1 + bi+1 ’ bi’1 )

6.3 Modeling ¬‚ows with

(6.36)

temperature-dependent

Figures 6.6 and 6.7 illustrate example solutions

viscosity

for 1D ice sheets of different sizes and with

varying bed topography. In the simple case of a

¬‚at bed, the ˜˜two-sided™™ solution predicts gen- Lava is the classic geological example of a non-

tler slopes near the divide compared to the ˜˜one- Newtonian ¬‚uid with a temperature-dependent

sided™™ parabolic solution. In the case of tilted viscosity. As lava cools, it becomes stiffer. The

beds, the divide position is solved for simultane- temperature-dependent viscosity of lava is partly

ously with the ice geometry using the ˜˜two-sided™™ responsible for the rich dynamics of lava ¬‚ows.

model. The divide position is located farther up- Pillow lavas, for example, form underwater

slope on steeper beds. through a periodic surge-type mechanism in

6.3 MODELING FLOWS WITH TEMPERATURE-DEPENDENTVISCOSITY 131

which pressure builds up beneath a cool, stiff parabolic vertical temperature distribution func-

tion that goes to zero at z = 0 and z = H :

¬‚ow margin until the margin fails, causing a

rapid surge of hot, rapidly ¬‚owing lava out in z z

T = 6T av 1’ (6.38)

front of the ¬‚ow. The surging lava then slows H H

rapidly as the lava cools in contact with the H

where T av = (1/H ) 0 T dz is the depth-averaged

water. Slow velocities are again maintained un-

temperature. The ¬‚ow itself is governed by the

til suf¬cient pressure builds up behind the ¬‚ow

radially symmetric Stokes equation in which the

margin to cause another surge.

pressure gradient is balanced by viscous stresses:

Lava ¬‚ows in the geologic record often oc-

‚H ‚ ‚v r

cur as low-relief deposits with a relatively ¬‚at

ρg = · (6.39)

‚r ‚z ‚z

top. When exhumed by uplift and erosion, such

deposits are often more resistant than the sur- where g is gravity, v r is the radial velocity, and

rounding rocks and create a topographic ˜˜inver- ρ is the density contrast between the ¬‚uid and

sion™™ in which the resistant lava ¬‚ow comprises its surrounding medium (e.g. air or water). The

the highest outcrops in the landscape, forming a vertically-averaged radial velocity is given by

mesa. Such lava-¬‚ow morphologies are not con-

ρg H 2

H

1

U= v r dz = ’

sistent with the domal morphology of isoviscous

3·h ·c

H 0

gravity ¬‚ows. Motivated by the observation that

‚H

(·h + (·c ’ ·h )T av )

lava-¬‚ow deposits often lack a domal morphol- (6.40)

‚r

ogy analogous to ice sheets, Bercovici (1994) in-

Conservation of mass for this system can be writ-

vestigated a radially-symmetric 1D model of the

ten as

coupled evolution of mass and heat ¬‚ow in a

‚H 1 ‚(rU H )

¬‚uid spreading by gravity on a ¬‚at plane. In this + =0 (6.41)

‚t ‚r

r

model, strongly temperature-dependent viscosi-

ties cause most of the relief of the gravity ¬‚ow and heat transport is given by

to be accommodated along a cool, stiff gravity- ‚ T av ‚ T av 12κ

+U = ’ 2 T av (6.42)

¬‚ow margin, while the hot, central portion of the ‚t ‚r H

¬‚ow is able to ¬‚ow at very low gradients. Hence,

where κ is the thermal diffusivity. Substituting

the temperature-dependent viscosity of lava is re-

Eq. (6.40) into Eqs. (6.41) and (6.42), Bercovici ob-

sponsible for the wide, tabular nature of some

tained

lava ¬‚ows according to this model. In this sec-

‚H 1‚ ‚h

tion we review Bercovici™s model as an example = r (1 + νT av )H 3 (6.43)

‚t r ‚r ‚r

of modeling the coupled evolution of ¬‚uid ¬‚ow

and rheology in gravity ¬‚ows.

‚ T av ‚h ‚ T av T av

Bercovici (1994) assumed that the lava viscos- = (1 + νT av )H 2 ’2 (6.44)

‚t ‚r ‚r

ity, ·, was inversely proportional to its tempera- H

ture, T : where H has been nondimensionalized by

(3·c Q /2π ρg) (where Q is the ¬‚ux of lava

·h ·c to the surface), t by H0 /12κ, and r by

2

·(T ) = (6.37)

·h + (·c ’ ·h )T ( ρg H0 /12C κνc ). The only free parameter in

5

Eqs. (6.43) and (6.44) is the relative viscosity con-

where ·h and ·c are the viscosities of the hottest trast ν = (·c ’ ·h )/·h .

and coldest parts of the lava, respectively. His Here we consider the case of a constant rate

model also assumed that cooling takes place effusion, Q to the surface through a cylindrical

along the upper and lower boundaries of the ¬‚ow conduit of radius ri (Figure 6.8). This case corre-

where the temperature is kept at a constant ref- sponds to a constant-¬‚ux boundary condition at

erence value T = 0. This thermal boundary con- r = ri :

dition, together with the assumption of a uni- ‚H

r (1 + νT av )H 3 = 1 at r = ri (6.45)

form heat ¬‚ux with depth z, is consistent with a ‚r

132 NON-NEWTONIAN FLOW EQUATIONS

mer Laurentide Ice Sheet and alpine glaciers of

H

T=0 the Alaskan Brooks Range.

H z Before considering the extension of Eq. (6.33)

0 to 2D, we will investigate the spatial variability

T

constant

U

flux

of basal shear stress beneath modern ice sheets.

Q

Basal shear stresses are not uniform. However,

T=1

by characterizing the spatial variability observed

in real ice sheets, we can extend the threshold-

ri r

T=0 sliding model to include the observed variations

as a further development of the model. Recently

Fig 6.8 Schematic diagram of the lava ¬‚ow model with compiled datasets for the ice-surface and bed

coupled mass and heat ¬‚ow.

topography of Greenland (Bamber et al., 2001)