<< . .

. 26
( : 51)



. . >>

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
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)

<< . .

. 26
( : 51)



. . >>