Bispectral Fingerprint Identifies 100 kyr Glacial Cycle:

Orbital Inclination

Gordon J. MacDonald, University of California, San Diego

Richard A. Muller, University of California, Berkeley

Abstract

Introduction

    Much progress has been made in the understanding of past climates by using linear spectral analysis, often Blackman-Tukey, to identify periods of interest in proxy records of climate (Hays et al., 1976). Spectral analysis establishes that the dominant period in many proxy climate records is 100,000 years (100 kyr), but the origin of this peak remains enigmatic (Imbrie et al., 1993). Linear spectral analysis techniques are obviously of limited value when various spectral components interact with one another due to some nonlinear or parametric process. With respect to climate, the orbital equations of motion for the Earth are highly nonlinear. An analysis technique that describes nonlinearities in climate records may allow identification of which orbital parameters control climate. In addition, a variety of proposals involving nonlinear behavior of the Earth have been advanced to explain paleoclimatic variations. An early nonlinear model was described by Imbrie and Imbrie (1981), involving an on-off response of ice sheets. Since then, dozens of models have been proposed to alter the observed astronomical forcing function so as to come into concordance with observations. These are reviewed in Imbrie et al. (1993). Curiously, nonlinear statistical techniques have not been applied to climate records.

    The technique we use in this paper is bispectral analysis. The first empirical bispectrum was calculated by Hasselman et al. (1963) using ocean height data. They observed a strong interaction of a dominant spectral component with itself (harmonic generation) and appreciable interaction with other spectral components. These interactions resulted in the peaking of the crests of the shallow water wave as the ocean waves moved onshore. Since this example was published, hundreds of applications have been described in the literature in such varied fields as turbulent fluid flow, electrical brain activity, plasma interactions, dynamical astronomy, and astrophysics (Tyron, 1981; Nikias and Petropulu, 1993).

    We begin by defining the bispectrum and the steps involved in calculating the bispectrum from a digital time series. We then present a brief tutorial example based on the suggestion of Imbrie and Imbrie (1980). We next discuss an analysis of the bispectrum of the theoretical variation of the inclination of the Earth's orbital plane (Muller, 1994) and compare it with the bispectrum of a proxy climate record, the revised SPECMAP (Muller and MacDonald, 1994). The bispectrum of SPECMAP contains the readily identifiable fingerprint of orbital inclination as well as a number of other interesting features.

Bispectrum

    We consider a zero mean time series x(t) that may be represented by its Fourier components Xk

x(t) = S[Xk exp(-2 ¼ i fk t)]

where S stands for the "sum" over k from ­(infinity) to +(infinity). The power spectrum P(k) and the bispectrum B(k,l) are then defined by

P(k) = <Xk Xk*>

B(k,l) = <Xk Xl Xk+l*>

The operator < > is the expectation value or averaging operator. Xk* is the complex conjugate of Xk, Xk* = X-k.

    The power spectrum describes the contribution to the variance of x(t) (<x2(t)>) from the spectral component at fk of width Df (=1/T) where T is the length of the record. The bispectrum represents the contribution to the skewness <x3(t)>) of those spectral components at fk and fl of area (Df)2. With this interpretation we have

<x2(t)> = SkPk

<x3(t)> = Sk,lB(k,l)

