. 8
( 14)


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
h i+1 ’ 2h i + h i’1 + bi+1 ’ 2bi + bi’1
2 „0
2x h 0.8
= ’( x) 2 2 3
ρ g h i (h i+1 ’ h i’1 + bi+1 ’ bi’1 ) (km)

where x is the resolution cell size of the grid.
This equation is solved iteratively with 0.2
L i (h iold ) 20 40 100
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
„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

Figures 6.6 and 6.7 illustrate example solutions
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

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
U= v r dz = ’
sistent with the domal morphology of isoviscous
3·h ·c
H 0
gravity ¬‚ows. Motivated by the observation that
(·h + (·c ’ ·h )T av )
lava-¬‚ow deposits often lack a domal morphol- (6.40)
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
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
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
·(T ) = (6.37)
·h + (·c ’ ·h )T ( ρg H0 /12C κνc ). The only free parameter in

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

mer Laurentide Ice Sheet and alpine glaciers of

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
of basal shear stress beneath modern ice sheets.
Basal shear stresses are not uniform. However,
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)
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

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

1.0 H



1.5 (b) (e)




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



5.0 0
4.0 5.0
0 2.0
1.0 2.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-

(a) (b)




basal shear stress

t (bars)

basal shear stress
2“5 bars
0.2 0.2
0.2 0.2
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
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
|∇h + ∇b|0.45 =
Our approach is to apply the empirical results (6.48)
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

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

grid continues
100 grid sweeps go from
Sx lowest to highest


x terminus cells

t = 18S 0.55
Fig 6.12 Illustration of the “sandpile” method for solving
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
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

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

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

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

2 km
2 km

3 km
3 km

ice topography
ice topography depth-averaged
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
better match to the modern ice-surface topogra-
reconstruction of the ice-surface topography (with same
phy than we obtain by using the Greenland rela-
grayscale and shading as in (b)). The divide position and
elevations are a good match to (B) except that the divides are tionship, which yields a maximum elevation of
too peaked. (d) Grayscale map of depth-averaged velocities just under 4 km.
computed using the balance-velocity algorithm of Budd and
These reconstructions illustrate that the sand-
Warner (1996) assuming uniform accumulation. The grayscale
pile model is capable of accurately predicting
map is scaled logarithmically.
ice-sheet geometries if the threshold basal shear
stress is well constrained. However, the principal
either case are small. The ¬‚ow map for Green- objective of glaciological work within the con-
land is indicated in Figure 6.13d and is similar text of geomorphology is to reconstruct former
to the map provided by Bamber et al. (2000b). The ice sheets and glaciers.
grayscale map is scaled logarithmically in this Reconstructions of the Laurentide Ice Sheet
and all subsequent ¬‚ow maps to illustrate the (LIS) have had a long history that includes many
internal ¬‚ow distribution. A linear map, in con- classic contributions. However, despite decades of
trast, would show signi¬cant velocities only at effort there is still vigorous debate over the mean
the outlet areas. thickness of the LIS and the locations of its major
Input data for the East Antarctic Ice Sheet re- domes. These questions bear on many important
construction (Figure 6.14) includes the DEM of issues, including the magnitude of Late Quater-
bed topography provided by Lythe et al. (2001). nary sea-level change. Numerical reconstructions
Lythe et al.™s ice-thickness data were also used to of the Laurentide Ice Sheet at the Last Glacial
derive a mask of ice coverage. Rather than includ- Maximum have generally reproduced a relatively
ing all areas of ¬nite ice thickness in our mask, thick (≥ 2 km) ice sheet with a central east-west
however, we used a threshold ice thickness of 150 divide (e.g. Denton and Hughes, 1981). In con-
meters instead. This elevation was used to elim- trast, Peltier™s (1994) ICE4G relative sea-level (RSL)

(a) (d)

bed topography velocity

(b) (c)

4 km
4 km
3 km
3 km

2 km
ice topography
ice topography 2 km

Fig 6.14 Reconstruction of the modern Greenland Ice
to account for bed de¬‚ection (e.g. Tomkin and
Sheet from Eq. (6.46) with „ = 18S 0.55 . (a)“(d) are identical
Braun, 2002). The ratio of 0.27 comes from the
to the corresponding images for Greenland in Figure 6.13.
relative density of ice and mantle. This approach
is straightforward to incorporate in the sandpile
model by replacing h with (1 + ρi /ρm )h in the
inversion reproduces a thinner ice sheet with a
more complex ice-surface topography. threshold condition. It should be noted, however,
An inherent uncertainty in many ice-sheet that this approach is valid only for ice sheets
reconstructions is the amount of isostatic re- larger than the ¬‚exural wavelength. Smaller ice
bound experienced by the bed between the sheets and glaciers may be partially or fully ¬‚ex-
time of glaciation and the present. This uncer- urally compensated.
tainty is not signi¬cant for areas that have been Including the effects of isostasy is a much
unglaciated for several tens of thousands of years. greater challenge for LGM ice sheet reconstruc-
In these cases, we can be sure that the modern to- tions because postglacial rebound has not yet run
pography includes no dynamic component: post- its course. It is dif¬cult to estimate the amount
glacial rebound has run its course. It is common of postglacial rebound that has taken place or
in these cases to ignore the time delay of man- to estimate the amount left to occur because
tle response and add an instantaneous de¬‚ection the dynamics of postglacial rebound depend on
amplitude proportional to the local ice load: for the size and shape of the load through time.
every kilometer of ice thickness, approximately In particular, we cannot extrapolate the present
270 meters of ice are added to the thickness rate of sea-level rise along an exponential curve

