Research Article

The 100,000-Year Ice-Age Cycle Identified and Found to Lag Temperature, Carbon Dioxide, and Orbital Eccentricity

See allHide authors and affiliations

Science  15 Sep 2000:
Vol. 289, Issue 5486, pp. 1897-1902
DOI: 10.1126/science.289.5486.1897


The deep-sea sediment oxygen isotopic composition (δ18O) record is dominated by a 100,000-year cyclicity that is universally interpreted as the main ice-age rhythm. Here, the ice volume component of this δ18O signal was extracted by using the record of δ18O in atmospheric oxygen trapped in Antarctic ice at Vostok, precisely orbitally tuned. The benthic marine δ18O record is heavily contaminated by the effect of deep-water temperature variability, but by using the Vostok record, the δ18O signals of ice volume, deep-water temperature, and additional processes affecting air δ18O (that is, a varying Dole effect) were separated. At the 100,000-year period, atmospheric carbon dioxide, Vostok air temperature, and deep-water temperature are in phase with orbital eccentricity, whereas ice volume lags these three variables. Hence, the 100,000-year cycle does not arise from ice sheet dynamics; instead, it is probably the response of the global carbon cycle that generates the eccentricity signal by causing changes in atmospheric carbon dioxide concentration.

It has long been apparent that there is an ∼100,000-year (100-ky) cyclicity in δ18O records from deep-sea cores. The first evidence for the duration of these cycles came from the use of radiometric dates from uplifted corals on the island of Barbados as a tool for dating the marine δ18O record (1). The second piece of evidence arose from anchoring the older part of the sequence of cycles to the radiometrically dated Brunhes-Matuyama magnetic polarity reversal (2). It is generally accepted that this 100-ky cycle represents a major component of the record of changes in total Northern Hemisphere ice volume (3). It is difficult to explain this predominant cycle in terms of orbital eccentricity because “the 100,000-year radiation cycle (arising from eccentricity variations) is much too small in amplitude and too late in phase to produce the corresponding climatic cycle by direct forcing” (4, p. 700). Most published explanations for the large amount of ice volume variability at this wavelength draw on the long response time that is associated with large ice sheets (5), especially when the response time of the underlying bedrock is taken into account. Imbrie and Imbrie (6) devised a simple mathematical model in which the slow and nonlinear response of a major ice sheet could generate an ice volume response with the 100-ky period of orbital eccentricity when the precession-dominated 65°N summer insolation record is used as a forcing. The nonlinearity of the ice sheet response (rapid deglaciation but slow ice growth) is an essential contribution to the relative success of this model.

It is difficult to test the proposed mechanisms for explaining the 100-ky cycle for several reasons. In particular, the strong 100-ky cycle has only dominated paleoclimate variability over about the past million years, too short a period to rigorously investigate with time-series analysis. Alternative explanations that have been proposed include forcing by orbital inclination (7) (which has a single 100-ky period, in contrast to the mixture of 95 and 125 ky that characterizes orbital eccentricity) and a wide range of models [reviewed in (4)] that seek to generate this cycle with or without the possibility of free oscillations. None of these investigations have questioned the primary association of the 100-ky δ18O cycle with global ice volume.

The time scale for the records that are used for the investigation of the 100-ky cycle is based on the existence in these δ18O records of variability that is confidently associated with variations in axial obliquity (41 ky) and with the precession of the equinoxes (23 and 19 ky) (8). Within the SPECMAP project (9,10), a time scale was developed for the marine δ18O record that is based on locking the phase of the observed cycles to an orbital template. The validity of this time scale is supported by the convincing coherence between marine δ18O (and other climate-related proxies) and the assumed obliquity and precession forcings. With these astronomically based age models, the 100-ky portion of the variability is significantly coherent with orbital eccentricity, but the coherence is not overwhelming. Given the difficulty in successfully explaining the 100-ky signal (4), doubts remain regarding its association with orbital eccentricity.

Atmospheric oxygen exchanges isotopically with the oxygen in ocean water chiefly as a consequence of global photosynthesis. The δ18O of atmospheric oxygen differs from that of the ocean by ∼23.5 per mil (‰); this difference is known as the Dole effect, and its origin is reviewed in (11). The record of δ18O of atmospheric oxygen trapped in bubbles in the polar ice sheets must reflect changes in marine δ18O with a small delay that is controlled by the balance between the rate of total global photosynthesis and the size of the atmospheric oxygen reservoir (12). Bender et al. (11) evaluated this balance and derived a turnover time of ∼1 ky. Here, the 100-ky component of the atmospheric δ18O signal is investigated. This proves to shed light on the 100-ky signal in marine δ18O.

