<<

. 12
( 14)



>>

functional relationship between temperature change, T / t,
et al., 1981)
and global temperature, T , for each process. (a) A high global
temperature results in more outgoing long-wavelength
Tn
= ’c 1 T n + D ·n (8.39)
radiation (and relative cooling). This effect is proportional to
t
the temperature difference from equilibrium where the
where · is a white noise with standard deviation
constant of proportionality, c 1 , is equal to the radiation
damping time, constrained to be c 1 = 180 yr as described in equal to D . Short-term internal temperature vari-
the text. (b) Lower global temperatures result in increased ations generate signi¬cant low-frequency temper-
albedo (and relative cooling) with a function that depends on ature variations because the climate system inte-
the ice sheet distribution in Figure 8.25. (c) Long-wavelength
grates the short-term random variations through
outgoing radiation and the ice-albedo feedback taken
their effect on the outgoing radiation ¬‚ux. The
together yield a system with two ¬xed points, corresponding
result of this process is Brownian-walk behavior
to glacial and interglacial states. (d) Bedrock subsidence and
for global temperature. Equation (8.39) has the
rebound modify the effects of the ice-albedo feedback by
behavior of a damped Brownian walk, however, be-
shifting the curve towards greater relative warming by an
amount proportional to |v |. From Pelletier (2003). cause the equation also includes a negative feed-
back term.
The growth of ice sheets provides a positive
system equal to 180 yr (the inverse of c 1 ). The feedback that dominates the negative feedback
simple linear relationship between temperature of outgoing radiation if ice sheets grow nonlin-
change and mean global temperature in Eq. (8.38) early with global temperature. The effects of the
is plotted in Figure 8.24a. ice--albedo feedback on global climate may be
The negative feedback of Eq. (8.38) is balanced modeled by including a term in Eq. (8.38), a(T ),
by natural variability in Earth™s radiation balance which re¬‚ects the dependence of Earth™s albedo
driven by short-term (decadal and centennial) on global temperature, to obtain
¬‚uctuations in atmospheric temperature. These Tn
= ’c 1 T n ’ c 2 a(T n ) (8.40)
¬‚uctuations create random variability in long- t
214 STOCHASTIC PROCESSES


We assume that the modern interglacial state is Vostok data provide the independent data to de-
a steady state of the climate system. This provides termine a(T ). The albedo data were obtained by
multiplying each 1—¦ — 1—¦ grid square in Peltier™s
the constraint a(0) = 0 and ¬xes the temperature
of any other stable states of the system to be with data by its area and the amount of incident solar
respect to an origin at the modern interglacial radiation it receives to obtain the relative change
temperature. in albedo as a function of global temperature:
Stable states of Eq. (8.40) are determined by 180 90
setting the right hand side equal to zero. T = a(T ) = ai, j (T ) cos2 ( j) (8.41)
i=’180 j=’90
0, the modern interglacial, is a ¬xed point be-
cause we de¬ned a(0) = 0. Whether or not an- where ai, j = 0.25 if the grid square is ice-free and
ai, j = 0.85 if the grid square is ice-covered.
other ¬xed point exists depends on how a(T ) in-
creases with decreasing T . If the Earth™s albedo Equation (8.41) provides an estimate of the rel-
increases greater than linearly with T , the pos- ative change in Earth™s albedo due to the pres-
itive ice--albedo feedback dominates the system. ence of Late-Pleistocene continental ice sheets.
A dominant ice--albedo feedback leads to an ice- The albedo difference between modern and LGM
covered Earth if a(T ) increases greater than lin- conditions calculated in Eq. (8.41) equals 10% of
early for all T . However, if the ice--albedo feed- the modern albedo assuming that ice-free areas
back increases less than linearly for any T , the are also continuously snow free as well. 10% is
combination of outgoing radiation and the ice-- an overestimate, however, because many areas as-
albedo feedback will be negative or self-limiting sumed to be ice-free in Eq. (8.41) have seasonal
for those values of T . snow and ice cover. Accounting for modern sea-
Determining the form of a(T ) is complicated sonal snow and sea-ice cover, Oerlemans (1980a)
by uncertainties in the geographic distribution of estimated that the albedo difference between the
ice as a function of global temperature. Most nu- LGM and the present was approximately 5%. We
merical ice-sheet-climate models determine a(T ) have used this value to scale the albedo curve to
by ice-sheet dynamics in their full complexity, a maximum change of 5% between the present
including internal ¬‚ow, basal sliding, and mass and the LGM. In other words, we have used Oer-
wasting. The complexity of these models results lemans™ albedo estimate for the total change be-
in dozens of free parameters for the motion of tween the LGM and the present, and Peltier™s
continental ice sheets. These processes, especially reconstructions to estimate the relative change
basal sliding, are poorly constrained. The distri- through time. The result is plotted in Figure 8.25.
bution of sea ice as a function of global tempera- The data of Figure 8.25 are well approximated
ture is also highly uncertain. As an alternative to by a square-root function:
modeling ice-sheet dynamics, the extent of land 1
(T h ’ T ) 2 T < Th
if
ice through time during the last deglaciation a(T ) = (8.42)
T ≥ Th
0 if
may be inferred from the pattern of postglacial
where T h = ’1—¦ C (relative to the modern tem-
rebound (Peltier, 1994). Peltier™s reconstruction
provides ice extents at 1-kyr intervals from the perature) is the temperature at which subpolar
Last Glacial Maximum (LGM) to the present. His ice-sheet growth is initiated in Hudson Bay. The
reconstruction provides data on ice-sheet topog- value of the coef¬cient of a(T ), c 2 , is given by
raphy as well as extent, but only the data on ice-
Q a(T g ) ’ a(T h )
c2 =
sheet extent are necessary to constrain variations (8.43)
C (T h ’ T g ) 1
2
in albedo. We can estimate the relationship be-
where T g is the temperature of the full glacial
tween global albedo and temperature from the
climate state. Q , the incoming solar energy, is
LGM to the present using Peltier™s ice extents and
equal to 340 J/s/m2 (North, 1991). Equation (8.43)
the Vostok time series. First, a global-temperature
yields c 2 = 5.4 — 10’10 —¦ C1/2 /s. Scaled to the ra-
time series can be constructed by averaging the
diative time scale of 180 yr this is c 2 = 3.0 —¦ C1/2 .
Vostok time series in 1-kyr intervals centered
The relationship between incoming energy and
on each integer value of kyr. The time-averaged
8.9 COHERENCE RESONANCE AND THE TIMING OF ICE AGES 215