Fig 6.15 Reconstruction of the Laurentide Ice Sheet at Last between 0 and 27%. If we believe that postglacial
Glacial Maximum (18 ka) assuming (a) „ = 18S 0.55 (referred
rebound beneath the LIS is nearly complete, we
to as standard basal shear) and (b) „ = 9S 0.55 (one-half the
should choose a value close to 0.27. In contrast, if
standard basal shear). Bed and ice-surface topography are
little postglacial rebound has occurred, the value
illustrated with shaded-relief (100— vertical exaggeration)
should be closer to zero. For our Laurentide re-
grayscale maps. (c) Grayscale map of depth-averaged
constructions, we will assume a ratio of 0.2. This
velocities for the standard basal shear model assuming
value assumes that the majority of postglacial re-
uniform accumulation.
bound has taken place. The uncertainty of this
value, which could plausibly vary between 15 and
25%, introduces a 10% uncertainty into determi-
because the Maxwell-relaxation time of the man-
nation of the LIS thickness.
tle depends on the wavelength of the load (e.g.
Two reconstructions of the LIS and LGM are
Turcotte and Schubert, 2002). Therefore, post-
presented in Figure 6.15. Figures 6.15a and 6.15c
glacial rebound is not a simple function of the
illustrate the ice-surface topography and ¬‚ow us-
local ice thickness but is also controlled by the
ing the basal shear stress observed in East Antarc-
large-scale morphology of the ice sheet, particu-
tica („ = 18S 0.55 ). We refer to this shear--stress re-
larly at long time scales.
lationship as the ˜˜standard™™ basal shear model,
Nevertheless, we can take a similar approach
because it characterizes the observed trends in
to including isostatic de¬‚ection that is often
the Greenland and Antarctic Ice Sheets (Figure
taken in cases with no dynamic topography.
6.11). For comparison, we have also provided a
Rather than adding 27% of the ice thickness to
second reconstruction which uses a shear stress
account for de¬‚ection, we can add some fraction

equal to one half of the standard basal shear, or and south between Figures 6.15a and 6.15b, but
„ = 9S 0.55 . otherwise the ¬‚ow pattern is unchanged.
The DEM of bed topography was constructed Figures 6.16 and 6.17 illustrate the ice-surface
by projecting the ETOPO5 DEM (Loughridge, topography and ¬‚ow of the Laurentide Ice Sheet
1986) for North America to a Lambert Equal- at 14, 13, 12, 11, 9, and 8 ka. These reconstruc-
Area projection. The ETOPO5 data include tions do not include ice coverage between north-
both bathymetry and topography data. The ern Canada and Greenland (Dyke et al., 2002), and
bathymetry is important in this case because are therefore only approximate in their northern-
the LIS covered many areas that are now sub- most regions.
merged. The ice margins were digitized from The reconstructions of Figures 6.16a--6.16f are
Dyke and Prest (1987) by Eric Grimm of the characterized by a thick central region that
Illinois State Museum (personal communication, changes little as it shrinks in size. By 11 ka the
2002) and recti¬ed to the DEM by the author. The ice sheet has shrunk to less than half its orig-
bed and ice-surface topography are illustrated inal area (including coverage in the Cordillera)
with shaded-relief (100— vertical exaggeration) but still maintains a signi¬cant central region
grayscale maps in Figure 6.15. The bathymetry in with elevations above 3 km. Only by 9 ka has the
these images has been removed so that the coast- overall shape of the ice sheet and the altitude of
lines could be indicated instead. Use of the stan- its central region changed signi¬cantly. At this
dard basal shear model yields an ice sheet with point, the central region of the ice sheet has col-
a central divide striking east-west. The average el- lapsed into three distinct domes. This collapse is
evation and thickness of the ice sheet in Figure followed by the rapid deglaciation of Hudson Bay
6.15a (the standard basal shear case) are 2211 m between 9 and 8 ka.
and 2011 m. Decreasing the shear stress preserves An important caveat should be given regard-
the overall shape of the ice sheet but results in a ing the reconstruction of the ice in the Cana-
migration of the divide to the southwest where dian Cordillera. In alpine terrain it is important
the ice is propped up by the topographic front to use DEM data that enable the model to es-
of the Canadian Rockies. The average elevation timate bed slopes accurately. It is unlikely that
and thickness of Figure 6.15b are 1482 m and the ETOPO5 dataset, with a resolution of 5 km,
1134 m, respectively. Although the average eleva- is an adequate representation of the bed topog-
tion and thickness of the standard reconstruc- raphy in the Cordillera. Coarse DEMs such as
tion are a relatively modest 2.2 and 2.0 km, re- ETOPO5 may underestimate bed slopes, especially
spectively, much of the central core of the ice in high-relief terrain where closely spaced ridges
sheet is quite high, with average elevations above and valleys may not be resolved. This problem
3 km and greater than 3.5 km in thickness. As a can be minimized by using DEMs of the highest-
caveat, it should be noted that we did not re- available resolution. Nevertheless, some size re-
move the post-glacial sediment from Hudson Bay duction may be necessary to maintain reason-
or the Great Lakes prior to this reconstruction. able grid sizes and computing times. When size-
This will not affect the large-scale geometry of reduction must be done, it should be done by
the ice sheet signi¬cantly. Certain areas, however, subsampling the DEM rather than by averaging
including the Great Lakes, will have ice overes- neighboring pixels. Subsampling will tend to pre-
timated thicknesses due to the postglacial sedi- serve the small-scale topographic relief in the
ment ¬ll in the DEM. A more precise reconstruc- DEM better than averaging. In the Laurentide re-
tion would require that any post-glacial sediment construction of Figure 6.15, the thickness of the
be removed from these areas to re¬‚ect the topog- Cordilleran Ice Sheet may be overestimated be-
raphy that the ice sheet experienced. cause of the coarse scale of the DEM (5 km). Appli-
The results of Figure 6.15 indicate that even cations focused on extensive areas of high-relief
though the thickness of the LIS is sensitive to the terrain should be reproduced on multiple scales
basal shear stress, the overall shape is not sensi- to verify that the results are independent of the
tive. The principal divide migrates to the west DEM scale.

