<< . .

. 27
( : 51)



. . >>

and East Antarctica (Lythe et al., 2001) are avail-
able for this purpose. Figures 6.10a and 6.10b are
Also, the lava is assumed to erupt onto the sur-
face with its maximum temperature T av = 1 and grayscale images of the basal shear stresses in
r = ri . Greenland and East Antarctica calculated from
Eq. (6.4). The brightness scale for these images is
Figure 6.9 illustrates the ¬‚ow-thickness and
given in the lower left of the ¬gure. The scale
temperature pro¬les of a spreading lava ¬‚ow for a
varies from 0 (black) to 5 bars (white), with typ-
case of (a) no viscosity contrast and (b) a large vis-
ical values in the range of 0.5 to 3 bars. The av-
cosity contrast. For the case of isoviscous ¬‚ow, a
erage basal shear stress for Greenland and East
domal shape develops in which the surface slope
Antarctica, calculated from these maps, are al-
is steepest at the margin, but varies continuously
most identical: 1.41 bars for Greenland and 1.39
from the margin to the conduit radius, de¬ned
to be at r = 1. This solution is broadly similar to for East Antarctica. These ice sheets also share a
common pattern of spatial variation.
the parabolic geometry of a perfectly plastic ice
Figures 6.10a and 6.10b clearly show that basal
sheet on a ¬‚at bed. In the case of a strong viscos-
shear stresses increase with distance from di-
ity contrast, a more tabular ¬‚ow develops. In this
vides, with values close to 0.5 bars near divides
case, relief is accommodated along the cool, stiff
to an average of 3 bars near margins. Notably, the
¬‚ow margin, while the hot, ¬‚uid central portion
bed-elevation data for East Antarctica are interpo-
of the ¬‚ow lacks the stiffness to support signi¬-
lated from a much lower density of observations
cant surface slopes.
than for Greenland, especially south of 87—¦ . Nev-
ertheless, in both Greenland and East Antarctica,
6.4 Modeling of threshold-sliding where the data are most complete, this pattern
is clear.
ice sheets and glaciers over
Data from example pro¬les are plotted in Fig-
complex 3D topography ure 6.10c. Although there is signi¬cant scatter,
the basal shear stresses increase from 0.5 to an
Analytic solutions are generally restricted to 1D average of 3 bars in each of the pro¬les. The scat-
ice sheets and glaciers with relatively simple ter observed in these data is partly a function of
bed geometries. Real ice sheets and glaciers, of the high resolution of the source data. At these
course, are more complex. In this section we de- small scales (i.e. 5 km), variations in ice-surface
scribe a method for modeling the 2D geometry of slope occur that are not signi¬cant at the larger
an ice sheet or glacier above complex topography. scales of greatest interest in ice-sheet reconstruc-
The solution method relies upon a cellular algo- tions.
rithm that iteratively builds up discrete packets One means of incorporating the observed spa-
of ice until a threshold condition is achieved. The tial variability of basal shear stresses is to con-
model is illustrated using applications to modern sider them to be a function of ice-surface slope.
Greenland and Antarctic ice sheets and the for- In Figure 6.11a and 6.11b, the basal shear stresses
6.4 MODELING OF THRESHOLD-SLIDING ICE SHEETS AND GLACIERS OVER COMPLEX 3D TOPOGRAPHY 133



n = 100
1.5 (a) n = 0 t=1
t=1 (d)


H,T
1.0 H


T
0.5


0.0

t=3
t=3
1.5 (b) (e)


H,T
1.0


0.5


0.0

t = 10
t = 10
1.5 (c) (f)


H,T
1.0


0.5


0.0
5.0 0
4.0 5.0
4.0
0 2.0
1.0 2.0
1.0
3.0 3.0
r r
Fig 6.9 Plots of ¬‚ow thickness, H , and depth-averaged ues as slope values increase. Least-squares power-
temperature, T , as a function of time, t, for a ¬‚ow with function ¬ts to these data are indicated by the
(a)“(c) constant viscosity, and (d)“(f) temperature-dependent
solid lines in Figure 6.11. The least-squares ¬t for
viscosity with ν = 100.
Greenland is „ = 15S 0.55 , while the ¬t for East
Antarctica is slightly higher at „ = 18S 0.55 . Cer-
tainly a higher-order ¬t would characterize the
of Greenland and East Antarctica are plotted
data more precisely (because both data sets have a
as a function of ice-surface slope on logarith-
signi¬cant convex curvature in their dependence
mic scales. These plots illustrate that basal shear
on ice-surface slope), but a power-law function of
stresses increase from values as low as 0.5 bars
slope provides a good ¬rst-order approximation.
at low slope values to higher, more variable val-
134 NON-NEWTONIAN FLOW EQUATIONS



(a) (b)


3




1



2
4



basal shear stress