where, once again, the "S" stands for "Sum" (usually Greek"Sigma").

    The definition of the bispectrum clearly shows that the bispectrum B(k,l) measures the statistical dependence between three frequency bands centered at fk, fl , and f-k-l. If there is energy at these bands, as there would be in any noisy record, then each band is characterized by an amplitude and a phase. If the phases are statistically independent, the sum or difference phases will be randomly distributed over (-p,p). When a statistical averaging, denoted by the operator < > is carried out, the bispectrum will vanish due to random phase mixing. Alternatively, if the three frequency bands are coupled to each other, the total phase will not be random at all, although the phase of each band may be randomly changing. This situation is one of the phase coherence among the three bands. In this case, the statistical averaging will not lead to a vanishing bispectrum. Phase coherence among various bands in the spectrum of a time series is evidence that the series was produced by a nonlinear process (MacDonald, 1989).

    A simple mechanical example illustrates the concept of the bispectrum. Consider a machine in which three wheels are interconnected through new teeth gears that initially operate without chatter. Suppose that their angular frequencies are w1, w2, and w1 + w2. The machine will radiate acoustic energy at the corresponding frequencies f1, f2, and f1+f2. A power spectrum of the acoustic field radiated by the machine would show peaks at f1, f2, and f1 +f2 and possibly at their harmonics, 2f1, 2f2, etc., with a continuum of background noise. The bispectrum of the acoustic signal would have a peak at (f1,f2) and, if the harmonics are generated, at (f1,f1), (f2,f2), etc. As the gears wear, perfect phase lock is lost, and there is a certain chatter as the gears slip randomly. Frequency coupling is maintained, with the power spectrum showing peaks at f1, f2, and f1 +f2, but phase coupling is lost. As phase coupling is lost, the bispectral peak diminishes in amplitude. Indeed, gear wear in machinery is operationally monitored by continuously calculating the instantaneous bispectrum (Sato et al., 1977) and noting the amplitude of the bispectrum.

    We note from its definition that the bispectrum is highly symmetric,

B(k,l) = B(l,k) = B*(-k,-l) = B(-k-l,l) = B(k,-k-l).

Further, the bispectrum detects only the quadratic (symmetric) part of a nonlinearity. If a linear signal x(t) passes through a cubic device, y(t) = x3(t), then the bispectrum of y(t) vanishes.

Computation of the Bispectrum

    Numerous schemes have been advanced over the years to estimate the complex bispectrum of a time series. Hasselman et al. (1963), in the pre-FFT days, used narrowband pass filters to estimate Xk, Xl, and Xk+l. Alternatively, in analogy with Blackman-Tukey, the triple correlation function r3(t1,t2) is calculated

r3(t1,t2) = (1/N) Sk=11+N-1[x(k) x(k + t1) x(k + t2)]

and the Fourier transform is then taken with respect to t1 and t2 (Lohmann and Wirnitzer, 1984).

    Our procedure follows current practice (Nikias and Petropulu, 1993). We assume that the sampling interval Dt is sufficiently small that the Nyquist frequency fN = 1/(2Dt) is larger than any spectral component of interest present in x(t). We further assume that the record length, T = NDt, is large enough to assure adequate frequency resolution. If frequency bands fk and fl contain significant power, then we desire that the elementary bandwidth Df = 1/T be much smaller than |fk - fl|. This consideration determines how we carry out the phase averaging. One method would break the record up into M separate records and then average the computed bispectrum over the M records. This procedure reduces the effective resolution by a factor of M. Instead, from the full record we compute a moving average of ten nearest neighbor estimates of the bispectra in the (f1,f2) plane, making use of the full symmetry of the bispectrum. This procedure smudges the resulting peaks in the bispectrum but retains the full resolution for estimates of individual points of the bispectrum.

In summary, the computational procedure involves:

1. Subtracting the mean value for the record;

2. Windowing with a Hanning window to reduce leakage (Harris, 1978);

3. Computing the Fourier amplitudes and phases using the FFT techniques (Press et al., 1993); and

4. Estimating the bispectrum by averaging in the (f1,f2) plane.

A Nonlinear Model of Climate

    Imbrie and Imbrie (1980) proposed a simple nonlinear model of ice volume behavior. The model has recently been used by Shackleton et al. (1990) to devise a time scale for glacial epochs that is widely accepted (Imbrie et al., 1993; Muller and MacDonald, 1994). The Imbrie and Imbrie model relates the ice volume, y(t), to an approximation of the northern high latitude insolation x(t) where

x(t) = e + a e sin[w -f),

where e is the obliquity, e the eccentricity, e sin w the precessional parameter, and a and f are adjustable constants. The model equations relating ice volume y(t) to insolation x(t) are:

dy/dt = (1+b)(x-y)/Tm if x > y

dy/dt = (1-b)(x-y)/Tm if x < y

