10 December 1976, Volume 194, Number 4270 SCIENCE

the last interglacial on the basis of these

curves have ranged from 80,000 to
180,000 years ago (22).
The second and more critical problem
in testing the orbital theory has been the
uncertainty of geological chronology.
Until recently the inaccuracy of dating
Variations in the Earth's Orbit: methods limited the interval over which
a meaningful test could be made to the
Pacemaker of the Ice Ages last 150,000 years. Hence the most con-
vincing arguments advanced in support
of the orbital theory to date have been
For 500,000 years, major climatic changes have based on the ages of 80,000, 105,000, and
followed variations in obliquity and precession. 125,000 years obtained for coral terraces
first on Barbados (15) and later on New
Guinea (23) and Hawaii (24). These struc-
J. D. Hays, John Imbrie, N. J. Shackleton tures record episodes of high sea level
(and therefore low ice volume) at times
predicted by the Milankovitch theory.
Unfortunately, dates for older terraces
For more than a century the cause of hypothesis has been formulated so as to are too uncertain to yield a definitive test
fluctuations in the Pleistocene ice sheets predict the frequencies of major Pleisto- (25).
has remained an intriguing and unsolved cene glacial fluctuations. Thus it is the More climatic information is provided
scientific mystery. Interest in this prob- only explanation that can be tested geo- by the continuous records from deep-sea
lem has generated a number of possible logically by determining what these fre- cores, especially the oxygen isotope rec-
explanations (1, 2). One group of theo- quencies are. Our main purpose here is ord obtained by Emiliani (26). However,
ries invokes factors external to the cli- to make such a test. the quasi-periodic nature of both the iso-
mate system, including variations in the Previous work has provided strong topic and insolation curves, and the un-
output of the sun, or the amount of solar suggestive evidence that orbital changes certain chronology of the older geologic
energy reaching the earth caused by induced climatic change (13-20). How- records, have combined to render plau-
changing concentrations of interstellar ever, two primary obstacles have led to sible different astronomical interpre-
dust (3); the seasonal and latitudinal dis- continuing controversy. The first is the tations of the same geologic data (13, 14,
tribution of incoming radiation caused by uncertainty in identifying which aspects 17, 27).
changes in the earth's orbital geometry of the radiation budget are critical to
(4); the volcanic dust content of the atmo- climatic change. Depending on the lati-
sphere (5); and the earth's magnetic field tude and season considered most signifi- Strategy
(6). Other theories are based on internal cant, grossly different climatic records
elements of the system believed to have can be predicted from the same astro- All versions of the orbital hypothesis
response times sufficiently long to yield nomical data. Milankovitch (4) followed of climatic change predict that the obliq-
fluctuations in the range 104 to 106 years. Koppen and Wegener's (21) view that uity of the earth's axis (with a period of
Such features include the growth and the distribution of summer insolation (so- about 41,000 years) and the precession of
decay of ice sheets (7), the surging of the lar radiation received at the top of the the equinoxes (period of about 21,000
Antarctic ice sheet (8); the ice cover of atmosphere) at 65°N should be critical to years) are the underlying, controlling
the Arctic Ocean (9); the distribution of the growth and decay of ice sheets. variables that influence climate through
carbon dioxide between atmosphere and Hence the curve of summer insolation at their impact on planetary insolation.
ocean (10); and the deep circulation of this latitude has been taken by many as a Most of these hypotheses single out
the ocean (11). Additionally, it has been prediction of the world climate curve.
argued that as an almost intransitive Kukla (19) has pointed out weaknesses The authors are all members of the CLIMAP
system, climate could alternate between in Koppen and Wegener's proposal and project. J. D. Hays is professor of geology at Colum-
bia University, New York 10027, and is on the staff
different states on an appropriate time has suggested that the critical time may of the Lamont-Doherty Geological Observatory,
scale without the intervention of any ex- be September and October in both hemi- Palisades, New York 10964. John Imbrie is Henry L.
Doherty professor of oceanography, Brown Univer-
ternal stimulus or internal time constant spheres. However, several other curves sity, Providence, Rhode Island 02912. N. J. Shackle-
(12). have been supported by plausible argu- ton is on the staff of the Sub-department of Quaterna-
ry Research, Cambridge University, Cambridge,
Among these ideas, only the orbital ments. As a result, dates estimated for England.
10 DECEMBER 1976 1121
 mechanisms of climatic change which Methods reproducibility for independent samples
are presumed to respond to particular from the sediment is about + 0.1 per mil
elements in the insolation regime (28). In Core selection. From several hundred [1 standard deviation (S.D.)]. Tests of G.
our more generalized version of the hy- cores studied stratigraphically by the bulloides are formed primarily at some
pothesis we treat secular changes in the CLIMAP project, we selected two depth below the sea surface so that
orbit as a forcing function of a system (RC1 1-120 and E49-18) whose locations changes in surface temperature do not
whose output is the geological record of (Fig. 1 and Table 1) and properties make greatly affect the temperature at which
climate-without identifying or evaluat- them ideal for testing the orbital hypothe- the carbonate is secreted. Down-core
ing the mechanisms through which cli- sis. Most important, they contain togeth- variations in 8180 reflect changes in oce-
mate is modified by changes in the global er a climatic record that is continuous, anic isotopic composition, caused pri-
pattern of incoming radiation (29). Most long enough to be statistically useful marily by the waxing and waning of
of our climatic analysis is based on the (450,000 years), and characterized by ac- Northern Hemisphere ice sheets (31, 33).
simplifying assumption that the climate cumulation rates fast enough (> 3 cen- Thus, the 8180 in our subantarctic cores
system responds linearly to orbital forc- timeters per 1,000 years) to resolve cli- is a Northern Hemisphere climatic rec-
ing. The consequences of a more realis- matic fluctuations with periods well be- ord.
tic, nonlinear response are examined in a low 20,000 years. In addition, the cores Estimates of Ts were made by apply-
final section here. are centrally located between Africa, ing statistical techniques (18) to subant-
Our geological data comprise measure- Australia, and Antarctica, and therefore arctic radiolarians. The data base for
ments of three climatically sensitive pa- little influenced by variations of erosion- writing the transfer functions has been
rameters in two deep-sea sediment al detritus from these continents. Final- expanded from that of Lozano and Hays
cores. These cores were taken from an ly, as explained below, a Southern Hemi- (34) by G. Irving and J. Morley to in-
area where previous work shows that sphere location provides an opportunity clude cores in the vicinity of E49-18. The
sedimeht is accumulating fast enough to to monitor simultaneously both North- accuracy of Ts as an estimate of sea-
preserve information at the frequencies ern Hemisphere ice volume and South- surface temperature is + 1.5°C; its repro-
of interest. Measurements of one vari- ern Hemisphere temperature. No other ducibility as an index of faunal change is
able, the per mil enrichment of oxygen- cores known to us have all these attri- 0.32°C (1 S.D.).
18 (8180), make it possible to correlate butes. The percentage of C. davisiana rela-
these records with others throughout the Geological data. We have measured tive to all other radiolarians was deter-
world, and to establish that the sediment (i) 6180, the oxygen isotopic composition mined by techniques previously de-
studied accumulated without significant of planktonic foraminifera; (ii) Ts, an scribed (35). These counts have a pre-
hiatuses and at rates which show no estimate of summer sea-surface temper- cision of about + 0.74 percent (1 S.D.).
major fluctuations. atures at the core site, derived from a The recent distribution of the cosmopoli-
To be used in tests of the orbital hy- statistical analysis of radiolarian assem- tan species C. davisiana Petrushevskaya
pothesis, these data are transformed into blages; and (iii) percentage of Cy- shows no relation to present sea-surface
geological time series. In our first test we cladophora davisiana, the relative abun- temperature (34, 36), so that the remark-
make the simplest geochronological as- dance of a radiolarian species not used in able Pleistocene abundance variations of
sumption, that sediment accumulated in the estimation of Ts. Identical samples this species (which comprised up to 50
each core at a constant rate throughout were analyzed for the three variables at percent of the fauna during glacial maxi-
the period of study. Later we relax this 10-cm intervals through each core (30). ma in some areas) are also probably not
assumption and allow slight changes in Oxygen isotope analyses have been due to temperature (35). The unique high
accumulation rate, as indicated by addi- made in tests of Globigerina bulloides by abundance of this species in recent sedi-
tional geochronological data. well-established techniques (31, 32). The ments of the Sea of Okhotsk has been
Our frequency-domain tests use the related to the structure of summer sur-
numerical techniques of spectral analysis face waters, where a low-salinity surface
and are designed to seek evidence for a Table 1. Core locations and depths. layer is underlain by a strong temper-
concentration of spectral energy at the ature minimum (37). The high abundance
frequencies of variation in obliquity and Water of C. davisiana during glacial times in
precession. We consider that support for Coe Length
Core
(cm) dph
Lati-
tude
Longi-
td the Antarctic may be due to a similar
the hypothesis can be decisive if both surface water structure.
frequencies are detected and, to allow RC11-120 954 3135 43031'S 79052'E Measurements of these parameters
for geochronological uncertainties, if it E49-18 1459 3256 460031S 90009'E therefore reflect changes in three parts of
can be clearly demonstrated that the ra- the climate system: Northern Hemi-
tio of the two frequencies detected does sphere land ice, subantarctic sea-surface
not differ significantly from the predicted _ _ -, temperature, and Antarctic surface wa-
ratio (about 1.8). w ter structure.
30'
Finally, our time-domain tests are de- Orbital data. Since the pioneering
signed to examine the phase coherence work of Milankovitch (4) the chronolo-
between the three climatic records and
between each record and the postulated
0E49.18 gies for orbital and insolation changes
Antarctic Po_ar Front have been recalculated several times
forcing function. To this end we use the (38-40). These papers and those of Kuk-
numerical techniques of bandpass filter 60' la (19) should be consulted for an expla-
analysis. Such an approach makes it pos- nation and evaluation of the numerical
sible to examine separately the variance procedures used and for a discussion of
components of the geological records 30' 60' 90' 120' the manner in which orbital changes bear
that correspond in frequency to the varia- Fi;ig. 1. Locations of cores in the southern In- on the earth's insolation regime as a
tions of obliquity and precession. di,-an Ocean. function of latitude, season, and time.
1122 SCIENCE, VOL. 194
When Vernekar's (39) calculations are neity between the 8180 minima and Ts surface temperature are nearly synchro-
compared with those of Berger (40), the maxima. However, changes in Ts pre- nous. Between stage 10 and the upper
timing of inflection points in the obliquity cede changes in 8180 by a small amount; part of stage 11, however, the higher-
and precession curves do not differ by this is most evident at extreme Ts maxi- frequency fluctuations in Ts and 8180
more than 1000 years over the past ma. We conclude that over this time appear to be out of phase; and the low-
400,000 years. For intervals older than interval changes in Northern Hemi- temperature extremes are colder than
500,000 years ago, however, discrepan- sphere ice volume and subantarctic sea- they are at higher levels in the cores.
cies between the calculations become
significant, and the work of Berger is to
be preferred because it includes the ef-
fect of higher-order terms.

