. 7
( 14)


for narrower loads. The transient response of the We can use the symmetry of the problem to spec-
lithosphere to loading and unloading can gener- ify the boundary conditions. The value of w must
ally be neglected over time scales of interest in go to zero as x goes to ∞ and dw/dx must go to
mountain building. zero at x = 0. This implies that c 1 and c 2 are zero
The equation for the 1D ¬‚exural-isostatic dis- and c 3 = c 4 to give
placement of a uniformly rigid lithosphere to ver- x x
w(x) = c 3 e’x/± cos + sin (5.10)
tical loading and unloading is: ± ±
d4w The value of c 3 is proportional to the applied load
+ (ρm ’ ρc )gw = q(x)
D (5.6)
dx 4 V . From Eq. (5.3) and Eq. (5.10) we have
This equation combines ¬‚exural bending stresses d3 w
1 4D c 3
V=D = (5.11)
(¬rst term on left side), the upward buoyancy ±3
dx 3
2 x=0
force exerted on the bottom of the lithosphere
We used 1 V for the load in Eq. (5.11) because we
(second term on left side), and the applied load 2
are solving for only half of the pro¬le (i.e. the
q(x) (right side).
other half of the pro¬le supports the other half
of the load). Substituting Eq. (5.11) into Eq. (5.10)
5.2 Methods for 1D problems gives
V ± 3 ’x/± x x
w(x) = cos + sin
e (5.12)
± ±
5.2.1 Displacement under line loading 8D
for x ≥ 0 (Figure 5.2b).
In some cases the topographic load is nar-
row enough compared to the ¬‚exural wave- In the case of de¬‚ection of a diving board,
length of the lithosphere that the load can be we found that the solution is characterized by

thickness of 40 km, for example, the ¬‚exural pa-
rameter under crustal loading is ≈ 130 km, while
for ice loading it is approximately 90 km.
(km) 5.2.2 Series solutions and Fourier ¬ltering
The ¬‚exure equation is a linear equation. As such,
ice we can use Fourier series and superposition to de-
velop general solutions for the de¬‚ection of the
50 lithosphere under any topographic load, given
the solution for de¬‚ection under a periodic load.
Let™s assume a periodic topography with ampli-
20 40 80 tude h 0 and wavelength »:
T e (km)
2π x
h(x) = h 0 sin (5.16)
Fig 5.3 Plot of ¬‚exural parameter ± as a function of elastic »
thickness Te for crustal and ice loads, assuming E = 70 GPa
The vertical load corresponding to this topogra-
and ν = 0.25.
phy is given by
a cubic polynomial with no characteristic length 2π x
qa (x) = ρc gh 0 sin (5.17)
scale. In the case of ¬‚exural-isostasy, however, the
de¬‚ection is characterized by the length scale ±.
and hence the ¬‚exural-isostatic equation we must
Also, the zone of depression is accompanied by
solve is
a smaller zone of uplift ¬‚anking the margin of
d4 w
the load. This uplift zone is the forebulge. The + (ρm ’ ρc )gw = qa (x)
D (5.18)
dx 4
location of the forebulge is determined by dif-
ferentiating w with respect to x and setting the The displacement of the lithosphere will follow
result equal to zero: the same periodic form as the loading -- the mag-
nitude of the de¬‚ection will depend sensitively on
dw 2w 0 ’x/± x
=’ sin = 0
e (5.13) the wavelength of the load and ¬‚exural param-
± ±
eter of the lithosphere, but there will be no off-
set or shift between the load and the response.
xb = ± sin’1 (0) = π± (5.14) Therefore, we look for solutions of the form
The amplitude of the forebulge can be found by 2π x
w = w 0 sin (5.19)
substituting xb into Eq. (5.12) to obtain
w b = ’w 0 e’π = ’0.0432w 0 (5.15) where w 0 is yet to be determined. Substituting
Eq. (5.19) into Eq. (5.18) gives
Therefore, the height of the forebulge is equal
to about 4% of the maximum de¬‚ection of the
w0 = (5.20)
ρm 2π 4
lithosphere beneath the load. This calculation as- D
ρc ρc g »
sumes that the lithosphere is entirely rigid. In
For small wavelengths, i.e. where
some cases, the loading causes the lithosphere
to break under the stresses of the load. In these 1
D 4
cases, the sin term in Eq. (5.12) becomes zero and 2π (5.21)
ρc g
the width of the zone of subsidence becomes nar-
rower (Turcotte and Schubert, 2002). the denominator in Eq. (5.20) is very large
The relationship between ± and T e is plotted and hence w 0 h 0 . This means that small-
in Figure 5.3 for loading by crustal rocks (i.e. wavelength topographic variations cause very lit-
ρ ≈ 600 kg/m3 ) and ice (i.e. ρ ≈ 2300 kg/m3 ), as- tle displacement of the crust and we say that the
suming E = 70 GPa and ν = 0.25. For an elastic crust is rigid for loads of this scale. At large scales,

where (so called because it ¬lters out low wavenumbers)
before the inverse Fourier transform is applied to
D 4
bring the function back to real space.
» 2π (5.22)
ρc g Appendix 4 provides code that implements
the Fourier ¬ltering technique for any input to-
Eq. (5.20) gives
pographic load. Figure 5.4 illustrates example
w0 = h0 output of that program for an east--west pro¬le
ρm ’ ρ c
through the central Andes. The east--west topo-
graphic pro¬le of the central Andes is character-
which is equivalent to the case of pure isostatic
ized by a high-elevation plateau (the Altiplano-
balance given by Eq. (5.2). In the large-scale limit,
Puna basin), high peaks on the margins of the
therefore, ¬‚exure plays no role and the topo-
Altiplano-Puna, and an active fold thrust belt to
graphic load is said to be fully compensated. At
the east of the Altiplano-Puna. In this example,
smaller scales, where the de¬‚ection is smaller
the topographic load of the Andes is considered
than the Airy isostatic limit, the topography is
to be all topography above 1 km in elevation. Fig-
said to be only partially compensated.
ure 5.4b plots the modeled de¬‚ection of the crust
The solution for periodic topography is use-
for ρc = 2.7 g/cm3 , ρm = 3.3 g/cm3 , and three hy-
ful for determining the degree of compensation
pothetical values of ±. The depth of total de¬‚ec-
of a topographic load. The real power of the pe-
tion beneath the Andes (≈ 15 km) depends lin-
riodic solution, however, is that it can be used
early on the relative density contrast between the
in conjunction with the Fourier series approach
crust and mantle, but only weakly on the ¬‚ex-
to determine the response of the lithosphere to
ural parameter for wide topographic loads. Low
any load using superposition. If the topographic
values of the ¬‚exural parameter predict narrow
load within a domain of length L is expressed as
foreland sedimentary basins, while higher val-
a Fourier series:
ues predict broad basins (Figure 5.4b--5.4c). The

h(x) = height of the ¬‚exural forebulge (Figure 5.4c) also
an sin (5.24)
L increases with higher values of the ¬‚exural pa-

rameter. Gravity modeling has been used to con-
with coef¬cients
strain the width of the Andean foreland basin by
2 nπ x taking advantage of the density contrast between
an = h(x ) sin dx (5.25)
L L the basin in¬ll and the crust below. Those con-