4
(c)
3
t (bars)




basal shear stress
2
2“5 bars
1.5“2
1
1“1.5
0“1
0
1.0
0.2 0.2
0.2 0.2
0
x/L (normalized distance from divide)

Fig 6.10 Basal shear stresses beneath Greenland and East
S = ∇h + ∇b to obtain
Antarctica, calculated from Eq. (6.4). Legend gives brightness
scale. (a) Greenland, calculated using datasets of Bamber „ (S)
|∇h + ∇b| = (6.46)
et al. (2001), (b) East Antarctica, calculated using datasets of
ρgh
Lythe et al. (2001). (c) Pro¬les of „ versus distance from the
divide for pro¬les 1“4. In these pro¬les, „ varies from about where
0.5 bars at divides to an average of 3 bars near the margins.
„ (S) = 15S 0.55 (6.47)
Basal shear stresses are also more spatially variable near
margins, as indicated by the large scatter on the right side of
to incorporate the spatial variations in basal
the plot.
shear stress observed in Figure 6.10. Equation
(6.46) then becomes
15
|∇h + ∇b|0.45 =
Our approach is to apply the empirical results (6.48)
ρgh
of Figure 6.10 directly to the threshold-sliding
model by replacing „ with „ (S) in Eq. (6.31). In Equation (6.48) can be transformed to a second-
2D, Eq. (6.31) applies along the direction of ¬‚ow order differential equation analogous to Eq.
lines. Since ¬‚ow lines are parallel to the local (6.33). The resulting equation cannot be solved
ice-surface gradient, the 2D version of Eq. (6.33) by relaxation methods, however, because the
is obtained by replacing S = dh/dx + db/dx with slopes vary both in direction and magnitude. This
6.4 MODELING OF THRESHOLD-SLIDING ICE SHEETS AND GLACIERS OVER COMPLEX 3D TOPOGRAPHY 135



add cell if S < t/rgh
101 (a)
t = 15S 0.55
grid continues

t
(bars)
Sy
grid continues
100 grid sweeps go from
Sx lowest to highest




10


y
(b)
101
x terminus cells

t = 18S 0.55
t
Fig 6.12 Illustration of the “sandpile” method for solving
(bars)
Eq. (6.46). At each iteration of the algorithm, a unit of ice is
100 added to each grid point within the area of ice coverage if the
addition does not violate the condition that S> „/ρg h
(where S is the ice-surface slope, equal to S x + S x ). The
2 2

ordering of the grid points is important in this algorithm. In
order to avoid oversteepenings, the sweep through the grid
10 should be from lowest to highest elevations.
10 10 10
S
lated as S = (S x + S y )1/2 , where S x is the downhill
2 2

Fig 6.11 Relationship of basal shear stress, „ , to ice-surface slope in the x-direction and S y is the downhill
slope for (a) Greenland, and (b) East Antarctica. Each point slope in the y-direction. During each ˜˜sweep™™ of
represents a pixel from Figure 6.10. Solid lines are the grid, the algorithm attempts to add ice to all
least-squares power-law ¬ts.
of the allowed grid points. The grid is swept re-
peatedly until no additional blocks can be added.
additional degree-of-freedom renders the relax- The end result is an exact solution to Eq. (6.46)
for any bed topography. „ can be taken to be uni-
ation method unstable in 2D.
Reeh (1982, 1984) developed a method for ap- form in this model or it can be a function of
plying the perfectly-plastic model to 2D ice sheets the ice-surface slope S. In addition to mimicking
by solving the 1D model on a curvilinear coor- the accumulation of ice in a thickening ice sheet,
dinate system that follows ¬‚ow lines. Equation the method is also analogous to the thickening
(6.48) can be solved more simply, however, us- of a sandpile until a threshold angle of repose is
ing a simple algorithm based on the accumula- reached. A fundamental difference between this
tion of discrete ˜˜blocks™™ of ice on a grid. This model and a real sandpile, however, is that the
method mimics the accumulation of ice thick- angle of repose is not uniform but is lower in
ness and slope until a threshold is reached. The areas of thicker ice.
algorithm is illustrated schematically in Figure Two additional aspects need to be speci¬ed to
6.12. Before the algorithm begins, grid points fully de¬ne the algorithm. First, in what order
within the ice margin are identi¬ed. These are should blocks be added to the grid? The order
the points where ice will accumulate; other grid matters because the addition of a block at one
points will remain ice free. The fundamental ac- grid point may increase the slope of one of its
tion item for each allowed grid point is simple: neighboring grid points to an oversteepened con-
add a discrete unit of ice thickness if the result- dition (i.e. one with a larger basal shear stress
ing surface slope is less than a threshold value, than the threshold value). To avoid this problem,
given by „ (S)/ρgh. The ice-surface slope is calcu- it is best to move through the grid in order of
136 NON-NEWTONIAN FLOW EQUATIONS


