Odekon M. Encyclopedia of paleoclimatology and ancient environments

Подождите немного. Документ загружается.

Hysteresis behavior of the MOC

Sea water density is a function of not only temperature, but

also salinity. Furthermore, at low temperatures, such as those

that exist in the deep water formation regions, the density of

sea water is more sensitive to changes in salinity than to

changes in temperature. In addition, throughout the Quaternary,

large quantities of freshwater have been periodically stored

on land at middle to high latitudes in the form of continental

ice sheets. During the growth or melt of these continental ice

sheets, tens of meters of sea level equivalent freshwater have

been taken from or released back to the ocean. As a conse-

quence, numerous studies have been conducted in an attempt

to understand the role of freshwater perturbations on the stabi-

lity of the MOC.

The dependence of sea water density on salinity is a key

player in what is now known as a hysteresis behavior of the

AMO. Recently it has been discovered that the AMO behaves like

a “see-saw,” with relatively warm surface conditions in the

North Atlantic and relatively cool surface conditions in the South

Atlantic occurring when NADW formation is active, and conver-

sely, with relatively cool surface conditions in the North Atlantic

and relatively warm surface conditions in the South Atlantic

occurring when NADW formation is inactive (Figure T2). Such

a finding is supported by both modeling studies, proxy records,

and from combined model/proxy record studies. More recently

still, it has been realized that the “see-saw” nature of the AMO is

fundamentally linked to a coupling between the strength of AAIW

formation and the strength of NADW formation (Figure T3).

Figure T2 Atlantic meridional overturning circulation (Sv) for: (a) the stable active mode of NADW formation; (b) the stable inactive mode of

NADW formation. The inactive mode (b) was obtained by adding freshwater (100 Sv years) to region A. Notice the enhanced intrusion of AAIW

and the absence of NADW in (b). (c) Surface air temperature (



C) difference between the climatic states represented by (a) and (b). Region B in

(c) is where freshwater is added to get the hysteresis of Figure T3, although similar results are found when freshwater is added to region

C. Taken from Weaver et al. (2003).

THERMOHALINE CIRCULATION 945

The meridional overturning circulation and

Quaternary climate change

Early ice core records for the last glaciation have revealed large

amplitude variability on the millennial timescale characterized

by abrupt warming events (interstadials) lasting from several

hundred to several thousand years. These oscillations (Figure T4

and T5), known as Dansgaard-Oeschger (D-O) oscillations,

also appear in North Atlantic sediment records suggesting a

role or response of the ocean. Evidence from the Santa Barbara

basin and the Northeast Pacific suggests that a signature of

these D-O oscillations is also present in the Pacific Ocean,

while further recent sediment analyses suggest they may be

an inherent part of late and early Pleistocene climate.

In an attempt to provide a mechanism for the observed D-O

variability, it was initially proposed that a stable AMO was not

possible during glacial times, when the northern end of the

Atlantic Ocean was surrounded by ice sheets. Furthermore,

when the AMO was weakened or shut down and ice sheets

were growing, there was little oceanic salt export from the

Atlantic to the other ocean basins. Assuming a net evaporation

over the North Atlantic, its salinity continued to increase as

moisture was deposited on land as snow, thereby expanding

the ice sheets. Upon reaching a critical salinity, deep convec-

tion, and subsequently the AMO, turned on, transporting and

releasing heat to the North Atlantic and thereby melting back

the ice sheets. The fresh water flux into the North Atlantic from

the melting ice sheets (or enhanced iceberg calving) eventually

reduced or shut off the AMO, causing the process to begin

anew. This hypothesis fundamentally assumed that most of

the evaporation from the Atlantic basin was balanced by a net

freshwater transport into the Atlantic associated with an active

AMO. While this in itself has not yet been established, it could

be true for a climate with an active AMO. However, it is clear

that the weakening of the AMO would trigger a reorganization

of the circulation in the South Atlantic, but it is not clear

whether this reorganized circulation would be able to balance

the required Atlantic evaporation.

Recent models, including an interactive continental ice

sheet, showed that the earlier mechanism for D-O oscillations

was in fact opposite to what actually occurred within the

coupled system. During the cold, glacial climate, when the

AMO was still active, the mass balance over continental ice

sheets was positive. That is, rather than melting the ice sheets

when the AMO was active, they grew, since the warmer atmo-

sphere allowed greater precipitation (in the form of snow). The

models also found a mechanism for D-O oscillations that

involved an interaction between the AMO and the adjacent

continental ice sheets, by assuming that the rate of iceberg

calving into the North Atlantic was not constant in time, which

is consistent with observations of ice-rafted debris. It was

shown that after a delay of several hundred years following a

major surging event and a reduction of the iceberg calving

rate, the AMO intensified and the North Atlantic experienced

an abrupt shift to a much warmer climate. This mechanism

appears to be more intuitively appealing than a competing

hypothesis involving stochastic resonance of the AMO, which

relies upon the existence of an unknown 1,500 year external

periodic forcing.

Heinrich (1988), in analyzing marine sediments in three

cores from the North Atlantic, noted the presence of six anom-

alous concentrations of lithic fragments over the last glaciation.

Since the source for these fragments was the land (and in

particular Canada), he argued that this provided evidence for

six anomalous surges of icebergs into the North Atlantic

(Figure T5). The so-called Heinrich events (H) were even more

striking when expressed in terms of the ratio of lithic fragments