straints, coupled with forward modeling of ¬‚ex-
then, according to Eq. (5.20), the de¬‚ection
ural pro¬les, suggest that ± = 100 km is a good
caused by the load is given by
estimate for the ¬‚exural parameter of the central
∞ nπ x
w(x) = L
an (5.26)
ρm nπ 4
’1+ D
n=1 ρc ρc g L

Equation (5.26) can be used to obtain analytic ex-
5.3 Methods for 2D problems
pressions for the ¬‚exural-isostatic displacement
underneath loads with simple (rectangular, tri-
5.3.1 Integral method
angular, etc.) loads.
Most real topography, however, does not lend Thus far we have considered solutions in 1D only.
itself to simple analytic forms. Equation (5.26) These solutions are very useful for long moun-
can be implemented numerically for any topo- tain chains or other landforms where the loading
graphic load, however, using the Fourier ¬lter- is relatively constant in one direction. In many
ing technique. In this technique, the Fast Fourier cases, however, (e.g. an isolated seamount or curv-
Transform (FFT) is used to quickly calculate the ing mountain belt) a full 2D solution is needed.
Fourier coef¬cients an for an input load function. In such cases, 2D methods analogous to those we
Then, each coef¬cient is multiplied by the ¬lter developed for 1D are available.

(a) 5


(c) 0.4
1 (km)
(b) 0
± = 50 km
75 km
100 km
(km) ± = 50 km
800 1000 1400
75 km
100 km x (km)

600 800 1000 1200
x (km)
Solutions corresponding to a spatially dis-
Fig 5.4 Flexural model of the response of the lithosphere
to the distributed load of the central Andes. (a) Topographic tributed load can be obtained by representing
pro¬le of the central Andes. (b) De¬‚ection of the lithosphere the load as a collection of many point sources
beneath the central Andes for three hypothetical values of
and integrating the response to each point. For
± = 50, 75, 100 km. (c) Close-up plot of the ¬‚exural pro¬le
a spatially distributed topographic load given by
of the foreland basin and forebulge, illustrating the
h(x, y), the de¬‚ection is given by:
dependence of basin width and forebulge amplitude on ±.

∞ ∞
First we consider the solution to the 2D w(r ) =
ρm ’ ρc
¬‚exural-isostatic equation: 0 0

— kei (x ’ x )2 + (y ’ y )2
D ∇ 4 w + ρgw = q ±
— h(x , y )dx dy (5.29)
corresponding to a point load V at r = 0. The
solution is the Kelvin function kei:

V ±2 This approach is broadly analogous to the series
w(r ) = kei (5.28)
πD ± solution method of Eq. (5.26).
Few analytic solutions can be obtained using
The shape of the kei function is broadly similar Eq. (5.29), both because topographic loads rarely
to Eq. (5.12), i.e. it has a deep zone of subsidence conform to analytic expressions and because the
beneath the load and a smaller forebulge ¬‚ank- kei function is dif¬cult to integrate. Appendix 4
ing the load. provides code for series approximations to the kei

ral parameter, the de¬‚ection follows the disk-like
shape of the load except at the edges where ¬‚ex-
0 0.2
ural bending stresses limit the curvature of the
w/h de¬‚ection pro¬le.
1 R/a =
5.3.2 Fourier ¬ltering
2 In the previous section we showed how the de¬‚ec-
tion beneath spatially distributed loads could be
calculated by dividing the load into many small
point sources and summing the contributions of
each point. This is a powerful approach, but it
can be computationally challenging when calcu-
’1 0 1
lating the response at high spatial resolution over
a large region. Such calculations can require eval-
uating the kei function 104 --106 times for each of
Fig 5.5 De¬‚ection of disk-shaped loads for various values
104 --106 points. The Fourier-¬ltering method pro-
of R /±.
vides a faster way by taking advantage of the
computational speed of the Fast Fourier Trans-
function that can be used to integrate Eq. (5.29)
form algorithm.
numerically for any topographic load.
In order to implement the Fourier ¬ltering al-
One analytic expression is available, how-
gorithm in 2D, the FFT of the 2D load h(x, y) over
ever, that nicely illustrates the response of the
a domain of width L x and length L y is computed
lithosphere to loads of different size. The re-
to obtain a set of Fourier coef¬cients an,m . Each
sponse of the lithosphere to a disk-shaped load
Fourier coef¬cient is then multiplied by the fac-
depends on the ratio of the disk radius R to the
¬‚exural parameter ±. The solution is given by
⎛ ⎛ ⎞4 ⎞’1
(Hetenyi, 1979): 2 2
⎝ ρm ’ 1 + D ⎝π n m ⎠⎠
hρc R R r
ρc ρc g
w(r ) = Lx Ly
ker ber
ρm ’ ρ c ± ± ±
R R r
’ kei +1
bei (5.30)
± ± ± Then, the inverse Fourier transform is applied to
the ¬ltered data set.
for r < R and
Appendix 4 provides a code for implementing
hρc R R r
w(r ) = the Fourier ¬ltering method in 2D. Figure 5.6 il-
ber ker
ρm ’ ρ c ± ± ±
lustrates example output from this program for
R R r
the central Andes for ± = 50, 75, and 100 km. Fig-
’ bei kei (5.31)
± ± ±
ures 5.6b--5.6d are grayscale plots of the absolute
for r > R . value of the de¬‚ection on a logarithmic scale. As
the value of ± increases, the pattern of de¬‚ection
Figure 5.5 plots Eq. (5.31) for loads of differ-
ent radii, expressed as a ratio of the ¬‚exural pa- becomes smoother. In each grayscale image the
rameter ±. For cases in which the load radius forebulge is represented by the ¬rst ring around
is much less than the ¬‚exural parameter, the the central zone of subsidence. Subsequent rings
load is only partially compensated. For cases in are the backbulge and other smaller zones of al-
which the radius is equal to the ¬‚exural parame- ternating uplift and subsidence with increasing
ter, the maximum de¬‚ection is equal to the Airy distance from the load. Figure 5.6e illustrates the
isostatic limit (i.e. Eq. (5.2)), but the maximum de- variations in forebulge amplitude in the vicinity
¬‚ection is attained only beneath the center of the of the Bolivian Orocline using a contour map. The
load. For load radii much larger than the ¬‚exu- forebulge amplitude varies from 160 m to 400 m