continues only as long as ice sheets expand non-
6
linearly with global temperature. The fact that
5 ice sheets grow preferentially on the continents
eventually weakens the positive feedback as the
4 North American ice sheets expand into a taper-
”a
ing continent. The result is a second stable state
(%) 3
characterized by an equilibrium between the ice-
-albedo feedback and the negative feedback of
2
outgoing radiation. The square-root form may be
approximately understood by considering the
1
constraints on the North American ice sheets (in-
0 cluding the Laurentide and former ice sheets) as
’10 8 6 4 2 0
they expand from their initial location in Hudson
”T (°C)
Bay. When the North American ice sheet is ex-
panding southward and westward across Canada
Fig 8.25 Observed relationship between albedo and global
during its initial growth, the Earth™s albedo is
temperature based on Peltier™s (1994) gravitationally
self-consistent sea-level inversion and Vostok temperatures sensitive to small changes in global temperature.
for the same time periods. The data closely approximate a The result is a strong positive feedback between
square-root function between albedo and global temperature cooler temperatures and larger ice sheets. The ice
1
(a(T ) = c 2 (Th ’ T ) 2 with Th = ’1). The shape of this curve
sheet eventually reaches its maximum westward
largely re¬‚ects the sensitivity of the Laurentide ice sheet
extent in northern Canada as the global temper-
coverage to global temperatures as the ice sheet retreated
ature further drops, however. If we assume that
westward and southward from Hudson Bay during the later
ice sheets expand preferentially on land, the ice-
stages of retreat. Earlier stages of retreat were not as
-albedo feedback weakens at this stage because a
sensitive because ablation resulted in limited retreat due to
further decrease in global temperature results in
the tapered ice sheet geometry in southern Canada. From
Pelletier (2003). smaller changes in ice-sheet extent because ex-
pansion is limited to a southerly direction and
a tapering North America. The combined effect
global temperature is plotted in Figure 8.24b. The of long-wavelength outgoing radiation and the
combination of the ice--albedo feedback and long- ice--albedo feedback is positive and of large mag-
wavelength outgoing radiation are illustrated in nitude for small ice sheets, but becomes progres-
Figure 8.24c. The combined effect of these pro- sively weaker as ice sheets grow in extent and are
cesses results in two ¬xed points in the model eventually geographically limited in where they
system corresponding to the full glacial and can expand. As global temperature continues to
interglacial states. drop, a balance is eventually achieved between
Bistability in the climate system is based on the long-wavelength outgoing radiation and the
a balance between the negative feedback of long- increase in Earth™s albedo. This is the full-glacial
wavelength outgoing radiation and the positive climate.
ice--albedo feedback. The negative feedback of Figure 8.24c illustrates the radiation balance
outgoing radiation is dominant in an Earth with- of the Earth as a function of global temperature
out extensive subpolar ice sheets. When large for the combination of long-wavelength outgoing
continental ice sheets expand in North Amer- radiation and the ice--albedo feedback. The ¬xed
ica and Northern Europe, however, temperature points are obtained by setting Eq. (8.43) equal to
drops result in an albedo increase large enough zero to obtain
to dominate the negative feedback of outgoing 1
T = c 2 (’1 ’ T ) 2 (8.44)
radiation. The net positive feedback driven by ice-
with c 2 = 3.0 —¦ C1/2 . This quadratic equation has
sheet expansion sends the climate system into a
two roots at T = ’1.1 and T = ’8.0 —¦ C. The root
self-enhancing feedback of colder temperatures
at T = ’8.0 —¦ C is a stable ¬xed point; cooling
and larger ice sheets. This feedback, however,
216 STOCHASTIC PROCESSES


below this point leads to a net warming and a re- 2.0
turn to the ¬xed point. The root at T = ’1.1 —¦ C is
an unstable ¬xed point. Equation (8.44) predicts
h
two stable states of the Late-Pleistocene climate
(km)
system: an interglacial climate with T = 0 —¦ C and
a full-glacial climate with T = ’8 —¦ C. This analy- 0.0
sis is independent of the system dynamics and
depends only on equilibria of Eq. (8.40).
The ¬‚exural isostatic de¬‚ection of the litho- slope = 0.001
sphere beneath an ice sheet may have at least
slope = 0.002
two effects on the geometry of the ice sheet.
First, lithospheric subsidence may lower accumu- 0 2000
1000
lation rates by decreasing the ice-sheet elevation. x (km)
This ˜˜load-accumulation™™ feedback is perhaps the
Fig 8.26 Steady-state ice-surface pro¬les of a threshold-
most well-studied effect of lithospheric de¬‚ection
sliding ice sheet with different scenarios of bedrock
on climate (North, 1991). The strength of this
topography, calculated from Eqs. (8.45) and (8.46). The
feedback, however, is limited by another effect
elevation of the ice surface and bedrock topography are
of lithospheric subsidence: changes in lateral ice-
shown for ¬‚at bedrock (solid line), and uniformly sloping
¬‚ow velocities. An ice sheet may be considered bedrock (dipping away from the margin) with a slope of
to be a conveyor belt in which ice ¬‚ows from ± = 0.002 (small-dashed line), and ± = 0.004 (long-dashed
zones of accumulation to zones of ablation. In line). As the bedrock depression deepens, the ice surface
steady-state, this ¬‚ow occurs without a change elevation decreases, but not as much as the bedrock. The
in surface topography or ice thickness. In fact, result is a counter¬‚ow within the ice sheet that causes
thickening or thinning of the ice sheet depending on whether
an increase in accumulation may be completely
the lithosphere is subsiding or rebounding. From Pelletier
accommodated by an increase in lateral ice-¬‚ow
(2003).
velocities such that no change in ice volume oc-
curs. This possibility highlights the dif¬culty of
modeling ice-sheet response to the small changes ing on its bed. On a ¬‚at bed, the pro¬le of a
in insolation experienced during the Pleistocene. threshold-sliding ice sheet is given by
In particular, a small change in basal-sliding pa-
h= 2h o (L ’ x) (8.45)
rameters may greatly change the sensitivity of ice
sheets to external forcing. Basal-sliding parame- and the thickness for an ice sheet on a sloping
ters are poorly known, introducing great uncer- base of gradient ± is
tainty into forward models of climate based on
h ho
ice-sheet dynamics. x’L = + 2 ln |±h ’ h o | (8.46)
± ±
Although the absolute values of basal-sliding
where L is the length of the ice sheet and h o is
parameters are poorly known, the relative effects
the maximum thickness of the ice sheet on a ¬‚at
of lithospheric de¬‚ection on ice-sheet sliding may
surface. h o is a function of the sliding shear stress
be determined more precisely. Lithospheric sub-
if we assume a sliding rather than a plastic ice
sidence perturbs a steady-state ice sheet by in-
sheet. To illustrate the effects of lithospheric de-
creasing the basal slope of the ice sheet, re-
¬‚ection on the ice-sheet thickness, we have plot-
ducing the gravitational force available to drive
ted ice-sheet pro¬les from Eqs. (8.45) and (8.46)
lateral ice ¬‚ow. This ¬‚ow reduction thickens the
for a ¬‚at surface and for surfaces with ± = 0.002
ice even under conditions of uniform mass bal-
and 0.004 in Figure 8.26. These pro¬les illustrate
ance. Conversely, lithospheric rebound enhances
that as the lithosphere subsides, an ice sheet
lateral ice-¬‚ow velocities and results in ice-sheet
will thicken because lateral ice ¬‚ow is inhibited
thinning. The effects of lithospheric subsidence
by the increase in basal slope. Conversely, litho-
and rebound may be illustrated by considering
spheric rebound leads to ice-sheet thinning.
the steady-state geometry of an ice sheet slid-
8.9 COHERENCE RESONANCE AND THE TIMING OF ICE AGES 217



glaciation deglaciation Fig 8.27 Schematic illustration of
the coupling between ice sheet
(d) flow-driven thinning only
(a) flow-driven thickening only
adjustment by lithospheric
de¬‚ection and ice-sheet advance or
retreat. The initial bedrock and
ice-sheet geometry is shown
(dashed line) along with the ¬nal
geometry (solid line) in each case.
v During glaciation: (a) Bedrock
v
subsidence considered alone leads
to ice-sheet thickening. (b)
viscoelastic rebound
viscoelastic subsidence Accumulation growth considered
alone leads to uniform thickening
(b) accumulation growth only (e) ablation decay only (assuming uniform accumulation
rate along the ice sheet) with a
corresponding advance. (c) Both
”agrowth ” adecay processes taken together act
destructively because some
accumulation growth is canceled out
by subsidence, reducing the effective
albedo change. During deglaciation:
(d) Bedrock rebound alone leads to
(c) thickening plus growth (f) thinning plus decay ice-sheet thinning. (e) Ablation alone
leads to uniform thinning of the ice
sheet with a corresponding retreat.
” anet ” anet (f) Together these processes act
more retreat constructively to modify albedo
less advance
than decay only
than growth only because ablation results in a greater
retreat for a thinner ice sheet.
v
Modi¬ed from Pelletier (2003).
v