to the sum of lithic fragments and foraminiferal shells, due to

the low foraminiferal counts in the Heinrich sediment layers.

A simple model illustrates the mechanism for Laurentide Ice

Sheet instability, which ultimately gives rise to an ice-sheet

surge and a Heinrich event. It has been argued that the Heinrich

Figure T3 The UVic Earth System Climate Model (Weaver et al., 2001)

hysteresis of NADW (blue), AAIW (red), and AABW (green), showing: (a)

surface densities in the source regions of these water masses; (b)

maximum (or minimum) MOC cells associated with these water masses.

The model was started from a mode with no NADW formation with

zero externally imposed freshwater flux. A freshwater perturbation was

then applied to the AAIW formation region with a linearly-increasing

rate of 0.2 Sv / 1,000 years. This rate was small enough to ensure a

nearly steady-state system, except for abrupt transitional events. Upon

transition to a mode with active NADW formation, the freshwater

perturbation was reduced at the same rate back to zero, and then to

negative values, in order to trigger a transition from the active to the

inactive mode of NADW formation. To complete the whole hysteresis

loop, the model was integrated for more than 11,000 years. Transitions

between the active and inactive modes of NADW formation occur

when the surface densities at the site of AAIW and NADW formation

cross. Taken from Saenko et al. (2003).

946 THERMOHALINE CIRCULATION

Figure T4 Variations in local Antarctic atmospheric temperature, as derived from oxygen isotope data, as well as concentrations of atmospheric

and CH

from Vostok, Antarctica ice core records over the last four glacial cycles. The rapid oscillations represent Dansgaard-Oeschger

cycles. Figure redrawn from Petit et al. (1999).

Figure T5 Sediment foraminiferal records from Bond et al. (1993) and ice core oxygen isotope records from Dansgaard et al. (1993) and GRIP

(1993) as correlated by Bond et al. (1993). The Younger Dryas (YD) and six Heinrich events (H1–H6) are indicated. Also shown at the bottom is the

saw tooth Bond Cycle pattern of successively weaker interstadials associated with a sequence of D-O oscillations following a Heinrich event (from

Bond et al., 1993).

THERMOHALINE CIRCULATION 947

cycle consists of two phases. In the growth phase, the Lauren-

tide Ice Sheet grew through snow accumulation while remain-

ing attached to the bedrock. In the purge phase, the high

pressures caused by the deep ice sheet caused thawing near

the base of the ice sheet thereby allowing the ice sheet to surge

seaward over a lubricated base through Hudson Strait. The

resulting freshwater discharge into the North Atlantic would

be of the order of 0.1–0.2 Sv over a period as short as

250–500 years.

Heinrich events, appearing about every 10,000 years, occur

at the end of a sequence of D-O cycles during a prolonged cold

period. Bond et al. (1993) noted that the sequences of D-O

oscillations tend to follow a saw-tooth cycle (now termed a

Bond Cycle) with successive D-O oscillations involving pro-

gressively cooler interstadials (Figure T5). They argued that

this Bond Cycle was terminated by a Heinrich event, after

which a rapid warming occurred and the process began anew.

The last cold event, known as the Younger Dryas (YD) (also

known as Heinrich event zero, H0) took place between 12,700

and 11,650 years

BP and terminated abruptly within a few

decades. Proxy data suggest that during the YD, even the

AMO was weakened. Further analysis of the proxy record

indicates that the YD was preceded by a warm period in the

North Atlantic climate, known as the Bølling-Allerød (B-A)

warm interval (Figure T5). The B-A interval was accompanied

by a strong global sea-level rise, apparently due to the melting

of continental ice sheets. Until recently, it was not clear how to

reconcile the warm North Atlantic climate during the B-A with

melting of adjacent ice sheets. The latter would have caused a

reduction of AMO, leading to a cooling, rather than warming

of the North Atlantic. A possible explanation lies in the hypoth-

esis that a significant portion of the observed rapid sea level

rise during the B-A could have originated from the melting of

Antarctic ice sheets, rather than solely from melting of the ice

sheets surrounding the North Atlantic. A meltwater pulse into

the Southern Ocean, originating from the partial collapse of

an Antarctic ice sheet, could trigger an enhanced formation of

NADW and an onset of warming in the North Atlantic (Weaver

et al., 2003).

Despite these fundamental advances over the last few years,

many challenges rem ain. In particular, a complete explanation

including the essential role of the MOC, for the existence of

millennial timescale (D-O) variability and its packaging into

Bond Cycles in cold climates, its association with Heinrich

events, and its dependence on the mean climatic state remains

elusive to the pa leoclimate community. Another major climate

mystery is the cause of atmospheric CO

and CH

changes

during the glacial cycles (Figure T4) and their possible

connection to long -term changes of str ucture and strength of

the GOC.

Andrew J. Weaver and Oleg A. Saenko

Bibliography

Bond, G., Broecker, W., Johnsen, S., McManus, J., Labeyrie, L., Jouzel, J.,

and Bonani, G., 1993. Correlations between climate records from North

Atlantic sediments and Greenland ice. Nature, 365, 143–147.

Dansgaard, W., Johnsen, S.J., Clausen, H.B., Dahl-Jensen, D., Gundestrup,

N.S., Hammer, C.U., Hvidberg, C.S., Steffensen, J.P., Sveinbjörnsdottir,

A.E., Jouzel, J., and Bond, G., 1993. Evidence for general instability of