0
Geological Time Series Go

Stratigraphic sequence. Because the

8180 record in deep-sea sediments pri-
marily reflects global ice volume, it is
globally synchronous in open-ocean
cores and provides (along with standard
biostratigraphical techniques) a basic I-

stratigraphy for the last million years (33, LO

41, 42). This stratigraphy has a resolu-

tion limited only by ocean mixing (about
1000 years) and bioturbation. The 8180 a
record was divided by Emiliani (26) into 4- .1
numbered stratigraphic stages, which are C
tb/C9-el,.r!201e '\t .I. 'i2\ l/ rm l n /o\\
used here.
3:
Ia 2e3 j*
i
The oxygen isotope record in core 8- ./~~~~~~~~~~~~~~~-
ar I.* I.
RC 1 1-120 is complete and typical back to
late stage 9 (Fig. 2). Stage 5 shows its 12-
characteristic three peaks (26, 42); and V ~~ ~ ~ ~ ~ ~ ~
stage 7 is interrupted by a sharp positive 16-
excursion, which is typical (31, 42). II
In core E49-18 (Fig. 3) the high per- 0 1 2 3 5 6 4 7 8 9
centage of C. davisiana at the top in- Depth (mi)
dicates that the Holocene is missing (35); Fig. 2. Depth plots of three parameters measured in core RC I 1-120: 6180 (solid line), Ts (dashed
the top may be as old as 60,000 years. In line), and percentage of C. davisiana (dash-dot line). Letter designations of peaks on the latter
curve are informal designations of various parts of the record.
addition, visual inspection of the core
shows that it has been mechanically
stretched between 300 and 400 cm during
the coring process. Consequently, the
core was analyzed for 8180 only from the
lower part of stage 5 to the base. The 0
record of stages 6 to 9 is similar to the Go
equivalent part of core RC 11-120; the
remainder can be compared stage by
stage with other cores (31, 42). As in 03
other cores studied, stage 12 is bracketed 2-
by the extinction of the coccolith Pseu-
doemiliania lacunosa below (43) and the
radiolarian Stylatractus universus above
(44). This confirms the presence of the
entire oxygen isotope sequence from
stage 6 to stage 13.
The Ts curve in both cores differs 0
0co
markedly in shape from the 8180 curve,
being characterized by abrupt increases
of estimated temperature of up to 6°C in
stage 1 and at the base of stages 5, 7, 9, 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
and 11. Elsewhere in the cores Ts fluctua- Depth (m)
tions do not exceed 3°C. The spiky char- Fig. 3. Depth plots of four parameters measured in core E49-18: 8180 (solid line at the top), Ts
acter of this record is similar to that of (dashed line), percentage of C. davisiana (dash-dot line), and percentage of CaCO3 (solid line at
certain Atlantic cores from high northern the bottom). The technique used for CaCO3 measurement is that of Hulsemann (81). A com-
parison of the lettered intervals of the C. davisiana curve for this core with those for core RC 11-
latitudes (45). 120 (Fig. 2) shows that the time represented by the top 1.5 m of RC1 1-120 is not present in E49-
Above stage 10 there is near synchro- 18.
10 DECEMBER 1976 1123
Table 2. Chronologic assumptions of age models. Interpolation within and exttrapolation Although the percentage C. davisiana
beyond control points shown is linear. For the combined PATCH ELBOW and PATC H TUNE- curve has a character distinct from the
UP records, data from 0 to 785 cm in RC1 1-120 were combined with data below 825 cm in
E49-18. other two, its maxima are generally cor-
related in timing, but not in amplitude,
Sedimen- with Ts minima and 8180 maxima, except
Age C Depth Age tation rate in stages 8 and 9.
ore (CM/1O3
model (cm) (x 103 years) Time control. A basic chronological
framework for these sequences is estab-
SIMPLEX RC 11- 120 0 0 3.46 lished by determining the absolute ages
440 127*
E49-18 490 127* 2.92 of certain horizons. In RC1 1-120, car-
1405 440t bon-14 dating at the 36- to 39-cm level
ELBOW RC11-120 0 0 yields an age of 9400 + 600 years (35).
39 9.4t 4.14 This level marks the most recent Ts maxi-
440 127* 3.40 mum and substantially precedes the
785 251t 2.78 Northern Hemisphere hypsithermal,
E49-18 490 127*
825 251t 2.70 which at many sites has been dated at
1405 440t 3.06 about 6000 years ago (1).
TUNE-UP RC11-120 0 0
39 9.4t 4.14 The age of the boundary between
440 127* 340 stage 12 and stage 11 was taken from
785 247 2.87 Shackleton and Opdyke (31), who esti-
E49-18 490 127* mated it at 440,000 years in an equatorial
825 247 2.79 Pacific core (V28-238) by assuming uni-
1405 425 3.26 form accumulation between the core top
*Age of isotopic stage 6-5 boundary (17). tAge of isotopic stage 8-7 boundary (251,000 years) and bound- and the magnetic reversal marking the
ary 12-11 (440,000 years) (31). tCarbon-14 determination (3S). Brunhes-Matuyama boundary. Extinc-
tion of S. universus occurred globally on
this stage boundary (44), so that the esti-
Table 3. Frequency-domain test of orbital theory based on SIMPLEX chronology for two deep- mates of 400,000 years for the age of this
sea cores. Values are mean periods (in thousand years per cycle) of peaks in unprtewhitened extinction (46) in the North Pacific and
geologic and orbital spectra.
Antarctic constitute independent deter-
Geologic data Orbital d ata minations for the age of the 12-11 bound-
Time Fre- C. ary. The range of these figures expresses
Core interval quency Ts 8180 davisiana Spectral E lement the current age uncertainty of this bound-
(x 103 B.P.) band (%) estimate ary (47).
a 87 91 106 In many areas there is evidence for a
RC11-120 0-273* b 38 38 37 40.8 Oblicquity (e) change in accumulation rate around
c 21 23 22.6 Prec ession (P) stage 6 (48) so that we have used an
b/c 1.8 1.7 1.8 e/P independent estimate of 251,000 years
a 94 109 119
E49-18 127-489t b 43 47 41.1 Oblicquity (E) for the stage 8-7- boundary. This esti-
c 24 24 21.9 Prec ession (P) mate, like that for the 12-11 boundary,
b/c 1.8 1.9 1.9 e/P was taken from Pacific core V28-238 (31).
*Geologic and orbital spectra for this interval were calculated with n = 91 and m = 40 (57). tGeologxc The age of the stage 6-5 boundary is
and orbital spectra for this interval were calculated with n = 121 and m = 50 (S7). within the range of several radiometric
dating techniques, and has been the cor-
nerstone of some previous attempts to
Table 4. Frequency-domain test of orbital theory using ELBOW chronology for PAT'CH core. support specific versions of the orbital
Values are mean periods (in thousand years per cycle), of peaks and subpeaks in gec)logic and hypothesis. We have used an age of
orbital spectra [n = 163 and m = 50 (57)]. Orbital data calculations cover the passt 468,000 127,000 years which has an analytical
years. error estimated at + 6000 years (17). At
Geologic data Orbital data this point in our analysis we do not exper-
iment with other published ages (49) be-
Fre- Unprewhitened spectra Prewhitened spectra cause it is very difficult to reconcile all
quency C. C. Spectr domain Ele- the terrace coral ages from Barbados,
band Ts 8180 davisiana Ts 8180 davisiana estimate estimate ment New Guinea, and Hawaii (15, 23, 24, 50)
(%) (%o) and data from deep-sea cores (51) with
a 94 106 122 105 97 Eccentri- any substantially different age (52).
city (e) Chronological models. To test the or-
b 40 43 43 42 43* 42* 41.1 40.6 0bliquity bital theory, age models (Table 2) must
(E) be developed to express each geological
cl 23 24 24 24* 24 24 23.1 P
sicens(1) variable as a function of time. We do this
cm 22* 22 22 21.8 21.6 P reces- by assuming that sediment accumulated
sion (Pm) at a constant rate between the horizons
C2 19.5 19.5 19.5 19.5 18.8 Pteces- for which we have independent esti-
sion (P2)
b/c1 1.7 1.8 1.8 1.8 1.8 1.7 1.78 dE mates of age. Although uniform sedimen-
tation is an ideal which is unlikely to
*Peaks in prewhitened spectra are significant at P = .05. prevail precisely anywhere, the fact that
1124 SCIENCE, VOL. 194
the characteristics of the oxygen isotope 30 - A Orbital data B f-23 K Insolation C Insolation
record are present throughout the cores 25 41 K (past 468,000 15 550S 60°N
25- (23K yea rs) Winter Summer
suggests that there can be no substantial (past 468,000 10t (past 468,000
lacunae, while the striking resemblance 20- -( years) years)
Bandwidth 10
to records from distant areas shows that c ,19 K (f) R(f)
15-
there can be no gross distortions of accu-
mulation rate. 411
10 H,g9
10 19 K
5
SIMPLEX models assume uniform Obliquity
sedimentation rates through each core 5 --Precession
and fix the core top of RC ll- 120 at zero, 0- O0 0-
the 6-5 stage boundary at 127,000 years f 0 .033 .067 .100 .133 .167 0 .033 .067 .100 .133 .167 0 .033 .067 .100 .133 .167 f
(cores RC I 1-120 and E49-18), and the 12- 1/f 100 30 15 10 7.5 6 100 30 15 10 7.5 6 100 30 15 10 7.5 6 1/f
11 stage boundary at 440,000 years (core Frequency (cycles/1000 years)
E49-18; see Table 2). It is important to Fig. 4. High-resolution spectra of orbital and insolation variations over the past 468,000 years.
note that although the SIMPLEX model Variance (as percentage of total variance per unit frequency band) is plotted as a function of
presents each variable as a function of frequency (cycles per thousand years). Arrows indicate weighted mean cycle lengths of spectral
peaks (in thousands of years). (A) Spectra for obliquity and precession (AesinfI). (B) Spectrum
age, the age estimated is an exact linear for winter insolation at 55'S. (C) Spectrum for summer insolation at 60'N. [All data are from
function of depth in the core. Vemekar (39)]
ELBOW models use more chronologi-
cal information and no longer represent
the records as a simple function of depth quencies according to a gain function considered unimportant (59). Obliquity
(Table 2). Sedimentation rates change (55). Therefore, whatever frequencies (E, the angle between the equatorial and
slightly at three control points and are characterize the orbital signals, we will ecliptic planes) ranges from 22.10 to
assumed constant between them. expect to find them emphasized in pa- 24.50, with an average period of about
PATCH is a time series in which we leoclimatic spectra (except for fre- 41,000 years (39).
have joined the records of the two cores quencies so high that they would be The climatic effect of precession is a
at the 8-7 stage boundary, thereby pro- greatly attenuated by the time constants function of HI, the longitude of perihelion
viding a longer and statistically more of response). based on the moving equinox, and e (60).
useful record. Another model in which Numerical procedures. The statistical Specifically, the intensity of incoming
the cores were joined at the 6-5 stage techniques we use to calculate both orbit- solar radiation at a particular latitude and
boundary was found to be statistically al and climatic spectra are taken without season varies as esinfl. To express these
almost identical and is not considered modification from the work of Blackman variations as a time series, the value of
here. and Tukey (56) with elaborations by Jen- esinfl for June 1950 A.D. is subtracted
kins and Watts (55). Our procedure in- from the same quantity calculated for
volves ten sequential steps: selection of particular times in the past. The resulting
Frequency-Domain Tests an absolute chronology, calculation of a precessional index (Aesinfl) is approxi-
time series, detrending, prewhitening mately equal to the deviation from its
Assumptions. From the viewpoint of (optional), lagging, calculation of the 1950 value of the earth-sun distance in
linear-systems modeling (53) the astro- autocovariance function, smoothing with June, expressed as a fraction of the semi-
nomical theory of Milankovitch postu- a Hamming lag window, spectral estima- major axis of the earth's orbit (61). Over
lates two systems operating in series. tion, scaling, and statistical evaluation the past 4 million years, this index ranges
The first is a radiation system which (57, 58). In all calculations the sampling from about +0.03 to -0.07 and has an
transforms orbital signals (obliquity and interval is fixed at 3000 years; hence the average period of about 21,000 years
precession) into a set of insolation sig- spectral estimates cover a frequency (39).
nals (one for each combination of lati- band ranging up to the Nyquist fre- We have used spectral techniques. to
tude and season). The insolation signals quency of 0.167 cycle per thousand analyze secular variations in E, e, sinIl,
are transformed by a second, explicitly years. and Aesinll (62). Because these varia-
formulated climate-response system into Frequency analysis of astronomical tions are quasi-periodic, it.is necessary
(predicted) climate curves. In contrast, data. For any particular latitude the in- to specify the interval over which fre-
we postulate a single, radiation-climate tensity of solar radiation received at the quencies are to be analyzed (63); our
system which transforms orbital inputs top of the atmosphere varies quasi-peri- analyses (Tables 3 and 4 and Fig. 4) treat
into climatic outputs. We can therefore odically with three elements of the time intervals corresponding to the core
avoid the obligation of identifying the earth's orbit: eccentricity, obliquity, and records.
physical mechanisms of climatic re- the longitude of perihelion based on the Spectra calculated for variations in ec-
sponse and specify the behavior of the moving equinox. Over the last 4 million centricity and obliquity over the past
system Qnly in general terms (54). The years the eccentricity (e, the ratio of the 468,000 years (Table 4) are both uni-
dynamics of our model are fixed by as- focal distance to the length of the major modal; that is, they have spectral peaks
suming that the system is a time-in- axis) ranges from a value near zero to a dominated by a single frequency (64).
variant, linear system-that is, that its maximum of about 0.06, and exhibits an These spectral peaks correspond to cy-
behavior in the time domain can be de- average period of 93,000 years (39). Vari- cles of 105,000 years (eccentricity) and
scribed by a linear differential equation ations in e, unlike variations in the other 41,000 years (obliquity). Spectra calcu-
with constant coefficients. The response orbital elements, also slightly affect the lated for sinfl (not tabulated) and for
of such a system in the frequency do- total annual insolation received by the Aesinrl are more complex. Both are bi-
main is well known: frequencies in the earth. Because the extreme range of this modal, with subpeaks corresponding to
output match those of the input, but their effect over the past 500,000 years is cycles of 23,000 years (P1) and 19,000
amplitudes are modulated at different fre- about 0.1 percent, it has generally been years (P2).
10 DECEMBER 1976 1125
 We have also calculated spectra for estimates of their constituent fre- series for the combined (PATCH) time
two time series recording variations in quencies. Nevertheless, five of the six series. As before, the climatic variance is
insolation (Fig. 4), one for 55SS and the spectra calculated are characterized by distributed mainly in three discrete spec-
other for 60°N (39). To the nearest 1000 three discrete peaks, which occupy the tral peaks (Table 4 and Fig. 5). More
years, the three dominant cycles in these same parts of the frequency range in than half of the total variance is con-
spectra (41,000, 23,000, and 19,000 each spectrum (Table 3). Those corre- tained in the low-frequency a peak (62
years) correspond to those observed in sponding to periods from 87,000 to percent for Ts, 58 percent for 8180, and
the spectra for obliquity and precession. 119,000 years are labeled a; 37,000 to 51 percent for C. davisiana). All three
This result, although expected, under- 47,000 years b; and 21,000 to 24,000 peaks are unimodal. Estimates of the
scores two important points. First, in- years c. This suggests that the b and c dominant cycle in the a peaks are 94,000,
solation spectra are characterized by fre- peaks represent a response to obliquity 106,000, and 122,000 years for Ts, 8180,
quencies reflecting obliquity and pre- and precession variation, respectively. and percentage of C. davisiana, respec-
cession, but not eccentricity. Second, The ratios of obliquity period to pre- tively. Because these estimates are de-
the relative importance of the insolation cession period calculated for the inter- rived from the low-frequency end of the
components due to obliquity and pre- vals of time covered by cores RC11-120 spectrum, they should be regarded as
cession varies with latitude and season. and E49-18 are 1.8 and 1.9, respectively rough estimates. However, there can be
Frequency analysis of geological (Table 3). The observed ratios of the no doubt that a spectral peak centered
data. Using techniques identical to those dominant period in peak b to the domi- near a 100,000-year cycle is a major fea-
applied to the astronomical data, we nant period in peak c are i.7 + 0.1 and ture of the climatic record.
have calculated spectra for each of the 1.9 ± 0.1 for RC11-120 and E49-18, re- A substantial fraction of the variance
three geological variables: Ts, 8180, and spectively. Because the observed ratios is contained in the b peaks: 19 percent
percentage of C. davisiana (Fig. 5). The from these geological data are indepen- for Ts, 27 percent for 8180, and 30 per-
SIMPLEX chronology is used to analyze dent of the ages used to calibrate the cent for C. davisiana (Fig. 5). Each of
the individual cores, and the ELBOW SIMPLEX time series and closely match the three peaks is unimodal. Estimates of
chronology is applied to the combined ratios derived for orbital data, we con- this dominant cycle in the peaks range
(PATCH) core. clude that the ratio test supports the from 40,000 to 43,000 years (Table 3).
Because the SIMPLEX time series are hypothesis of orbital control. A smaller fraction of the total variance
undesirably short and are based on limit- To obtain accurate estimates of the is contained in the c peaks: 11 percent
ed chronological control, we do not climatic frequencies, we now examine for Ts, 9 percent for 8180, and 7 percent
place much reliance on the accuracy of spectra calculated from ELBOW time for C. davisiana. The Ts and C. davisia-
na peaks, both unimodal, correspond to
cycles of 23,000 and 24,000 years, respec-
Ts 8180 % C. davisiana tively (Table 4). The c peak in the 8180
LoJ
-J
9-
R(f ) Bandwidth g(f) R(f) spectrum is bimodal, with subpeaks (cl
and c2) corresponding to cycles of 24,000
0-
(%
6- and 19,500 years (Table 4).
Although the statistical significance of
a 3
ML the a peaks can hardly be in doubt, for
peaks b and c a more detailed examina-
tion and statistical evaluation are clearly
12- |
desirable. These objectives are achieved
x '
(f) ~ R(f) R(f)
by expressing the variance on a log scale

I
-J
cn (Fig. 6, top) and by prewhitening the
4, signal to reduce distortions in the higher-
u
L -
0 - frequency part of the spectrum caused
by variance transfer from the a peak
during analysis. The resulting spectra,
-J a. which have the flat trend desired (Fig. 6,
0 6-
bottom), should be used in the higher-
12 | R(f) R(f) R(f) frequency part of the spectrum, not only
3-
to conduct statistical tests but also to
0
I
.. 9- obtain more accurate frequency esti-
km
n0*- mates.
Estimates of b peak frequencies in the
I-
C.) prewhitened signals are changed very
little; their dominant cycles now range
v-
f 0
'~ ~ ~ ~~T'
.033 .067 .100 .133 .167 0
l'
.033 .067 .100 .133 .167
IL
0 .033 .067 .100 .133 .167 f
from 42,000 to 43,000 years (Table 4).
The c peaks differ mainly in the appear-
1/f 100 30 15 10 7.5 6 100 30 15 10 7.5 6 100 30 15 10 7.5 6 1/f ance of a small subsidiary peak in the Ts
Frequency (cycles/1000 years) spectrum. The midpoints of the three cl
Fig. 5. High-resolution spectra of climatic variations in Ts, 8180, and percentage of C. davis- peaks now correspond to cycles 24,000
iana. Variance (as percentage of total variance per unit frequency band) is plotted as a function years long; and the midpoints of the c2
of frequency (cycles per thousand years). Prominent spectral peaks are labeled a, b, and c. Ar- peaks in the Ts and 8180 spectra corre-
rows indicate weighted mean cycle lengths (in thousands of years). The age models used in the
calculations are given in Table 2. (A) Spectra for core RC1 1-120 are calculated for the SIM- spond to a 19,500-year cycle (Fig. 6 and
PLEX age model. (B) Spectra for core E49-18 are calculated for the SIMPLEX age model. (C) Table 4).
Spectra of the combined (PATCH) record are calculated for the ELBOW age model. Statistical evaluation. Based on the
1126 SCIENCE, VOL. 194
 null hypothesis that the data are a sample the assumption of constant accumulation tions differs from ours by about a factor
ofa random signal having a general distri- rates. of 2. The explanation of this paradox is
bution of variance like that in the ob- Having found the frequencies of obliq- to be found in the dominant climatic
served low-resolution spectrum, con- uity and precession in our geological rec- periodicity, which in all of the cores is
fidence intervals are calculated as a ords-as predicted by a linear version of the 100,000-year cycle, and not, as ex-
guide to the statistical significance of the theory of orbital control-we should pected (2, 26), the geologic response to
spectral estimates in the high-resolution consider again the lower-frequency cli- the 41,000-year obliquity cycle (Table 5).
spectrum (Fig. 6). A particular peak in matic components which, although not The spectral peak identified by Kemp
that spectrum is judged significant if it predicted, actually contain about half of and Eger (16) and by Chappell (20) as
extends above the low-resolution spec- the observed variance. These com- due to precession, corresponds to our b
trum by an amount that exceeds the ap- ponents form the a peaks in our spectra. peak, and therefore (we argue) is ac-
propriate one-sided confidence interval. Concentrated at periods near 100,000 tually an effect of obliquity. The spectral
Of the three b peaks, one is significant at years (Figs. 5 and 6), they are close to peak identified by van den Huevel (14) as
P = .02 (C. davisiana) and one at the 105,000-year period in the eccentric- due to the precession half-cycle corre-
P = .05 (8180). Of the three c1 peaks, ity spectrum (Table 4). This similarity of sponds to our c peak, and therefore is
one (Ts) is significant at P = .05 (65). the dominant frequencies in late Quater- now to be understood as the effect of a
We carefully considered aliasing and nary records of climate and eccentricity full precession cycle. Only with the ad-
harmonic problems (66) and conclude has been noted before (13, 17, 18) and vent of chronologies based on the
that our spectral peaks are not an artifact demands an explanation. One hypothesis Brunhes-Matuyama magnetic reversal
of procedure. This conclusion was sup- (developed further below) is that the radi- (17, 18, 31) was the dominant climatic
ported by examining more detailed ation-climate system responds nonlinear- period in the 8180 record identified as
spectra calculated by the maximum en- ly (68) to changes in the geographic and approximately 100,000 years (Table 5).
tropy technique (67). seasonal distribution of insolation. An-
Discussion. Having examined the con- other is that the small control eccentric-
fidence intervals for individual climatic ity exerts on the total annual insolation is Time-Domain Tests
spectra, and having eliminated the alias- significant climatically (59).
ing and harmonic problems, we can now An apparently independent confirma- Assumptions. As with the frequency-
ask if the frequencies found in the three tion of our conclusions about spectral domain test, we start here with the as-
spectra are those predicted by our linear peaks b and c can be found in reports of sumption that the radiation-climate sys-
version of orbital theory. In making this climatic periodicities in the 8180 record tem is time-invariant and linear. One
frequency-domain test we must note that closely matching those of obliquity and well-known characteristic of such a sys-
the geologic spectra contain substantial precession (14, 16, 20, 26, 69). However, tem forms the basis for our time-domain
variance components at many fre- the time scale used in these investiga- tests: any frequency component of the
quencies in the range of interest, and not
simply ask whether there are statistically
significant frequencies which match
those predicted but whether the spectra
observed show sufficient emphasis at the
frequencies predicted to be accepted as
nonrandom results (65). Our answer is
"yes" for the following reasons: (i) Us-
ing a chronology that is completely inde- o

pendent of the astronomical theory, we U 1.

find that modal frequencies observed in -J c

the geological record match the obliquity to

and precession frequencies to within 5 b-3

oo
percent. The coherence of these results
across different variables measured in
ZCI
two cores we regard as very powerful 0
support for the theory. (ii) Two geologi-
cal spectra (8180 and percentage of C. CIL-

davisiana) have peaks corresponding to

the predicted obliquity frequency that
are significant at P = .05. One geological
spectrum (Ts) has a peak corresponding
to the dominant precession frequency -3 -D-.0 I ; I ' .I 21 i i
which is also significant at P = .05. (iii) .

f 0 .033 .067 .100 .133 .167 0 .033 .067 .100 .133 .167 0 .033 .067 .100 .133 .167 f
In addition, the predicted ratios of obliq- l/f 100 30 15 10 7.5 6 100 30 15 10 7.5 6 100 30 15 10 7.5 6 1/f
uity to precession frequencies (calcu- Frequency (cycles/1000 years)
lated for the time intervals represented Fig. 6. Spectra of climatic variations (in Ts, 8180, and percentage of C. davisiana) in the com-
bined (PATCH) record from two subantarctic deep-sea cores. Calculations are based on the
by cores RCl 1-120 and E49-18) match ELBOW age model (Table 2). Arrows without crossbars indicate weighted mean cycle lengths
the ratios of measured climatic fre- of spectral peaks (in thousands of years). Arrows with crossbars show one-sided confidence
quencies in the two cores (using the SIM- intervals (C.I.) attached to estimates in the high-resolution spectrum. Prominent spectral peaks
PLEX chronology) to within 5 percent- are labeled a, b, and c. (Top row) High-resolution spectra from Fig. 5C expressed as the natural
log of the variance as a function of frequency (cycles per thousand years). (Bottom row) High-
a result that is independent of absolute resolution spectra (solid line) and low-resolution spectra (dashed line) after prewhitening with a
age specifications and depends only on first-difference ifiter.
10 DECEMBER 1976 1127
input signal will appear at the same fre- Although our predictive model is non- pect that curve to lag Ts and percentage
quency as a component of the output, specific, in the sense that it does not say of C. davisiana at all frequencies.
but exhibit there a certain phase shift. what the shape of a particular climatic Numerical procedures. Discrete spec-
How much that phase shift will be we record should be, it does predict that tral peaks in our geological time series
cannot say, for at any particular fre- whatever orbital-geological phase shift is have already been observed and identi-
quency the magnitude of the phase lag observed at a particular (obliquity or pre- fied as variance components related to
will depend not only on the time constant cession) frequency, that shift should be obliquity and precession. Using a band-
of the system, but also on the exact form constant (19). pass digital filter (Tukey filter), we will
of the system's linear response. If, for Furthermore, we can assume that each now extract these components from the
example, the system has a single-expo- climatic index reflects a response of one signal and display them in the time do-
nential response without delay, then the physical part of the climate system and is main (70). To use this phase-free filter
output will lag behind the input by no characterized by a certain time constant. the investigator chooses f0, the center
more than a quarter of a cycle. But if the Therefore, whatever phase shift is ob- frequency in the band to be studied, and
system does exhibit delay, then the lag served (at a particular frequency) be- fixes the bandwidth by determining r, the
can exceed a quarter-cycle. An addition- tween a pair of climatic indices should resolution of the filter (71). The impact of
al source of uncertainty attaches to the also be constant, and the subsystem with a particular ifiter on the frequency do-
23,000-year (but not the 40,000-year) the larger time constant should lag the main is described by its frequency-re-
component of the orbital input, for the other. Because the 8180 curve reflects sponse function, H(f). We used two fil-
phase of the precession index itself is changes in the cryosphere (a part of the ters, one centered at a frequency of 0.25
arbitrarily defined with respect to a par- climate system that must have longer cycle per 1000 years and one at 0.043
ticular time of the year. time constants than the ocean), we ex- cycle per 1000 years (Fig. 7). These will
be referred to as the 40K (40,000-year)
and 23K (23,000-year) filters, respective-
Table 5. Comparison of published 8180 chronologies and spectra. Dominant periods (in thou- ly. Their bandwidths have been chosen
sands of years) as calculated by investigators cited are correlated with spectral peaks a, b, and c
documented in this article. so that the adjacent half-power points
nearly coincide, and the variance in the
Chronological models Correlation of dominant ifitered signal is approximately the same
periods identified as the variance in the corresponding
Age of 12-11 spectral peak.
boundary Reference a b c Because filters of this type yield
(x 10 years) smooth curves no matter what data are
220* Emiliani (26) 40 processed (72), the objective in using
220 van den Heuvel (14)t 40 13 them is simply to determine in the time
240t Kemp and Eger (16) - 52 - 27 domain, and for each frequency com-
270§ Emiliani (27) 50 ponent of interest, the position of in-
270 Chappell (20)11 - 46 - 25 flection points in the filtered record. Al-
350 Imbrie and Kipp (18)1 80 30
375 Broecker and van Donk (17) 90 though the average interval between in-
380 Pisias (82)# 23 flection points must match some
This article: frequency within the passband of the ifl-
425 TUNE-UP age 99 42 22 ter, information in the data controls the
440 ELBOW age 106 43 22
exact timing of each inflection. Phase
*Age calculated on the assumption that 8180 stages reflect obliquity. tData from Emiliani (68).
selected to match spectral peaks calculated from Emiliani's (27) 8'8O data with orbital frequencies.
*Age shifts can be estimated visually by com-
§Age
calculated by extrapolation beyond a section of the curve estimated to be younger than 150,000 years by early paring two geological variables filtered at
results of the Parh technique. IlData from Emiliani (27). 1Data for 80 from Broecker and van Donk the same frequency, or by comparing an
(17); an unconformity was later recognized in this record (78). #Data for 8180 from Shackleton and Op-
dyke (31). orbital curve with a geological curve fil-
tered at the same frequency. Cross-spec-
tral techniques are employed on unfil-
1.2 H(f) 40K filter A H(f) 23K filter B tered data to obtain for the frequency of
r 18 r 15 interest a quantitative estimate of the av-
fo. 025 fo=.043 erage phase shift over the study interval
0.9
(73).
Patterns in the geologic record. Two
sets of curves (Fig. 8, A and C) show the
results of applying filters to the ELBOW
0.6
time series. The three geological curves
filtered at 40K are approximately in
phase throughout their length. Cross-
spectral analysis shows 8180 lagging Ts
by 2000 years and percentage of C. davis-
iana by 1000 years. The data filtered at
23K show a nearly constant phase rela-
0 .033 .067 .100 .133 .167 0 .033 .067 .100 .133 .167 tionship between 8180 and percentage of
Frequency (cycles/1000 years) C. davisiana throughout the record, with
Fig. 7. Gain function for band-pass filters used to calculate curves in Figs. 8 and 9: (A) 40K ifiter 8180 lagging percentage of C. davisiana
centered on a frequency of 0.025 cycle per thousand years; (B) 23K filter centered on a fre- by an average of 4000 years. After
quency of 0.043 cycle per thousand years. Tukey filters were used (70). 350,000 years ago 8180 systematically
1128 SCIENCE, VOL. 194
 lags Ts by an average of about 2000 Discussion. We regard the results of magnitude of the phase shift (7000 to
years. Before that time, however, 8180 the time-domain test as strong evidence 9000 years) is less than a quarter of the
and Ts are clearly out of phase. oforbital control of major Pleistocene cli- cycle length. (ii) Averaged over the same
The general pattern of these relation- matic changes, for the following reasons: interval, the 23K components of the geo-
ships can be observed on the strati- (i) Over the past 300,000 years, each of logical curves exhibit in-phase relation-
graphic diagrams (Figs. 2 and 3) and is in- the 40K components of the geological ships with precession which are as con-
dependent of assumptions about chronol- records exhibits a phase relationship stant as could be expected from a geo-
ogy. Therefore, the climate changes of with obliquity which is as constant as logical record. Again, the observed
the two hemispheres are nearly in phase could be expected from a geological rec- regularity is too great to be explained as
during the last 300,000 years with ord. Monte Carlo tests we have con- a random result. When the chronology is
changes in the Southern Hemisphere ducted with our filters show that the de- most accurate for the filtered record
ocean appearing to lead changes in the gree of regularity observed would be (50,000 to 150,000 years ago), the 23K
Northern Hemisphere ice sheets by up to highly unlikely as a random result. The components of the geological curves lag
a few thousand years. However, before
300,000 years ago Ts and 8180 appear to
be out of phase in the 23K but not the I_
40K frequency band. 23K
Orbital-geologic phase relationships.
Displayed on the ELBOW chronology,
phase relationships between the filtered -0.15._
geological curves and orbital curves are 0 E
systematic over the past 300,000 years
(Fig. 8, A and C); that is, back to the up- +0.1 X
per part of stage 9. Over that interval, .2 1- +0.5T-
U) r
times of low temperature, high 8180 ra- C.)
L-n4
tios, and abundant C. davisiana follow o '--
W
times of low obliquity. The 40K com-
ponents of the geological curves (Fig.
8C) lag obliquity by about a quarter- - 1.0
cycle. Measurements made at the maxi-
ma and minima of 12 obliquity half-cy-
cles covering the interval 70,000 to
300,000 years ago show that 8180, Ts, Go
and percentage of C. davisiana lag obliq- -1 %
uity by 9000 + 3000, 8000 + 3000, and
7000 + 4000 years, respectively. Below
the top of stage 9, however, the geologi-
cal curves on the average lead obliquity
by several thousand years.
Over the interval 0 to 150,000 years -0.1 --
0
.
ago on the ELBOW time scale, the inter- E
val that has the most certain chronology,
+0.1 ff
the 23,000-year components of the geo-
logical curves systematically lag pre- ._
+0.5 ,c
cession by about 3000 years (Fig. 8A). 0
._
However, when averaged over the inter- 0 -0.5
val 0 to 300,000 years they are approxi-
mately in phase with precession. Mea-
surements made at the maxima and mini-
ma of 22 precession half-cycles covering
the interval 50,000 to 300,000 years ago
show that 8180, Ts, and percentage
of C. davisiana lag precession by
1500 ± 3500, 0 ± 3000, and-500 + 4500
years, respectively. Times of low temper-
ature, high 8180 ratios, and abundant C.
davisiana are associated with times of 0 50 100 150 200 250 300 350 400 450 500
high positive values of the precession in- Years x 103 B.P.
dex-that is, times when the earth-sun Fig. 8. Variations in obliquity, precession, and the corresponding frequency components of cli-
distance in June is greater than normal. mate over the past 500,000 years. Orbital data are from Vernekar (39). Climatic curves are varia-
The systematic orbital-geologic phase tions in 8180, Ts, and percentage of C. davisiana plotted against alternate geological time scales
relationships just described do not, how- (ELBOW and TUNE-UP) as defined in Table 2. The variations shown are frequency com-
ever, exist below the 300,000-year hori- ponents extracted from the raw-data curves by means of digital band-pass filters (Fig. 7). The
two sets of curves in (A) and (B) include the precession curve and the 23,000-year frequency
zon. There, a confusing pattern of leads components of climate based on the ELBOW (A) and TUNE-UP (B) time scales. The two sets
and lags among all of the curves is dis- of curves in (C) and (D) include the obliquity curve and the 40,000-year frequency components
played. of climate based on the ELBOW (C) and TUNE-UP (D) time scales.
10 DECEMBER 1976 1129
precession by about 3000 years. (iii) The Implications of Test Results time domain (Fig. 8, B and D) is to ex-
expectation that the 23K and 40K com- tend the systematic phase relationships
ponents of the 8180 curve should show a Quaternary time scale. Takern together previously noted with obliquity back
constant lag against the corresponding with the results of the frequenc y-domain over the whole range of the filtered rec-
components of Ts and percentage of C. test, the systematic phase relaLtionships ord-some 425,000 years. There is also a
davisiana is confirmed by the geological just described suggest that a srnall error definite improvement in the phase rela-
record of the last 350,000 years. The out- occurs in the older portions off the EL- tionship with precession. For the 8180
of-phase relationship between Ts and BOW chronology. To explore this hy- record, this relationship is generally regu-
8180 that appears in stage 10 and the up- pothesis we have made the minLimum ad- lar and includes 13 cycles extending back
per part of stage 11 may be related to the justments in the ELBOW chironology 340,000 years. For Ts the relationship ex-
low amplitudes of precession variation which would extend farther bacAk in time tends back 325,000 years, and for the per-
resulting from near-zero eccentricity at the systematic phase relationw ships ob- centage of C. davisiana back 280,000
this time. This may have allowed a de- served in the younger parts off the EL- years (with one exception near 240,000
coupling of Ts variation from precession BOW record. These adjustmenits (Table years). We have experimented with oth-
over this short interval. (iv) The fact that 2) are easily within the absolut(e error of er ages for our chronological control
the large irregularities which do occur in the radiometric dates on whichi the EL- points and find that further adjustments
the observed orbital-geologic phase rela- BOW time scale is based. On thiis revised result in a deterioration of the match be-
tionships at both frequencies are strati- chronology (called TUNE-UP)>,the age tween orbital parameters and the geologi-
graphically concentrated in the early part of the isotope stage 11-12 boundlary is re- cal records filtered at the same fre-
of the record (before 300,000 years ago) duced by 3 percent (25,000 years), and quencies (74). In the frequency domain,
where the chronology is least accurate, the age of isotope stage 7-8 boiundary is estimates for spectral peaks calculated
rather than randomly distributed, sug- reduced by 2 percent (4000 yearrs) (Table from the TUNE-UP time series match
gests that these irregularities result from 6). those predicted for obliquity and pre-
errors in the older chronology. The impact of these adjustme nts in the cession within 1000 years (Table 5).
The 100,000-year cycle. The dominant
cycles in all of our spectra (Fig. 5C and
Table 4) are about 100,000 years long-
3 an observation which merely confirms a
geological opinion now widespread (17,
18, 75-77). Yet this cycle would not arise
as a linear response of the climate sys-
tem to variations in obliquity and pre-
a, .04 >, cession. Mesolella et al. (15) and Broeck-
h..
F , er and van Donk (17) account for this
0 -.02 c cycle by relating it to eccentricity or to
0
LO
F the phasing of obliquity with precession.
° Our data, displayed on a time scale de-
rived without reference to eccentricity,
dramatically confirm the empirical asso-
ciation of glacial times with intervals of
0 low eccentricity (Fig. 9).
0 +1.5 Because this time-domain match
VI,(A TUNE-UP Ts '

O- ,\ ", ,,
,II agrees with independent evidence in the
A.,~~~~~~~~~~~~~~~~~~~~~~~~~~~~.
frequency domain (8180 and Ts in Table
-1.5 -
4), we conclude, as others have (15, 18),
10.5- . .! .04 >' that the 100,000-year climate cycle is
/
- A 2 driven in some way by changes in orbital
8.0 ,'
-.02 @ eccentricity. As before, we avoid the ob-
5.5- 0 (, ligation of identifying the physical mech-
3.0- anism of this response, and instead char-
+1 .0 A
acterize the behavior of the system only
,

o- \\ ,
A, ,w
~ ~ ,
\ in general terms. Specifically, we aban-
-1.0- A~ ~ ~ ~ ~ ~ ~ ~ ~~~\
don the assumption of linearity (78). In
such a nonlinear system (68) there are
0 50 1oo 150 200 250 300 350 400 450 500 many ways in which the modulating ef-
Years x 103 B.P. fect of eccentricity on the precession in-
dex could generate 100,000-year vari-
Fig. 9. Variations in eccentricity and climate over the past 500,000 years. Climatic 4curves are ance components in the geological rec-
obtained from the combined (PATCH) record of two subantarctic deep-sea cores and
the TUNE-UP time scale (Table 2). (A) Solid line in center shows variations in 8180. 1)otted line ord. Among these we are attracted to a
is a plot of orbital eccentricity (39). Upper curve is the 23,000-year frequency comlponent ex- hypothesis that ice sheets waste faster
tracted from 8180 by a band-pass digital filter (Fig. 6). Lower curve is the 40,000-year frequency than they grow; that is, that two different
component extracted from 8180 by a band-pass digital filter (Fig. 6). (B) Dashed lin e in center time constants of the cryosphere's re-
shows variations in estimated sea-surface temperature (Ts). Dotted line is a plot of orb
tricity data from Vernekar (39). Upper curve is the 23,000-year frequency component
from Ts by a band-pass digital filter (Fig. 6). Lower curve is the 40,000-year frequc corn-
eitraclteend
tncy
sponse to orbital forcing are involved.
This concept (79) can be deduced glacio-
ponent extracted from Ts by a band-pass digital filter (Fig. 6). logically (7) or arrived at inductively
1130 SCIENCE, VOL. 194
 from the climatic record (17). Because Table 6. Estimates of the ages of stage bound- 12. E. N. Lorenz, Meteorol. Monogr. 8, 1 (1968).
13. W. S. Broecker, Science 151, 299 (1966);
the amplitude of each precession-index aries based on the TUNE-UP chronology. Ex- D. L. Thurber, J. Goddard, T.-L. Ku,
cycle is proportional to eccentricity, cept for the lowest two boundaries, the esti- R. K. Matthews, K. J. Mesolella, ibid. 1S9, 297
mates are considered accurate within a range (1968).
such a nonlinear response of the ice of +5000 to -1000 years (74). 14. E. P. J. van den Heuvel, Geophys. J. R. Astron.
sheets would bring out the 100,000-year Soc. 11, 323 (1966).
15. K. J. Mesolella, R. K. Matthews, W. S. Broeck
eccentricity signal in the geologic record tIopic Depth in core (cm) Age er, D. L. Thurber, J. Geol. 77, 250 (1969).
16. W. C. Kemp and D. T. Eger, J. Geophys. Res.
by forcing the mean of the 23,000-year stage (x 103 72, 739 (1967).
climatic cycle to approach values direct- bound- RCl 1-120 E49-18 years) 17. W. S. Broecker and J. van Donk, Rev. Geophys.
Space Phys. 8, 169 (1970).
ly proportional to eccentricity. dary 18. J. Imbrie and N. G. Kipp, in The Late Cenozoic
Future climate. Having presented evi- Glacial Ages, K. K. Turekian, Ed. (Yale Univ.
2-1 40 10 Press, New Haven, Conn., 1971), pp. 71-182.
dence that major changes in past climate 3-2 105 29 19. G. Kukla, Boreas 1, 63 (1972); Nature (London)
were associated with variations in the ge- 2S3, 600 (1975).
4-3 215 61 20. J. Chappell, Quat. Res. (N. Y.) 3, 221 (1973).
ometry of the earth's orbit, we should be 5-4 255 73 21. W. Koppen and A. Wegener, Die Klimate der
6-5 440 490 127 Geologischen Vorzeit (Berlin, 1924).
able to predict the trend of future cli- 7-6 620 640 190 22. N. J. Shackleton, The Phanerozoic Time-Scale
mate. Such forecasts must be qualified in (Geological Society, London, 1971), part 1, pp.
8-7 785 825 247 35-38.
two ways. First, they apply only to the 9-8 900 920 276 23. A. Bloom, W. S. Broecker, J. Chappell, R. K.
natural component of future climatic 10-9 1115 336 Matthews, K. J. Mesolella, Quat. Res. (N. Y.) 4,
185 (1974).
trends-and not to such anthropogenic 11-10 1180 356 24. T.-L. Ku, M. N. Kimmel, W. H. Easton, and T.
effects as those due to the burning of fos- 12-11 1405 - 425 J. O'Neil [Science 183, 959 (1974)] document
13-12 1510 - 457 only the oldest of the terraces.
sil fuels. Second, they describe only the 25. M. L. Bender, F. T. Taylor, R. K. Matthews,
Quat. Res. (N.Y.) 3, 142 (1973).
long-term trends, because they are 26. C. Emiliani, J. Geol. 63, 538 (1955).
linked to orbital variations with periods 27. , ibid. 74, 109 (1966).
ic component has an average period 28. More than a century ago, J. Croll [Climate and
of 20,000 years and longer. Climatic os- Time (Appleton, New York, 1875)] of Scotland
cillations at higher frequencies are not close to, and is in phase with, orbital ec- employed the same basic strategy. He compared
astronomical calculations of orbltal history with
predicted. centricity. Unlike the correlations be- the geologic record of climate, and pointed to
tween climate and the higher-frequency evidence of multiple glaciations as confirming
One approach to forecasting the natu- the astronomical theory of the ice ages.
ral long-term climate trend is to estimate orbital variations (which can be ex- 29. W. L. Donn and D. M. Shaw, Science 157, 722
plained on the assumption that the cli- (1966); J. M. Suarez and I. M. Held, inProceed-
the time constants of response necessary ings of the WMOIIAMAP Symposium on Long
to explain the observed phase relation- mate system responds linearly to orbital Term Climatic Fluctuations (World Meteor-
forcing), an explanation of the correla- ological Organization, Geneva, Switzerland,
ships between orbital variation and cli- 1975).
matic change, and then to use those time tion between climate and eccentricity 30. The raw data on which this study is based will be
published (J. D. Hays, A. D. Vemekar, J. Im-
constapnts in an exponential-response probably requires an assumption of non- brie, N. J. Shackleton, in preparation).
model. When such a model is applied to linearity. 31. N. J. Shackleton and N. D. Opdyke, Quat. Res.
(N. Y.) 3, 39 (1973).
Vernekar's (39) astronomical projec- 6) It is concluded that changes in the 32. N. J. Shackleton, Colloq. Int. CNRS 219, 203
earth's orbital geometry are the funda- (1974).
tions, the results indicate that the long- 33. , Nature (London) 215, 15 (1967); W.
term trend over the next 20,000 years is mental cause ofthe succession of Quater- Dansgaard and H. Tauber, Science 166, 449
nary ice ages. (1969).
toward extensive Northern Hemisphere 34. J. Lozano and J. D. Hays, Geol. Soc. Am. Mem.
glaciation and cooler climate (80). 7) A model of future climate based on 145 (1976), pp. 303-336.
35. J. D. Hays, J. Lozano, N. Shackleton, G. Ir-
the observed orbital-climate relation- ving, ibid., pp. 337-372.
ships, but ignoring anthropogenic ef- 36. M. G. Petrushevskaya, in Biological Reports of
the Soviet Antarctic Expedition (1955-1958), A.
Summary fects, predicts that the long-term trend P. Andriyashev and P. V. Ushakov, Eds. (Lenin-
over the next sevem thousand years is gradskoe Otdelenie, Leningrad, 1967; Israel Pro-
gram for Scientific Translations, Jerusalem,
1) Three indices of global climate have toward extensive Northern Hemisphere 1968), p. 2.
glaciation. 37. J. H. Robertson, thesis, Columbia University
been monitored in the record of the past (1975).
450,000 years in Southern Hemisphere 38. A. J. van Woerkom, Climatic Change (Harvard
References and Notes Univ. Press, Cambridge, Mass., 1953), pp. 147-
ocean-floor sediments. 157; S. G. Sharaf and N. A. Budnikova, Tr. Inst.
1. R. F. Flint, Glacial and Quaternary Geology Teor. Astron. 11, 231 (1967) (translated by the
2) Over the frequency range 10-4 to (Wiley, New York, 1971); J. M. Mitchell, Jr., in Clearinghouse for Federal Scientific and Techni-
10-5 cycle per year, climatic variance of The Quaternary of the United States, H. E.
Wright, Jr., and D. G. Frey, Eds. (Princeton
cal Information, Springfield, Va.).
39. A. D. Vernekar, Meteorol. Monogr. 12 (1972).
these records is concentrated in three dis- Univ. Press, Princeton, N.J., 1965), pp. 881- 40. A. L. Berger, Astron. Astrophys., in press.
crete spectral peaks at periods of 23,000, 901. 41. N. J. Shackleton and N. D. Opdyke, Geol. Soc.
2. C. Emiliani and J. Geiss, Geol. Rundsch. 46, 576 Am. Mem. 145 (1976), pp. 449464; CLIMAP
42,000, and approximately 100,000 (1957). Project Members, Science 191, 1131 (1976).
years. These peaks correspond to the 3. B. Dennison and V. N. Mansfield, Nature (Lon- 42. D. Ninkovich and N. J. Shackleton, Earth Plan-
don) 261, 32 (1976); F. Hoyle and R. A. Lyttle- et. Sci. Lett. 27, 20 (1975).
dominant periods of the earth's solar or- ton, Proc. Cambridge Philos. Soc. 35, 405 43. H. Thierstein, K. Geitzenauer, B. Molfino, N. J.
bit, and contain respectively about 10, (1939). Shackleton, in preparation; S. Gartner,
4. M. Milankovitch, K. Serb. Akad. Beogr. Spec. Palaeogeogr. Palaeoclimatol. Palaeoecol. 12,
25, and 50 percent of the climatic vari- Publ. 132 (1941) (translated by the Israel Pro-
gram for Scientific Translations, Jerusalem,
169 (1972).
44. J. D. Hays and N. J. Shackleton, in preparation.
ance. 1969). 45. T. Kellogg, in Climate of the Arctic, G. Weller
3) The 42,000-year climatic component 5. J. P. Kennett and R. C. Thunell, Science 187, and S. A. Bolliny, Eds. (Geophysical Institute,
497 (1975). Univ. of Alaska, Fairbanks, 1975), pp. 3-36;
has the same period as variations in the 6. G. Woilin, D. B. Ericson, W. B. F. Ryan, J. H. Geol. Soc. Am. Mem. 145 (1976), pp. 77-110.
obliquity of the earth's axis and retains a Foster, Earth Planet. Sci. Lett. 12, 175 (1971); 46. J. D. Hays and D. Ninkovich, Geol. Soc. Am.
G. Wollin, D. B. Ericson, W. B. F. Ryan, Na- Mem. 1Xe6 (1970), p. 263; J. D. Hays and N.
constant phase relationship with it. ture (London) 232, 549 (1971). Opdyke, Science 158, 1001 (1967).
4) The 23,000-year portion of the vari- 7. J. Weertman,J. Glaciol. 5, 145 (1964). 47. We have used only those ages for the stage 12-11
8. A.T. Wilson,Nature (London) 201,147(1964). boundary based on interpolations between the
ance displays the same periods (about 9. M. Ewing and W. L. Donn, in Polar Wandering core top and the Bnunhes-Matuyama magnetic
and Continental Drift (Society of Economic Pa- reversal. Although Hays et al. [J. D. Hays, T.
23,000 and 19,000 years) as the quasi-pe- leontologists and Mineralogists, Tulsa, Okla., Saito, N. D. Opdyke, L. H. Burckle, Geol. Soc.
riodic precession index. 1956), pp. 94-99. Am. Bull. 30, 1481 (1969)] estimated the time of
10. G. N. Plass, Tellus 8, 140 (1956). extinction of Stylatractus universus by this
5) The dominant, 100,000-year climat- 11. R. E. Newell, Quat. Res. (N.Y.) 4, 117 (1974). method and obtained an estimate of 341,000
10 DECEMBER 1976 1131
 years, Shackleton and Opdyke (31) argued con- BMDO2T program (58) modified by Y. Yera- domain, r data points are lost at the end of the