The two time scales in the problem-rate of ice accumulation and rate of ice melting-are incorporated in b and Tm is a mean time constant for the problem.

    The nonlinearity is introduced in the model by the different functional form of y(t) depending on the relative scaled values of ice volume and insolation. The nonlinearity has a quadratic component, and sum and difference frequencies should be generated by the model. The three leading terms in the precessional parameter e sin w have frequencies:

f1 = 0.042166 cpky

f2 = 0.044587 cpky

f3 = 0.052698 cpky

(Berger, 1978). The sum frequencies give lines in the spectra at periods corresponding to 11.5, 10.5, and 10.3 kyr. These peaks are not observed in proxy climate records (MacDonald, 1990). The difference frequencies and corresponding periods are: f2 -f1 = 0.00242 413 kyr

f3 -f1 = 0.01053    95 kyr

f3 -f2 = 0.00811    123 kyr.

The strongest peak would be at 413 kyr, but as is widely recognized, proxy climate records do not show this difference frequency (Imbrie et al., 1993; MacDonald and Muller, 1994). The split eccentricity peaks at 123 and 95 kyr are also not observed in high-resolution spectra of isotopic records (MacDonald and Muller, 1994). Further bispectral analysis fails to detect any interaction involving f1, f2, and f3. On the basis of these observations, it is unlikely that the Imbrie and Imbrie model is applicable, and the analysis raises further questions about the validity of the Shackleton et al. (1990) time scale (Muller and MacDonald, 1994).

Bispectrum of SPECMAP

    Muller (1994) has suggested that the dominant peak in proxy climate records at about 100 kyr is due to the variation in inclination of the Earth's orbit relative to the invariable plane (the plane perpendicular to the angular momentum vector of the solar system). Quinn et al. (1990) have carried out a detailed numerical integration of the dynamics of the solar system for 3.05 million years into the past. They find significant differences with Berger's (1978) secular solutions, particularly for times earlier than one million years ago (Ma). The variation of inclination as obtained by Quinn et al. (1990) and transformed with respect to the invariable plane is shown in figure 1. The strong 100 kyr oscillation is clearly evident, as is a longer-term modulation.

 

    Figure 2a gives the power spectrum of the calculated inclination.

 

The principal peaks are at

f1T = 0.0102 cpky

f2T = 0.0086 cpky

but both of these peaks are split with

f3T = 0.091 cpky

f4T = 0.006 cpky

The record length is insufficient to determine the peak frequencies of the split peaks accurately. We note, however, the possibility of frequency coupling, since

f2T + f3T = 0.0996 Å f1T

The power spectrum of the revised SPECMAP is shown for comparison in figure 2b.

    A contour plot of the bispectrum of the theoretical inclination is shown in figure . A statistically significant peak shows at (0.009, 0.001), indicating that the frequency bands near 0.001, 0.01, and 0.009 cpky are frequency and phase locked. The bispectrum detects the split peak near 0.01 cpky and relates it to the peak near 0.001. The calculation has two unsatisfactory features. The record length and averaging broaden the bispectral peaks. The bispectrum of the full 3 million year record shows a much sharper peak. Since the record is only 900 kyr long, the peak at 0.001 cpky has less than a complete cycle, but still maintains a phase lock with peaks at 0.009 and 0.01, which go through nine cycles.

 

Figure 3. Bispectra of eccentricity, of the revised Specmap oxygen isotope data, and of orbital inclination. The bispectra all cover the period 0-894 ka. The axes are frequency, in units of cycles per kiloyear. Note the similarity between Specmap and inclination bispectra, and the disagreement between Specmap and eccentricity bispectra.

    Figure 3 also shows the bispectrum of revised SPECMAP. The dominant peak is at (0.009, 0.001) cpkyr. The peak is significantly broader than the theoretical bispectrum peak, as might be expected for the noisier observational record. We suggest that the close coincidence of the dominant peak in the bispectrum of SPECMAP with that of inclination argues strongly that the 100 kyr cycle observed in climate records is associated with variations in the inclination of the Earth's orbital plane. Also note the disagreement between the bispectrum of Specmap and that of eccentricity.