(a) (b)

2 km
2 km
3 km
3 km

13 ka
14 ka

(c) (d)

2 km

3 km
3 km

11 ka
12 ka

(e) (f)

2 km

2 km
3 km

8 ka
9 ka
Fig 6.16 Reconstruction of the ice-surface topography of the Laurentide Ice Sheet from 14“8 ka assuming the standard basal
shear model and using the margin positions of Dyke and Prest (1987). Bed and ice-surface topography are illustrated with
shaded-relief (100— vertical exaggeration) grayscale maps. Ice coverage between Greenland and northern Canada has not been
included in this reconstruction.

(a) (b)

13 ka
14 ka

(c) (d)

11 ka
12 ka

(e) (f)

8 ka
9 ka

Fig 6.17 Depth-averaged ice velocities corresponding to
the reconstructions of Figures 6.16a“6.16f.

ment in which only a portion of the ice sheet
can be modeled at one time. In addition, the Fin-
Next we illustrate the model application to
ger Lakes are an example where both the bed
a high-resolution reconstruction of the Lauren-
topography and the distance from the ice mar-
tide Ice Sheet in the Finger Lakes Region of New
gin have an important in¬‚uence on ice-surface
York State. This example will illustrate the appli-
topography and ¬‚ow. In contrast, large ice sheets
cation of the model to an ice-marginal environ-

and alpine glaciers are predominantly in¬‚uenced topography. We will not model the full dynam-
by only one of these variables. For example, the ics of Finger Lakes formation in this chapter but
local topography and ¬‚ow in large ice sheets will rather determine what the ice-surface topog-
is primarily a function of distance from the raphy and ¬‚ow looked like after the Finger Lakes
margin, while the thickness and ¬‚ow in alpine had formed.
glaciers are almost entirely controlled by bed The reconstruction of this region requires the
topography. following inputs: bed topography (Figure 6.18a),
The Finger Lakes Region has been the focus a mask grid incorporating the position of the
of many glacial-geologic studies, including glacio- Valley Heads Moraine, and the threshold basal
logical reconstructions (Ridky and Bindschadler, shear stress. The topography of the region was
1990), stratigraphic studies of meltwater produc- obtained by joining several 90-m resolution USGS
tion and lake in¬lling (Mullins and Hinchey, DEMs. The position of the Valley Heads Moraine
1989), and geomorphic studies of the Finger Lakes was used to create the mask grid by overlay-
themselves (von Engeln, 1956). Geomorphically, ing the sur¬cial geologic map of the area (New
the Finger Lakes Region is dominated by subpar- York State Geological Survey, 1999) onto the bed
allel, glacially-scoured troughs with their south- topography to delineate the ice margin. At the
ernmost extents in Seneca and Cayuga Lakes, the spatial scale of the Finger Lakes, the relation-
two largest of the Finger Lakes. The ¬ve largest ship between basal shear stress and ice-surface
troughs of the region comprise the Finger Lakes slope observed in Greenland and East Antarc-
proper, but there are numerous other troughs cut tica is not useful because ice-surface slopes this
into the Allegheny Plateau of smaller size that close to the margin are uniformly steep. In ad-
are not deep enough to be enclosed depressions. dition, the average basal shear stress observed
Figure 6.18a is a shaded-relief image of the topog- for modern ice-sheet margins (3 bars) may also
raphy of the region. The troughs vary in spacing be inapplicable because shear stresses near ice-
from 10 to 30 km along strike, with the greatest sheet margins are spatially variable and may be
spacing between Seneca and Cayuga Lakes. The much smaller than 3 bars locally. In lieu of a bet-
southern tips of the Finger Lakes coincide with ter constraint on basal shear stresses in the re-
the Valley Heads Moraine (14 ka). As such, the gion, we have used the traditional end-member
scouring of the Finger Lakes Region most likely values of 0.5 and 1.5 bars in two alternative re-
took place when the ice margin was coincident constructions. In addition, we have modi¬ed the
with this moraine, although several phases of ice-¬‚ow portion of the reconstruction so that all
glaciation may have contributed to their forma- of the ¬‚ow originates at the top boundary of the
tion. grid. This modi¬cation represents the incoming
Little is known about the processes and dy- ice ¬‚ow from the Laurentide Ice Sheet, which we
namics of the scouring of glacial troughs in gen- will assume dominates over any locally-derived
eral or the Finger Lakes in particular. The regular accumulation.
spacing of the Finger Lakes suggests a positive- The results for the ice-surface topography and
feedback or instability mechanism in which in- depth-averaged ¬‚ow are given in Figures 6.18b--
cipient depressions in the bed topography focus 6.18c and 6.18d--6.18e for 1.5 and 0.5 bars, respec-
ice ¬‚ow, resulting in enhanced deepening and tively. These values bracket an important tran-
focusing of ice ¬‚ow. In order for this model to sition in ice-sheet behavior. If the basal shear
be valid, the ice-surface topography in the region stress is 1.5 bars or larger, the topography and
must re¬‚ect variations in bed topography in or- ¬‚ow are only weakly dependent on bed topogra-
der for ice ¬‚ow to be focused into troughs. The phy. Figure 6.18b illustrates the ice-surface topog-
¬rst step towards testing this hypothesis is to re- raphy in this case. The reconstructed ice sheet
construct the ice-surface topography and ¬‚ow in in this ¬gure is dominated by a forked central
the region to determine whether the ice-¬‚ow pat- divide and distributary ice ¬‚ow. In contrast, for
terns were likely to have been focused by the bed basal shear stresses of 0.5 bars or smaller, the