vincingly that this estimate was too young, be- caris to avoid recoloring the spectrum when the filtered signal. All of our filter calculations were
cause of sedimentation rate changes between prewhitening option is used. carried out with a BMDOIT program (58).
the lower and upper Brunhes magnetic series. 58. W. J. Dixon, Ed., BMD Biomedical Computer 72. E. Slutsky, Econometria 5, 105 (1937).
We have not used estimates based on Io/Th Programs (Los Angeles School of Medicine, 73. Cross-spectral techniques (55) were first applied
techniques because we believe the intrinsic in- Univ. of California, Los Angeles, 1965), p. 620. to deep-sea cores by Moore et al. [T. C. Moore,
accuracy of these estimates is large. 59. A. L. Berger [in Symposium on Long-Term Cli- Jr., N. Pisias, G. R. Heath, in The Fate of Fossil
48. C. Emiliani and N. J. Shackleton, Science 183, matic Fluctuations, Norwich, August, 1975 Fuel C02 (Plenum, New York, in press)].
511 (1974). (World Meteorological Organization, Geneva, 74. In arriving at the TUNE-UP chronology, we
49. E. Rona and C. Emiliani, ibid. 163, 66 (1969). 1975), pp. 65-72] concludes that variations in first expenmented with ages of the stage 12-11
50. H. H. Veeh and J. M. A. Chappell, ibid. 167, 862 eccentricity are positively correlated with plan- boundary ranging from 380,000 to 440,000 years
(1970). etary insolation receipt, and that, although ago, while retaining the stage 6-5 and younger
51. T.-L. Ku and W. S. Broecker, ibid. 131, 448 small, this effect might be climatically signifi- ELBOW control points (Table 2). The phase
(1966). cant. coherences between the orbital and filtered geo-
52. C. Sancetta, J. Imbrie, N. G. Kipp, Quat. Res. 60. Changes in n reflect the interaction of pre- logical signals reached a maximum in these ex-
(N.Y.) 3, 110 (1973). cession with the changing orientation of the periments when the 12-11 boundary was fixed at
53. L. A. Zadeh and C. A. Desoer, Linear Sys- orbital ellipse. 425,000 years. With that boundary determined,
tem Theory (McGraw-Hill, New York, 1963), 61. This interpretation of the precession curve, due we then experimented with different ages of the
p.628. to Broecker and van Donk (17), was confirmed 8-7 and 6-5 boundaries to arrive at the TUNE-
54. In this section of the article we have benefited by A. L. Berger (personal communication) as UP chronology. To evaluate the uncertainty of
from discussions with W. A. Wolovich. accurate to the first order of the earth's eccen- this chronology over the interval 50,000 to
55. G. M. Jenkins and D. G. Watts, Spectral Analy- tricity. 350,000 years ago (where we can check phase
sis and Its Applications (Holden-Day, San Fran- 62. Digital data were provided by A. D. Vernekar. coherence at both frequencies) we next con-
cisco, 1968). 63. We have estimated astronomical frequencies by ducted trials in which the three oldest TUNE-
56. R. B. Blackman and J. W. Tukey, The Measure- two methods: averaging the number of peaks in UP control points were all shifted by the same
ment of Power Spectrafrom the Point of View of the time domain, and spectral analysis (Tables 3 experimental interval, while holding constant
Communication Engineering (Dover, New and 4). The averaging gives satisfactory results the 9400-year level dated by carbon-14 and the
York, 1958), p. 190. for simple sinusoidal curves such as obliquity, zero-age core top. With a shift of either plus or
57. The procedures are as follows: (i) Selection of but cannot resolve the component frequencies in minus 22,000 years (that is, an interval equal to
an absolute chronology. For astronomical vari- more complex precession, insolation, or climate an average precession cycle and approximately
ables, we adopt the chronology of Vernekar curves. For uniformity, spectral estimates are equal to half an obliquity cycle), the phase coher-
(39). For geological variables, we use our chron- used exclusively in this article. ences are drastically disturbed in the younger
ological models discussed above. (ii) Calculation 64. Higher-resolution spectra of much longer eccen- part of the record. A final experimental series in
of a time series. A uniform sampling interval (t) tricity records actually have two discrete peaks which the three oldest TUNE-UP control points
approximately equal to the average interval be- (corresponding to penods of 96,000 and 125,000 were shifted together over a range of t 6000
tween sample points is chosen (in this study, years). years shows that the only chronologies worth
3000 years), and a time series of n values is 65. Discussions with J. W. Tukey were helpful in serious consideration (under the hypothesis of
calculated by linear interpolation. For the framing statistical hypotheses and in properly phase coherence at the frequencies of obliquity
PATCH time series, n ranges from 157 to 163 prewhitening the signals. and precession) show deviations from the
according to the chronology used. (iii) Detrend- 66. Questions of sampling variance apart, there are TUNE-UP values ranging between +5000 and
ing. Any long-term trend is removed by calculat- two ways in which our results might be artifacts -1000 years, with the optimum phase matches
ing residuals from linear regression. (iv) Pre- of procedure. First, the b and c peaks could be from 50,000 to 350,000 years ago being very
whitening. If the spectrum is dominated by low- aliases of cycles present in the cores but not close to that fixed by the TUNE-UP chronology
frequency components the signal may be pre- visible in our data because their frequencies are (Table 2).
whitened by means of a first-difference filter in higher than the Nyquist frequency. However, 75. G. Kukla, Geol. Foeren. Stockholm Foerh. 92,
which the intensity of the prewhitening is con- the mixing of deep-sea sediments by bottom- 148 (1970).
trolled by a constant, C. In our study, C = 0.998 living animals rules out the possibility that sig- 76. W. F. Ruddiman and A. McIntyre, Geol. Soc.
(56). (v) Lagging. The number of lags (m) to be nals with sufficient power at these frequencies Am. Mem. 145 (1976), pp. 111-146.
used in the next step is chosen according to the could remain in the record. Second, a significant 77. U.S. Committee for the Global Atmospheric
bandwidth (amount of detail) desired in the esti- part of the variance in the data for 8flhS and Research Program, National Research Council,
mated spectrum. (vi) Calculation of autocovari- particularly for Ts has the form of a DiMrac Understanding Climatic Change, A Program for
ance function. The covariance between the origi- comb-that is, has equally spaced spikes (Figs. Action (National Academy of Sciences, Wash-
nal time series and each of m lagged versions is 2 and 3). As the Fourier transform of a comb ington, D.C., 1975), figure A-14.
calculated. (vii) Smoothing. The autocovariance yields a comblike spectrum with frequencies 78. J. Imbrie, J. van Donk, N. G. Kipp, Quat. Res.
function, smoothed by a Hamming lag window, that are harmonics of the fundamental, it might (N. Y.) 3, 10 (1973).
yields a bandwidth of approximately 1.258/mAt be supposed that our b and c peaks (42,000, 79. N. Calder, Nature (London) 252, 216 (1974).
cycles per thousand years. To establish the gen- 24,000, and 19,500 years) could be statistically 80. J. Imbrie and J. Z. Imbrie, in preparation.
eral trend of the spectrum, we fix m as 8; the blurred versions of the third, fifth, and sixth 81. J. Hfilsemann, J. Sediment. Petrol. 36, 622
resulting spectrum has high accuracy but low harmonics of a fundamental 122,000-year cycle (1966).
resolution. For more detail, m is fixed at 40 to (- 41,000, - 24,000, and - 20,000 years). But 82. N. G. Pisias, Geol. Soc. Am. Mem. 145 (1976),
50, depending on n. (viii) Spectral estimation. several arguments can be advanced which elimi- pp. 375-392.
Variance contributions to each of m + 1 fre- nate this possibility: (i) No discrete spectral 83. Supported by NSF IDOE grants IDO 71-04204
quency bands are calculated by a Fourier trans- peaks occur at frequencies corresponding to the and OCE 75-19627 to Lamont-Doherty Geologi-
form of the smoothed autocovariance function. second and fourth harmonics (61,000 and cal Observatory of Columbia University; NSF
The resulting spectrum has a frequency range - 30,000 years) of the presumed fundamental. grant OCD 75-14934 to Brown University; and
from zero to the Nyquist frequency (fn = 1/ (ii) The dominant frequencies of peaks b and c NERC grant GR 3/1762 to Cambridge Universi-
2At), the highest frequency that can be deter- are not harmonics of the dominant cycles we ty, England. We are grateful to those who mate-
mined with a given At. Frequencies in cycles per estimate in a (94,000 years for Ts and 106,000 rially contributed to this article in the following
thousand years are plotted on a linear scale; a years for 8180). (iii) Although the percentage of ways: data generation, G. Irving, K. Jare (radio-
reciprocal scale indicates cycle length in thou- C. davisiana a peak actually has a dommant larian counts), M. A. Hall (operation of V.G.
sands of years. (ix) Scaling. The initial estimates period near 122,000 years, and its 42,000-year Micromass mass spectrometer), A. Vernekar
of an unprewhitened signal specify variance spectral peak is strongly developed, this time (orbital and insolation data in digital form), and
(" wer") as a function of frequency, that is, series does not approach a Dirac comb. A. L. Berger (digital data on seasonal in-
Pp. For prewhitened spectra, the correspond-
ing symbol is Pc(f). To simplify the use of con-
67. This method [R. T. Lacoss, Geophysica 36, 661
(1971)] was applied to our data by T. E. Land-
solations); data processing, Y. Yeracaris, J.
Morley, N. Kipp, D. Kirkpatrick, and J. Z.
fidence intervals, variance is expressed on a log ers, Lincoln Laboratory, Massachusetts Insti- Imbrie; production of manuscript, R. M. Cline,
scale. To facilitate the comparison of spectra tute of Technology. M. Perry, and R. Mellor. We are also grateful to
calculated from variables with different scales, 68. H. S. Tsien, Engineering Cybernetics (McGraw- all those who freely contributed their ideas and
values of P(f) are rescaled as A(f) so that the Hill, New York, 1952), p. 164. critically read the manuscript; in particular we
area under the curve equals unity. (x) Statistical 69. C. Emiliani, Ann. N.Y. Acad. Sci. 95, 521 wish to thank J. M. Mitchell, Jr., A. Cox, T.
evaluation. A one-sided confidence interval at- (1961). Hughes, J. W. Tukey, G. Kukla, A. Berger, and
tached to each estimate in the high-resolution 70. N. R. Goodman, J. Franklin Inst. 270, 437 S. Sachs. We are grateful to all our CLIMAP
spectrum is calculated from the chi-square distri- (1960). colleagues who encouraged and criticized this
bution, using 2n/m degrees of freedom, and ex- 71. The larger the value of r, the narrower the work during the past 2 years. This article is
pressed on a log scale (55, 56). Calculations in passband. Because the filter is an array of 2r + I Lamont-Doherty Geological Observatory Con-
steps iii through viii were carried out with a weights for a moving average applied in the time tribution 2434.

1 132 SCIENCE, VOL. 194