Pacemaking the Ice Ages by Frequency Modulation of Earth's Orbital Eccentricity

See allHide authors and affiliations

Science  23 Jul 1999:
Vol. 285, Issue 5427, pp. 564-568
DOI: 10.1126/science.285.5427.564


Evidence from power spectra of deep-sea oxygen isotope time series suggests that the climate system of Earth responds nonlinearly to astronomical forcing by frequency modulating eccentricity-related variations in insolation. With the help of a simple model, it is shown that frequency modulation of the approximate 100,000-year eccentricity cycles by the 413,000-year component accounts for the variable duration of the ice ages, the multiple-peak character of the time series spectra, and the notorious absence of significant spectral amplitude at the 413,000-year period. The observed spectra are consistent with the classic Milankovitch theories of insolation, so that climate forcing by 100,000-year variations in orbital inclination that cause periodic dust accretion appear unnecessary.

Understanding the climate of the recent geological past is of importance, because finding out how Earth's environment has transformed the heat input from the sun into climate variations can help predict future global climate change. Most of our knowledge about climatic variations in the Plio-Pleistocene (the last 5.2 million years) is extracted from time series of oxygen isotope ratios, a proxy for global ice volume generally known as the δ18O records (Fig. 1A). The records show that during the last million years, Earth has experienced at least 10 major glaciations, which according to the astronomical theory of the ice ages (1) are the consequence of secular variations in insolation caused by changes in Earth's orbital eccentricity, axial tilt, and longitude of perihelion (Earth's closest approach to the sun). The theory finds support in the fact that the spectra of the δ18O records contain some of the same frequencies as the astronomical variations (2–4), but a satisfactory explanation of how the changes in orbital eccentricity are transformed into the ∼100-ky (1 ky = 1000 years) quasi-periodic fluctuations in global ice volume indicated by the data has not yet been found (5).

Figure 1

(A) Oxygen isotope time series from ODP Site 806 shows the fluctuations in ice volume over the last 2.1 million years. (B) Time-frequency spectrogram showing the distribution of signal power in the isotope record. The horizontal band around 0.025 cycles/ky is the linear, single-frequency, time-invariant response of the climate system to the 41-ky variation in axial tilt, enhanced by orbital tuning. The much wider band of power centered at 0.01 cycles/ky is the response to variations in orbital eccentricity. In sharp contrast to the response to axial tilt, this is an age- and frequency-dependent response, likely to be nonlinear. The dashed sinusoid overlay indicates that the response appears to switch from one predominant low frequency to a predominant high frequency in cycles of ∼413 ky (the period of the sinusoid), strongly suggesting the presence of frequency modulation (FM). The spectrogram is computed with 350-ky-long windows for the Fast Fourier Transform and 10-ky lag between overlapping windows. (C) Around the 95-ky eccentricity component the power spectrum shows a characteristic multi-peaked power distribution that can be easily predicted assuming a 95-ky “carrier” frequency-modulated by a 413-ky modulating signal and a 826-ky subharmonic. The positions of the peaks are essentially insensitive to orbital tuning of the original time series.

For instance, a fundamental difficulty is to understand the notable absence in the δ18O data of a significant response to the 413-ky component of the orbital eccentricity, whose signal power is of the same order of magnitude as the 95-ky component. This has been called the “400-ky problem” by Imbrie and Imbrie (6).

One of the well-known features of the δ18O time series from deep-sea sediments is the slight but evident change in the duration of consecutive glacial periods (7), which oscillates from about 120 to 80 ky. For instance, the interval between the last two interglacials is 120 ky, while 400 thousand years ago (ka) the interglacial interval was nearly 80 ky, with three successive interglacials occurring in less than 200 ky. This is reflected clearly in the spectrogram of Fig. 1B (moving window Fourier transform) of Ocean Drilling Program (ODP) Site 806 (8), a representative paleoclimate record. The power of the eccentricity band (∼0.01 cycles/ky) in the spectrogram occurs along an approximately sinusoidal curve (as indicated by the dashed curve overlay), apparently because of periodic switching of main power from ∼120 ky to ∼80 ky every 400 ky or so. This fact suggests frequency modulation (FM) of the ∼100-ky eccentricity signals by the longer-period, 413-ky component. Such frequency modulation is entirely similar to electronic modulation of a high-frequency carrier by a low-frequency modulating signal, as used in FM radio and television broadcasting (9).

In addition to the characteristic sinusoidal shape of the spectrogram, frequency modulation of a single frequency carrier widens the spectral band by the addition of sidebands. Sidebands are spectral peaks distributed symmetrically on both sides of a carrier's peak at intervals equal to integer multiples of the modulating frequency (10). The appearance of orderly distributed sidebands in the power spectra of all the δ18O data analyzed here further supports the presence of frequency modulation. For instance, the power spectrum of Site 806 (Fig. 1C) shows evidence of frequency modulation as the observed 75- and 123-ky peaks coincide with the predicted sidebands of a 95-ky carrier frequency-modulated by a 413-ky signal, and the 85-, 107-, and 143-ky peaks coincide with the predicted sidebands of a 95-ky carrier frequency-modulated by the 826-ky subharmonic of the 413-ky signal (11).

Another important prediction of the FM hypothesis is that the power spectrum of an FM signal does not contain the modulating signal as a spectral peak, but as the separation interval between sidebands (9), a fact that is entirely consistent with the well-known absence of a distinctive 413-ky spectral peak in paleoclimate data.

It is important to point out that the spectra in Fig. 1are obtained from an orbitally tuned (12) time series. Orbital tuning is routinely applied to raw δ18O time series in order to establish a common time scale for all records regardless of location and sedimentation rate differences. However, the “clock” used to develop the chronology is the rhythm of any of the orbital periodicities, including eccentricity. Hence, it is possible that orbital tuning introduces bias or even spectral peaks not present in the data, and thus it is useful and prudent to analyze and compare both tuned and untuned records. In their original form the orbitally untuned records are of δ18O variation versus depth; the depth scale is converted to time using known published ages (determined by biostratigraphic or magnetic methods, or both) at one or two anchor points (13), and assuming a constant sedimentation rate for the rest of the sediment core. In spite of its obvious drawbacks, this assumption produces spectra similar enough to the tuned data to permit easy cross-identification of spectral peaks. As shown in Figs. 2 through 4, sidebands separated from the main spectral power as predicted by the FM hypothesis are a common feature in both tuned and untuned spectra from different sites, and as discussed below (Fig. 4), the effect of tuning appears to reduce, rather than to increase, the number of recognizable spectral peaks in the eccentricity band.

Since what follows hinges on the interpretation of power spectra, it is important to note that because of the strong age dependence of the eccentricity band (Fig. 1B), resolving individual spectral peaks depends critically on the window length chosen to calculate the spectrum. This becomes an important issue, because the apparent presence of a single power peak at 100 ky (Fig. 3) in several δ18O recordings has been interpreted as due to variations in Earth's orbital inclination, which force the planet to periodically cross a ring of extraterrestrial dust, thus causing the ice ages of the last million years (14). This conclusion is, however, based on power spectra computed from time windows spanning just the last 600 ky, an interval half the duration of typical eccentricity signals, which have significant power from 0 to about 1200 ka (Fig. 1B). As shown in Fig. 3, a 600-ky window limits the resolving power in the eccentricity band so much that important details are lost, and it is thus plausible that the presence of the 100-ky single spectral peak is due simply to merging of the 95-, 107-, and 125-ky eccentricity-related spectral peaks because the time window is too short, rather than from periodic changes in orbital inclination (15, 16).

If frequency modulation indeed occurs, it should be possible to replicate the observed variation in ice age duration seen in the δ18O records using an appropriate FM formula. A simple mathematical form of a frequency-modulated signal isEmbedded Image(1)where Ω = 2πf c, ω = 2πf m and f cand f m are the frequencies of the carrier and the modulating signals respectively, ɛ = Δf c/f m is the modulation index (17), and Δf c is the frequency deviation, a measure of the bandwidth produced by the modulation (9, 10). Since the predominant frequency shifts from ∼120 to ∼80 ky, which is also the observed shift in ice age duration, ɛ ≈ 0.7 (moderate modulation) for f m = 1/413 cycle/ky.

Trigonometric functions of the form in Eq. 1 can be represented by an infinite series of Bessel functions of increasing order, such asEmbedded Image Embedded Image Embedded Image(2)which clearly indicates that the Fourier spectrum of a FM signal such as given by Eq. 1 is composed of a central spectral peak at the carrier frequencyf c and additional spectral peaks (sidebands) at the frequencies f c ±nf m for n = 1, 2, 3 … , so that the modulating frequency f m = ω/2π appears in the power spectrum only as the interval between spectral peaks, but not as a spectral peak itself (9). A plausible model to simulate the δ18O records should thus contain the 125-, 100-, and 95-ky harmonics as carrier signals and the 413-ky and, as explained below, a 826-ky subharmonic as the modulating signals. The following model is possibly the simplestEmbedded Image Embedded Image Embedded Image Embedded Image(3)where t is the time in kiloyears. The constants a through c are adjustable parameters, with a typically three times greater than b or c, reflecting the relative strength of the eccentricty components. The term ɛ′ is the modulating index for the subharmonic, and since ɛ ∼ 1 then ɛ′ ∼ 2, sincef m is halved. In spite of the simplicity of Eq. 3, its spectrum (Fig. 4), reproduces every spectral peak in the power spectra of the tuned and untuned records, including the prominent combination tone at 107 ky (Figs. 2 and 4) which occurs only if the additional modulation by a 826-ky subharmonic is included in the model, as in Eq. 3. Frequency modulation as described by Eq. 3 is nonlinear in the sense that the effect produced by two modulating frequencies is not equal to the sum of the signals produced separately by each (18).

Figure 2

The multiple-peak character of the δ18O spectra is common to both the orbitally tuned and orbitally untuned records. Orbital tuning may, however, change the relative amplitudes of some peaks, as shown in (C). Modeling these differences provides further support for the FM hypothesis. The strong 41-ky peak in the orbitally tuned spectra is caused by tuning at that frequency. Data are from ODP Site 677 (A), Site 846 (B), and Site 849 (C). Bandwidths are 0.0011 ky–1 (A), 0.0014 ky–1 (B), and 0.00093 ky–1 (C). Age models are described in (13).

The 107-ky spectral peak is here interpreted as a combination tone (1/107 = 1/95 – 1/826) produced by nonlinear response of the climate system. Although the precise mechanism by which this may occur in the real climate system is not resolved here, an explanation is attempted below. On the other hand, Eq. 3 uses superposition of each modulated carrier to construct the total signal, and so in this sense the response is assumed linear.

Figure 4 shows that even when a spectrum is distorted by tuning, the sidebands can still be identified and the general structure of the spectrum understood. The distortion occurs because tuning in this case targeted the eccentricity frequencies (19), thus enhancing the amplitude of the 95-ky eccentricity carrier to the expense of its sidebands, as becomes apparent when comparing Fig. 4B with Fig. 4A. Accordingly, to compensate for tuning and simulate the untuned spectra required just a 20% increase in ɛ′ in Eq. 3.

It should be noted that the model spectra shown in Fig. 4 account not only for all the periods present, but also for those not present. As mentioned before, FM theory predicts that the modulating frequencies appear only as intervals between sidebands, and not as spectral peaks. Hence, if the role of the 413-ky eccentricity component is indeed to frequency-modulate the higher frequency components, we can expect to find signal power at 413 ky in all the δ18O records, though not in the form of a spectral peak, but rather as the interval between sidebands (clear examples are Figs. 1C, 3A, and 4). Thus, the apparent absence of a 413-ky spectral peak in Earth's climate response, also known as “the 400-ky problem” (6), may not be a problem after all, but the logical consequence of the fact that the 1/413 ky−1 signal acts as a frequency modulator.

Figure 3

The effect of time window length is crucial to the interpretation of power spectra of long-period δ18O data. Time windows here are (A) 0 to 2140 ky B.P., (B) 0 to 1070 ky B.P., and (C) 0 to 600 ky B.P. From (A) to (C), the spectra depict the strong loss in resolution when the δ18O time window is shorter than about 1000 ky. As shown, the worst spectral resolution occurs for a 600-ky window (C), which obscures the multiple-peak character of the spectrum. Both tuned and untuned records were studied, and identical loss of resolution and a similar single peak at 100 ky emerges when the short 600-ky-long window is used in both. This short window length was used, however, by Muller and MacDonald (14) to determine the power spectra of nine δ18O records, all of which look very much like (C) above (Fig. 3) (14), a feature which they interpreted as an indication of climate forcing by 100-ky oscillations of Earth's orbital plane. Bandwidths are 0.00093 ky–1(A), 0.0019 ky–1 (B) and 0.0033 ky–1 (C). The data for Site 849 is from (18).

Figure 4

(A) High resolution spectra (signal bandwidth 0.00093 cycles/ky) of the tuned Site 849 record (seeFig. 3) is shown compared to the Fourier transform of the model FM signal given by Eq. 3. As predicted by FM theory, the modulating power does not appear as a spectral peak, but as separation between each carrier and its sidebands. (B) High-resolution spectra (signal bandwidth 0.00093 cycles/ky) of the untuned Site 849 record. The comparison is made again with the Fourier transform of the model signal in Eq. 3. In order to optimize the fit, the modulating parameters ɛ and ɛ′ need to be increased by 20%. This suggests stronger modulation intensity, which causes the peak of the 95-ky carrier to be reduced in size relative to the sidebands, as predicted by FM theory. Thus, the FM mechanism appears robust for the untuned record. All the peaks predicted by the FM theory are present in both the tuned and the untuned spectra. A 95% confidence interval of 0.0005 cycles/ky was calculated from the mismatch between model peaks and the observed spectra. This is about 40% of the sideband separation of 1/826 ky−1.

The simulated FM signal ΦFM(t), shown in Fig. 5 compared to selected δ18O time series, successfully reproduces the subtle change in time interval between consecutive ice ages. For times before 1 Ma, the fit deteriorates, because the period of the ice ages becomes dominated by the 41-ky obliquity signal and the eccentricity carriers gradually disappear (Fig. 1, A and B).

Figure 5

Comparison of the FM signal constructed with Eq. 3and selected δ18O records. The synthetic FM signal closely reproduces the variable duration of the ice ages. The fit with the data deteriorates for times earlier than ∼900 ky B.P., because the 41-ky cycle becomes predominant. The best-fitting signal was obtained by trial-and-error varying the values of the adjustable parameters a, b, and c.

How does the climate system frequency-modulate the astronomical forcing? A natural place to search for an analogy is the mechanism of the resonant circuit of an electrical oscillator that generates FM signals as its capacitance or inductance are made to vary (9). The simplest analogue paleoclimatic mechanism is a conceptual representation of the climate system as a (potentially) resonant oscillator (20), where the linear extent of the ice sheet, L, satisfies an equation of the formEmbedded Image(4)where Ω = (1/C L C T)1/2is the natural angular frequency of the climate model,C T the ocean's thermal overturning time, andC L the ice cap's relaxation time. There are at least two ways to include external forcing in Eq. 4: with an additive term of the form F cos ω0 t, which introduces the possibility of resonance when ω0 → Ω, and/or with a time-varying parameter (parametric excitation), assuming that Ω changes with the dimension of the ice sheet (21) through CL, which is proportional to a critical dimension of the ice sheet (22).

If C L is of the formC L(1 + ΔC L/C L cos ωt) where ΔC LC L and ω is the modulating frequency, Eq. 4becomesEmbedded Image(5)which is the Mathieu equation (23), whose periodic solutions are frequency-modulated sinusoids (24) of the formEmbedded Image(6)where A and φ are integration constants. The solution to Eq. 6 is of the same form as each of the terms of Eq. 3used to simulate the δ18O records. Nonlinear forms of Eq. 5 are known to develop instabilities, bifurcations, and period-doubling cascades into chaos (25).


View Abstract

Stay Connected to Science

Navigate This Article