(b) 3 km (c)

2 km

1 km


ice topography
1.5 bars

(d) (e)

1 km

bed topography
30 km

0.5 bars ice topography

Fig 6.18 Reconstruction of the ice-surface topography and
depth-averaged velocities in the Finger Lakes Region. Two
Finger Lakes be removed to more accurately re-
reconstructions are presented with uniform basal shear
¬‚ect the bed topography experienced by the ice
stresses of (a) and (c) 1.5 bars and (d) and (e) 0.5 bars. (a)
Shaded-relief image of bed topography. (b) Shaded-relief image sheet. Also, drumlin orientations can be used to
of ice-surface topography and (c) depth-average velocity provide a constraint on the ¬‚ow pattern and,
assuming 1.5 bars. For this case, representing a relatively rigid in conjunction with reconstructions of different
base, ice topography and ¬‚ow are only weakly controlled by
shear stresses, provide an indirect constraint on
subglacial topography. (d) Shaded-relief image of ice-surface
the appropriate value of shear stress. Ridky and
topography and (e) depth-average velocity assuming 0.5 bars.
Bindschadler (1990), for example, used drumlin
orientations to guide their 2D ¬‚ow line recon-
ice-surface is strongly controlled by the subglacial structions in this area.
topography. Ice ¬‚ow is focused into the troughs One additional point should be made regard-
of the Finger Lakes in this case. If the Finger Lakes ing this example. Although the solution of par-
were formed by focused ice ¬‚ow, as the positive- tial differential equations such as Eq. (6.46) typi-
feedback model requires, these reconstructions cally require the application of boundary condi-
suggest that the basal shear stress was probably tions at the edges of the grid, the sandpile algo-
closer to 0.5 bars than to 1.5 bars. Future work rithm does not require boundary conditions at
will require that the postglacial sediment in the grid boundaries in upslope directions. The reason

Fig 6.19 Reconstruction of the ice
(a) 30 km thickness for the Wisconsin Brooks
Range glaciation assuming uniform
basal shear stresses of 1 bar. (a) Bed
topography and ice margins mapped
by the Alaska PaleoGlacier Atlas
Group (data available at
data intro.html). (b) Grayscale map
of ice thickness, with brightness
scale at lower left.
bed topography and ice margin


ice thickness
1 bar

300“400 m

boundary conditions are not required is that the margins. We have assumed a uniform threshold
algorithm refers only to downslope grid points. shear stress of 1 bar for this reconstruction.
This is useful for applications in ice-marginal ar- In the alpine environment, the bed-slope com-
eas because only a portion of the ice sheet need ponent of the basal shear stress is dominant. As
be considered. such, variations in ice thickness are small and are
One of the bene¬ts of the threshold-sliding associated with areas of bed curvature: thicken-
model is that it can be broadly applied to the re- ing occurs in areas of bed concavity and thinning
construction of any ice sheet or glacier in the con- occurs in areas of bed convexity. Rapid changes
tinuum from large ice sheets to alpine glaciers. in glacier thickness are, therefore, restricted to
As a robust and broadly applicable model, it can slope breaks and glacier snouts.
provide preliminary reconstructions for any area Figure 6.19b presents the results of our re-
and requires the constraint of just a single input construction based on the input data in Figure
variable, the basal shear stress. 6.19a. Ice is up to several hundred meters thick
As an example of the model application in glacial valleys and is generally thicker where
in an alpine environment, we will recon- bed slopes are gentler, except for areas close to
struct the Wisconsin glaciation of the central glacier snouts where thicknesses taper. Modify-
Brooks Range of Alaska using input data com- ing the assumed basal shear stress from 1 bar
piled by the Alaska PaleoGlacier Group (http:// to a greater or lesser value leads to thicker or
instaar.colorado.edu/QGISL/data intro.html). This thinner glaciers but does not signi¬cantly change
data includes a 100-m resolution DEM of bed to- the relative thicknesses illustrated in Figure
pography and a vector dataset of Wisconsin ice 6.19b.

In this reconstruction, steep areas in the high- ground surface
lands of the Brooks Range are predicted to have
thin ice, as indicated by the dark tones in the
grayscale map in Figure 6.19b. However, it is likely
that below some minimum thickness, ice cannot
form and locally the area will be ice free. decollement
For example, in steep terrain the snow depth
may be insuf¬cient to insulate deeper layers from
summer melting or it could be too thin to initi-
ate compaction. To account for this, the alpine re- (b)
constructions could be post-processed to remove H
ice from areas that have insuf¬cient thickness q
for ice formation. In such cases, the reconstruc-
tion should then be recomputed with the revised
mask to provide an improved reconstruction.

