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