basin environments are also attractive environ-
Fig 5.6 2D ¬‚exural-isostatic response of the lithosphere to
ments for studying the linkages of climate, tec-
the topographic load of the central Andes (a) for ± = 50 (b),
tonics, and erosion, because the tectonic param-
75 (c), and 100 km (d). In (b)“(d), the grayscale images
eters controlling foreland basin geometry can be
represent the absolute value of the de¬‚ection on a
unusually well constrained. Speci¬cally, thrust-
logarithmic scale. This color scale clearly shows that the
width of the forebulge increases with increasing ±. (e) A belt migration rates can be determined using
contour map of the forebulge in the vicinity of the Bolivian a mass-balance framework that relates crustal
Orocline. This example shows that the height of the shortening to thrust-belt migration (DeCelles and
forebulge changes dramatically (from a minimum of 160 m to
DeCelles, 2001) and basin subsidence is often
a maximum of 400 m) in the vicinity of a curved load.
well characterized using ¬‚exural modeling, grav-
ity measurements, and seismic stratigraphy (e.g.
within a distance of ≈ 100 km, illustrating the
Watts et al., 1995). Sediment transport in ¬‚uvial
importance of 2D effects in the vicinity of a bend
channels has a linear relationship with channel
in the topographic load.
slope, resulting in a diffusion model for the basin
topographic pro¬le (Begin et al., 1981; Paola et al.,
5.4 Modeling of foreland basin
Foreland basin evolution can be modeled as
geometry a combination of thrust-belt migration, ¬‚exural
subsidence, and diffusive ¬‚uvial sediment trans-
port (Flemings and Jordan, 1989). Modeling this
Flexure plays an important role in large-scale
system requires solving the diffusion equation
landform evolution because, as erosion and de-
with a moving boundary (the thrust front, where
position change the topographic load, ¬‚exural
sediment enters the basin) and moving sink
isostasy causes uplift and subsidence responses
term (the ¬‚exural depression). Typically, mov-
that vary as a function of spatial scale and litho-
ing boundaries and sinks in diffusion problems
spheric rigidity. In this section we explore the
require numerical approaches. However, in this
interaction between ¬‚exural isostasy and sedi-
case the problem can be solved analytically by
mentation in a foreland basin adjacent to a mi-
working in the frame of reference of the migrat-
grating thrust belt. Foreland basins provide a po-
ing thrust front and assuming steady-state con-
tentially useful record of tectonics and erosion
ditions. Conceptually, basins are in steady state
of the adjacent range, but our ability to inter-
when sediment progradation occurs at the same
pret that record uniquely is limited. Foreland

rate as thrust-belt migration. Basins naturally where w(x) is the ¬‚exural pro¬le, given by Tur-
evolve toward this condition because when sed- cotte and Schubert (2002) as
iment progradation lags behind thrust-belt mi- x
w(x) = ’w 0 e’x/± cos (5.36)
gration at any point along its pro¬le, the basin
steepens locally until the two are in balance. where w 0 is the basin depth beneath the thrust
Conversely, if sediment progradation outstrips front and ± is the ¬‚exural parameter. Equa-
thrust-belt migration, slope decreases and progra- tion (5.36) corresponds to the ¬‚exural pro¬le cor-
dation slows. This approach is analogous to that responding to a line load. This is only an approx-
of the Kenyon and Turcotte model of deltaic sed- imation to the Andes, where loading is spatially
imentation considered in Chapter 2. In this case, distributed.
however, we will also consider the pattern of Combining Eqs. (5.33) and (5.35) and as-
¬‚exural subsidence that controls accommodation suming steady-state conditions gives the time-
space in the basin. independent equation:
Basin evolution in the model is governed by
d2 h dh dw
κ +F =F
the diffusion equation: (5.37)
dx dx dx
‚h ‚ 2h which can be written as a ¬rst-order differential
=κ 2 (5.33)
‚t ‚x equation for the slope S = dh/dx with solution
where h is the elevation, t is time, κ is the diffu-
’e’F x/κ x
e’F x /κ
S(x) = Qs ’ F dx
sivity, and x is the distance from the thrust front.
κ dx
Boundary conditions are needed at both ends of
the basin to solve Eq. (5.33). At the thrust front,
Substituting Eq. (5.36) into Eq. (5.38) gives
a constant sediment ¬‚ux Q s (in m2 /yr) enters the

’e’F x/κ ±F ±F
F x x x
e’(F /κ+1/±)x
S(x) = Qs ’ ’ sin ’ 2sin ’
cos (5.39)
κ (± F /κ)2 + 2± F /κ + 2 κ ± ± ± κ

The basin topographic pro¬le is obtained by inte-
basin, providing a ¬‚ux boundary condition
grating Eq. (5.39) analytically or numerically and
‚h Qs enforcing the basin-outlet boundary condition to
= (5.34)
‚x κ constrain the integration constant.

Model solutions are plotted in Figure 5.7, rep-
In this 2D framework, the sediment ¬‚ux is given
resenting the effects of varying Q s on basin pro-
by ELd , where E is the basin-averaged erosion rate
¬les around a reference case with parameter val-
and L d is the upstream drainage basin length.
ues F = 10 mm/yr, ± = 150 km, w 0 = 5 km, κ =
At the basin outlet, a constant base-level eleva-
5000 m2 /yr, Q s = 10 m2 /yr, and assuming an in-
tion boundary condition is prescribed. Two cases
¬nite basin. Figure 5.7a illustrates three impor-
can be considered. For an in¬nite basin, h = 0 as
tant points. First, the basin pro¬le is highly sen-
x = ∞. For a basin of length L b and base-level ele-
sitive to incoming sediment supply: the topo-
vation h b (determined by sea level or a valley-¬‚oor
graphic relief of the basin doubles when Q s is
channel), the boundary condition is h(L b ) = h b . In
increased by 50% from 8 to 12 m2 /yr. This suggests
the moving reference frame, the foreland basin
that basin pro¬les may provide useful constraints
moves towards the thrust front and enters the
on upstream sediment supply through forward
¬‚exural depression with velocity F , the thrust-
modeling and comparison with observed pro¬les.
migration velocity. This motion is represented by
Second, as sediment supply decreases, the basin
an advection equation:
undergoes a transition from over¬lled to under-
‚h ‚w ‚h ¬lled (closed) basins. As such, the model provides
=F ’ (5.35)
quantitative criteria for topographic closure in
‚t ‚x ‚x

reference case:
e’Fx/k F = 10 mm/yr
a = 150 km
Qs = 12 m2/yr
0.6 w0 = 5 km 10 000 m2/yr
Qs = 10 m2/yr k = 5000 m2/yr
5000 m2/yr
Qs = 8 m2/yr Qs = 10 m2/yr
h (km)

1500 —
0.4 2500 m2/yr
vert. exag.
(closed basin)

0.2 forebulge k=
1000 m2/yr
sensitivity to Qs sensitivity to k
0.0 infinite basin
infinite basin
(a) no forebulge erosion
no forebulge erosion

0.8 Ld
sensitivity to Qs
finite basin
(c) (d)
no forebulge erosion
0.6 Qs = 12 m2/yr
Qs = 10 m2/yr flux boundary
h (km)

Qs = 8 m2/yr
dh/dx = Qs /k sediment flux s
F k
0.4 diffusivity
basin-fill profile
0.2 slope
flexural profile
800 1000
400 600
x (km)
cause sediment back¬lling at the thrust front to
Fig 5.7 (a) and (b): Model solutions for the topographic
produce steep, short basins. Surprisingly, higher
pro¬le of an in¬nite basin assuming no forebulge erosion,
diffusivities (e.g. more humid conditions) pro-
showing sensitivity to variations in Q s and κ, respectively,
around a reference solution with F = 10 mm/yr, ± = 150 km, mote basin closure if sediment supply is held
w 0 = 5 km, κ = 5000 m2 /yr, and Q s = 10 m2 /yr. (c) Solutions constant. Naively, we might expect more humid
for a ¬nite basin of length L b = 1000 km for the same model conditions to result in steeper, more over¬lled
parameters as in (a). (d) Schematic diagram of the model.
basins. However, to understand the in¬‚uence of
Modi¬ed from Pelletier (2007b).
climate on basin geometry we must consider
the effects of climate on sediment supply (i.e.
foreland basins. Third, the distal basin pro¬le is weathering) and basin transport rates indepen-
approximated by a simple exponential function dently. Greater sediment supply puts more total
with length scale κ/F . This result suggests that sediment into the basin and therefore promotes
distal basin pro¬les can be used to uniquely and steeper and more over¬lled basins. Diffusivity val-
easily constrain κ values if thrust-belt migration ues, on the other hand, control how uniformly
rates are known. that sediment is spread across the basin. Higher
Figure 5.7b illustrates the sensitivity of model diffusivity values promote closed basins because
solutions to variations in κ around the reference a larger fraction of the total sediment is trans-
value of κ = 5000 m2 /yr. Diffusivity values are ported past the foredeep, leaving behind a depres-
sion if the ratio Q s /κ is suf¬ciently small. These
controlled by upstream basin length, precipita-
tion, and sediment texture (with longer basins, results emphasize that the key to understanding
higher precipitation rates, and ¬ner textures pro- basin geometry is not simply the sediment sup-
moting greater diffusivity values) (Paola et al., ply to the basin, but the ratio of sediment supply
1992). Low diffusivity values (e.g. κ = 1000 m2 /yr) to transport capacity (e.g. as emphasized in the

geomorphic context by Bull (1979)). Figure 5.7c on crustal shortening rates. Flexural modeling
illustrates model solutions for a ¬nite basin us- coupled with geophysical and geomorphic obser-
ing the same model parameters as in Figure 5.7a. vations constrain the ¬‚exural parameter to be
Shorter basins have lower basin relief, all else be- approximately 150 km (resulting in a foredeep ap-
ing equal. proximately 250 km wide and a forebulge at a
Equation (5.40) assumes no forebulge erosion distance 400 km from the thrust front) (Coud-
(i.e. that the forebulge is buried or, if exposed, ert et al., 1995; Watts et al., 1995; Horton and
that the erosion rate is low). In some cases, it may DeCelles, 1997; Ussami et al., 1999). These data
be more appropriate to assume that the forebulge provide relatively ¬rm constraints on three of
the ¬ve model parameters: F = 10 mm/yr, ± =
erodes completely to sea level. In this alternative
150 km, w 0 = 4 km.
end-member scenario, Eq. (5.36) is amended to
w(x) = 0 for x > 3π±/4 (i.e. the point at which Figure 5.9a illustrates a map of the topog-
the forebulge rises above sea level). In this case, raphy of the central Andes and the adjacent
foreland basin from 17—¦ S to 31—¦ S. This region
the slope is given by Eq. (5.36) for x < 3π±/4,
and by

’e’F x/κ ±F ±F
F 1 3π ±
1.414e’( κ + ± ) 4
S(x) = Qs ’ 1+ + (5.40)
κ /κ)2 + 2± F /κ + 2 κ κ
(± F

for x > 3π ±/4. Model solutions assuming fore- corresponds to a humid-to-arid transition in the
Andes based upon mean winter precipitation val-
bulge erosion are given in Figure 5.8.
ues above 1500 m elevation. Associated with this
This model framework also provides a simple
climatic transition is a geomorphic transition
means to estimate the basin sediment delivery
from strongly over¬lled (700 m in relief) to weakly
ratio (the ratio of the sediment yield to erosion
over¬lled (100 m in relief) basins. Closed basins
rate). The sediment ¬‚ux at the entrance to the
(e.g. L. Ambargasta and L. Mar Chiquita) also oc-
foreland basin is Q s . For a basin of length L b ,
cur near the southernmost transect, but they are
the ¬‚ux leaving the basin is computed by evaluat-
ing Eq. (5.40) at x = L b and multiplying by κ. For located in wedgetop depozones of a broken fore-
land and hence are not well described by the
basins that have a base-level control located be-
yond the forebulge distance, Eq. (5.40) reduces to
Three topographic pro¬les were extracted
a simple exponential function with length scale
κ/F . In the reference case, for example, a 500-km from the Shuttle Radar Topographic Mission
(SRTM) data and are plotted in Figure 5.9b. To
wide basin has a sediment delivery ratio of 1/e
= 0.37, and a 1000-km wide basin has a value model these pro¬les, it is most accurate to use
of 1/e2 = 0.14. This latter value is broadly consis- ¬nite-basin solutions constrained by values of h b
and L b for each pro¬le. The base level eleva-
tent with measured sediment delivery ratios for
tion h b of the Andean foreland is controlled by
large drainage basins in areas of active tectonism
Rio Paraguay, which ranges in elevation from 80
(Schumm, 1977).
to 10 m (decreasing from north to south) at dis-
As a test of the model, we consider the mod-
tances L b varying from 350 to 720 km from the
ern foreland basin of the central Andes. The
thrust front. Model results in Figure 5.7a sug-
latest phase of Andean deformation responsible
gest that the distal basin pro¬le is an exponential
for Eastern Cordilleran and Subandean uplift is
function with length scale κ/F . The inset plot in
generally considered to be late Miocene in age
Figure 5.9b showing hh b (note vertical logarith-
(Gubbels et al., 1993). The foreland basin is ap-
mic scale) as a function of x con¬rms this pre-
proximately 3.5--4.5 km thick beneath the thrust
diction for the central Andean foreland. Best-¬t
front based on seismic re¬‚ection data (Horton
exponential pro¬les (straight lines on this loga-
and DeCelles, 1997). DeCelles and DeCelles (2001)
rithmic scale) provide estimates of κ = 1500 m2 /yr
estimated the thrust-migration rate to be approx-
for the northern pro¬le and κ = 4000 m2 /yr for
imately 10 mm/yr for the central Andes based

result of assuming uniform κ values. Down-
(a) 0.8 reference case:
stream ¬ning occurs in all basins, and therefore
F = 10 mm/yr
± = 150 km
κ values are likely to increase somewhat down-
0.6 w0 = 5 km
h stream, resulting in gentler distal-basin slopes
Qs = 12 m2/yr κ = 5000 m2/yr
Qs = 10 m2/yr
Qs = 10 m2/yr
(km) compared to the model solutions of Figure 5.9b.
Qs = 8 m 2
/yr 1500x
0.4 vert. exag.
(closed basin)
5.5 Flexural-isostatic response to
0.2 forebulge
glacial erosion in the western
sensitivity to Qs
0.0 US
infinite basin
forebulge erosion

(b) 0.8
The pace of alpine glacial erosion increased sig-
ni¬cantly starting 2--4 Myr ago when global cool-
h ing initiated the Plio-Quaternary era of large con-
10,000 m2/yr
(km) tinental ice sheet advances and retreats (Zhang
0.4 et al., 2001). This erosion must have triggered sig-
5000 m2/yr
ni¬cant ¬‚exural-isostatic unloading of glacially
2500 m2/yr carved mountain belts. Since glacial erosion is
0.2 κ=
strongly focused in the landscape, scouring out
1000 m2/yr
sensitivity to κ
valley ¬‚oors where ice ¬‚ow is focused but erod-
0.0 infinite basin
forebulge erosion
ing the highest points in the landscapes rela-
800 1000
400 600
0 tively slowly, it is likely that glaciated mountain
x (km) belts have increased in relief during the Plio-
Quaternary era due to isostatic rock uplift, even
Fig 5.8 (a) and (b) Model solutions for the topographic
in areas that are tectonically inactive. Determin-
pro¬le of an in¬nite basin assuming complete forebulge
ing precisely how much glacial erosion has taken
erosion, showing sensitivity to variations in Q s and κ,
place is dif¬cult to do in most cases. Nevertheless,
respectively, for the same model parameters as Figures 5.7a
it is still useful to ask how much rock uplift oc-
and 5.7b. These results show that model results are broadly
similar to those obtained by assuming no forebulge erosion; curred for a given amount of glacial erosion. By
the only signi¬cant difference is a minor reduction in basin using the geographic extent of glacial cover as a
relief when forebulge erosion is assumed.
proxy for the spatial distribution of erosion, we
can use the 2D ¬‚exural-isostatic model to com-
pute the ratio of late Cenozoic erosionally driven
the two southern pro¬les using F = 10 mm/yr. rock uplift to glacial erosion. In this section we
Given observed or inferred values for h b , L b , and compute that ratio for the western United States,
κ, as well as estimated values for F , ±, and w 0 following the approach of Pelletier (2004a). Ar-
from the published literature, a family of solu- eas where this ratio is large have likely under-
tions corresponding to a range of Q s values can gone the greatest relief production and should
be generated for each pro¬le and compared to be the focus of future efforts to identify signa-
the observed data. Figure 5.9c illustrates the ob- tures of erosionally driven rock uplift and re-
served pro¬les plotted with their corresponding lief production. In this application we will use
best-¬t pro¬les (κ values for the two southern- the Alternating-Direction-Implicit technique, tak-
most pro¬les were also varied to ¬nd the optimal ing into account spatial variations in ¬‚exural
¬t for both κ and Q s ). Model solutions match rigidity.
the observed data very closely for the strongly The ratio of rock uplift to erosion depends on
over¬lled pro¬le, and they match the ¬rst-order the area and shape of the eroded region and on
trends of the weakly over¬lled pro¬le. Discrep- the local ¬‚exural wavelength of the lithosphere.
ancies in the southernmost pro¬les may be the If the ¬‚exural wavelength is large compared to

68° W 57° W 0.8 1
k = 1500 m2/yr
17° S h ’ hb
(km) k = 4000 m2/yr
sediment-rich (km)
0 200 600
x (km)

0 200 400
x (km)
Rio Paraguay

k = 1500 m2/yr
Qs = 9 m2/yr
h ’ hb
sediment-starved 0.4 k = 2000 m2/yr
Qs = 6.5 m2/yr

k = 2000 m2/yr
Qs = 5 m2/yr

100 km (a)
31° S L. Mar Chiquita
L. Ambargasta
0 200 400
100 400 2000 6000 m
0 x (km)
spatial distribution of glaciated regions in the
Fig 5.9 (a) Grayscale map of topography in the central
western United States. Porter et al. (1983) mapped
Andes indicating sediment-rich, strongly over¬lled basins in
the northern central Andes (near 20—¦ S) and areas of Pliocene--Quaternary glaciation in the
western United States, based on the elevations
sediment-starved, under¬lled and weakly over¬lled basins
(near 30—¦ S) corresponding to a humid-to-arid transition. of the lowest cirques. Their map identi¬es the
(b) Basin topographic pro¬les extracted from SRTM data. approximate local elevation of vigorous glacial
Also plotted using logarithmic scale (with base-level elevation
erosion. To construct a map of the extent of late
subtracted) to illustrate exponential trend and associated
Cenozoic glacial erosion, the contours of glacial-
best-¬t κ values. (c) Best-¬t model pro¬les with associated
limit elevation from Porter et al. (1983) were ¬rst
Q s values. Sediment supply decreases from north to south
digitized, projected to Lambert Equal-Area, and
with precipitation rates and as deformation becomes more
then interpolated to obtain a grid of glacial-
spatially distributed. Note that 100 m has been added to the
limit elevations. Second, a 1 km resolution DEM
northernmost pro¬le so that the plots can more easily be
(digital elevation model) of the western United
distinguished. Modi¬ed from Pelletier (2007). Reproduced
with permission of the Geological Society of America. States was recti¬ed to the glacial-limit-elevation
grid. Third, the glacial-limit-elevation and DEM
grids were compared pixel-by-pixel to determine
the width of the glaciated region, uplift will be whether local elevations were above or below
limited by ¬‚exural bending stresses. In contrast, the glacial-limit elevation. Pixels with elevations
if the ¬‚exural wavelength and the width of the above the glacial limit were considered areas of
eroded region are comparable, bending stresses glacial erosion and were given a value of 1. Con-
will be smaller and uplift will approach the Airy versely, pixels with elevations below the glacial
isostatic limit in which 80% of the eroded mate- limit were considered not to have been subjected
rial is replaced through rock uplift. to glacial erosion and were given a value of 0.
To map the compensation resulting from The resulting binary grid is a map of relative
glacial erosion, we ¬rst must determine the glacial erosion; an underlying assumption is that

to gravity, and q(x, y) is the vertical load (Watts,
117° W 103° W
49° N
2001). In order to account for spatial variations
300 km in elastic thickness we retained D (x, y) inside the
parentheses in Eq. (5.41). This results in addi-
tional terms involving gradients of D that are
not present when D is assumed to be spatially
42° N
The binary grids of glacial erosion were di-
rectly input into the ¬‚exure equation as the
loading term, so that q(x, y) = 1 where erosion
occurred and q(x, y) = 0 where erosion did not
occur. The solution to the ¬‚exure equation with
these loading terms is the ratio of the positive de-
¬‚ection of the lithosphere to the depth of glacial
SJ erosion. This is the compensation resulting from
glacial erosion. Here we use the alternating dif-
ference implicit (ADI) technique to solve Eq. (5.41),
following the work of van Wees and Cloetingh
109° W
Te (1994). Previous techniques we have discussed
0 100 km
(e.g. Fourier ¬ltering) cannot be used when D
Fig 5.10 Grayscale map of elastic thickness in the western varies spatially.
United States. Data from Lowry et al. (2000) and Watts For the analysis, we will use Lowry et al.™s
(2001). Boundaries of glaciated areas are shown in white.
(2000) data set for elastic thickness (T e ) in the
Labeled ranges are: SN (Sierra Nevada), SJ (San Juan), and
western United States. These data are based on
YP-A (Yellowstone Plateau-Absaroka). Modi¬ed from
the coherence method and provide a direct mea-
Pelletier (2004a).
sure of lithospheric strength that includes the
weakening effects of faults. These data have the
erosion is uniform within the glaciated regions. highest resolution available for the region, but do
The boundaries of the glaciated regions obtained not completely cover the western United States.
in this way are mapped in white in Figure 5.10 To cover the entire area, the Lowry et al. data was
(which also shows state boundaries for reference). augmented with the global compilation of Watts
Also shown in Figure 5.10 is a grayscale map of (2001), and the data sets joined smoothly with a
the elastic thickness in the western United States. low-pass ¬lter. The resulting T e map was recti¬ed
Ranges along the Paci¬c coast were not included to the glacial erosion maps and is shown as a
in Porter et al. (1983) and have also been ex- grayscale map in Figure 5.10. T e values were con-
cluded from this analysis. In addition, the north- verted to a grid of ¬‚exural rigidity, D (x, y), for
ern Cascades were excluded because they were use in Eq. (5.41), by means of the relationship
completely covered by the Cordilleran Ice Sheet
E T e3
and hence subject to a different style of glacia- D= (5.42)
12(1 ’ ν 2 )
tion than other ranges of the western United
States. assuming E = 70 GPa and ν = 0.25.
The de¬‚ection of a thin elastic plate with vari- A map of glacial-erosion compensation in the
able elastic thickness under a vertical load is western United States is given in Figure 5.11a.
given by Also shown in Figure 5.11b is the compensa-
tion determined by assuming that glacial ero-
∇ 2 (D (x, y)∇ 2 w) + ρgw = q(x, y) (5.41)
sion was concentrated near the equilibrium line,
where w is the de¬‚ection, D (x, y) is the ¬‚exural not uniformly within the glacial region as in Fig-
rigidity, ρ is the density contrast between the ure 5.11a. Although compensation is strictly de-
crust and the mantle, g is the acceleration due ¬ned to be between 0 and 1, we have mapped

103° W
103 W 117° W
117° W 49 N
49 N
(a) (b)
300 km

42 N

uniform erosion above glacial-limit erosion concentrated near glacial-limit
109° W
109° W
0 compensation 1
Fig 5.11 Map of compensation of glacial erosion assuming Plateau--Absaroka Range. The high values in this
(a) uniform erosion within glaciated areas and (b) erosion region re¬‚ect the extensiveness of the glaciated
concentrated within 500 m altitude of glacial limit. For color
area, which is nearly 200 km in width. In ad-
version, see plate section. Modi¬ed from Pelletier (2004a).
dition, however, nearby ranges are also being
eroded, reducing the bending stresses that would
otherwise be concentrated at the edges of the
the full solution to the ¬‚exure equation in Fig-
Yellowstone area, limiting its potential uplift. In-
ure 5.11, including areas of minor subsidence (in
stead, nearby ranges, including the Beartooth,
black). All of the major ranges in the western
Gallatin, Wyoming, and Wind River Ranges, ac-
United States have compensation values greater
commodate the bending stresses to produce a
than or equal to 0.4, indicating that nearly half
˜˜bulls-eye™™ over the central part of the area -- the
of all glacially eroded mass is returned to the
Yellowstone Plateau--Absaroka Range.
region™s ranges as erosionally driven rock uplift.
The pattern of compensation varies little between
Figures 5.11a and 5.11b because ¬‚exure is insen-
sitive to small-scale variations in loading.
The greatest difference between Figures 5.11a
and 5.11b is the mean values: values of compen- 5.1 Calculate the degree of isostatic compensation
for continental topography with a wavelength of
sation obtained by assuming concentrated ero-
20 km and an elastic thickness of 5 km, assum-
sion near the glacial limit are 20% lower than
ing E = 70 GPa, nu = 0.25, ρm = 3300 kg/m3 , and
for uniform erosion. These lower values re¬‚ect
ρc = 2700 kg/m3 .
the additional bending stresses present in the
5.2 In addition to driving rock uplift, the isostatic re-
case of localized erosion; i.e., for the same eroded
sponse to erosion also causes subsidence in some
mass, a circular area (i.e., uniform erosion within
locations. Consider a 1D model of periodic line un-
the range) results in less intense bending stresses loading (Figure 5.12). Each line unloading might
than a ˜˜ring-shaped™™ area (i.e., localized erosion represent glacial erosion concentrated in the high
near the glacial limit). peaks of a Basin and Range-style topography. Calcu-
The highest compensation values in the west- late the distance between unloading, L , that pro-
ern United States are found in the Yellowstone duces the largest value of subsidence in spaces

L LIS (9 ka)

Fig 5.12 Schematic diagram of the 1D periodic line loading
calculation used in Exercise 5.2.

Lake Agassiz
between the locations where unloading is concen-
5.3 Calculate the lithospheric de¬‚ection associated
with a narrow seamount of radius 10 km and
height 1 km using the Kelvin function solution.
5.4 Use the Fourier series approach to construct an
expression for the ¬‚exural response of the litho- Fig 5.13 Outline of the Laurentide Ice Sheet at 9 ka,
sphere to a periodic topographic load of trian- showing the extents of proglacial lakes at the ice sheet
gular mountain ranges with width L and height margins.
5.5 Proglacial lakes formed at the margins of the
this model, assuming that the lake ¬lls up to sea
Laurentide Ice Sheet. The largest of these, Lake
Agassiz, had a width of several hundred kilome-
5.6 Following the Andean example, download topo-
ters (e.g. Figure 5.13). Model the 1D bathymetric
graphic data for a portion of the Main Central
pro¬le of Lake Agassiz perpendicular to the ice
Thrust of the Himalaya. Using published estimates
margin assuming that the lake formed as a ¬‚ex-
for the elastic thickness of the lithosphere, model
ural depression. Assume that the proglacial to-
the 1D ¬‚exural pro¬le and determine the ap-
pography was initially ¬‚at and at sea level, and
proximate location of the Himalayan forebulge in
that the ice sheet acts as a rectangular load with
height 2 km. How wide is the lake according to
Chapter 6

Non-Newtonian ¬‚ow equations

extend in¬nitely far into and out of the cross sec-
6.1 Introduction tion. The horizontal force F exerted on the left
side of a thin vertical slice of the ¬‚uid of width
x is equal to the total hydrostatic pressure in
Chapter 1 introduced the concept of non-
the ¬‚uid column. The hydrostatic pressure at a
Newtonian ¬‚uids and highlighted their impor-
depth y below the surface is given by ρgy. The
tance for our understanding of the mechanics
total hydrostatic pressure is calculated by inte-
of glaciers. In this chapter we explore methods
grating ρgy over the ¬‚uid column from y = 0 to
for modeling non-Newtonian ¬‚uids with applica-
y = h:
tions to ice sheets, glaciers, lava ¬‚ows, and thrust
sheets. In most cases we will make the simplify-
ρgh 2
F (x) = ρgydy =
ing assumption that the rheology of the ¬‚ow is (6.1)
constant in time. In nature, however, many types
of non-Newtonian ¬‚ows are complicated by the On the right side of the vertical slice (i.e. at po-
sition x + x), the pressure is lower because the
fact that the rheology of the ¬‚ow depends on
parameters that evolve with the ¬‚ow itself. The ¬‚ow is thinner:
viscosity of lava, for example, depends on its tem-
ρg(h ’ ρg 2
F (x + x) = = (h ’ 2h h
perature, which, in turn, depends on the ¬‚ow
2 2
behavior (e.g. slower moving lava will cool more ρg 2
+ ( h)2 ) ≈ (h ’ 2h h) (6.2)
quickly, increasing its viscosity and further slow- 2
ing ¬‚ow in a positive feedback). Here we focus
The ( h)2 term in Eq. (6.2) is much smaller than
primarily on modeling ¬‚ows with a constant, pre-
the 2h h term for ¬‚ows with small (i.e. < 0.1) sur-
scribed rheology. In cases where the rheology de-
face slopes, and hence the second-order term can
pends strongly on the state of the ¬‚ow the meth-
be neglected in such cases. This ˜˜small-angle™™ ap-
ods introduced in this chapter will serve as the
proximation is commonly used when modeling
foundation for more complex models which fully
ice sheets, glaciers, and lava ¬‚ows on Earth. The
couple the rheology with the ¬‚ow.
difference between the horizontal force exerted
on the opposite sides of the ¬‚uid slice must be
equal to the shear stress „0 exerted over the width
6.2 Modeling non-Newtonian and of the slice, x:
perfectly plastic ¬‚ows
’„0 = F (x + x) ’ F (x) = ρgh h (6.3)

First we consider the forces acting on a ¬‚uid Rearranging Eq. (6.3) gives
spreading over a ¬‚at surface (Figure 6.1). We be-
„0 = ’ρghS
gin with a 1D example, so the ¬‚uid is assumed to (6.4)

”x (a)
h4 t0 = 1.5 bar
1.0 bar
F 3
F+ ”F
0.5 bar
Fig 6.1 Schematic diagram of the forces acting on a gravity
¬‚ow atop a ¬‚at surface. 0
1.5 bar
where S is equal to the local ice-surface gradient, 3
h/ x. Field studies of modern ice sheets and 1.0 bar
glaciers have found the basal shear stress calcu- 2
lated with Eq. (6.4) to be typically between 0.5
and 1.5 bars (Nye, 1952), a surprisingly narrow 1
0.5 bar
range considering the spatial variability of ob-
served basal conditions, including gradients in 0
temperature, meltwater content, basal debris, till
rheology, and other variables.
200 300
100 400 500
6.2.1 Analytic solutions x (km)
Equation (6.4) is a general result for any ˜˜slowly™™
moving gravity ¬‚ow (i.e. slow enough that acceler- Fig 6.2 Plots of ice sheet pro¬les predicted by the perfectly
plastic model with (b) and without (a) isostatic adjustment
ations can be neglected). In cases where the ¬‚uid
included, for a range of threshold basal shear stresses.
moves rapidly when a threshold shear stress is
exceeded, Eq. (6.4) can be used to model the ¬‚ow
geometry by treating „0 as a constant. Taking the the right side of Eq. (6.7) has units of length
limit as x goes to zero and expressing Eq. (6.4) and is often written as h 0 (not to be confused
as a differential equation gives with an initial, or maximum thickness). For an
„0 ice sheet with „0 = 1 bar (105 Pa), ρ = 920 kg/m3 ,
h (6.5)
and g = 9.8 m2 /s, h 0 ≈ 12 m. It should be noted
that this approach breaks down near the divide,
Multiplying both sides of Eq. (6.5) by dx and in-
where additional stress components become im-
tegrating yields
h2 = x+C (6.6) The parabolic solution of Eq. (6.7) is usually
associated with perfectly-plastic deformation of
The integration constant C is equal to the termi- near-basal ice. Perfectly plastic deformation, how-
nus position. If the ¬‚ow has a half-width of L , for ever, is an idealization of just one of the com-
example, then C = L , and Eq. (6.6) becomes plex basal ¬‚ow mechanisms we recognize today,
including debris-controlled frictional sliding and
h= (L ’ x) (6.7) till deformation. So, how can Eq. (6.7) provide
even an approximate model for real ice sheets
Figure 6.2a illustrates several examples of this and glaciers? The reason is that many deforma-
simple parabolic ice-sheet model solution. Equa- tion mechanisms mimic the threshold-controlled
tion (6.7) is usually referred to as a ˜˜perfectly ¬‚ow of the perfectly plastic model. Many ¬‚ow
plastic™™ model for ice sheets. The ratio „0 /ρg on mechanisms are threshold controlled and result

in low effective friction when a threshold shear Plots of Eqs. (6.12) and (6.13) are shown in Figure
6.2b for the same cases as Figure 6.2a with δ =
stress is exceeded. Threshold behavior, for exam-
ple, characterizes the plastic deformation of sub- 0.4. Equation (6.13) indicates that the volume of
glacial till (e.g. Tulaczyk et al., 2000) and the shear an ice sheet increases when isostatic effects are
stress necessary to overcome kinetic friction dur- taken into account. In essence, the ˜˜hole™™ created
ing basal sliding (e.g. Lliboutry, 1979). The per- by the weight of the ice provides an additional
fectly plastic model, therefore, can be generally con¬ning force that enables the ice to be thicker
applied to model the geometry of ice sheets and than it would be on a ¬‚at bed.
glaciers that move at low effective friction when Equation (6.4) can be generalized to cases with
a threshold shear stress is exceeded, whether rel- variable bed topography. In such cases, the form
ative motion occurs internally or by sliding over of Eq. (6.4) does not change, but the surface slope
deforming till or directly over the bed. To em- S is no longer simply equal to dh/dx. Instead, S
is equal to dz/dx = dh/dx + db/dx where b(x) is
phasize this broader applicability, we will also
use the term ˜˜threshold-sliding™™ to refer to the the bed elevation. Equation (6.4) then becomes
perfectly plastic model.
dh db
+ =
The results of Figure 6.2a are unrealistic in h (6.14)
dx dx
at least one major respect: the subsidence of the
crust beneath the ice sheet was neglected. For
Equation (6.14) can be used to derive an ana-
ice sheets larger than a few hundred kilometers,
lytical expression for the approximate shapes of
the weight of the overlying ice will cause signif-
alpine glaciers with a uniform bed slope. In such
icant bedrock subsidence beneath the ice sheet
cases, db/dx = Sb and Eq. (6.14) becomes
in order to achieve isostatic equilibrium. This ef-
fect can be quanti¬ed with a simple adjustment dh
+ Sb = h0
h (6.15)
to the threshold-sliding model. If the bed eleva- dx
tion is given by b(x), then isostatic equilibrium
where h 0 = „0 /ρg. Dividing both sides of Eq. (6.15)
requires that
by h gives
ρi z = ’ρm b (6.8)
dh h0
= Sb ’
where ρi is the density of ice, z = b + h is the (6.16)
dx h
ice-surface elevation, and ρm is the density of the
which integrates to
mantle. Substituting z = b + h, Eq. (6.8) becomes
ρi h h0
b=’ z = ’δz x’C = + 2 ln(h 0 ’ Sb h)
(6.9) (6.17)
ρm ’ ρi Sb Sb
where where C is the integration constant. Figure 6.3a
ρi illustrates solutions to Eq. (6.17) for C = 0, Sb =
δ= (6.10)
ρm ’ ρi 0.05 and 0.1, and h 0 = 12 m (i.e. „0 = 1 bar). In or-
der to better visualize the results, it is helpful to
The ice thickness can be written as
plot z = ’Sb x + h instead of simply h. Figure 6.3b
h = (1 + δ)z (6.11)
presents these results. The value of C can be var-
ied to shift the glacier terminus to the right or
Substituting Eq. (6.11) into Eq. (6.7) yields the
left, but otherwise the value of C does not affect
following equations for the ice-sheet elevation
the solution.
and thickness:
Another solution relevant to alpine glaciers
z= (L ’ x)
2 (6.12) involves the variations in thickness caused by mi-
ρg(1 + δ)
nor variations in valley slope. If the bed topog-
raphy is slowly-varying and the region close to
„0 (1 + δ) the glacier terminus is not considered, the dh/dx
h= (L ’ x)
2 (6.13)
ρg term can be neglected relative to the db/dx term

250 300
h 200 (km) h
Sb = 0.05
(m) 2.0 200

100 100
Sb = 0.1
10 30
x (km)
Fig 6.4 Plots of ice thickness variations in an alpine glacier
with slowly varying bed topography, according to Eq. (6.18),
z 400 with h 0 = 12m.

0 scribed, however, the geometry of the ice sheet
or glacier does not explicitly depend on mass bal-
ance. When considering ice bodies obeying Glen™s
Sb = 0.1
Flow Law, however, it is necessary to prescribe a
Sb = 0.05 mass balance rate M . Here we will assume M to
be uniform and constant, but in general it will
10 15
0 5
x (km) vary in both space and time.
The starting point for modeling ice-sheet or
Fig 6.3 Plots of (a) thickness and (b) ice-surface elevation
glacier geometries with a power-law ¬‚uid rheol-
for Sb = 0.05 and 0.1 for h 0 = 12 m (i.e. „0 = 1 bar).
ogy is the relationship between shear stress and
strain rate, given in Chapter 1 for the case of sim-
to give ple shear:
db du
h = h0 = A„ n
(6.18) (6.19)
dx dy
Figure 6.4 plots solutions to Eq. (6.18) for h 0 = where u is the ¬‚uid velocity, y is the depth, „ is
12 m, a mean valley slope of 0.05, and an ideal- the shear stress, and A and n are rheological pa-
ized sinusoidal bed topography. The ice thickness rameters. Substituting Eq. (6.4) into Eq. (6.19) and
is inversely proportional to the bed slope, thick- integrating from 0 to h gives the velocity pro¬le
ening in areas of gentle slopes and thinning in
areas of steep slopes. u(y) = ’ (ρgS)n (h ’ y)n+1 ’ h n+1 (6.20)
2(n + 1)
When considering the geometry of ice sheets
and glaciers within the framework of the per-
Vialov (1958) was the ¬rst to obtain an analytic
fectly plastic model, it is not necessary to pre-
solution for the pro¬le of an ice sheet on a ¬‚at
scribe the accumulation or ablation rates because
bed evolving according to a power-law ¬‚uid rhe-
the shape of the ice sheet or glacier is con-
ology. Vialov assumed a radially symmetric ice
trolled entirely by the threshold shear stress „0
sheet with radius R . Equation (6.20) must be in-
(given standard values for ice density and gravity).
tegrated from 0 to h to obtain the total horizontal
Changes in the spatial distribution of accumula-
¬‚ux of ice:
tion may cause the ice terminus to migrate. As
such, the position of the terminus is controlled dh
U = ’C h 4 (6.21)
by mass balance. If the terminus position is pre- dr

where C is a single constant that combines pa-
rameters for viscosity, density, and gravity. Con- perfectly plastic
servation of mass is given by h
‚h ‚(r hU ) Vialov
= ’∇ · (H U ) + M = ’ +M (6.22)
‚t r dr 2.0
For an ice sheet in steady state,
d(r hU ) = M r dr (6.23)
which integrates to
U= 0.0
2h 1000
400 600 800
r (km)
Combining Eqs. (6.21) and (6.24) yields
M Fig 6.5 Comparison of the perfectly plastic model solution
dh = ’
r 1/n dr
h (6.25)
(with „0 = 1 bar) to the Vialov solution with Glen™s Flow Law
with a maximum thickness of 2.5 km.
which can be integrated to obtain
n n M
(h 2+2/n ’h 2+2/n ) = ’ r 1+1/n (6.26) Figure 6.5 prescribes h c = 2.5 km. The Vialov solu-
2n+2 n+1 2C
tion predicts gentler slopes in the ice sheet inte-
where h c is the thickness at the center where
rior and steeper slopes near the margin relative
r = 0. Equation (6.26) can be rearranged to give
to the perfectly plastic model.
’ = ’2
h 2+2/n r 1+1/n
h (6.27)
2C 6.2.2 Numerical solutions
Thus far we have considered only ˜˜one-sided™™
The value of h c can be constrained by setting a
ice sheets and glaciers, i.e. those with elevations
boundary condition of vanishing ice thickness at
the margin, i.e. h(R ) = 0, to give that decrease monotonically from a known di-
vide position or increase from a known terminus
h 2+2/n R 1+1/n position. Equation (6.14) is limited to one-sided
ice sheets because dh/dx + db/dx must be posi-
Substituting Eq. (6.28) into Eq. (6.27) gives the full tive throughout the domain, increasing from the
margin (at x = 0) to the divide (at x = L ). In the
1/(2n+2) more general case that includes variable bed to-
h=2 n/(2n+2)
pography, only the positions of the two margins
are known and the topography on both sides of
R 1+1/n ’ r 1’1/n (6.29)
the ice sheet must be solved for. In such cases, the
For n = 3, Eq. (6.29) becomes divide position cannot be prescribed but must
be solved for along with the geometry. For such
M 3/8
h=2 R 4/3 ’ r 4/3
(6.30) cases, Eq. (6.14) must be modi¬ed to include both
positive and negative ice-surface slopes to obtain
It is dif¬cult to compare the Vialov solution „
dh db
+ = (6.31)
precisely to the perfectly plastic solution because ρgh
dx dx
one was derived for a 1D ice-sheet pro¬le and
The boundary conditions for Eq. (6.31) are h = 0
the other for a radially symmetric ice sheet. Nev-
ertheless, Figure 6.5 compares solutions obtained at the two ice margins.
for the two models. The ratio M/C acts as a free Equation (6.31) can be solved with stan-
parameter in the Vialov solution that relates to dard methods for boundary-value problems. The
the maximum height h c . The solution plotted in equation must ¬rst be modi¬ed to remove the

absolute-value sign by squaring both sides: 2.5
h 2.0
2 2
dh db
+ = “one-sided”
ρgh solution
dx dx
and then by differentiating both sides of Eq.
(6.32) to obtain:

„2 0.5
d2 h d2 b 1
+ 2 = ’ 20 2 (6.33)
ρ g h3 +
dh db
dx 2 dx dx dx
200 800 1000
400 600
x (km)
Equation (6.33) is nonlinear and singular at
divides (i.e. where dh/dx + db/dx = 0). As such, Fig 6.6 Comparison of analytic and numerical perfectly
care must be taken to solve it. One approach is plastic solutions for an ice sheet over a ¬‚at bed with „0 = 1


. 7
( 14)