ascending elevations. To do this, a table is cre- ties in Greenland (Bamber et al., 2000) and Aust-
ated at the start of each grid sweep that indexes fonna, in eastern Svalbard (Dowdeswell et al.,
grid points from lowest to highest elevation. Ef¬- 1999). These studies have generally emphasized
cient methods based on the quicksort algorithm the important of the distance from the ice mar-
are available for this purpose (e.g. Press et al., gin in controlling the spatial distribution of ¬‚ow.
1992). Oversteepenings do not occur with this or- The interiors of thick ice sheets are generally
dering because slopes are dependent only on the characterized by relatively slow, distributed ¬‚ow
values of downslope points. Therefore, adding ice while outlet glaciers and ice streams are charac-
to higher elevations does not change the slopes terized by much faster and more focused ¬‚ow.
of grid points that have already been swept, and The ¬‚ow-distribution algorithm at the heart of
oversteepenings are avoided. The reconstructions Budd and Warner (1996) is similar to ¬‚uvial ¬‚ow-
in this chapter were obtained using this ordering distribution algorithms (e.g. Murray and Paola,
scheme. Another ordering scheme that works is 1994; Coulthard et al., 2002), suggesting that
to proceed from the outermost to the innermost glacial drainage networks may be modeled using
pixels. similar approaches developed for ¬‚uvial drainage
Second, how thick should the discrete ˜˜unit™™ networks (e.g. Willgoose et al., 1991).
of ice be for each iteration? In order to achieve Balance velocities are those required to main-
a smooth, precise ice-surface topography, the tain a time-independent ice-sheet topography.
unit should generally be less than or equal to The analogy of an ice sheet or glacier as a
1 m thick. However, for an ice sheet that is sev- conveyor belt is useful here; balance veloci-
eral kilometers thick, this would require adding ties are those that transport exactly the right
blocks of ice to each grid point thousands of amount of ice at each point to prevent thick-
times over, resulting in an inef¬cient method. ening or thinning of the ice at any point. For
The optimal approach is to start with a coarse threshold-sliding ice sheets and glaciers, balance
block size (e.g. 100 m) and reduce the size by a velocities are equivalent to instantaneous veloc-
constant factor (e.g. 3) whenever no more blocks ities because in a threshold-sliding ice sheet or
can be added to the grid. The model gradually re- glacier the entire ice mass continuously adjusts
duces the block size until a prescribed minimum to a change in accumulation or ablation. Al-
thickness (e.g. 10 cm) is reached. This approach though Budd and Warner™s method computes
builds up a coarse solution quickly in the early depth-averaged velocities rather than basal veloci-
phase of the run, gradually smoothing the rough ties, the two are nearly identical for threshold-
edges as the block size is reduced. sliding ice sheet and glaciers in which the ice-
For many glacial-geologic applications, ice ve- surface slope is much less than one. Nye (1951)
locities are an important variable to constrain provided equations for both depth-averaged and
and map. In glacial bedrock erosion, for ex- basal velocities in threshold-sliding ice sheets and
ample, basal-sliding velocities are likely to be glaciers. He found that the magnitude of the
the key variable controlling bedrock-erosion rates correction term relating depth-averaged veloci-
(Andrews, 1972; Hallet, 1979, 1996). Therefore, a ties to basal velocities is on the order of the ice-
better understanding of the formation of ¬nger surface slope, and hence is very small for most
lakes, U-shaped valleys, cirques, and other glacial- applications.
erosional landforms will likely require a method In the Budd and Warner algorithm, the ice-
for calculating the velocity ¬eld of former ice- surface topography, ice thickness, and rates of ac-
sheets and glaciers in 2D. cumulation and ablation are used as input. Local
Budd and Warner (1996) developed an algo- accumulation rates are added to each grid point
rithm for computing the depth-averaged balance and incoming ice is distributed to the downslope
velocities of ice sheets. These authors illustrated nearest neighbors. The amount of ice distributed
their method by computing a map of balance ve- to each downslope grid point is proportional
locities for East Antarctica. This important work to the ice-surface slope in the direction of that
quickly inspired calculations of balance veloci- grid point. The input and output through each
6.4 MODELING OF THRESHOLD-SLIDING ICE SHEETS AND GLACIERS OVER COMPLEX 3D TOPOGRAPHY 137