past climate from a 250-kyr ice-core record. Nature, 364, 218–220.

GRIP Project Members, 1993. Climate instability during the last intergla-

cial period recorded in the GRIP ice core. Nature, 364, 203–207.

Heinrich, H., 1988. Origin and consequences of cyclic ice rafting in the

northeast Atlantic Ocean during the past 130,000 years. Quaternary

Res., 29, 143–152.

Killworth, P.D., 1983. Deep convection in the world ocean. Rev. Geophys.

Space Phys., 21,1–26.

Petit, J.R., Jouzel, J., Raynaud, D., Barkov, N.I., Barnola, J.-M., Basile, I.,

Benders, M., Chappellaz, J., Davis, M., Delayque, G., Delmotte, M.,

Kotlyakov, V.M., Legrand, M., Lipenkov, V.Y., Lorius, C., Pépin, L.,

Ritz, C., Saltzman, E., and Stievenard, M., 1999. Climate and atmo-

spheric history of the past 420,000 years from the Vostok ice core,

Antarctica. Nature, 399, 429–436.

Rintoul, S.R., 1991. South Atlantic interbasin exchange. J. Geophys. Res.,

96, 2675–2692.

Saenko, O.A., Weaver, A.J., and Gregory, J.M., 2003. On the link between

the two modes of the ocean thermohaline circulation and the formation

of global-scale water masses. J. Clim., 16, 2797–2801.

Schmitz, W.J. Jr., and McCartney, M.S., 1993. On the North Atlantic circu-

lation. Rev. Geophys., 31,29–49.

Weaver, A.J., Eby, M., Wiebe, E.C., Bitz, C.M., Duffy, P.B., Ewen, T.L.,

Fanning, A.F., Holland, M.M., MacFadyen, A., Matthews, H.D.,

Meissner, K.J., Saenko, O., Schmittner, A., Wang, H., and Yoshimori,

M., 2001. The UVic Earth System Climate Model: Model description,

climatology and application to past, present and future climates.

Atmosphere-Ocean, 39, 361–428.

Weaver, A.J., Saenko, O.A., Clark, P.U., and Mitrovica, J.X., 2003.

Meltwater pulse 1A from Antarctica as a trigger of the Bølling-Allerød

warm interval. Science, 299, 1709–1713.

Wunsch, C., 2002. What is the thermohaline circulation? Science

, 298,

1179–1181.

Cross-references

Antarctic Bottom Water and Climate Change

Binge-Purge Cycles of Ice Sheet Dynamics

Bølling-Allerød Interstadial

Dansgaard-Oeschger Cycles

Heinrich Events

Millennial Climate Variability

North Atlantic Deep Water and Climate Change

Ocean Paleocirculation

Younger Dryas

TILLS AND TILLITES

Tills, and their lithified equivalents, tillites, are sediments

deposited by glacier ice. In the past, this definition has been

employed in a very broad sense, to include unsorted sediments

deposited in subglacial, supraglacial, ice-marginal, and subaqu-

eous environments. Such a broad definition is difficult to apply

in practice, because it may hard to differentiate sediments

deposited directly from ice and those later reworked by non-

glacial processes. Many researchers now prefer to restrict the

term till to sediments deposited directly by grounded glacier

ice, in subglacial and supraglacial environments, and which

have not undergone significant reworking by fluvial, gravita-

tional, or other processes (Dreimanis, 1989; Benn and Evans,

1998). The non-genetic terms diamict and diamictite have

now generally replaced the older, broad usage of till and tillite.

The correct identification of tills and tillites is important,

because they provide direct evidence for the former distribu-

tion of glaciers, and hence have wide-reaching paleoclimatic

implications.

The main processes of formation of tills are: lodgement,

deformation, melt-out, and sublimation. Lodgement refers to

frictional retardation of clasts or slabs of debris below sliding

948 TILLS AND TILLITES

ice, and deformation is the re-molding of pre-existing sedi-

ments or debris in a subglacial deforming layer (Boulton,

1982). Melt-out and sublimation remove interstitial ice from

englacial debris (Paul and Eyles, 1990). Tills rarely form by a

single process (Hicock, 1990). Many tills exhibit evidence of

subglacial deformation, indicating shear of unlithified sedi-

ments below glaciers and ice sheets although this is contested

by some researchers. Unmodified melt-out tills are rare,

because high pore-water pressures during de-icing encourage

re-molding by gravitational processes. Sublimation tills are

restricted to arid polar environments, where temperatures never

exceed 0



C and melting does not occur.

The identification of tills and tillites should be based on

multiple criteria, including: (a) faceted and striated clasts, with

stoss-lee forms; (b) far-traveled, erratic clasts; (c) moderate to

strong a-axis fabrics, which show regional consistency; (d) evi-

dence of shear (e.g., low-angle faults, attenuated stringers, and

overfolds) at macro- and micro-scales; (e) tabular geometry;

(f) massive, matrix-supported character; (g) multi-modal

grain-size distributions, often with a significant silt fraction

(‘rock flour’); and (h) association with folded and thrust sedi-

ments (glacitectonites). Not all tills will exhibit all of these

characteristics, and some non-glacial sediments may exhibit

some of them. Tills form a very variable group of sediments,

and their recognition requires careful study.

Douglas I. Benn

Bibliography

Benn, D.I., and Evans, D.J.A., 1998. Glaciers and Glaciation. London:

Arnold.

Boulton, G.S., 1982. Subglacial processes and the development of glacial