A time scale for the atmospheric gas record from Vostok.Petit et al. (13) published 400-ky-long records of three signals embedded in the air bubbles of the Vostok ice cores. The concentrations of CO2 and CH4 are particularly interesting because these are both important greenhouse gases. The third record is of the δ18O of atmospheric oxygen, which varies in response to the δ18O of ocean water as noted above and to changes in the conditions of global terrestrial photosynthesis and other factors that change the magnitude of the Dole effect (11). Petit et al.(13) illustrated the obvious similarity between the δ18O of atmospheric oxygen on their time scale and 65°N summer insolation (Fig. 1A) but did not generate a new time scale based on this comparison.

Figure 1

(A) Vostok air δ18O record of (13), published time scale. (B) The “classic” Milankovitch forcing, June insolation at 65°N. (C) Benthic δ18O record of core V19-30 (18), published time scale.

Sowers et al. (12) generated a time scale for the upper part of the Vostok core by correlating their record of air δ18O (at that time extending to ∼130 ky) to the SPECMAP time scale (10). Instead of taking this approach, the similarity between Fig. 1A and 1B [noted in (13)] was used here to generate a more refined time scale that is independent of SPECMAP (10). The advantage of this approach is that the strong amplitude modulation of the precession forcing provides a very powerful criterion for judging the success of the time scale where there is a precession-related signal (14). None of the records used in constructing the SPECMAP time scale (10) contain such a “clean” precession signal; therefore, it appears likely that this procedure will lead to a more accurate time scale than could be achieved by correlating to SPECMAP. In order to determine the phase of the record relative to precession, an independent age is required for at least one event. The upper part of the Greenland Ice Sheet Project 2 ice core has an age scale that is based on annual layer counting and is accurate to better than ±2% (15); the record of δ18O of atmospheric oxygen in this core (16) shows that the age of the midpoint of the most recent glacial-to-interglacial transition is 12,000 years ago (12 ka). The phase lags in the tuning target, 5 ky for precession and 8 ky for obliquity, are determined so as to be consistent with this age (once the precession phase is set, the phase relative to obliquity is determined by the record itself). In the final tuning target, the amplitudes of the precession and obliquity components are matched to their coherent amplitudes in the record. Steps in the processing of the time series are available as supplemental Web material (17). To minimize “spikes” in the differences between pairs of records, tuning was carried out by assigning ages to the midpoints of transitions; Web table 1 (17) gives the published and tuned transition ages. The data on the new time scale, together with the tuning target, are shown in Fig. 2.

Figure 2

(A) Vostok air δ18O record (13), on the tuned time scale of this paper. (B) Tuning target for the Vostok air δ18O record. (C) Benthic δ18O record of cores V19-30 and V19-28 on the tuned time scale of this paper. (D) Tuning target for the V19-30/V19-28 benthic δ18O record. (E) Linear variance amplitude spectrum for the air δ18O record, obtained by cross-spectral analysis versus its tuning target; the shaded area illustrates the amplitude of coherent variance in the obliquity and precession bands (the tuning targets do not include eccentricity). (F) Linear variance amplitude spectrum for the benthic δ18O record, obtained by cross-spectral analysis versus its tuning target; the shaded area illustrates the amplitude of coherent variance in the obliquity and precession bands (the tuning targets do not include eccentricity). Arrows labeled “e,” “t,” and “p” identify the frequencies associated with eccentricity, obliquity (tilt), and precession.

A new time scale for the benthic marine δ18O record. In order to construct a strictly comparable deep-sea chronology, the well-known benthic foraminiferal δ18O record of core V19-30 (18) was retuned by using a procedure analogous to that described above. The V19-30 record was slightly extended by using the data from nearby core V19-28 (19) to make a high-quality record of the same length as the Vostok record. In designing the final tuning target for V19-30, an account was taken of the fact that the target (and hence the tuned data) should lead rather than lag the Vostok target, because the Vostok air δ18O record must incorporate the history of the changing marine δ18O record with an ∼1-ky lag (11). In order to satisfy this requirement, the modeled midpoint of the last deglaciation must be 13 ka, yielding lags of 3 ky with respect to precession and 7 ky with respect to obliquity. The V19-30 (extended) record on this time scale is shown in Fig. 2C, together with its tuning target. This time scale is very similar (20) to that already published but differs in detail because the record was tuned directly to a specially designed orbital target instead of being tuned with SPECMAP age controls (10) and also because only transition midpoints were used for age controls, rather than a mixture of controls. The amplitudes and coherent amplitudes and the phase lags in the obliquity and precession bands of the two records after this tuning are given in Table 1.