grid point is calculated using a highest-to-lowest Figure 6.13 illustrates the bed and ice-surface
ordering scheme. This scheme ensures that all topography for modern Greenland (Bamber
of the contributing discharge from upstream et al., 2001). Areas below sea level are indicated in
sources is accumulated at each grid point before black. This and all other reconstructions in this
it is distributed to downslope neighbors. In the ¬- chapter require three inputs: a DEM of bed to-
nal step of the algorithm, the local ice discharge pography, a ˜˜mask™™ grid that de¬nes the areas of
is divided by the local ice thickness to obtain the the ice margin, and a function de¬ning the basal
shear stress, such as „ = 15S 0.55 . Appendix 5 pro-
depth-averaged velocity for each grid point.
As a complement to our reconstructions of vides a C code that implements the sandpile and
ice-surface geometry, we have used the Budd balance velocity method for the Greenland and
and Warner method to calculate relative depth- Antarctic ice sheets.
averaged velocities for each of our reconstruc- The bed topography is input directly into the
tions. For these calculations we have assumed model. The ice-thickness data, derived from the
that accumulation is uniform and that there difference between the ice-surface and bed topog-
is no ablation. Although these assumptions are raphy, is used to provide a binary ˜˜mask™™ grid
clearly invalid, the resulting maps are useful be- de¬ning the grid points that are allowed to accu-
cause ¬‚ow is so strongly controlled by the dis- mulate ice in the sandpile model. The mask grid
tance from divides. Therefore, even if accumu- has values of 1 where the ice thickness is greater
lation rates differ by a factor of ten or more than 0 (i.e. areas covered by ice), and values of
from one side of an ice sheet to the other, the 0 where the thickness is 0 (i.e. no ice coverage).
basic pattern of ¬‚ow is essentially unchanged. For this reconstruction we used the shear-stress
relationship „ = 15S 0.55 observed in Figure 6.11.
This insensitivity is analogous to the spatial varia-
tions in discharge within a river basin. Discharge Figure 6.13c is a shaded-relief and contour
in a river network is predominantly a function map of the solution to Equation (6.46) obtained
of drainage area; variations in precipitation and with the sandpile model. Figure 6.13c has been
runoff within a basin play a much less signi¬cant constructed with the same vertical exaggeration
role. The ¬‚ow maps we have created also serve (30—) as Figure 6.13b to provide a direct, side-by-
to indicate the extent to which the ice ¬‚ow is side comparison. The similarity of the location
controlled by the subglacial topography. For ice and shape of the contours indicate that the over-
sheets in which the thickness is much greater all solution is in good agreement with the ob-
that the subglacial topographic relief, ice ¬‚ow is served topography. The location of the major di-
only weakly controlled by the topography. In con- vide, offset from center by approximately 100 km,
trast, if the subglacial-topographic relief is com- is also reproduced in the model. One major differ-
parable to the ice thickness, ¬‚ow can be strongly ence between the observed and modeled topog-
focused into troughs and depressions beneath raphy, however, is the steepness of divides in the
the ice. model solution. The model topography has a sig-
In this section we present reconstructions of ni¬cantly more angular appearance than the ob-
the ice-surface topography and ¬‚ow of the mod- served topography. This angularity can be traced
ern Greenland and Antarctic ice sheets using the to the poorness of ¬t between the power-function
sandpile and Budd and Warner (1996) algorithms. and the data of Figure 6.11a. In the lower-left cor-
Greenland is perhaps the most important case ner of the plot, representing areas of low slope,
study in this chapter because the bed topogra- the data fall far below the power function, in-
phy is both well constrained and has a signif- dicating that the threshold basal shear stress
icant in¬‚uence on the morphology of the ice near divides is signi¬cantly lower than the val-
sheet. The principal divide in Greenland is off- ues predicted by the power function. As a result,
set from central Greenland by approximately 100 the model solution overestimates the ice-surface
kilometers. This asymmetric pro¬le re¬‚ects a bed- slopes near divides by as much as a factor of 2.
topographic control associated with higher bed This does not introduce signi¬cant errors in the
elevations in eastern Greenland. elevations of these regions because the slopes in
138 NON-NEWTONIAN FLOW EQUATIONS



(a) (b) (c) (d)


2 km
2 km

3 km
3 km




ice topography
ice topography depth-averaged
(modeled)
(measured)
bed topography velocity

Fig 6.13 Reconstruction of the modern Greenland Ice
Sheet using Eq. (6.46) and the observed correlation inate ice shelves and areas of sea ice from the
„ = 15S 0.55 from Figure 6.11. Bed and ice-surface topography reconstruction.
are illustrated with shaded-relief (30 — vertical exaggeration)
Here we used the shear--stress relationship
grayscale maps. (a) Input bed topography from Bamber et al.
„ = 18S 0.55 observed in Figure 6.11b. For com-
(2001). Darkest areas are below sea level. (b) Ice topography
parison, we also performed a reconstruction us-
observed from radar interferometry and given by Bamber
ing the relationship of „ = 15S 0.55 observed for
et al. (2001). Contours are 1 km in spacing (1 km contour is
Greenland. Using „ = 18S 0.55 results in a slightly
too close to the margin to be easily seen). (c) Numerical

<< . .

. 27
( : 51)



. . >>