Ice-sheet thickening and thinning, in turn, effects of subsidence-induced thickening and ac-
modify the effects of changes in global tempera- cumulation growth taken together are illustrated
ture to retard glaciations and accelerate deglacia- in Figure 8.27c. The combination results in a
tions. This coupling may be called the ˜˜load- smaller ice-sheet advance relative to the case
advance™™ feedback and is illustrated in Figure with accumulation alone. In cases in which ac-
8.27. Figures 8.27a--8.27c illustrate how a de- cumulation is not uniform (i.e. accumulation in-
crease in lateral ice-¬‚ow velocities results in a creases with distance from the margin) the lim-
reduced rate of ice-sheet advance. Conversely, Fig- iting effect of subsidence on ice-sheet advance
ures 8.27d--8.27f illustrate how ice-sheet retreat is is even greater. Figures 8.27d--8.27f illustrate the
enhanced by an increase in lateral ¬‚ow driven by load-advance feedback during deglaciation. Fig-
lithospheric rebound. In both cases, lithospheric ure 8.27d illustrates lithospheric rebound and re-
adjustment acts to increase the albedo relative sultant ice-sheet thinning. Figure 8.27e illustrates
to cases with a ¬xed lithosphere. Figure 8.27a the ice-sheet retreat resulting from ablation. Fi-
illustrates the effect of subsidence on ice-sheet nally, the combination of thinning and ablation
geometry: a reduction in lateral ice-¬‚ow veloci- are shown in Figure 8.27f. This ¬gure illustrates
ties results in ice-sheet thickening as shown in how ice-sheet thinning enhances retreat.
Figure 8.26. In Figure 8.27b the effects of ac- The effects of lithospheric subsidence and
cumulation growth are shown in the absence rebound on the motion of an ice margin are de-
of subsidence. Accumulation is assumed to oc- pendent on the vertical velocity of the lithosphere,
|v |, not on its elevation. This is a key difference
cur uniformly along the ice-sheet pro¬le. The
218 STOCHASTIC PROCESSES



Table 8.1 Values of the parameters in the reference model.

c 1 (s’1 ) c 2 (—¦ C 2 /s) c 4 (—¦ C/s)
1
„ (kyr)
c3

1.8 — 10’10 5.4 — 10’10
original N/A N/A 10
scaled to c 1 1 3.0 25.0 0.9 56


between the load-accumulation feedback and Eq. (8.37). The ¬rst two terms on the right side of
the load-advance feedback. In the load-advance (8.50) represent the bistability in the system. The
feedback, the retarding and amplifying effects of third and fourth terms represent the delay feed-
Figures 8.27c and 8.27f are observed only when back and random variability, respectively. The val-
ues of c 1 through c 4 and „ used in the reference
the lithosphere is actively subsiding or rebound-
ing. With the load-accumulation feedback, in model run of Eq. (8.50) are given in Table 8.1.
contrast, the reduction in accumulation rate The values of c 1 and c 2 have been determined
that accompanies subsidence is a function of the from empirical data in the preceding sections.
bedrock elevation, not its velocity. This is an im- The value of c 4 may be constrained using the
portant distinction, because asymmetric glacial- Vostok time series. To do this, we ¬rst averaged
interglacial transitions are only reproduced in a the Vostok time series in equal intervals of the
model with velocity-dependent feedbacks. model scaling time of 180 years. Second, we com-
To model the load-advance feedback, we in- puted the standard deviation of T n from the
Vostok data. The result is c 4 = 0.9 —¦ C and quan-
clude a term in the model equation (Eq. (8.40))
proportional to the magnitude of the velocity, |v |, ti¬es the magnitude of short-term global tem-
given by perature changes. This value is consistent with
other available estimates of natural century-scale
Bn
vn = (8.47) global-temperature variations. The Maxwell relax-
t
ation time of the mantle was chosen to be 10 kyr,
where a value consistent with postglacial rebound. The
n =n viscous response time of the mantle actually
T h ’ T ’(n’n )/„
Bn = e (8.48) varies with the size of the load so a constant

n =0
value will only be an approximation. The remain-
and „ is the Maxwell-relaxation time of the vis- ing model coef¬cient, c 3 , is the only free pa-
coelastic crust (Watts, 2001) and B is the bedrock rameter in the model. This coef¬cient character-
depth. izes the strength of the load-advance feedback,
The model equation, including long- which depends on the geometries and sliding
wavelength outgoing radiation, the ice--albedo rates of the ice sheets. We have not attempted
feedback, and lithospheric de¬‚ection, is to constrain this parameter with empirical data
because of the complexity of the processes in-
Tn 1
= ’c 1 T n + c 2 (T h ’ T n ) 2 + c 3 |v |n + c 4 ·n
volved. Instead, the model was run with a range
t
(8.49) of parameters to qualitatively identify the sys-
tem behavior. Coherence resonance occurs over
To simplify, we rescale Eq. (8.49) to make
a broad range of values for c 3 . c 3 = 25 was cho-
c 1 = 1:
sen as the reference parameter value for the nu-
Tn 1
merical model experiments because it leads to
= ’T n + c 2 (T h ’ T n ) 2 + c 3 |v |n + c 4 ·n
t both coherence resonance and behavior similar
(8.50)
to the Vostok time series, including asymmetric
transitions.
Model time steps are equal to the radiative-
The behavior of Eq. (8.50) is illustrated in Fig-
damping time of 180 yr in Eq. (8.50). Equation
ure 8.28b with the reference values of Table 8.1.
(8.50) is analogous to the coherence resonance
8.9 COHERENCE RESONANCE AND THE TIMING OF ICE AGES 219



2
2 (a) (b)
0
0
T
T
o
o
( C)
( C)




400
100 200 300
0
0
400 300 200 100
t (kyr)
age (ka)
10 (c)

|v|, B
5


0
400
100 200 300
0
t (kyr)
strength of the ice--albedo feedback in this early
Fig 8.28 (a) Time series of atmospheric temperature
stage. The glaciation slows as it proceeds, re¬‚ect-
variations inferred from deuterium concentrations at Vostok,
ing a decrease in the ice--albedo feedback that
Antarctica (Petit, 1999). (b) Time series of model
temperature variations with c 2 = 3.0, c 3 = 25.0, c 4 = 0.9, limits ice advance. Deglaciation is triggered by
and „ = 10 kyr. (c) Time series of bedrock elevation and a warming trend that is long enough to ini-
velocity in the model example of (b). From Pelletier (2003). tiate lithospheric rebound. The length of the
warming trend is important because lithospheric
motion lags behind ice loading. Therefore, a
Variations in bedrock depth B and the magni-
tude of bedrock velocity |v | are shown in Figure sustained warming trend is required to initiate
deglaciation through the combined effects of ice
8.28c for this run. The Vostok time series is plot-
sheet melting and a lithospheric rebound. Dur-
ted in Figure 8.28a for comparison. Note that the
ing deglaciation, the rebounding lithosphere acts
Vostok data is plotted here with time, not age,
to thin the ice sheet and amplify the ice--albedo
increasing toward the right in order to simplify
feedback in an additional feedback mechanism.
comparison with the model time series.
In this additional mechanism, faster ice-sheet
The model time series reproduces behavior
melting results in faster lithospheric rebound
very similar to the Vostok data, including the
and an enhanced ice--albedo feedback. The ice--
same temperature range (with full glacial cli-
mates 8 —¦ C cooler than interglacials), a 100-kyr albedo feedback, in turn, drives rapid ice-sheet
melting in a positive feedback. It should be noted
periodicity, and asymmetric glacial--interglacial
that even though the relaxation time of the man-
transitions. More speci¬cally, glaciations in both
tle is 10 kyr, the dominant oscillation is close to
the model and observed data are relatively rapid
100 kyr.
at ¬rst but slow prior to a very rapid deglacia-
The similarity between the model behavior
tion. This pattern may be understood in terms
and the Vostok time series may be documented
of the strength of radiation feedbacks and the
more thoroughly by comparing the histograms
variations in bedrock velocity through time.
and power spectra of both data sets. The power
Glaciations begin rapidly at ¬rst, re¬‚ecting the
220 STOCHASTIC PROCESSES



1.0
(a) (a)
10
S(f ) ∝ f ’2 0.8
f (T)
S( f )10
0.6
S( f ) = C
(oC2 kyr)
10
S(f ) ∝ f ’1/2 0.4

10
0.2
10
0.0
(b)
1.0
10’1 (b)
0.8
’2
S( f ) 10
f (T)
(oC2 kyr)
0.6
10’3

0.4
10’4