bedforms. In Davidson-Arnott, R., Nickling, W., and Fahey, B.D.

(eds.), Research in Glacial, Glacio-fluvial, and Glacio-lacustrine Sys-

tems. Norwich: Geobooks, pp. 1–31.

Dreimanis, A., 1989. Tills: Their genetic terminology and classification. In

Goldthwait, R.P., and Matsch, C.L. (eds.), Genetic Classification of

Glacigenic Deposits. Rotterdam: Balkema, pp. 17–84.

Hicock, S.R., 1990. Genetic till prism. Geology, 18, 517–519.

Paul, M.A., and Eyles, N., 1990. Constraints on the preservation of diamict

facies (melt-out tills) at the margins of stagnant glaciers. Quaternary

Sci. Rev., 9,51–69.

Cross-references

Diamicton

Glacial Geomorphology

Glacial Sediments

Glaciomarine Sediments

Ice-Rafted Debris (IRD)

Moraines

TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA

Introduction

This article is a brief synopsis of time-series analysis methods

that have proven to be useful in the analysis of paleoclimate

data. Space limitations precludes all but the briefest outline of

the actual mathematics involved, but references are provided.

The following sections define some basic functions such as

spectrum and autocorrelations, outline the basic multitaper

methods that are recommended to estimate them, describe

periodic components and Cramér–Rao bounds for frequency

estimation, and coherence estimates. Other sections below discuss

jackknife confidence limits for spectra and similar functions, spec-

trogram methods useful with nonstationary processes, discusses

wavelets, Lomb–Scargle and Blackman–Tukey estimates which

should be avoided for these applications. Finally, the last section

summarizes the points made.

One should not expect the analysis of paleoclimate data to

be easy. If one considers the effort that has been made over

the last three or four decades in making general circulation

models reproduce even fairly gross features of the climate sys-

tem, one should not expect that the analysis of data produced

by the real, and obviously far more complicated, climate system

can be done simply. Added to this intrinsic difficulty, paleocli-

mate data is usually a proxy, the geographic distribution of

samples is far from ideal, and from a data analysis viewpoint,

the fact that the timescale is uncertain can be a major problem.

Worse than these, however, is the problem that the climate sys-

tem is nonstationary. Stationarity in this context implies that the

statistics of the process do not depend on time.

Theoretical quantities

In this section we define some functions of a time series. Here

we assume that the process x( t) is stationary, that is its statisti-

cal properties do not vary with time; that the expected value of

x(t) is zero, denoted E{x(t)} = 0; and that x(t) has a Gaussian

distribution. The samples are assumed to be equally spaced in

time at a unit sampling rate. Thus the Nyquist frequency,

2Dt

(where Dt is the time between samples) is

. Generally, the sin-

gle most useful description of a time series is its power spec-

trum, also known as the power spectral density, or plain

spectrum. The spectrum is a decomposition of the variance of

the process as a function of frequency, and functionals of the

spectrum give most of the interesting descriptive statistics of

a time series. Thus, if we denote the spectrum at frequency f

by S( f ), with 

 f <

, one has the process variance



Sðf Þdf ð1Þ

The units of the spectrum are (data  units)

per (cycle/time-

unit). As an example, if one has d

O data with samples measured

at a rate of one per thousand years, the frequency range will be

from 

21000

and the spectrum will be in units of

per (cycle/1000years). Most of the series encountered

in paleoclimatology are real-valued so the spectrum is sym-

metric, that is the power at  f is the same as the power at + f,

so it is sufficient to plot the spectrum for positive frequencies

only. However, if one has a series of directional data, say the

north-south and east-west components of a magnetic field or

flow, it makes sense to treat them as a complex-valued time-ser-

ies, and in this case, the spectrum S(f ) is different from S( f ).

Possibly as important as the variance, the one-step ahead

prediction variance is given by:

¼ exp



ln Sðf Þdf

ð2Þ

Climate studies are often associated with prediction and this

formula shows that the overall shape of the spectrum on a

logarithmic power scale is the critical parameter, not just the

largest peak.

TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA 949

Similarly, the single-sample interpolation variance given by



Sðf Þ

ð3Þ

shows that a missing sample may be interpolated with an accu-

racy approaching that of the measurement error. Together,

equations (1), (2), and (3) are the arithmetic, geometric, and

harmonic means of the spectrum, respectively.

The autocovariance function,

RðtÞ¼EfxðtÞxðt þ tÞg ð4Þ

is the Fourier transform of the spectrum

RðtÞ¼



i2pf t

Sðf Þdf ð5Þ

so, comparing (1) and (5), R(0) = s

. The autocorrelations,

r(t) are just the autocovariances normalized by the variance

r(t)=R(t)/s

and are, as the name implies, the correlation of

the process with delayed versions of itself. The Fourier trans-

form of the autocovariance is the spectrum. Experience shows

that the autocovariance is useful for prediction and interpola-

tion. Even though the autocovariance and spectrum carry

equivalent information, plots of the autocorrelation are usually

remarkably uninformative. (John Tukey remarked that “I only

knew one man who could look at an autocorrelation function

and make sense of it, and he’s dead.”) One useful formula is

that for the zero-crossing rate, the expected number of zeroes

per sample (Kedem, 1994),

EfR

g¼

arccosðrð1ÞÞ ð6Þ

As an example, the Lag–1 correlation of the data shown in

Figure T6 is 0.978 with (6) evaluating to 0.0668, so one