Discussion

    In a series of previous papers (Muller, 1994; MacDonald and Muller, 1994; Muller and MacDonald, 1994), we have argued on a number of grounds that variations in the Earth's orbital inclination have played a key role in determining past climate, at least over the past 900 kyr. We believe that the present paper adds weight to the argument, even while recognizing the deficiencies in our treatment. The agreement between the theoretical and observational bispectra are remarkably close, but are not perfect. Limitations on length of record proscribe seeing the split peaks and their interaction in detail. While we believe that revised SPECMAP represents an improved time scale over the past 900 kyr, there still may remain significant inaccuracies in dating.

    Extension of the time scale beyond 900 ka is essential if we are to make further progress in understanding past global climate change. Almost certainly this will require acquisition of reliable isotopic age determinations. We are concerned that relying on cycle counting may not be sufficient. If orbital inclination is important and dust accretion is a mechanism for cooling, there may be periods, for astronomical reasons, when the solar system is more or less dusty. In these circumstances, the relative magnitudes of various peaks will vary with possible confusion as to the identification of the appropriate dating cycle.

References

Berger, A. (1978), Long-term variation of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362-2367.

Harris, F. (1978), On the use of windows for harmonic analysis with discrete Fourier transform, Proceedings IEEE, 66, 51-56.

Hasselman, K., W. Munk, and G. MacDonald (1963), Bispectra of ocean waves, in Time Series Analysis, M. Rosenblatt, Ed., New York: Wiley, 125-139.

Hays, J., J. Imbrie, and N. Shackleton (1976), Variations in the earth's orbit: Pacemaker of the ice ages, Science, 194, 1121-1132.

Imbrie, J., and J. Z. Imbrie (1980), Modeling the climate response to orbital variations, Science, 207, 943-952.

Imbrie, J., A. Berger, E. Bogle, S. Clemens, A. Duffy, W. Howard, G. Kukla, J. Kutzbach, D. Martinson, A. McIntyre, A. Mir, B. Molfino, J. Morley, L. Peterson, N. Pisias, W. Prell, M. Raymo, N. Shackleton, and J. Toggweiler (1993), On the structure and origin of major glaciation cycles. 2. The 100,000 year cycle, Paleoceanography, 8, 699-735.

Lohmann, A., and B. Wirnitzer (1984), Triple correlations, Proceedings IEEE, 72, 889-901.

MacDonald, G. (1989), Spectral analysis of time series generated by nonlinear processes, Rev. Geophys., 27, 449-469.

MacDonald, G. (1990), Global climate change, in Global Climate and Ecosystem Change, G. MacDonald and L. Sertorio, Eds., New York: Plenum, 1-95.

MacDonald, G., and R. Muller (1994), High-resolution spectral analysis of glacial cycles: eccentricity is ruled out, Lawrence Berkeley Laboratory Report LBL-35660, May 1994 (submitted for publication).

Muller, R. (1994), Glacial cycles and orbital inclination, Lawrence Berkeley Laboratory Report LBL-35665, January 1994.

Muller, R., and G. MacDonald (1994), Glacial cycles: Orbital inclination dominated for the last 900,000 years, Lawrence Berkeley Laboratory Report LBL-35667, May 1994 .

Nikias, C., and A. Petropulu (1993), Higher-Order Spectra Analysis, Englewood Cliffs, NJ: Prentice-Hall.

Press, W., B. Flannery, S. Teukolsky, and W. Vetterling (1993), Numerical Recipes, Cambridge University Press.

Quinn, T., S. Tremaine, and M. Duncan (1991), A three million year integration of the earth's orbit, Astron. J., 101, 2287-2305.

Sato, T., K. Sasaki, and Y. Nakamura (1976), Real-time bispectral analysis of gear noise and its application to contactless diagnosis, J. Acoust. Soc. Am., 62, 382-387.

Shackleton, N., A. Berger, and W. Peltier (1990), An alternative astronomical calibration of lower Pleistocene time scale based on ODP Site 677, Trans. Roy. Soc. Edinburgh, 81, 251-261.

Tyron, P. (1981), The bispectrum and higher order spectra: A bibliography, National Bureau of Standards, Technical Note 1036, Washington, D.C.