10’5 0.2
10’2 10’1
10’3 100 101
f (1/kyr)
0.0
0 2
Fig 8.29 (a) Lomb periodogram, S( f ), of Vostok time
T (oC)
series data plotted as a function of frequency f . Note
logarithmic scales on both axes. Spectrum has a dominant Fig 8.30 Histograms of (a) Vostok time series data,
frequency corresponding to a period near 100 kyr, a constant uniformly sampled at 200 yr intervals, and (b) model time
background power spectrum between time scales of 1 Ma and series. These curves re¬‚ect the fraction of time the climate is
40 kyr, an f ’2 spectrum between 40 kyr and 2 kyr, and an in the interglacial state, the glacial state, or ¬‚uctuating in
f ’1/2 spectrum at higher frequencies. (b) Average power between. The time spent near the glacial state exceeds the
spectrum of 1000 samples of the model time series, with time spent near the interglacial state by a ratio of about 5:1.
same scale as (a). The spectrum has a dominant periodicity at From Pelletier (2003).
120 kyr superimposed on a background spectrum similar to
(a) (except at the highest frequencies). In addition,
small-amplitude periodicities appear above the background at
superimposed on a Lorenztian spectrum. Notably,
odd harmonics of the dominant periodicity. From Pelletier
both spectra also exhibit a high-frequency transi-
(2003).
tion to f ’1/2 behavior. At the highest frequencies,
power in the model spectrum drops as an arti-
fact of the simulation method. Speci¬cally, the
spectra of the Vostok and model time series are viscoelastic load was averaged at time scales less
plotted in Figures 8.29a and 8.29b, respectively. than 1 kyr to speed up the code for long runs.
We estimated the power spectrum of the Vostok Histograms of the Vostok and model data are
time series using the Lomb periodogram (Press plotted in Figures 8.30a and 8.30b, respectively.
et al., 1992). This technique is commonly used Histograms of 1000 model samples have been av-
for nonuniformly sampled time-series data. The eraged to obtain Figure 8.30b. Both histograms
power spectrum of the model time series was es- are bimodal, re¬‚ecting the stable states of the sys-
timated by averaging the spectra of 1000 indepen- tem. The asymmetry of the model and observed
dent samples to obtain a precise spectrum. Both data are re¬‚ected in the larger peak near the full
spectra exhibit a dominant 100-kyr periodicity glacial climate.
EXERCISES 221



Exercises Such a model, however, does not incorporate the
effects of random spatial variations in porosity
8.1 Using the code in Appendix 7, generate a log-
within the aquifer. As an alternative to the diffu-
normal fractional Gaussian noise with β = 1/2, a
sion model, construct a 2D random walk model
mean value of 1 acre-ft/s, and a coef¬cient of vari-
of contaminant transport in aquifers. De¬ne the
ation of 1.0. Use this series as a model for the
central grid point to be the contaminant source.
discharge into a reservoir. Using Excel, model the
Simulate the variation in concentration from the
volume of the reservoir as a function of time, as-
source using random walkers. Assuming that the
suming that each year the reservoir loses 1 acre-ft/s
contaminant migrates in a random direction 10 m
to evaporation. In a simulation of 1000 yr, deter-
each year, what is the most likely time following
mine the duration of the longest drought (de¬ned
the leak that the contaminant will ¬rst arrive a dis-
as consecutive years below 50% of normal volume).
tance of 10 km from the source? Consider the worst
8.2 Develop a code that implements Masek and Tur-
possible scenario. What is the minimum time re-
cotte™s DLA model for drainage network evolution.
quired for the contaminant to travel 10 km from
However, instead of growing the network each time
the source? Imagine that the aquifer has a net
a random walker lands on a site whose neighbor
transport direction to the south. Model this effect
is part of the network, introduce a threshold con-
by making the probability that the walker will mi-
dition that requires that N walkers visit a site be-
grate south twice as large as the probability that it
fore the site becomes part of the network. Quali-
will migrate north. What is the relative contamina-
tatively compare the drainage networks produced
tion risk north and south of the source, expressed
with N = 10 with those produced with N = 1.
as the ratio of the most likely arrival times 10 km
8.3 The spreading of contaminants in underground
north and south of the contaminant source?
aquifers can be modeled as a diffusion process.
Appendix 1




Codes for solving the diffusion equation

A1.1 Preliminaries

The appendices of this book provide sample code for implementing many of the methods we have
described. The codes are written in the standard C language and make extensive use of Numerical
Recipes (Press et al., 1992) routines. In addition to Numerical Recipes, the codes make use of the
standard C libraries malloc (memory allocation, useful for constructing vectors and matrices), math
(standard math functions), stdio (standard input and output of data), and stdlib (miscellaneous
standard functions). In order to inform the C compiler that we will be using these functions, the
¬rst few lines of every program that uses code from this book should include:

#include<malloc.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include"nr.h"
#include"nrutil.h"

Most C compilers come standard with the ¬rst four libraries. The ¬nal two libraries, nr.h and nrutil.h
(Numerical Recipes and its utilities), should be copied into the local directory where the code will
be compiled or appropriately linked from another location.
The codes in this book have been successfully compiled using the UNIX cc compiler and the
Microsoft Visual C++ compiler. Every effort has been made to ensure that they are error free, but
no guarantees are made. Moreover, these codes will generally need to be modi¬ed to work with
other data sets, boundary conditions, etc. In addition, codes sometimes need to be modi¬ed on
different systems and as systems change over time. As such, code used from this book should
be thoroughly tested by the user for accuracy and internal consistency when applied to speci¬c
problems.
The codes in this book use float data types, which provide accuracy to six decimal places for
individual operations. For greater accuracy, data structures of type double should be used. Numerical
Recipes includes a version of each of their routines for double precision. For example, if a user is
implementing a code from the book that calls the qsimp function and he or she requires optimal
accuracy, all float variables should be changed to type double and the routine dqsimp should be
called instead of qsimp.
Output from these codes is usually sent to an ASCII text ¬le in two-line format: x h <return> (for
1D models) and as a series of rows and columns for 2D models (with one grid point on each line).
A1.1 PRELIMINARIES 223