expects the curve to cross its average about every 15 samples

or so.

Since the true spectrum is not known, the best estimate pos-

sible of a spectrum must be made. In estimating a spectrum,

there are several competing criteria to optimize: bias, frequency

resolution, variance, efficiency, ability to distinguish periodic

from non-periodic components, minimum assumptions, and

robustness. The spectrum is, of course, not the only useful

function and one should always plot the raw data and ancillary

functions such as histograms, plots of x(t) vs x(t  1), Q-Q

plots, as well as basic data checks such as histograms of bit

and digit frequencies. Also, if the data is nonstationary (as

much climate data is) one plots a spectrogram, see below, basi-

cally a sequence of spectrum estimates, usually computed on

overlapping time segments and displayed as a function of both

time and frequency.

Estimating the spectrum

As a data example, one can use the first 3 Myr of Shackleton’s

composite benthic

O data from V19-30, ODP 677, and ODP

846 (see Shackleton et al., 1990, 1995a,b)

. The time step was

scaled by a factor of the ratio 41090.0/41234.9 to bring the

strong obliquity line frequency into exact agreement with

Laskar (1990). For the 3 Myr period being studied, this may

be the most accurate orbital frequency available, see Laskar

(1999). Assuming that a timescale has been established (see

Pälike and Shackleton, 2000), and N data samples, x(t), are avail-

able at a time-spacing of Dt years, the preferred method of esti-

mating the power spectrum is the multitaper method. While the

details of this method has been described in a series of papers

(Thomson, 1982, 1990a,b, 2001; Thomson et al., 2007), and

several books (Percival and Walden, 1993; Rao et al., 2003;

Weedon, 2003), the basic idea is simple and consists of the fol-

lowing steps:

1. Choose the time–bandwidth product controlling frequency

resolution and bias,

¼ NW ¼ BT ¼ B  NDt ð7Þ

so, at unit sampling rate the bandwidth, W,isW = c

/N or, in

physical units one has B = c

/T where T is the total duration

of the data sample, T = NDt. Usually one should choose

2  c

 10 and, for most paleoclimate data sets, a choice

of 4 or 5 is a reasonable starting point. In this example

= 5 was used.

Figure T6 Shackleton’s composite benthic d

O record, minus a linear trend, for the past 3 Myr.

The data was obtained from //delphi.esc.cam.ac.uk/ coredata/ v677846.

html and the few gaps interpolated; see Appendix A of Thomson et al.

(2001).

950 TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA

The time–bandwidth product, c

, is the dimension of

the space of duration T and bandwidth T and to represent

the information in such a space one needs c

bases functions.

It is convenient to work with real-valued bases functions,

with the result that the true bandwidth spans the range

from  B to + B and consequently one has a space of dimen-

sion 2BT. Thus one has a set of K  2BT bases sequences that,

in spectrum estimation problems, have historically been

known as data windows or data tapers. A typical choice is

K =2c

 2, or 2NW  2 so if, takes c

= 5 as we do here,

one would use K = 8 windows.

The choice of c

is fairly critical: if c

is too small the esti-

mate will have a high variance and relatively poor bias pro-

tection, but if c

is too large, one will have poor frequency

resolution. In this example values of the spectrum go from

aminimumof 15 to a maximum of  4.7  10

, a range

of  3,000, so extreme bias protection is not required.

2. Compute the first K lowest order discrete prolate spheroidal

sequences or Slepian sequences, v

(k)

(N,W) where 0  t 

N  1. These sequences (Slepian, 1978), have the greatest

possible energy concentration within the frequency band

(W,W) of any time-limited sequence and, in sampled data,

define the time–frequency uncertainty principle. They are

defined by the algebraic eigenvalue equation

ðN; WÞv

ðkÞ

ðN; WÞ¼

N1

m¼0

sin 2pWðt  mÞ

pðt  mÞ

ðkÞ

ðN; WÞð8Þ

The eigenvalues, l

(N,W), define the energy concentration in

(W,W). The first K b2NWcare close to one and are the ones

used in spectrum estimation. Although the Slepian sequences

are defined by (8), this is emphatically not the way to compute

them; see Appendix B of Thomson (1990a). For the rest of

this article, the explicit dependence on N and W in the

sequences and eigenvalues is dropped and l

and v

(k)

are used.

3. Subtract the mean, and possibly trend, from the data

denoted x(t).

4. For k = 0,1,...,K  1 compute the eigencoefficients

ðf Þ¼

N1

n¼0

xðnÞv

ðkÞ

i2pnf

ð9Þ

Normally the eigencoefficients are computed using a zero-

padded fast Fourier transform (FFT). The procedure is:

(a) Choose the number of transform points M  N and

usually M  N. The proper choice depends on how good

the signal-to-noise ratio of your data is, with better data

requiring larger M. It is generally unwise to choose

M < 2N. For paleoclimate data, the choice M = N is

almost always a poor one as it implies circular correla-

tions. (This means that the sample autocovariances

obtained by taking the Fourier transform of the spectrum

estimate, the sample equivalent of (5), would include both

“legitimate” terms such as x(n)x(n + t)and“illegitimate”

ones such as x(n)x(n + t  N) where the end of the sam-

ple wraps around and is correlated with the start.) The

efficiency of the FFT algorithm depends on M being fac-

torable into small primes, so choosing M as a product of

2’s, 3’s, 5’s, and other small primes up to 19 is necessary.