6.5 Thrust sheet mechanics
Fig 6.20 Schematic diagram of the 2D thrust sheet model.
Deformation within Earth™s crust takes place by
brittle fracture and slip along faults and joints. the shape of the fold depend on the yield stress
As such, it is not appropriate, in detail, to model and fault geometry? In this example we prescribe
the crust as a ¬‚uid. At large scales, however, plas- the tectonic geometry and slip rate, but it is also
tic ¬‚ow models do a surprisingly good job at re- possible to prescribe the stress state along the
producing the observed geometries of structural fault and allow the slip rate to respond to both
folds. The reason why plastic ¬‚ow models are suc- the tectonic stress and the overburden stress, al-
cessful is that they represent the collective ef- lowing the fault to slow down as overburden ac-
fects of slip along a complex system of faults and cumulates.
joints. As a rock is stressed, it deforms very little The thickness of rock above the fault surface
until the yield stress of the rock is exceeded and is governed by Eq. (6.46) with a uniform yield
a new fault or joint is formed. Subsequent defor- stress:
mation is accommodated by slip along faults and
h|∇h + ∇b| =
joints when the threshold frictional force is ex- (6.49)
ceeded. Rock fracture and slip along a network
of faults are both threshold-controlled motions. Equation (6.49) can be solved with the sandpile
As such, a perfectly plastic or threshold-sliding algorithm using a very similar approach to the
model is often a good model for the large-scale one used for modeling 2D ice sheets over com-
deformation of the crust. plex topography. The advance of the fault tip is
In this section we consider a simple exam- modeled by adding mass to the system directly
ple of thrust-sheet displacement using a per- above the fault tip at a rate equal to the vertical
fectly plastic model. We assume the existence of component of slip. Appendix 5 provides a C code
a crustal ramp of angle θ and width L separating that implements this approach.
the ground surface from a decollement surface at Figure 6.21 illustrates the fold evolution for
a case with „0 = 20 bars, ρ = 2700 kg/m3 , H =
depth H (Figure 6.20). Motion along the thrust
3 km, L = 300 km, θ = 2—¦ , and s = 1 cm/yr at
fault takes place with slip rate s . The key ques-
three time increments t = 10 Myr, 20 Myr, and
tions we will ask are: (1) what is the 2D shape of
the fold above the thrust fault, and (2) how does 30 Myr following the initiation of thrusting. The

(a) 3
100 km 100 km
(km) 2 (km) 1

t = 10 Myr t0 = 5 bar
4 3
(km) 2 (km) 1

t = 20 Myr L = 150 km
(c) (c)
(km) 2 (km) 1

t = 30 Myr q = 1° 100 200
0 300
200 400
x (km) x (km)

Fig 6.21 Shaded-relief images and cross-sectional plots of Fig 6.22 Sensitivity of the thrust sheet system to variations
in the (a) yield stress („0 = 5 bars compared to 20 bars for
the fold topography developed over a thrust fault for a
the reference case), (b) fault width (L = 150 km instead of
perfectly-plastic model of crustal deformation, assuming
„0 = 20 bar, ρ = 2700 kg/m3 , H = 3 km, L = 300 km, 300 km), and (c) fault angle (θ = 1—¦ ).
θ = 2—¦ , and s = 1 cm/yr at three time increments t = 10 Myr,
20 Myr, and 30 Myr following the initiation of thrusting.
total shortening equals the width of the thrust
sheet does the fault achieve a true steady-state
maximum structural relief in the model has also taper.
been restricted to be 2H (i.e. the fault becomes The sensitivity of fold geometry to model pa-
horizontal, forming a large syncline, once the rameters is shown in Figure 6.22. Each panel
fault tip has achieved an elevation equal to H shows the shaded relief and cross section of the
fold at t = 10 Myr for a case with (a) a low yield
above the ground surface). The shaded-relief im-
stress („0 = 5 bars), (b) a low yield stress („0 =
ages show the topography of the fold (given by
5 bars) and a narrow fault (L = 150 km), and (c)
the sum of the fault-surface elevation b and the
a low yield stress („0 = 5 bars) and thrust angle
rock thickness h). Early in the simulation, fault
(θ = 1—¦ ). Lower yield stresses lead to thinner rock,
displacement creates an asymmetric fold charac-
terized by a steeper backlimb than frontlimb. The as expected. Out in front of the fault, this thin-
forelimb progrades out in front of the fault to a ning has the effect of creating a broader zone of
maximum length controlled by the structural re- deformation and hence a higher degree of oro-
lief and the yield stress. Early in the simulation, genic curvature for faults with otherwise equal
the forelimb has a linear planform shape. Over geometries and shortening amounts.
time, the curvature migrates from the edge of the To date, most ¬‚uvial landform evolution
orogen into the center. Only when the amount of models have assumed relatively simple tectonic

