<< . .

. 42
( : 51)

. . >>

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

Table 8.1 Values of the parameters in the reference model.

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

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
vn = (8.47) global-temperature variations. The Maxwell relax-
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
(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
to the Vostok time series, including asymmetric
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

2 (a) (b)
( C)
( C)

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

|v|, B

100 200 300
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

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

10’1 (b)
S( f ) 10
f (T)
(oC2 kyr)


10’5 0.2
10’2 10’1
10’3 100 101
f (1/kyr)
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-
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 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:


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

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

if ((i>1)&&(i<lattice_size_x)&&(j>1)&&(j<lattice_size_y))

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;

for (i=1;i<=lattice_size_x;i++)
for (j=1;j<=lattice_size_y;j++)


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;

for (i=1;i<=lattice_size_x;i++)

<< . .

. 42
( : 51)

. . >>