Also note that, if one has a few hundred data points and

the estimated spectrum differs noticeably between the

choice of M =4000=2

 5

and M =4096=2

you are likely doing something wrong.

(b) Multiply the data by the k

taper. If one is using a stan-

dard (Singleton, 1969) FFT, say FFT(M,X,Y) where X

and Y are arrays of length M, X being the real part

and Y the imaginary, one sets

X ðnÞ¼

xðnÞv

ðkÞ

0  n  N – 1,

0N n  M – 1



ð10Þ

YðnÞ¼00 n  M  1 ð11Þ

The FFT computes

out

ðjÞþiY

out

ðjÞ¼

M1

n¼0

½X

ðnÞþiY

ðnÞe

i2p

ð12Þ

for j =0,...M  1. Because the data encountered in

paleoclimate work is usually real-valued, it is only

necessary to keep the outputs for positive frequencies,

that is for j =0,1,...,M/2. Taking the time at t he n

sample to be t

= nDt and the frequency at the j

FFT bin as f

= jdf and comparing exponents it is

evident that

df ¼

MDt

ð13Þ

so the effect of zero padding is simply to compute the

frequencies on a finer grid. Note that one occasionally

sees statements to the effect that “aFFTwasnotused

and the Fourier transform was computed exactly.” This

is a fallacy: not only is the FFT far quicker than naive

evaluation of the formula, but it also computes exactly

what it claims to compute and, because it uses far few

operations to do it, suffers less from roundoff error than

direct evaluation of the formula.

5. In practice one estimates the spectrum by solving (17)

below, but to understand the properties of a multitaper

estimate, it is easier to use the approximation



Sðf Þ¼

K1

k¼0

ðf Þj

ð14Þ

on the grid of frequencies f

= jdf, j = 0,1,...,M/2. This

estimate has the following properties:

(a) The bandwidth is

B ¼

¼

NDt

ð15Þ

Specifically, given a pure sinusoid at frequency f

, one-

half of the energy in the sinusoid would be distributed

across the band (f

 W, f

+ W), and the other half of

the energy across the symmetric band (f

 W, f

W ). As will be described in the next section, when sinu-

soids are known to be present, one can estimate their fre-

quencies much more accurately than B.

(b) Under reasonable conditions the estimate (14) will have

a chi-squared distribution with 2K degrees-of-freedom,

written w

. At frequencies where there is a periodic, or

other narrow-band component, the distribution will be a

non-central chi-square, while at frequencies between such

components i.e., at frequencies where the signal is noise-

like, the distribution is central chi-square. It is safer to

TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA 951

estimate the variance by using a jackknife than by relying

on classical distribution theory. See below for details.

As standard practice, both the classical and jackknife

estimates and are computed and compared.

band ( f  W, f + W) through the sidelobes of the win-

dows, can be bounded by

Bðf Þ <

ð1  l

K1

Þs

ð16Þ

where l

is the eigenvalue of the k

Slepian sequence

defined in (8), and s

the process variance.

Finally, while (14) is easy to describe, in practice one uses an

adaptively weighted estimate. Details are given in Thomson

(1982) or on p 546 of Thomson (1990a), but the critical step

is to solve

K1

k¼0

Sðf Þjx

ðf Þj

½l

Sðf Þþð1  l

Þs



¼ 0 ð17Þ

for

Sðf Þ. It can be shown that this solution is between

min

{jx

( f )j

} and max

{jx

( f ) j

} and, in practice, an iterative

solution starting at

Sðf Þ¼ðjx

ðf Þj

þjx

ðf Þj

Þ=2 usually

converges in 2 or 3 iterations. The properties are very similar

to those listed above except that, with (17), the effective

degrees-of-freedom will be frequency dependent and usually

lower near the minimum of the spectrum. This loss of

degrees-of-freedom can be reduced by prewhitening.

Lines and Crame

r–Rao bounds

A source of confusion in the analysis of paleoclimate data (and

elsewhere) is the distinction between the continuous and line

components of the spectrum. The lines are isolated periodic

components such as those occurring in the familiar, e.g., Berger

(1977, 1978), trigonometric expansions of eccentricity, obli-

quity, and precession, that is, sinusoids at discrete frequencies.

The continuous part of the spectrum is basically everything

else: measurement errors and related problems such as biotur-

bation, and non-deterministic parts of the climate process with

a short prediction horizon.

The basic multitaper estimate for periodic components can

be understood by considering the effects of a single sinusoid

in noise. Define

yðtÞ¼xðtÞþA cosð2pf

t þ fÞð18Þ

which can be written

yðtÞ¼xðtÞþ

i2pf

if

i2pf

; ð19Þ

where A2e

¼ m is the complex amplitude and x(t) is assumed

to be non-periodic. The eigencoefficients of y(t) are

ðf Þ¼x

ðf ÞþmV

ðf f

Þþm



ðf þf

Þð20Þ

where V

( f ) is the Slepian function

ðf Þ¼

N1

n¼0

ðkÞ

i2pfn

ð21Þ

Because of the energy concentration properties of the Slepian

sequences, the term m

∗

( f + f

) can usually be ignored if

> W. Because it is assumed that x(t) is a zero-mean process,

the expected value

Efy

ðf Þg ¼ mV

ðf f

Þð22Þ

One does a formal hypothesis test for the presence of a periodic

component. The null hypothesis, H

, is that the frequency band