Table 1

Spectral analyses of ETP (25) versus tuned air δ18O and tuned foraminiferal δ18O, carried out to evaluate the linearly forced components of the two records (41). Confidence limits for phase lags are at the 95% level.

View this table:

In the time scales developed under the inspiration of Imbrie (8,9), the assumption was made that the phase lags with respect to obliquity and precession ought to reflect the lag of a single-exponent system with a reasonable time constant for the response of a continental ice sheet. In the present work, this requirement was relaxed because it is now known that there are a number of different components contributing to the δ18O records, each with a different time constant. The approach whereby the tuning target is optimized to the data enables the difference (in thousands of years) between the lags with respect to the obliquity and precession components of midsummer insolation to be estimated; geological ages for the last deglaciation then give the absolute lags.

Comparing the Vostok air δ18O record and the benthic marine δ18O record. The linear variance spectra of the air and benthic δ18O records obtained by cross-spectral analysis versus their respective orbital tuning targets are also shown in Fig. 2 (21). In view of the relatively short (∼1 ky) turnover time of atmospheric oxygen with respect to marine photosynthesis (11), the amplitude of the marine δ18O signal must be fully captured by the air δ18O record, yet at the obliquity frequency, the foraminiferal signal is clearly larger. This demonstrates that a substantial part of the foraminiferal obliquity signal (∼0.12‰ coherent amplitude) must derive from changing deep-ocean temperature, which affects the isotopic fractionation between water and carbonate, and not from ice volume. On the other hand, in the precession band it is evident that almost half of the signal recorded in air δ18O must arise from changes in the Dole effect, because at this frequency the air δ18O signal is twice the amplitude of the foraminiferal δ18O signal (22). It is clear from the elegant and thorough review by Bender et al. (11) that there are many potential routes through which this could arise. However, the extremely high coherence with precession (see Table 1) suggests that there must be a relatively simple dominating mechanism that yields an almost linear response. Although part of the phase lag undoubtedly reflects the 1-ky response of the atmosphere with respect to photosynthesis, the remainder implies either a mechanism for changing the Dole effect that has a long time constant (which is unlikely) or control centered on a time of year that is later than midsummer. A relation with varying monsoon precipitation was suggested in (11), and both modeling and observation suggest that the monsoon precipitation lags the orbital forcing by a few thousand years (23), equivalent to an autumn response.

The precession component in air δ18O must be composed of an ice volume signal and a Dole effect signal. Because these two must have approximately the same absolute phase, the amplitudes may be subtracted without complication, yielding a maximum amplitude of 0.22‰ for ice volume and a minimum of 0.21‰ relating to the Dole effect. These values are adopted in the following discussion in the knowledge that if a higher value is assigned, the implication would be a reduced precession component allocated to ice volume variability, in which case a part of the precession-related variability in the V19-30 record would be reallocated to changing deep-water temperature. The best way to refine this partitioning will be to accurately determine the sea levels at ∼90 and 110 ka in comparison with levels at ∼80 and 100 ka, because these differences reflect the precession contribution to sea level variability.

It is impossible at present to analytically determine whether any of the 41-ky variance in the air δ18O signal is attributable to the Dole effect or whether it derives from the ocean water. However, if it is assumed that the precession variance in the Dole effect is related to monsoon precipitation, then there ought to be an obliquity component that is commensurate with the contribution of obliquity to the low-latitude insolation spectrum. In the region and season of the Indian monsoon, the proportion of variance arising from obliquity is ∼10% of the total, equivalent to an amplitude of ∼0.02‰, which is negligible in relation to the overall 0.27‰ coherent amplitude associated with obliquity.

The residual signal and eccentricity. It has been shown that almost all of the 41- and 21-ky variance in the air δ18O signal can be explained as a linear response to orbital forcing. The result of subtracting this linearly forced component (24) from the observed record is shown in Fig. 3A. The equivalent for the V19-30 foraminiferal δ18O record, obtained by subtracting the linearly forced component of the 41- and 21-ky variance from the full record is shown in Fig. 3B. Three features stand out. First, both residuals display the 100-ky cycle; second, the amplitude of this signal is far smaller in the Vostok atmospheric δ18O residual than in the V19-30 foraminiferal δ18O residual; and third, the residuals retain “jumps” in value at deglaciations, such as at the base of the Holocene (12 ka) and at the base of marine isotope stage (MIS) 5 (128 ka).