forcings. When more realistic tectonic forcings are calculated based on the reconstructed ice-
have been considered (e.g. Willett et al., 2000), surface topography. Our approach is to model the
the effects on mountain-belt geometry have been geometry and ¬‚ow in ice sheets and the resultant
shown to be very important. Although simpli¬ed, subglacial erosion iteratively through time. Fig-
the perfectly plastic model of topographic evolu- ure 6.23 illustrates the basic steps of the model
tion above a crustal-scale fault provides a start- and the input and output data at each step us-
ing point to quantify rock uplift rates above pre- ing the New York State Finger Lakes as the ex-
scribed structural elements. Coupling this model ample application. In this example, we use the
with bedrock erosion processes could help to fur- modern Finger Lakes topography as input to de-
ther unravel the relationships between tectonics termine the bedrock-erosion pattern that would
and erosion in thrust systems. occur if an ice sheet were to advance over the
area today. Input data for step 1 of the model
includes (1) the bedrock topography, (2) the po-
sition of the ice margin, and (3) the thresh-
6.6 Glacial erosion beneath ice
old basal shear stress for motion. These inputs
sheets are used to reconstruct the geometry of the ice
sheet using the threshold-sliding model and the
Bedrock erosion beneath ice sheets is the domi- sandpile method. Following the ice-surface recon-
nant agent of Quaternary erosion and landform struction, the lithospheric de¬‚ection under the
evolution in most of the Northern Hemisphere. ice load must be considered if the study area
Despite the dominance of subglacial erosion in is of large extent. If the study area is less than
≈ 100 km in width, the lithosphere is assumed
many parts of the world, our understanding of
subglacial erosional landforms and the feedback to be perfectly rigid. If the area is signi¬cantly
processes between ice ¬‚ow, bedrock erosion, and larger than 100 km the de¬‚ection of the litho-
bedrock topography beneath ice sheets is limited. sphere should be determined using the 2D ¬‚ex-
For example, what processes and feedbacks con- ure equation and estimates for the local elastic
trol the morphology of overdeepenings beneath thickness of the crust. If de¬‚ection is computed,
glaciers and ice sheets? How do these processes the model returns to step 1 where the ice-surface
differ on the margins of ice sheets from their in- topography is recomputed taking into account
teriors? What processes control the size and dis- the lithospheric de¬‚ection in the bed topogra-
tribution of regularly spaced glacial lakes such phy. Although 100 km is given as a rule-of-thumb,
as the New York State Finger Lakes? Although the precise threshold length scale where ¬‚exu-
considerable work has focused on modeling land- ral compensation becomes signi¬cant depends
form evolution in the alpine-glacial environment on the local elastic thickness.
(e.g. Harbor, 1992; Braun et al., 1999; MacGregor In step 2, basal-¬‚ow velocities are determined
et al., 2000), previous work has just touched on with the balance-velocity method of Budd and
landform evolution beneath ice sheets. In ad- Warner (1996). The balance-velocity method re-
dition to providing a better understanding of quires that the pattern of accumulation and abla-
glacial-landform evolution, determining the con- tion be speci¬ed. For many paleo ice sheets, this
trols on glacial-lake formation could help to con- pattern may not be well constrained. For these
strain paleo ice-sheet thickness, a critical param- cases, the accumulation/ablation pattern may be
eter that is dif¬cult to determine using exist- modeled using an equilibrium line altitude (ELA)
ing methods such as relative sea-level inversion and an accumulation/ablation gradient. For ex-
(Peltier, 1994). ample, assuming an ELA of 1 km and an accumu-
In this example application we use the lation gradient of 0.1 m/yr per km, the ablation
rate at an elevation of 500 m is ’0.5 m/yr. In step
threshold-sliding model as the ¬rst step in a
glacial-landform evolution model for erosion ice 3, the basal velocities are used to determine the
sheets. Following the ice-surface reconstruction, bedrock-erosion rate using Hallet™s model. Hal-
basal-sliding velocities and bedrock erosion rates let™s model estimates erosion rates as a power

(a) (c)


1 km
is area of study
> 100 km?

30 km


0.5 bars

0.5 bars

Fig 6.23 A ¬‚ow chart illustrating the model inputs and outputs with erosion of the modern New York State Finger Lakes as the
example application. (a) Inputs to the model include the initial bedrock topography (from US Geological Survey DEMs), the
ice-margin position (given by the position of the Valley Heads Moraine (New York State Geological Survey, 1999)), and the
threshold basal shear stress (assumed to be 0.5 bars). (b) In step 1 these inputs are used to reconstruct the ice-surface
topography. Following this step, the de¬‚ection of the lithosphere under the ice load should be determined if the study area is
equal to or larger than 100 km in extent. For this example we have neglected ¬‚exural subsidence. (c) In step 2 the ice-surface
topography and pattern of accumulation and ablation are used with the balance-velocity method of Budd and Warner (1996) to
determine the basal-sliding velocities. (d) In step 3, Hallet™s model, which relates the bedrock erosion rate to a power function of
the sliding velocity, is used to determine the spatial pattern of erosion. The model then returns to step 1 where the erosion
pattern is used to determine a new bed topography for the next iteration.

function of the basal-sliding velocity. In step 4, Any study of glacial erosion must contend
the pattern of erosion rate is used, together with with two complicating factors. First, the initial,
the model time step, to modify the bedrock to- preglacial surface is generally unknown. This
pography according to b = E t, where b is the complication applies to individual glacial ad-
bedrock topography, E is the erosion rate, and t vances as well as to erosion during the full dura-
is the time step. The new bedrock topography, re- tion of the Quaternary. As a speci¬c example, it
¬‚ecting the bedrock scour during the time step is commonly assumed that the ice ¬‚ow respon-
t, is then used as input for the next model it- sible for the carving of the New York State Fin-
eration, which begins back at step 1. ger Lakes took advantage of preexisting ¬‚uvial