In 2D cases, all grid output proceeds as rows ¬rst, columns second, beginning with the upper left
corner of the grid, moving from left to right and top to bottom. Many standard commercial and
freeware programs are available for visualizing 1D and 2D output. For 1D output, the results can be
imported from an ASCII text ¬le into Microsoft Excel (http://www.microsoft.com), for example. For
visualizing 2D grids, the RiverTools software package (http://www.rivix.com) provides a good com-
mercial option while 3DEM provides a good freeware option (http://www.visualizationsoftware.com).
Both RiverTools and 3DEM can import and export raster grids as ASCII text ¬les as well as a variety of
standard Geographic Information System (GIS) ¬le formats. Routines that create graphical output in
the PostScript language can also be written directly in the C language. Gershenfeld (1999) provides
a number of very useful visualization routines written in C.
The majority of methods in this book require that we refer to neighboring grid points when
performing grid operations. Many of the codes use a simple labeling scheme presented below. In
this scheme, the vectors iup[i] and idown[i] are de¬ned to be i+1 and i-1, respectively, for all
values of i except on the boundaries of the grid where i+1 and i-1 are not de¬ned. The values
of iup[i] and idown[i] at the boundaries depend on the type of boundary conditions. For peri-
odic boundary conditions, the value of iup[lattice size x] should be 1 and the value of idown[1]
should be lattice size x. For ¬xed boundary conditions, the value of iup[lattice size x] should be
lattice size x and the value of idown[1] should be 1. This labeling scheme has two advantages. First,
it enables the boundary conditions to be switched from periodic to ¬xed by changing just a couple
of lines of code. Second, it eliminates the need to include conditional statements. If, for example, we
used i+1 and i-1 to refer to grid neighbors, we would need to include conditional statements such
as

if ((i>1)&&(i<lattice_size_x)&&(j>1)&&(j<lattice_size_y))
h[i][j]+=timestep*D/(delta*delta)*
(hold[i+1][j]+hold[i-1][j]+hold[i][j+1]+hold[i][j-1]-4*hold[i][j]);

throughout the code in order to avoid a segmentation fault on the second line.
The following code is a standard routine for initializing grid-neighbor labels for a 2D grid with
¬xed boundary conditions. Many codes in this book use this function, but we will list it only once.
Some of the modules that appear in early appendices are used in later appendices. In subsequent
usage of setupgridneighbors,


int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void setupgridneighbors()
{ int i,j;

idown=ivector(1,lattice_size_x);
iup=ivector(1,lattice_size_x);
jup=ivector(1,lattice_size_y);
jdown=ivector(1,lattice_size_y);
for (i=1;i<=lattice_size_x;i++)
{idown[i]=i-1;
iup[i]=i+1;}
idown[1]=1;
iup[lattice_size_x]=lattice_size_x;
for (j=1;j<=lattice_size_y;j++)
{jdown[j]=j-1;
224 CODES FOR SOLVING THE DIFFUSION EQUATION


jup[j]=j+1;}
jdown[1]=1;
jup[lattice_size_y]=lattice_size_y;
}


If periodic boundary conditions are needed, setupgridneighbors can be modi¬ed as follows:


int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void setupgridneighborsperiodic()
{ int i,j;

idown=ivector(1,lattice_size_x);
iup=ivector(1,lattice_size_x);
jup=ivector(1,lattice_size_y);
jdown=ivector(1,lattice_size_y);
for (i=1;i<=lattice_size_x;i++)
{idown[i]=i-1;
iup[i]=i+1;}
idown[1]=lattice_size_x;
iup[lattice_size_x]=1;
for (j=1;j<=lattice_size_y;j++)
{jdown[j]=j-1;
jup[j]=j+1;}
jdown[1]=lattice_size_y;
jup[lattice_size_y]=1;
}

Another common trick we will use involves time stepping. In the case of the linear diffusion
equation with constant diffusivity, the stability criterion is t < ( x)2 /(2D ). In nonlinear and/or
coupled equations in which the diffusivity changes through time as the solution evolves, the stability
criterion is sometimes dif¬cult to determine. Stability can be enforced by varying the time step to
limit the maximum change in the variable or variables that are being modeled. This approach is
crude but it works. The downside to this approach is that it requires some trial and error in order
to determine how big a step is allowed before the solution becomes unstable. In the example code
below, maxchange is set equal to the largest change allowed in any value of topo[i][j] during a single
time step (determined by trial and error). If the largest change exceeds that maximum value, the
values of topo revert to the values at the last time step topoold and the time step is reduced by a
factor of 2. If the largest change is smaller than one tenth of the maximum value then the time step
is increased by a factor of 1.2. In this way, the time step is increased and decreased dynamically in
order to maintain similar maximum changes in the function or functions we are modeling. We will
use this scheme in a number of example codes in this book.
max=0;
for (i=1;i<lattice_size_x;i++)
{deltah=timestep/deltax*(flux[i]-flux[iup[i]]);
topo[i]+=deltah;
if (fabs(deltah)>max) max=fabs(deltah);}
elapsedtime+=timestep;
A1.2 FTCS METHOD FOR 2D DIFFUSION EQUATION 225


if (max>maxchange)
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<=lattice_size_x;i++)
topo[i]=topoold[i];}
else
{if (max<0.1*maxchange) timestep*=1.2;}



A1.2 FTCS method for 2D diffusion equation

The sample code below solved the 2D diffusion equation using the FTCS method. The code accepts an
input grid of type float, size 300 — 300, and resolution 10 m/pixel from the local directory, calculates
the solution for t = duration, and then outputs the result to an ASCII ¬le in the local directory.


int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

main()
{ FILE *fp1,*fp2;
float delta,**topo,**topoold,D,duration,timestep;
int i,j,t,nsteps;

fp1=fopen("inputdem","r");
fp2=fopen("outputdem","w");
lattice_size_x=300;
lattice_size_y=300;
delta=10.0; /* m */
D=1.0; /* m^2/kyr */
duration=100.0; /* kyr */
setupgridneighbors();
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f",&topo[i][j]);
topoold[i][j]=topo[i][j];}
timestep=0.5*delta*delta/(2*D);
nsteps=(int)(duration/timestep);
for (t=1;t<=nsteps;t++)
{for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topo[i][j]+=timestep*D/(delta*delta)*(topoold[iup[i]][j]+topoold[idown[i]][j]
+topoold[i][jup[j]]+topoold[i][jdown[j]]-4*topoold[i][j]);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
226 CODES FOR SOLVING THE DIFFUSION EQUATION


fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}

As a concrete example, we consider the case of terrace diffusion described in Section 2.3.4. The
following code can be used to reproduce Figure 2.23b, given the input ¬le rillmask, available from
the CUP website. This code also makes use of an ˜˜ef¬ciency mask™™ that identi¬es the grid points and
their neighbors that are actively undergoing erosion. This speeds up the code early in the simulation
when erosion is localized close to the rills.


#define PI 3.141592653589793

main()
{ FILE *fp1,*fp2;
int length_x,length_y,delta,lattice_size_x,lattice_size_y,**rills,**efficiencymask,i,j;
float deltah,simulationlength,elapsedtime,timestep,dissectionrelief,kappa,depositslope;
float depositaspect,northwestdatum,**topo,**topoold,**initialtopo;

fp1=fopen("rillmask","r");
fp2=fopen("terracediffusion2d","w");
length_x=300; /* m */
length_y=300;
delta=1.0; /* m */
dissectionrelief=10.0;
kappa=1.0; /* m^2/kyr */
depositslope=0.06; /* m/m */
depositaspect=300*PI/180; /* surface drains SE; degrees converted to radians */
northwestdatum=0.0; /* elevations are relative to NE corner = 0 m */
lattice_size_x=length_x/delta;
lattice_size_y=length_y/delta;
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
initialtopo=matrix(1,lattice_size_x,1,lattice_size_y);
rills=imatrix(1,lattice_size_x,1,lattice_size_y);
efficiencymask=imatrix(1,lattice_size_x,1,lattice_size_y);
timestep=0.5*delta*delta/(2*kappa);
simulationlength=50; /* kyr */
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%d",&rills[i][j]);
if (rills[i][j]>0) rills[i][j]=0; else rills[i][j]=1;
efficiencymask[i][j]=1;
if (i<lattice_size_x) if (rills[i+1][j]==1) efficiencymask[i][j]=0;
if (i>1) if (rills[i-1][j]==1) efficiencymask[i][j]=0;
if (j<lattice_size_y) if (rills[i][j+1]==1) efficiencymask[i][j]=0;
if (j>1) if (rills[i][j-1]==1) efficiencymask[i][j]=0;
topo[i][j]=northwestdatum-(i-1)*delta*depositslope*cos(depositaspect)+
(j-1)*delta*depositslope*sin(depositaspect);
A1.3 FTCS METHOD FOR 1D NONLINEAR DIFFUSION EQUATION 227


if (rills[i][j]==1) topo[i][j]-=dissectionrelief;
topoold[i][j]=topo[i][j];
initialtopo[i][j]=topo[i][j];}
elapsedtime=0.0;
while (elapsedtime<=simulationlength)
{printf("%f %f\n",elapsedtime,timestep);
for (j=2;j<=lattice_size_y-1;j++)
for (i=2;i<=lattice_size_x-1;i++)
if (efficiencymask[i][j]==0)
{if (rills[i][j]==0)
deltah=timestep*kappa/(delta*delta)*(topoold[i+1][j]+topoold[i-1][j]+
topoold[i][j+1]+topoold[i][j-1]-4*topoold[i][j]);
else
deltah=0;
topo[i][j]+=deltah;
if (fabs(topo[i][j]-initialtopo[i][j])>0.01*dissectionrelief)
{efficiencymask[i+1][j]=0;
efficiencymask[i-1][j]=0;
efficiencymask[i][j+1]=0;
efficiencymask[i][j-1]=0;}}
elapsedtime+=timestep;
i=1;
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=topo[i+1][j]+depositslope*cos(depositaspect)*delta;
i=lattice_size_x;
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=topo[i-1][j]-depositslope*cos(depositaspect)*delta;
j=1;
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=topo[i][j+1]-depositslope*sin(depositaspect)*delta;
j=lattice_size_y;
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=topo[i][j-1]+depositslope*sin(depositaspect)*delta;
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topoold[i][j]=topo[i][j];}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}