Figure 3

Components of the δ18O records not explained by linear forcing by insolation. (A) The residual obtained by subtracting the Vostok air δ18O tuning target from the measured record, attributed to ice volume. (B) The residual obtained by subtracting the V19-30 benthic δ18O tuning target from the measured record, attributed to both ice volume and deep-ocean temperature. foram, foraminiferal. (C) The difference between these two residuals, attributed to deep-ocean temperature.

With regard to the first observation, cross-spectral analysis versus eccentricity, tilt, and precession (ETP) parameters (25) shows that, in both cases, the residual is coherent with eccentricity with small, but different, phase lags. However, the records are too short to attach much importance to this observation on its own.

It is difficult to conceive of a mechanism that would prevent the atmospheric δ18O recording the full amplitude of a 100-ky cycle in ocean water δ18O; thus, the implication of the second observation must be that a substantial portion of the marine 100-ky cycle that has been the object of so much attention over the past quarter of a century is, in reality, a deep-water temperature signal and not an ice volume signal.

The third observation confirms the supposition that the rapid deglacial “terminations” cannot be explained as part of the linear response of the climate system (4). This observation has the important implication that the midpoint of the transition in a linear target cannot be expected to predict the exact age of the midpoint of one of these deglaciations (26). On the other hand, the rapid nonlinear part of the transition must have been close to this midpoint, or the residual would show an anomalous overshoot. Experiments were conducted to determine the oldest allowable base for the last interglacial that is consistent with the assumptions involved in the tuning (17). When Martinson et al. (10) explored the uncertainties in their age model, they were working with geological data sets containing a much larger and irreducible proportion of variance that could not be explained by linear forcing; therefore, they were obliged to attach a rather large uncertainty (an average of 5 ky) to their age estimates. The age model for Vostok air δ18O could not be varied by more than ∼2 ky without substantially reducing the amount of variance representing a linear response.

Exploring the 100-ky signal of the Vostok atmospheric CO2 records. A long record of atmospheric CO2 has been obtained from the Vostok core, and it has been noted that its variance is dominated by eccentricity (13). Because the new time scale refers to the air bubbles, it can be directly applied to the CO2 record (Fig. 4F). With this time scale, the earlier observation is confirmed: The record has a coherency of 0.92 with the 100-ky component of eccentricity and is essentially in phase, and the coherent amplitude of the 100-ky CO2 signal is equivalent to 30 parts per million by volume (ppmv). At obliquity, the coherence is 0.94, and the coherent amplitude of the signal is equivalent to 16 ppmv, with about the same phase as the δ18O signals.

Figure 4

(A) The reconstructed ocean water δ18O record attributable to ice volume and sea level changes. (B) The reconstructed ocean water δ18O record (150 ky), scaled as sea level and compared with sea level observations (32). (C) The reconstructed record of changes in the Dole effect. (D) The reconstructed component of foraminiferal δ18O attributable to changes in deep-Pacific water temperature. (E) The record of atmospheric CO2(13) on the time scale of this paper. (F) Vostok D/H (13) on the time scale of this paper.

The residual from the Vostok air δ18O signal, after subtracting the linearly forced part of the record, is interpreted in terms of change in the δ18O of the ocean. To a first approximation, it should represent the total of that part of the ocean δ18O (“ice volume”) variability that is not linearly forced by obliquity or precession (Fig. 3A). At a period of 100 ky, this residual signal is highly coherent (0.99) with the Vostok CO2 signal but lags CO2.

The residual in the V19-30 foraminiferal δ18O signal, after subtracting the linearly forced part (Fig. 3B), must contain contributions from global ice volume and from changing deep-ocean temperature. Thus, to obtain the component of the ocean temperature record that is not linearly forced, I subtracted the Vostok residual (ice volume) from the V19-30 residual (ice volume and temperature), taking account of the 1-ky lag of the air δ18O record (Fig. 3C). The coherent amplitude of this record at the 100-ky period is 0.26‰, equivalent to ∼0.7°C in temperature. Deep-water temperature is also highly coherent (0.97) with the Vostok CO2 signal in the 100-ky band. The coherences and phases in the 100-ky band strongly suggest that atmospheric CO2 has a direct and immediate control on deep-water temperature (presumably with high-latitude air temperature as an intermediary). In addition, temperature has a direct control on ice volume, mediated by the long response times of ice sheets.