(f  W, f + W) has an approximately flat (called “locally

white”) spectrum of unknown power. The alternative hypoth-

esis, H

, is that in addition to the locally white noise, there is

a sinusoid at frequency f. Estimate its amplitude m by least-

squares, that is minimize the error at frequency f between the

data, the y

( f )’s, and the model, (22). For an assumed complex

amplitude

mðf Þ, the error is

ðf ;

mðf ÞÞ ¼

K1

k¼0

ðf Þ

mðf ÞV

ð0Þj

; ð23Þ

where the assumption that the line frequency f

= f results in the

zero argument of the Slepian function. Taking the partial deri-

vative of E with respect to

mðf Þ



, one obtains

mðf Þ¼

K1

k¼0

ðf ÞV

ð0Þ

K1

k¼0

ð0Þj

; ð24Þ

the regression of the eigencoefficients on their expected values.

One tests this regression in the usual way by using an F–test.

This is the ratio of the energy explained by the assumed

periodic component at frequency f to that remaining in the

residuals,

Fðf Þ¼

mðf Þj

K1

k¼0

ð0Þj

2K2

ðf ;

mðf ÞÞ

: ð25Þ

To understand this equation, recall that an F–statistic is defined

as (w

/p)/(w

/q) where w

and w

are independent chi-squared

variables of p and q degrees-of-freedom, respectively. If there

is no periodic component present, just “noise”, then, because

mðf Þ is complex, j

mðf Þj

is chi-squared with two degrees-of-

freedom. The residual, E

, has 2K  2 degrees-of-freedom

and is independent of the estimated mean, so F( f ) has a central

F-distribution under H

It can be shown by some tedious algebra that if one uses the

frequency where F( f ) is maximized as an estimate of f

,itis

usually within a few percent

of the Cramér–Rao bound on

the achievable accuracy of a frequency estimate.

Notes on the har monic F–test

When interpreting F -tests for lines, some precautions should be

observed.

1. Because the F-test is computed as a function of frequency,

there are approximately N independent tests in the fre-

quency domain. Thus one should expect some false detec-

tions. For example, if one adopts a conventional 95%

level, one should expect about 0.05N “detections” in cases

when there are no lines present.

2. When there are strong lines present, the single-line F-test

given above will suppress false detections within  W of

these lines.

The actual factor is the estimate efficiency X defined in Thomson (1982)

}VII.

952 TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA

3. The test is sensitive to departures in assumptions. For exam-

ple, an amplitude or phase modulated sinusoid will usually

be rejected if the amplitude varies by a factor of two over

the sample interval, or if the frequency,

f, varies by

about one-third of a Rayleigh resolution, that is, 1/(3T ).

4. The test usually rejects multiple lines within the band.

5. When one observes a roughly rectangular peak of total

width ⪆ 2W in the spectrum without a peak in the F-test

at the center frequency of the peak, one should test for mul-

tiple lines or for a modulated periodic component. (The

lower precession band in Figure T8 is typical; this band,

42–45 c/Myr, contains four strong lines of which two are

detected.) Tests are described in Thomson (1990a) and

Johnson et al. (1996).

6. A peak in the F-test corresponding to a local minimum in

the spectrum is usually a false detect. The line marked

“20,850” in Figure T8 is almost certainly in this category.

Similarly, the lines marked “13,397,” and “12,540” may be

false detections, but could be a nonlinear interaction

between obliquity and precession. Similarly, maxima in

F on the sides of peaks, or on very steep slopes, are often

false detections. The line marked “37,430” in Figure T8 is

on the shoulder of the obliquity band. The reason for these

false detections is that the “locally white” condition is

violated, causing the y

’s to be correlated.

7. When a sharp peak in the F–test occurs near the center of a

roughly rectangular peak of total width  2W,itisusually

reliable. The 41,090-year obliquity line and 18,957-year pre-

cession line are prime examples, but again, the upper preces-

sion band contains lines at several closely-spaced periods.

The Crame

r–Rao bound

The Cramér–Rao bound is a fundamental statistical inequality

that gives lower limits on the accuracy with which various

quantities can be estimated; here the frequency. For a single,

isolated line, i.e., one that is isolated from other lines by at least

max(2,c

)/T, the bound is approximately

SDf

f g

2pT

ﬃﬃﬃ

ð26Þ

where r is the effective signal-to-noise power ratio,

r ¼

ðf Þ

ð27Þ

where S

( f ) is the noise spectrum. This is the spectrum without

the periodic components; in Figure T7 the noise spectrum near

the 41,000 year obliquity line is about 500. Combining (26)

and (27) it can be seen that, other things being equal, the var-

iance of a frequency estimate decreases as T

3

. Thus, when a

choice can be made, effort should be made obtaining longer

series, not in making more precise measurements in a shorter

series. As an example, if one wanted to double the accuracy

of a frequency estimate one could consider either making more

precise measurements or increasing the length of the series.

The first choice, implies that s

must be reduced by a factor

of four, usually corresponding to a four-fold increase in mea-

surement effort. The second requires increasing T by

ﬃﬃﬃ

,or

Figure T7 A multitaper estimate of the spectrum of the 3 Myr record

shown in Figure T6. In making this estimate, c

= NW = 5.0 with K =8

windows were used with adaptive weighting, (17). The dashed curve

near 25 c /Myr (the 41090-year obliquity line) is the “raw” spectrum

estimate, while the sharp black curve is the “reshaped” estimate.