valleys of the Allegheny Plateau formed by north- the abrading particles, embedded in the ice, is
¬‚owing rivers (von Engeln, 1965). It is dif¬cult the most important controlling factor for ero-
to test this hypothesis, however, since we have sion. Hallet™s model can be summarized as
no constraints on the depth of pre-glacial ¬‚u-
= ’c V n (6.50)
vial valleys. If valleys were deep, it may still be
possible that ice ¬‚ow and bedrock erosion were
where ‚b/‚t is the erosion rate, V is the basal
controlled principally by the ice-surface topogra-
sliding velocity, and n and c are empirical co-
phy of the Ontario Lobe rather than by the un-
ef¬cients. n may taken to be equal to 2 (e.g.
derlying bed topography. Numerical modeling of
MacGregor et al., 2000), although this value is
glacial erosion is ideal for coping with this un-
not well-constrained empirically. The value of
certainty because multiple working hypotheses
c may be constrained using sediment-¬‚ux esti-
can be tested and evaluated for their impact on
mates from modern glaciers (e.g. Hallet et al.,
model results. For example, a suite of initial to-
1996). In this chapter we use Hallet™s model, re-
pographies can be used to determine a suite of
¬‚ecting the consensus that sliding velocity is the
corresponding model results to assess the sensi-
dominant control on bedrock erosion rates. Nu-
tivity of the model results to the initial topogra-
merical modeling in glacial geomorphology is
phy. Results that are insensitive to variations in
still in its infancy. Three of the most important
initial topography can be considered to be robust.
contributions have been made by Harbor (1992),
A second complicating factor in studying ero-
MacGregor et al. (2000), and Braun and his col-
sion beneath ice sheets is the dif¬culty of separat-
leagues (Braun et al., 1999; Tomkin and Braun,
ing variations in the intensity of glacial erosion
2002). All of these contributions focused on the
from variations in the erodibility of bedrock. In
alpine glacial environment, but they represent
some cases, structural control is easy to recog-
the previous work most similar to the work in
nize. For example, the topography of the Finger
this chapter. Harbor™s model illustrated how U-
Lakes Region is dominated by a weathering scarp
shaped glacial valleys may form from pre-existing
separating the ¬ne-grained shale of the Ontario
¬‚uvial valleys. However, because he assumed a
Lowlands from the coarser-grained sandstone and
simple form for the glacier (i.e. a constant-area
shale of the Allegheny Plateau. In other areas,
cross-sectional ¬ll with a level surface), his model
however, minor differences in deformation may
cannot readily be generalized to more complex
have in¬‚uenced erosion more greatly than varia-
2D cases in which a level surface cannot be as-
tions in ice thickness or basal sliding. This is not
sumed. MacGregor et al. (2000) developed a 1D
a great concern for distinctive landforms such as
model to study the formation of overdeepenings
¬nger lakes. However, if bedrock geologic maps
using a model based on mass-balance, ice-creep,
and structural data are available for a region,
and stress-controlled basal sliding. They showed
structural control should be evaluated through
that several classic alpine-glacial landforms, in-
correlations between topography and structure.
cluding overdeepenings and hanging tributaries,
Determining the controls on glacial bedrock-
could be reproduced by their model. Their results
erosion rates is an active area of research. Boul-
were limited to 1D, however.
ton (1974) emphasized the role of ice thickness
Braun and colleagues (Braun et al., 1999;
in determining erosion rates, while Hallet (1979,
Tomkin and Braun, 2002) constructed the ¬rst
1996) emphasized the importance of sliding ve-
2D glacial-landform evolution model. These au-
locity. In Boulton™s model, the hydrostatic pres-
thors used a continuity equation to model local
sure induced by the ice overburden pressure in-
changes in ice thickness:
creases the normal stress acting on subglacial de-
bris, in turn increasing their ability to abrade the
=∇·F +M
bed. In Hallet™s model, ice thickness is not a fac- (6.51)
tor because hydrostatic pressures are assumed to
where the ¬‚ux F was dependent on creep and
act equally on the tops and bottoms of rounded
sliding velocity and M is the mass balance or rate
subglacial particles. In this case, the velocity of