A1.3 FTCS method for 1D nonlinear diffusion equation

The code presented below is used to solve for the evolution of a 1D hillslope according to the
nonlinear diffusion equation (described in Section 2.3.3) subject to instantaneous base level drop.
228 CODES FOR SOLVING THE DIFFUSION EQUATION


main()
{ FILE *fp1;
int *iup,*idown,printed,outputinterval,lattice_size_x,i;
float kappa,Sc,slope,den,*flux,duration,deltax,deltah,elapsedtime,max,timestep;
float *topo,*topoold,maxchange;

fp1=fopen("nonlinbasedrop","w");
lattice_size_x=102;
deltax=1.0; /* m */
kappa=1.0; /* m^2/kyr */
Sc=1.0; /* m/m */
iup=ivector(1,lattice_size_x);
idown=ivector(1,lattice_size_x);
topo=vector(1,lattice_size_x);
topoold=vector(1,lattice_size_x);
flux=vector(1,lattice_size_x);
duration=1000.0; /* kyr */
outputinterval=1000.0;
maxchange=0.001;
for (i=1;i<=lattice_size_x;i++)
{topo[i]=1.0;
if ((i==lattice_size_x)||(i==lattice_size_x-1)||(i==lattice_size_x-2)) topo[i]=0.0;
topoold[i]=topo[i];
flux[i]=0;
iup[i]=i+1;idown[i]=i-1;}
iup[lattice_size_x]=lattice_size_x;idown[1]=1;
timestep=0.0001*deltax*deltax;
elapsedtime=0.0;printed=0;
while (elapsedtime<duration)
{printed=elapsedtime/outputinterval;
for (i=1;i<=lattice_size_x;i++)
{slope=(topoold[idown[i]]-topoold[i])/deltax;
den=1-fabs(slope)*fabs(slope)/(Sc*Sc);
if (den<0.01) den=0.01;
flux[i]=kappa*slope/den;
if (i==1) flux[i]=0.0;
if (i==lattice_size_x) flux[i]=flux[idown[i]];
if (flux[i]<flux[idown[i]]) flux[i]=flux[idown[i]];}
max=0;
for (i=1;i<lattice_size_x;i++)
{deltah=timestep/deltax*(flux[i]-flux[iup[i]]);
topo[i]+=deltah;
if (fabs(deltah)>max) max=fabs(deltah);}
elapsedtime+=timestep;
if (max>maxchange)
{elapsedtime-=timestep;
timestep/=2.0;
for (i=1;i<lattice_size_x;i++)
topo[i]=topoold[i];}
A1.4 ADI METHOD FOR SOLVING THE 2D DIFFUSION EQUATION 229


else
{if (max<0.1*maxchange) timestep*=1.2;}
for (i=1;i<lattice_size_x;i++)
topoold[i]=topo[i];
if ((int)(elapsedtime/outputinterval)>printed)
{printf("%f %f\n",elapsedtime,timestep);
for (i=1;i<=lattice_size_x;i++)
fprintf(fp1,"%d %f\n",i,topo[i]);
printed=(int)(elapsedtime/outputinterval);}}
fclose(fp1);
}


A1.4 ADI method for solving the 2D diffusion equation

The following code accepts a 30 m/pixel US Geological Survey DEM covering Hanaupah Canyon from
¬le hanaupahtopo and uses the ADI technique to solve the diffusion equation on hillslopes. In the
hillslopediffusion subroutine, the diffusion equation is solved for each row and each column in
an alternating fashion using the tridag routine of Numerical Recipes, which inverts a tridiagonal
matrix. Hillslopes are de¬ned to be pixels with a contributing area less than 0.1 km2 . Channel pixels
act as ¬xed boundary conditions for the channel pixels. In this example, 10 Myr of diffusive evo-
lution is simulated with a diffusivity value of D = 10 m2 /kyr using time steps of 1 Myr. The ¬le
hanaupahtopo can be downloaded from the CUP website. The code below makes use of two functions:
fillinpitsandflats and mfdflowroute introduced in Appendix 2.

float D,delta,thresholdarea,*ax,*bx,*cx,*ux,*rx,*ay,*by,*cy,*uy,*ry;
float **topoold,*topovec,**topo,**flow,**flow1,**flow2,**flow3,**flow4,**flow5;
float **flow6,**flow7,**flow8;
int *topovecind,lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void hillslopediffusion()
{
int i,j,count;
float term1;

count=0;
while (count<5)
{count++;
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
topoold[i][j]=topo[i][j];
for (i=1;i<=lattice_size_x;i++)
{for (j=1;j<=lattice_size_y;j++)
{term1=D*timestep/(delta*delta);
if (flow[i][j]<thresholdarea)
{ay[j]=-term1;
cy[j]=-term1;
by[j]=4*term1+1;
ry[j]=term1*(topo[iup[i]][j]+topo[idown[i]][j])+topoold[i][j];}
230 CODES FOR SOLVING THE DIFFUSION EQUATION


else
{by[j]=1;
ay[j]=0;
cy[j]=0;
ry[j]=topoold[i][j];}
if (j==1)
{by[j]=1;
cy[j]=0;
ry[j]=topoold[i][j];}
if (j==lattice_size_y)
{by[j]=1;
ay[j]=0;
ry[j]=topoold[i][j];}}
tridag(ay,by,cy,ry,uy,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
topo[i][j]=uy[j];}
for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)
topoold[i][j]=topo[i][j];
for (j=1;j<=lattice_size_y;j++)
{for (i=1;i<=lattice_size_x;i++)
{term1=D*timestep/(delta*delta);
if (flow[i][j]<thresholdarea)
{ax[i]=-term1;
cx[i]=-term1;
bx[i]=4*term1+1;
rx[i]=term1*(topo[i][jup[j]]+topo[i][jdown[j]])+topoold[i][j];}
else
{bx[i]=1;
ax[i]=0;
cx[i]=0;
rx[i]=topoold[i][j];}
if (i==1)
{bx[i]=1;
cx[i]=0;
rx[i]=topoold[i][j];}
if (i==lattice_size_x)
{bx[i]=1;
ax[i]=0;
rx[i]=topoold[i][j];}}
tridag(ax,bx,cx,rx,ux,lattice_size_x);
for (i=1;i<=lattice_size_x;i++)
topo[i][j]=ux[i];}}
}