Because the significance level of the F–test exceeded the 99.999% level,

the raw estimate was replaced with a peak at the estimated frequency.

The width of the peak is given by the estimated frequency uncertainty,

and the height is scaled to preserve power. The precession peaks (near

43 and 52 c / Myr) were not reshaped because they both consist of

multiple lines and the eccentricity lines (7.8 – 10.5 c /Myr) are not

resolved.

Figure T8 The single-line F –test for periodic components, (25). The

data is the 3 Myr record shown in Figure T6, with parameters are as

in Figure T7, so the test has 2 and 14 degrees-of-freedom. The

low-frequency part is shown in Figure T10 and the horizontal dashed

line shows the theoretical 95% level of the F

2,14

distribution.

TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA 953

about 26% more samples. Note that this calculation ignores the

tuning problem.

If one is interested in estimating the period of a component,

as is common in paleoclimate data, a Taylor’s series argument

shows that (26) becomes

SDf

Pg

2pT

ﬃﬃﬃ

ð28Þ

The F–test on the 41,000 year obliquity line was  30, and can

be used as an estimate of r. Thus, if one oversimplifies things

to ignore fine structure, possible frequency-modulation from

changes in the earth’s moments of inertia (see }11 of Thomson

(1990a) and Figure T14 of Laskar et al., 1993) the accuracy of

any estimate of the period from the 3 Myr record is bounded by

SDf

Pg

41000

2p 3  10

ﬃﬃﬃﬃﬃ

 38 years ð29Þ

Because the period estimated from the original time scale was

41,234.9 years, 3.8 standard deviations from the nominal per-

iod of 41,090 years, the time was scaled by a factor of

41090.0/41234.9 to bring the two into agreement.

Equation (28) has further implications for data analysis.

Because the Cramér–Rao bound is a fundamental statistical

inequality, it does not apply to one particular estimation method

or another, it simply says that one cannot estimate parameters

with arbitrary accuracy. Arguments based on wavelet analysis

of paleoclimate data claim that the precession line is frequency

modulated. The usual trigonometric expansions used in celes-

tial mechanics are, in fact, the Bessel function expansions of

a frequency-modulated sinusoid. Since wavelets sequentially

use data blocks half the length of the preceding one, the effec-

tive T in (28) at a given stage is much less than the length of

the series. Thus (28) implies that the standard deviation will

be proportionally larger, negating the argument.

Equation (26) also determines the amount of zero-padding

necessary. A minimum requirement is that the frequency mesh,

(13), should be less than about one-half the standard deviation

implying

M ⪆ 4p

ﬃﬃﬃ

N ð30Þ

For the 3 Myr series this implies that the 1,000 samples should

be zero-padded to 2M ⪆ 28,000 points and, in fact, these

examples were zero-padded to 2M =2

= 131,072.

When lines are not isolated and, in particular, when they are

separated by less than twice the Rayleigh resolution, computing

the bounds is more difficult. The bounds depend on the amplitudes

and phases of the lines in addition to the frequency separation.

Looking at the eccentricity expansions, there are four major lines

with periods near 100,000 years. They are, explicitly, 94,782,

98,715, 123,818, and 130,615 years. Note that these come in

two pairs; the frequency difference between the two frequencies

in the first pair and those in the second both correspond to a period

of 2.4 Myr. This means that estimating these frequencies with less

than 4.8 Myr of data will be difficult even in utopian circum-

stances. Figure T9 shows the bounds for the four-line case com-

puted using the frequencies, phases, and amplitudes from Berger

(1977) as a function of the observation span. Note that, when the

data span is less than 800,000 years, the bound on the standard

deviation is 100,000 years, about the same as the period. Note also

that these bounds are for Gaussian noise with a known timescale;

when the timescale is estimated the bounds will be worse. More-

over, four pure sinusoids are assumed, whereas the eccentricity

expansion includes dozens of lines in this band so, again, the

bounds will be even worse than shown here. Thus, for example,

attempts to discriminate between frequencies of 77,000 years

and 95,000 lines with less than a few million years of data are

futile. (It may be futile even with a few million years of data

because 4.8 Myr is long enough for continental drift to signifi-

cantly change climate dynamics.)

The multitaper setup for multiple lines is discussed in

Thomson (1990a), and an example using the above data is

shown here in Figure T10. The two-line test used here assumed

a known frequency offset between the two lines, 1/2.2379 

cycle/year, but even with this assumption, the test is only

significant at about the 97% level.

Coherence

Many analysis problems in paleoclimate data involve studying

several different series. The simplest case is the bivariate one,

that is where one has two time series, x(t) and y(t). These are

assumed to be on the same timescale, are sampled at the same

time-step, and have the same number of samples

. Keeping the

Figure T9 The Crame

r–Rao lower bound for the standard deviation of

the 123,000-year eccentricity line as a function of sample duration. A

uniform sampling rate of Dt = 1, 000 years is assumed. The upper curve

corresponds to a signal-to-noise ratio r = 5, the lower to an

unrealistically high value of 100. For a perfectly-tuned 800 Kyr record

the standard deviation of the best period estimate is about the same as

the period.

If the time-steps and number of samples are different one can still do the

analysis as long as the total time span is the same and, obviously, one can

always truncate a longer series. The frequencies used must be the same, so it

may be necessary to use a “slow” Fourier transform (see e.g., Gentleman,

1969).

954 TIME–SERIES ANALYSIS OF PALEOCLIMATE DATA