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

112 FLEXURAL ISOSTASY

thickness of 40 km, for example, the ¬‚exural pa-

200

rameter under crustal loading is ≈ 130 km, while

for ice loading it is approximately 90 km.

crust

150

(km) 5.2.2 Series solutions and Fourier ¬ltering

The ¬‚exure equation is a linear equation. As such,

100

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-

0

60

20 40 80 tude h 0 and wavelength »:

0

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-

± ±

dx

eter of the lithosphere, but there will be no off-

Therefore,

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

h0

to about 4% of the maximum de¬‚ection of the

w0 = (5.20)

ρm 2π 4

’1+

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,

5.3 METHODS FOR 2D PROBLEMS 113

where (so called because it ¬lters out low wavenumbers)

before the inverse Fourier transform is applied to

1

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

ρc

w0 = h0 output of that program for an east--west pro¬le

(5.23)

ρ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

∞

πx

h(x) = height of the ¬‚exural forebulge (Figure 5.4c) also

an sin (5.24)

L increases with higher values of the ¬‚exural pa-

n=1

rameter. Gravity modeling has been used to con-

with coef¬cients

strain the width of the Andean foreland basin by

L

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-

0

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

Andes.

sin

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.

114 FLEXURAL ISOSTASY

(a) 5

4

h3

(km)

(c) 0.4

2

w

0.2

1 (km)

0

0.0

(b) 0

± = 50 km

75 km

0.2

100 km

w5

0.4

(km) ± = 50 km

1200

800 1000 1400

75 km

100 km x (km)

10

15

600 800 1000 1200

400

200

0

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 ±.

∞ ∞

ρc

First we consider the solution to the 2D w(r ) =

ρm ’ ρc

¬‚exural-isostatic equation: 0 0

√

2

— kei (x ’ x )2 + (y ’ y )2

D ∇ 4 w + ρgw = q ±

(5.27)

— 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

2r

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

5.3 METHODS FOR 2D PROBLEMS 115

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 =

1.0

4.0

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

3

point sources and summing the contributions of

each point. This is a powerful approach, but it

4

can be computationally challenging when calcu-

’1 0 1

lating the response at high spatial resolution over

x/a

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

tor

¬‚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 ± ± ±

(5.32)

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

116 FLEXURAL ISOSTASY

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.,

1992).

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