main()
{ FILE *fp1,*fp2;
float time,timestep,duration;
int i,j,t,dum;
A1.4 ADI METHOD FOR SOLVING THE 2D DIFFUSION EQUATION 231


fp1=fopen("hanaupahtopo","r");
fp2=fopen("hanaupahtopodiffuse","w");
lattice_size_x=640;
lattice_size_y=335;
delta=30.0; /* m */
D=10.0; /* m^2/kyr */
thresholdarea=0.1; /* km */
duration=10000.0; /* kyr */
timestep=1000.0; /* kyr */
setupgridneighbors();
/* grids needed for matrix inversion */
ax=vector(1,lattice_size_x);
ay=vector(1,lattice_size_y);
bx=vector(1,lattice_size_x);
by=vector(1,lattice_size_y);
cx=vector(1,lattice_size_x);
cy=vector(1,lattice_size_y);
ux=vector(1,lattice_size_x);
uy=vector(1,lattice_size_y);
rx=vector(1,lattice_size_x);
ry=vector(1,lattice_size_y);
/* other grids */
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
topoold=matrix(1,lattice_size_x,1,lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{flow[i][j]=(delta/1000)*(delta/1000); /* km */
fscanf(fp1,"%d",&dum);
topo[i][j]=dum;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
232 CODES FOR SOLVING THE DIFFUSION EQUATION


{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
time=0;
while (time<duration)
{time+=timestep;
hillslopediffusion();}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}


A1.5 Numerical integration of Fourier“Bessel terms

The following code accepts a list of the ¬rst 300 zeros of the Bessel function J 0 (r ) from the local
directory and performs numerical integration to calculate the value of the terms in Eq. (2.57) for the
case of a volcanic cone with r f = 200 m, rr = 500 m, and rc = 1500 m.

int nsum;
float *zeros;

float func(x)
float x;
{
return x*x*bessj0(zeros[nsum]*x);
}

main()
{ FILE *fp0,*fp1;
double a,zero=0.0,s1,s2,s3,rf,rr,rc;

fp0=fopen("besselzeros","r");
fp1=fopen("besselzerosintegral","w");
zeros=vector(1,300);
a=3000.0;
rf=200.0/a;
rr=500.0/a;
rc=1500.0/a;
for (nsum=1;nsum<=300;nsum++)
{fscanf(fp0,"%f",&zeros[nsum]);
s1=qsimp(func,zero,rf);
s2=qsimp(func,zero,rr);
s3=qsimp(func,zero,rc);
A1.5 NUMERICAL INTEGRATION OF FOURIER“BESSEL TERMS 233


fprintf(fp1,"%f %12.10f %12.10f %12.10f\n",zeros[nsum],s1,s2,s3);}
fclose(fp1);
fclose(fp2);
}

The following code accepts output from the code above, computes the results of the Bessel series
in Eq. (2.57), and prints the results to an output ¬le.
main()
{
FILE *fp1,*fp2;
float rf,rr,rc,*f,kappaT,deltax,a,*prefact,*topo,*J0rr,*J0rc,*J1rr,*J1rc;
float *J0rf,*J1rf,*Intrf,*J0r0,*J1r0,*Intrr,*Intrc;
int lattice_size_x,i,nsum;

fp1=fopen("besselzerosintegral","r");
fp2=fopen("volcanicconeoutput","w");
kappaT=0.0001; /* dimensionless scaled to a^2 */
lattice_size_x=300;
deltax=10.0;
a=lattice_size_x*deltax;
rf=200.0/a;
rr=500.0/a;
rc=1500.0/a;
topo=vector(1,lattice_size_x);
f=vector(1,lattice_size_x);
zeros=vector(1,300);
J0rc=vector(1,300);
J1rc=vector(1,300);
Intrc=vector(1,300);
J0rr=vector(1,300);
J1rr=vector(1,300);
Intrr=vector(1,300);
J0rf=vector(1,300);
J1rf=vector(1,300);
Intrf=vector(1,300);
prefact=vector(1,300);
for (nsum=1;nsum<=300;nsum++)
{fscanf(fp1,"%f %f %f %f",&zeros[nsum],&Intrf[nsum],&Intrr[nsum],&Intrc[nsum]);
J0rc[nsum]=bessj0(rc*zeros[nsum]);
J1rc[nsum]=bessj1(rc*zeros[nsum]);
J0rr[nsum]=bessj0(rr*zeros[nsum]);
J1rr[nsum]=bessj1(rr*zeros[nsum]);
J0rf[nsum]=bessj0(rf*zeros[nsum]);
J1rf[nsum]=bessj1(rf*zeros[nsum]);}
for (i=1;i<=lattice_size_x;i++)
{if (i*deltax/a<rf)
f[i]=1+(rf-rr)/(rc-rr);
else if (i*deltax/a<rr)
f[i]=1+(i*deltax/a-rr)/(rc-rr);
234 CODES FOR SOLVING THE DIFFUSION EQUATION


else if (i*deltax/a<rc)
f[i]=1+(rr-i*deltax/a)/(rc-rr);
else f[i]=0;
/* print initial condition first */
fprintf(fp2,"%f %12.10f\n",i*deltax,f[i]);}
for (i=1;i<=lattice_size_x;i++)
topo[i]=0;
for (nsum=1;nsum<=300;nsum++)
prefact[nsum]=2.0/zeros[nsum]*exp(-kappaT*zeros[nsum]*zeros[nsum])/
(bessj1(zeros[nsum])*bessj1(zeros[nsum]));
for (i=1;i<=lattice_size_x;i++)
{for (nsum=1;nsum<=300;nsum++)
topo[i]+=prefact[nsum]*bessj0(zeros[nsum]*i*deltax/a)*(rf*rf/(rc-rr)*
J1rf[nsum]-2*rr*rr/(rc-rr)*J1rr[nsum]+(1+rr/(rc-rr))*
rc*J1rc[nsum]+2*zeros[nsum]/(rc-rr)*Intrr[nsum]-zeros[nsum]/
(rc-rr)*(Intrf[nsum]+Intrc[nsum]));
fprintf(fp2,"%f %12.10f\n",i*deltax,topo[i]);}
fclose(fp1);
fclose(fp2);
}
Appendix 2




Codes for ¬‚ow routing


A2.1 Filling in pits and ¬‚ats in a DEM

The following code is used to eliminate pits and ¬‚ats from a DEM. The output from this type of
processing is sometimes called a hydrologically corrected DEM. The code is currently constructed to
accept a DEM from the local directory named inputdem that has 300 columns and 300 rows. In this
routine, all of the grid points are tested to determine whether they have a downslope neighbor. If
they do not have a downslope neighbor then they must be a ¬‚at or the bottom of a pit. In such
cases, the local elevation value is incremented by a small value fillincrement and the neighbors
are searched, recursively, for additional pits and ¬‚ats. Any additional pits and ¬‚ats that are found
also have their elevation values incremented. This procedure stops when all of the pixels in the pit
or ¬‚at have a downslope neighbor (i.e. the areas drains). For very large pits or ¬‚ats, it is possible
for this routine to run out of memory by calling itself too many times. In such cases, the value
of fillincrement can be raised, and/or a counter can be set up fillinpitsandflats that limits the
number of times that the routine can call itself before returning to the main program.
#define fillincrement 0.01

float **topo,
int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void fillinpitsandflats(i,j)
int i,j;
{ float min;

min=topo[i][j];
if (topo[iup[i]][j]<min) min=topo[iup[i]][j];
if (topo[idown[i]][j]<min) min=topo[idown[i]][j];
if (topo[i][jup[j]]<min) min=topo[i][jup[j]];
if (topo[i][jdown[j]]<min) min=topo[i][jdown[j]];
if (topo[iup[i]][jup[j]]<min) min=topo[iup[i]][jup[j]];
if (topo[idown[i]][jup[j]]<min) min=topo[idown[i]][jup[j]];
if (topo[idown[i]][jdown[j]]<min) min=topo[idown[i]][jdown[j]];
if (topo[iup[i]][jdown[j]]<min) min=topo[iup[i]][jdown[j]];
if ((topo[i][j]<=min)&&(i>1)&&(j>1)&&(i<lattice_size_x)&&(j<lattice_size_y))
236 CODES FOR FLOW ROUTING


{topo[i][j]=min+fillincrement;
fillinpitsandflats(i,j);
fillinpitsandflats(iup[i],j);
fillinpitsandflats(idown[i],j);
fillinpitsandflats(i,jup[j]);
fillinpitsandflats(i,jdown[j]);
fillinpitsandflats(iup[i],jup[j]);
fillinpitsandflats(idown[i],jup[j]);
fillinpitsandflats(idown[i],jdown[j]);
fillinpitsandflats(iup[i],jdown[j]);}
}

main ()
{ FILE *fp1,*fp2;
int i,j;

fp1=fopen("inputdem","r");
fp2=fopen("filleddem","w");
setupgridneighbors();
lattice_size_x=300;
lattice_size_y=300;
topo=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fscanf(fp1,"%f ",&topo[i][j]);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",topo[i][j]);
fclose(fp1);
fclose(fp2);
}




A2.2 MFD ¬‚ow routing method

The following code calculates the contributing area of a 30 m/pixel DEM using MFD ¬‚ow routing
assuming a local input ¬le named inputdem.


#define oneoversqrt2 0.707106781187

float **topo,**flow,**flow1,**flow2,**flow3,**flow4,**flow5,**flow6,**flow7,**flow8;
int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;
A2.2 MFD FLOW ROUTING METHOD 237


void mfdflowroute(i,j)
int i,j;
{ float tot;

tot=0;
if (topo[i][j]>topo[iup[i]][j])
tot+=pow(topo[i][j]-topo[iup[i]][j],1.1);
if (topo[i][j]>topo[idown[i]][j])
tot+=pow(topo[i][j]-topo[idown[i]][j],1.1);
if (topo[i][j]>topo[i][jup[j]])
tot+=pow(topo[i][j]-topo[i][jup[j]],1.1);
if (topo[i][j]>topo[i][jdown[j]])
tot+=pow(topo[i][j]-topo[i][jdown[j]],1.1);
if (topo[i][j]>topo[iup[i]][jup[j]])
tot+=pow((topo[i][j]-topo[iup[i]][jup[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[iup[i]][jdown[j]])
tot+=pow((topo[i][j]-topo[iup[i]][jdown[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[idown[i]][jup[j]])
tot+=pow((topo[i][j]-topo[idown[i]][jup[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[idown[i]][jdown[j]])
tot+=pow((topo[i][j]-topo[idown[i]][jdown[j]])*oneoversqrt2,1.1);
if (topo[i][j]>topo[iup[i]][j])
flow1[i][j]=pow(topo[i][j]-topo[iup[i]][j],1.1)/tot;
else flow1[i][j]=0;
if (topo[i][j]>topo[idown[i]][j])
flow2[i][j]=pow(topo[i][j]-topo[idown[i]][j],1.1)/tot;
else flow2[i][j]=0;
if (topo[i][j]>topo[i][jup[j]])
flow3[i][j]=pow(topo[i][j]-topo[i][jup[j]],1.1)/tot;
else flow3[i][j]=0;
if (topo[i][j]>topo[i][jdown[j]])
flow4[i][j]=pow(topo[i][j]-topo[i][jdown[j]],1.1)/tot;
else flow4[i][j]=0;
if (topo[i][j]>topo[iup[i]][jup[j]])
flow5[i][j]=pow((topo[i][j]-topo[iup[i]][jup[j]])*oneoversqrt2,1.1)/tot;
else flow5[i][j]=0;
if (topo[i][j]>topo[iup[i]][jdown[j]])
flow6[i][j]=pow((topo[i][j]-topo[iup[i]][jdown[j]])*oneoversqrt2,1.1)/tot;
else flow6[i][j]=0;
if (topo[i][j]>topo[idown[i]][jup[j]])
flow7[i][j]=pow((topo[i][j]-topo[idown[i]][jup[j]])*oneoversqrt2,1.1)/tot;
else flow7[i][j]=0;
if (topo[i][j]>topo[idown[i]][jdown[j]])
flow8[i][j]=pow((topo[i][j]-topo[idown[i]][jdown[j]])*oneoversqrt2,1.1)/tot;
else flow8[i][j]=0;
flow[iup[i]][j]+=flow[i][j]*flow1[i][j];
flow[idown[i]][j]+=flow[i][j]*flow2[i][j];
flow[i][jup[j]]+=flow[i][j]*flow3[i][j];
flow[i][jdown[j]]+=flow[i][j]*flow4[i][j];
238 CODES FOR FLOW ROUTING


flow[iup[i]][jup[j]]+=flow[i][j]*flow5[i][j];
flow[iup[i]][jdown[j]]+=flow[i][j]*flow6[i][j];
flow[idown[i]][jup[j]]+=flow[i][j]*flow7[i][j];
flow[idown[i]][jdown[j]]+=flow[i][j]*flow8[i][j];
}

main ()
{ FILE *fp1,*fp2;
int i,j,t,*topovecind;
float *topovec,delta;

fp1=fopen("inputdem","r");
fp2=fopen("MFDcontribarea","w");
setupgridneighbors();
lattice_size_x=300;
lattice_size_y=300;
delta=30.0; /* m */
topo=matrix(1,lattice_size_x,1,lattice_size_y);
topovec=vector(1,lattice_size_x*lattice_size_y);
topovecind=ivector(1,lattice_size_x*lattice_size_y);
flow=matrix(1,lattice_size_x,1,lattice_size_y);
flow1=matrix(1,lattice_size_x,1,lattice_size_y);
flow2=matrix(1,lattice_size_x,1,lattice_size_y);
flow3=matrix(1,lattice_size_x,1,lattice_size_y);
flow4=matrix(1,lattice_size_x,1,lattice_size_y);
flow5=matrix(1,lattice_size_x,1,lattice_size_y);
flow6=matrix(1,lattice_size_x,1,lattice_size_y);
flow7=matrix(1,lattice_size_x,1,lattice_size_y);
flow8=matrix(1,lattice_size_x,1,lattice_size_y);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
{fscanf(fp1,"%f ",&topo[i][j]);
flow[i][j]=delta*delta;}
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
fillinpitsandflats(i,j);
for (j=1;j<=lattice_size_y;j++)
for (i=1;i<=lattice_size_x;i++)
topovec[(j-1)*lattice_size_x+i]=topo[i][j];
indexx(lattice_size_x*lattice_size_y,topovec,topovecind);
t=lattice_size_x*lattice_size_y+1;
while (t>1)
{t--;
i=(topovecind[t])%lattice_size_x;
if (i==0) i=lattice_size_x;
j=(topovecind[t])/lattice_size_x+1;
if (i==lattice_size_x) j--;
mfdflowroute(i,j);}
for (j=1;j<=lattice_size_y;j++)
A2.3 SUCCESSIVE FLOW ROUTING WITH MFD METHOD 239


for (i=1;i<=lattice_size_x;i++)
fprintf(fp2,"%f\n",flow[i][j]);
fclose(fp1);
fclose(fp2);
}


A2.3 Successive ¬‚ow routing with MFD method

The following code performs a successive ¬‚ow routing analysis on an input DEM. The user is prompted
for the names of the input and output DEM ¬les, the size of the DEM, and the prescribed runoff
depth. The routine uses Manning™s equation to convert a prescribed discharge into an equivalent
depth, assuming steady ¬‚ow. As such, the local ¬‚ow depths depend sensitively on the local bed
slope. For some applications of this code it may be necessary to ¬rst smooth the input DEM in order
to remove topographic irregularities. In the example below it is assumed that the smallest bed slope
was 1% (any smaller value is made equal to 1%).

#define niterations 100
#define oneoversqrt2 0.707106781186

float deltax,**topo,**toposave,*topovec,**depth,**slope;
float **flow,**flow1,**flow2,**flow3,**flow4,**flow5,**flow6,**flow7,**flow8;
int lattice_size_x,lattice_size_y,*iup,*idown,*jup,*jdown;

void calculatealongchannelslope(i,j)
int i,j;
{ float down;

down=0;
if (topo[iup[i]][j]-topo[i][j]<down) down=topo[iup[i]][j]-topo[i][j];
if (topo[idown[i]][j]-topo[i][j]<down) down=topo[idown[i]][j]-topo[i][j];
if (topo[i][jup[j]]-topo[i][j]<down) down=topo[i][jup[j]]-topo[i][j];
if (topo[i][jdown[j]]-topo[i][j]<down) down=topo[i][jdown[j]]-topo[i][j];
if ((topo[iup[i]][jup[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[iup[i]][jup[j]]-topo[i][j])*oneoversqrt2;
if ((topo[idown[i]][jup[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[idown[i]][jup[j]]-topo[i][j])*oneoversqrt2;
if ((topo[iup[i]][jdown[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[iup[i]][jdown[j]]-topo[i][j])*oneoversqrt2;
if ((topo[idown[i]][jdown[j]]-topo[i][j])*oneoversqrt2<down)
down=(topo[idown[i]][jdown[j]]-topo[i][j])*oneoversqrt2;
slope[i][j]=fabs(down)/deltax;
}

main()
{ FILE *fp1,*fp2;
float runoff,manningsn,*topovec;
char inputfile[100],outputfile[100];
int i,j,k,t,*topovecind;
240 CODES FOR FLOW ROUTING


scanf("%s",&inputfile);
scanf("%s",&outputfile);
scanf("%d",&lattice_size_x);
scanf("%d",&lattice_size_y);
scanf("%f",&deltax); /* m */
scanf("%f",&runoff); /* m */
scanf("%f",&manningsn);
fp1=fopen(inputfile,"r");
fp2=fopen(outputfile,"w");

<<

. 12
( 14)



>>