It could be argued that the data equally support a model whereby deep temperature has a direct and immediate control over the concentration of CO2 in the atmosphere. However, it is hard to envisage a mechanism for generating a 100-ky deep-water temperature cyclicity directly, whereas the global carbon system that controls atmospheric CO2 includes components with long time constants, a necessary requirement for generating a 100-ky signal from the orbital forcing. In either case, the important conclusion that may be drawn from the examination of the 100-ky part of the variance in the several signals is that the prevalent 100-ky signal does not arise from the long time constant of ice volume response. Instead, ice volume responds, with an appropriate lag that is consistent with the value derived by Imbrie and Imbrie (5), to forcing at 100 ky by temperature (it can be assumed that the temperature of the deep ocean also reflects in some way the air temperature at high latitudes).

Important support for this interpretation comes from the Vostok ice D/H record, which is thought to reflect high-latitude temperature in the air over Vostok. Assigning a time scale to the ice itself poses additional problems that are not addressed here, because the ice-air age difference varies from ∼2 ky at times of higher temperature and accumulation rate to 6 ky or more at cold times when snow accumulation is low. The usual procedure is to derive a time scale for the ice and then compute a time-dependent air-ice age difference that takes account of the temperature and accumulation rates as a function of time. Here, the primary time scale has been developed for the trapped air, so a complex iteration would be required to formally derive an ice time scale. Instead, the ice time scale (13) is simply offset by the same amounts as the gas time scale was changed (Fig. 4E). Although this will slightly affect the high-frequency portion of the variance spectrum, it will not substantially affect the 100-ky signal. This signal has a phase similar to that of the deep temperature record and a coherent D/H amplitude at the 100-ky period of 15‰, equivalent to an air temperature of ∼3°C.

Testing the partitioning of the δ18O signal.On the basis of the argument above, the δ18O record of global ice volume may be reconstructed from three components: the precession component, from the coherent precession component of the V19-30 foraminiferal δ18O tuning target; the obliquity component, from the coherent obliquity component of the Vostok air δ18O tuning target; and the nonlinear residual that is obtained by subtracting the Vostok air δ18O tuning target from the measured record (27). This reconstruction is illustrated in Fig. 4A. Two aspects of this reconstruction support its being a better reconstruction of ice volume and sea level from δ18O than has been previously achieved. First, the ice volume contribution to the difference between the Last Glacial Maximum and the Holocene is ∼1.0‰, close to the reconstruction of 1.0 ± 0.1‰ by Schrag et al. (28) based on direct measurement of the δ18O of relict glacial pore water [also compare estimates in (29) rather than the 1.2‰ (2), 1.26‰ (30), or 1.4‰ (31) estimates based on earlier interpretations of marine isotope records]. Second, the record is consistent with recent sea level estimates (32), as shown in Fig. 4B.

With the derived history of the ice volume component of the foraminiferal δ18O record, the total temperature component may be obtained by subtraction as mentioned above. The resulting record is illustrated in Fig. 4D. Figure 5 shows the total amplitude and coherent amplitude spectra of ice volume (Fig. 4A), Pacific temperature (Fig. 4D), Antarctic air temperature (Fig. 4E), and atmospheric CO2 (Fig. 4F), together with their estimated phases relative to ETP (25). Table 2summarizes the statistical description of the 100-ky signal in the four records. The CO2 and deep-ocean temperature records are dominated by variance with a 100-ky period and vary in phase with eccentricity. Contrary to previous interpretations, global ice volume does lag changes in orbital eccentricity, but eccentricity is not the dominant peak in the ice volume spectrum (33). The effect of orbital eccentricity probably enters the paleoclimatic record through an influence on the concentration of atmospheric CO2.

Figure 5

Linear variance spectra obtained by cross-spectral analysis versus ETP (25) of (A) deep-Pacific temperature, (B) ocean δ18O (ice volume, sea level), (C) Vostok D/H (Antarctic air temperature), and (D) Vostok atmospheric CO2. In each panel, the top section shows the amplitude and (shaded area) coherent amplitude spectrum, and the bottom shows phases (with 95% confidence limits) for each orbital band (using the convention where a positive phase angle is a measure of the lag in degrees with reference to ETP, that is, with reference to insolation in the Northern Hemisphere midsummer). Arrows labeled “e,” “t,” and “p” identify the frequencies associated with eccentricity, obliquity (tilt), and precession. The bar labeled “BW” indicates the bandwidth.

Table 2

Spectral analyses (41) to investigate the 100-ky signal in geological records (Fig. 5). Confidence limits for phase lags are at the 95% level. For an explanation of the ETP used, see (25).

View this table:


View Abstract

Navigate This Article