of accumulation and ablation. Sliding velocities the basal shear stress can also be varied to simu-
in their model were given by late ˜˜sticky™™ spots on the bed. Both assumptions
yield nearly identical results.
2A s βc (ρgh)n
us = |∇(h + b)|n’1 ∇(h + b) (6.52) The ice-surface topography and basal-sliding
velocities above a ¬‚at surface with white noise
where h is the ice thickness, b is the bedrock variations (shown in Figure 6.24c) are illustrated
elevation, and A s , βc , N , n, and P are empiri- in Figures 6.24a and 6.24b for a basal shear stress
of 1 bar. The grid for this example is 50 km —
cal parameters. Importantly, basal sliding in this
model is dependent only on the local ice depth, 50 km in size with a pixel resolution of 100 m.
ice slope, and bed slope. In real glaciers, however, To construct this surface, we allowed all pixels
sliding velocity is controlled by the discharge of in the grid to accumulate ice to simulate an ice
ice, increasing with distance from the divide and sheet with no margin. The slopes in the y direc-
reaching a maximum beneath the ELA, even if tion, S y , were restricted to point down the grid,
the glacier and bed morphology remains con- rather than in either the up or down direction,
stant along the glacier pro¬le. to simulate an ice-sheet ¬‚owing from top to bot-
Here we illustrate the behavior of our glacial- tom. The slope in the x direction, S x , was uncon-
landform evolution model with several numer- strained. The sandpile algorithm for this example
ical experiments. We distinguish three types of was initiated with a uniform ice-surface topog-
depressions: (1) lakes in ice sheet interiors that raphy of 2 km. From this initial state, the sand-
are ¬‚exurally compensated, (2) large lakes that pile method was used to solve for the threshold-
are uncompensated (˜˜great™™ lakes), and (3) elon- sliding ice-surface topography. The resulting ice-
gate ice-marginal depressions (˜˜¬nger™™ lakes). In surface topography dips gently from the top of
each case, topographic analyses are introduced the grid to the bottom (Figure 6.24a). In addition,
for comparison with the model predictions. the ice surface has small-scale variations that re-
¬‚ect the white-noise variations in bed topogra-
6.6.1 Length scales < ¬‚exural wavelength phy. Equation (6.46) requires that the ice-surface
First we consider bedrock erosion beneath a slope and ice thickness be inversely proportional.
gently-sloping ice sheet at scales less than the In the model reconstruction, therefore, pixels
¬‚exural wavelength. In this example we assume with slightly lower bed elevations have lower ice-
that the lithosphere is perfectly rigid. The ini- surface slopes. In this way, variations in the bed
tial bed topography is ¬‚at with random micro- topography are re¬‚ected in the ice-surface topog-
topography characterized by a Gaussian white raphy.
noise with a standard deviation of 10 m. This The pattern of basal-sliding velocities in Fig-
initial topography is certainly an idealization; ure 6.24b was produced assuming uniform ac-
no initial condition would look like this. How- cumulation within the model domain. Although
ever, all topography has some small-scale rough- the ice-surface topography has only minor small-
ness, and it is essential to include these varia- scale variations, the basal sliding is strongly lo-
tions in the initial bed topography because the calized into bedrock depressions in a braided pat-
glacial-erosion instability ampli¬es initial topog- tern. This behavior results from an additive effect
raphy over time. White noise variations in ini- of the ice-surface topography on the basal sliding.
tial conditions are commonly used to initiate The input of ice to the system as accumulation is
¬‚uvial-geomorphic instabilities such as drainage- uniform throughout the domain. At each pixel,
network evolution (e.g. Willgoose et al., 1991) and however, slightly more ¬‚ow is diverted to steeper
hillslope rilling (Hairsine and Rose, 1992), for ex- downslope pixels than to pixels with a gentler
ample. A key question in this section will be slope. As ice ¬‚ows from pixel to pixel downice,
whether the glacial-erosion instability ampli¬es small variations in surface slope act cumulatively
initial topography at all wavelengths or whether to divert ice into bedrock depressions. The re-
certain wavelengths are preferred. As an alterna- sulting areas of intense sliding velocities result
tive to assuming variations in initial topography, in enhanced scouring of the bed. This, in turn,

2020 m iteration 1 iteration 1 iteration 1
(b) (c)


2000 m
ice-surface topography basal sliding velocity bed topography
(d) (e) (f)

iteration 10 iteration 10 iteration 10

Fig 6.24 Erosion and basal sliding beneath ice-sheet they depress the ice-surface topography over a
interiors. (a) Shaded-relief image of ice-surface topography,
larger area, focusing more ¬‚ow into the depres-
(b) grayscale image of basal-sliding velocities, and (c) shaded
sion. The dynamics of the model, therefore, can
relief image of bed topography in the ¬rst iteration of the
be characterized as depressions of different size
model. (d) Ice-surface topography, (e) basal-sliding velocities,
competing for ice drainage. The larger a depres-
and (f) bed topography after ten model iterations. Bedrock
sion is, the faster it expands to ˜˜capture™™ neigh-
depressions deepen and expand over time, focusing basal
boring lakes. As a result, Figure 6.24f contains
sliding and erosion in a positive feedback.
lakes or pits of many sizes, while Figure 6.24c
contains only small pits. The braided patterns of
Figures 6.24b and 6.24e meander and shift as the
results in localized bedrock erosion in depres-
bed is eroded and depressions expand.
sions. This is the fundamental glacial-erosion in-
In Figure 6.25, the lakes of Figure 6.24f are
stability. In the case of a gently-sloping surface
compared with the lakes of glaciated topography
of an ice sheet with no margin, the instability
southwest of Hudson Bay (Figures 6.25d--6.25e). In
results in lakes with little or no elongation and
Figure 6.25a, a shaded-relief image of the model
a wide range of sizes. Near the ice margin, in con-
topography from Figure 6.24f is given, with all
trast, lakes are elongated and regularly spaced.
closed depressions ¬lled. The binary image of Fig-
The ice-surface topography, basal sliding ve-
ure 6.25b includes all of the ¬lled areas, or lakes,
locities, and bed topography of this model af-
of Figure 6.25a.
ter 10 iterations are given in Figures 6.24d--6.24f.
The frequency-size distribution was intro-
These ¬gures illustrate that ¬‚ow becomes more
duced in Chapter 1 as a tool for characteriz-
focused into bedrock depressions as the relief of
ing the population of lakes in a given area. The
the bed topography increases through time. Ini-
frequency-size distribution is the number of lakes
tially, all depressions are limited to a few pixels in


. 8
( 14)