The wavelet time-frequency spectral analysis is applied to geological records that are proxies of paleoclimatic variations: δ18O in sedimentary cores, atmospheric CO2 concentration, and a loess magnetic susceptibility stratigraphy within the past million years. These spectra are compared with those for the astronomically predicted variations of the earth’s orbital eccentricity, obliquity, precession, and their resultant variations of the incoming insolation. The latter has been known to be unable to explain the characteristics of the observed 100-kyr paleoclimatic cycles. Based on similarities between the wavelet spectra of the orbital variations and paleoclimatic cycles, the authors introduce a signal–noise resonance theory to understand the dynamics of climate response to the orbital forcing. It is shown that the observed 100-kyr cycles are mainly caused by the period variation in the obliquity, which amplifies the small orbital forcing. But the observed flickers within these cycles are induced by the amplitude variation of obliquity and precession, which are two major components of the Milankovitch insolation deviations.
It is known that conventional astronomical theory for climate change cannot explain the magnitude of the 100- kyr paleoclimate cycles found in geological records. This quandary has called into question our fundamental understanding of how the climate responds to external or internal forcing (Broecker 1992; Ludwig et al. 1992;Winograd et al. 1992; Imbrie et al. 1993; Beauford 1994). Recently, new astronomical concepts (Liu 1992, 1995) and wavelet spectral methods (Chao and Naito 1995; Bolton et al. 1995; Lau and Weng 1995; Weng and Lau 1994) have been proposed for orbital and climate time series studies. In this paper, we compute the wavelet spectrum of the earth’s obliquity to explain the relationships between the observed paleoclimatic cycles and the earth’s orbital variations. Our purpose is to see how much of the observed climate variability can be accounted for by considering the temporal modulation in total insolation due to lengthening or shortening of the obliquity cyclic period.
2. Wavelet time-frequency spectral analysis
The wavelet transform of the time signal f(t) is defined as
where ψ(t) is the basic wavelet, a is a dilation/compression-scale factor that determines the characteristics frequency, and b represents the translation in time. Here we follow the scheme of Chao and Naito (1995) and select the Morlet wavelet (Morlet et al. 1982), which is a zero-mean, normalized, Gaussian-enveloped complex sinusoid:
The parameter ω, which determines the number of oscillations of the wavelet, is chosen to be ω = π(2/ln2)1/2 ≅ 5.336.
The wavelet spectral amplitude is customarily displayed in the time-frequency domain, that is, varying b and a. We shall examine only the real part of the amplitude, which, because of its symmetric nature, gives the amplitude undulation with the appropriate polarity and phase with respect to time. In contrast, the imaginary part, being antisymmetric, gives the amplitude undulation as well but imparts a 90° phase shift time. Alternatively, one can show the modulus of the wavelet transform; but the polarity and phase information becomes absent. Thus, amplitude peaks and troughs inhorizontal successions in high contrast shades (colors) indicate the presence of strong oscillations in the data, relative to the weaker, and hence less significant, “background” undulations of low contrast. The vertical location of the peak–trough succession gives the characteristic frequency (or period) of the oscillation in question; the corresponding horizontal location indicates the time when that oscillation occurs most strongly. In this manner the wavelet spectrum decomposes a time series into its time-frequency content.
Some limitations of the wavelet spectrum should be pointed out. Because of the temporal localization of the wavelet, the frequency resolution is limited to no better than about a quarter octave. Further, the lack of data outside our time span introduces the edge effect of the spectrum, which is more severe with longer periods.
3. Interpretation of paleoclimatic cycles
Advanced spectral tools such as Blackman–Tukey, maximum entropy, and the highly efficient Thomson technique have been applied to investigate the paleoclimatic variability using a set of four deep-sea core data (Berger et al. 1991). Dominant features are the presence of periods centered at 117.7, 43.6, and 24.9 or 19.3 kyr. One of the best paleoclimate records is the δ18O curve as shown in Fig. 1 (Shackleton et al. 1990). The geological timescale of the oxygen isotopic δ18O variations in ocean sediment cores is generated by astronomical tuning, which is more accurate than any other means of age–depth determination. The wavelet spectrum of δ18O (Bolton et al. 1995; Lau and Weng 1995) has shown three major cycles, their mean values centered at 100, 41, and 23 or 19 kyr. These periodicities correspond to the major periodicities of the earth’s orbital parameters: 1) orbital eccentricity (Fig. 2a), 2) obliquity (Fig. 2b), and 3) precession of the equinoxes (Fig. 2c). However, according to the conventional astronomical forcing for climate change (Milankovitch 1930; Vernekar 1972), the Milankovitch insolation deviations (Fig. 3) could explain only the 41 kyr and 23 or 19 kyr as an orbital response. Mathematically, Milankovitch insolation deviations consist of two components: 1) deviations due to the 41-kyr obliquity variation and 2) deviations due to the 23- or 19-kyr precession variation whose amplitude is modulated by the 100-kyr forcing of the orbital eccentricity. As shown in Fig. 3, the Milankovitch type of astronomical forcing provides an almost negligible contribution from eccentricity. For the past several decades, this has been a scientific puzzle. Geophysicists have speculated that the dominant 100-kyr cycle in Fig. 1 may be largely stochastic in nature because of internal forcing within the earth’s system (Imbrie et al. 1993).
4. Period variation of the obliquity
The problem concerns an intrinsic description of astronomical forcing for climate change in terms of periodvariation of the obliquity (Liu 1995). Figure 4 shows the wavelet spectrum of the period variations in the obliquity for the past million years, which is caused by solar and planetary gravitational torques on the equatorial bulge of the earth. We define this mechanism as astrodynamical forcing to distinguish from the Milankovitch astronomical forcing. It should be noted that the periods as illustrated in Fig. 4 (centered at 80, 120, and 160 kyr) are common multiples of the obliquity and the precession periods.
In the astronomical (Milankovitch) theory for climate change, the dynamical behavior of orbital oscillations has been limited to amplitude variations. According to phase theory in physics (Berry 1984; Berry and Robnik 1986; Monteoliva et al. 1994), however, there is another important factor that may be associated with the climate system. That is the phase function of orbital oscillations. The time derivative of this phase function is closely related to the observed climate change events (Liu 1992). Therefore, Milankovitch theory for climate change should be modified and extended to include the period variation effect of orbital forcing.
5. Climate response to orbital signal and noise
Based on Weertman’s (1976) ice sheet model, Liu (1995) has calculated ice sheet size changes in response to period variations of the obliquity, in addition to amplitude variations of the obliquity and precession. The half-width of the ice sheet size model of the Northern Hemisphere is defined by L(t). Temporal modulation of insolation due to lengthening or shortening of the obliquity cyclic period could cause the ice sheet size to fluctuate from glacial to interglacial conditions. Wavelet spectrum of the modeled ice sheet size in terms of half- width L(t) is shown in Fig. 5, which displays a similar characteristic time-frequency structure to that of Fig. 1, leading support to the validity of the astrogeodynamical forcing for the 100-kyr glacial cycles.
The similarity between the wavelet spectra of the observed and modeled time series suggests a possible way to understand the 100-kyr glacial cycles in a bistable climate system. It concerns the problem of detecting hidden signals in noisy environments. Climate response to orbital signals is intrinsically nonlinear in a noisy environment of orbital motion. The orbital amplitude signal alone, which causes insolation deviations, is unable to reproduce major spectral peaks in the observed climate records. The noise effect of the obliquity phase (or period), however, may produce a signal–noise resonance that amplifies the small orbital forcing. This new mechanism could enhance the probability of jumping between the two observable climate states.
6. Orbital signal–noise resonance
For the bistable climate system subject to the 100- kyr orbital eccentricity signal so weak as to be normallyundetectable, noise-induced cooperative phenomena can often be set up, leading to a stochastic resonance between the weak deterministic external signal and the stochastic internal noise (Benzi et al. 1981; Benzi et al. 1982; Sutera 1980, 1981; Matteucci 1989). Such a resonance can actually match the deterministic and stochastic timescales, making the signal apparent. This resonance theory of climate is based upon internal mechanisms within the bio–cryo–ocean–atmosphere–lithosphere system believed to have response time on the order of the observed mean periodicity (Hyde and Peltier 1985, 1987; Imbrie et al. 1993). However, internal mechanisms or feedbacks have difficulties in explaining the apparently episodic and irregular timing appearance of the 100-kyr cycle (Beaufort 1994).
Recently, a group of signal–noise resonance theories has appeared in the physics literature (Bulsara and Garamaitoni 1996; Wiesenfeld and Moss 1995). They yield the same resonance behavior of the output signal for subthreshold signal and noise inputs. It should be pointed out that the amplitude signal of the Milankovitch solar radiation has been embedded in the noise of its period variations (Liu 1992, 1995). A threshold system (Gammaitoni 1995) for the orbital signal and noise is shown in Fig. 6. The following statements can be made.
A bistable climate system, usually modulated by a weak periodic signal of the orbital eccentricity, is unable to scale the potential barrier in the absense of the obliquity period noise.
The input 100-kyr eccentricity signal is below threshold (dashed line), and its quantization would result in a constant output (red line) with a clear loss of signal detail.
Illustration of the signal–noise resonance consists of a threshold (shown by two straight dashed lines) and a subthreshold signal of the Milankovitch insolation at 65°N with added period noise of the obliquity. The addition of the period noise of appropriate intensity enables the signal to cross the thresholds in coincidence with its maximum amplitude. Each time the signal plus the noise increases across the threshold, a pulse can be written to the time series.
The quantized output then reproduces more closely the characteristics of the analog input signal, indicating the signal–noise resonance of orbital forcing as a paradigm for the mode (phase or period) locking (Ecke 1991; Hilborn 1994; Lau and Weng 1995) in driving the bistable climate system.
The presence of orbital noise alone is sufficient to induce irregular switching between the two wells shown in Fig. 6a. The dynamics of the bistable climate system can be modeled by the differential equation dx/dt = −d∪(x)/dx + S(t) + N(t) where ∪(x) is the barrier potential, S(t) is the amplitude variation signal of the Milankovitch insolation, and N(t) is the period variation noise of the obliquity. Since the signal–noise resonance is a nonlinear cooperative effect whereby the small signal entrainsthe noise-induced hopping (shown in Fig. 6c), we can ideally reconstruct the modulation function as a pulse train (Rowe 1965, chap. V). The orbital parameter of the pulse train consists of pulse amplitude, duration, and position. As expected, we have found that the noise power at the climate output is a smooth function of the noise power added to the input signal; in fact, for sufficiently large input noise, the output signal is dominated by the input noise. The wavelet spectrum of the pulse- modulated insolation time series exhibits characteristics similar to those shown in Figs. 1 and 5. In Figs. 1 and 5, each bistable element can be in either “red” or “purple” state and is subjected to a weak amplitude signal of insolation buried in the noise of period variation of the obliquity. Using the phenomenon of signal–noise resonance, the arrays are tuned to oscillate between states with 100-kyr cycles.
7. Climate variability induced by period and amplitude modulations of orbital parameters
Similarity between wavelet spectral characteristics in Figs. 1 and 5 provides a major opportunity to investigatehow much of the observed climate variability can be accounted for by considering the period and amplitude modulations of the orbital parameters. Model experiments indicate that period modulation of the obliquity (Liu 1992) associated with the Milankovitch solar radiation deviations (Verneker 1972) could have elicited sporadic rapid decay and prolonged regrowth in the history of the Northern Hemisphere ice sheet. These findings are best illustrated in Figs. 7a, b, which show clearly the similarity between the simulated and observed results. Figures 7a,b provide a reasonable physical framework for explaining the geological record (Shackleton et al. 1990) as the output signal of a bistable climate system responding to amplitude and period modulation of orbital forcing.
For the case of rapid melting, Fig. 7a shows that the 42 upward sawteeth in the half-width L(t) of ice sheet correspond to the 42 downward sawteeth in the δ18O record and that the minimum values of L(t) indicate rapid melting as seen in δ18O record. As for the case of prolonged glaciation, Fig. 7b shows that the 39 downward sawteeth in the half-width L(t) of ice sheet correspond to the 39 upward sawteeth in the δ18O record and that the times when the maximum values of L(t) proceed to occur correspond with the observed prolonged glaciation without a misfit event. Visual similarities in position and magnitude between the total 81 corresponding sawteeth are scientifically significant. The missing prediction of the initiation of large ice ages at −320 and −680 kyr in the previous work (Liu 1992) is thus corrected. Thus, model simulations have illustrated a kind of behavior that is consistent with the main features of the observed climate behavior, which contains flickers within cycles. The main cycles are mainly caused by the period variations of the obliquity. However, the flickers within cycles are induced by the amplitude variations of the obliquity and precession, which are two major components of the Milankovitch radiation deviations.
8. Comparison with other paleoclimate proxies
Geophysical records of atmospheric CO2 concentration (Shackleton and Pisias 1985) and loess magnetic susceptibility (Kukla et al. 1988) are proxies of paleoclimatic variations. Their wavelet spectra are shown in Figs. 8a,b. A good match of the spectra of these paleoclimate proxies (Figs. 8a,b) with the spectrum of the modeled ice sheet size (Fig. 5) can be found. Therefore, it seems conceivable that this new type of astrodynamical forcing, which includes both amplitude and period variations in the earth’s orbital parameters, may play an important role in explaining the 100-kyr cycle in the wavelet spectra of paleoclimate records.
In this paper, we have employed wavelet time-frequency spectrum analysis on two independent datasetsand compared them: the geologically observed δ18O variations and the astronomically derived orbital variations. Wavelet spectrum of the obliquity period in orbital variations represents a strong orbital noise that exhibits cyclic behavior centered at 80, 120, and 164 kyr. These periods of orbital noise are common multiples of the obliquity and the precession periods. Our analysis relies on linking this orbital noise to the orbital signal of the Milankovitch solar radiation and illustrates the significance of this noise, with which the substantive feature of the 100-kyr cycle in paleoclimate proxies can exist.
Similarity between the wavelet spectra of the observed and modeled time series of climate change reveals an orbital signal–noise resonance effect on climate response that could enhance the probability of jumping between the two observable climate states. Simulations have illustrated a kind of behavior that is consistent with the main features of the climate behavior, which contains flickers within cycles. The main cycles are mainly caused by the period variations of the obliquity, but the flickers within cycles are induced by the amplitude variations of the obliquity and precession, which are two important components of the Milankovitch insolation deviations.
We thank C. Wade Jr. at the Laboratory for Astronomy and Solar Physics for generating the astronomical time series. This study is supported by the NASA Geophysics Program.
Corresponding author address: Dr. Han-Shou Liu, Geodynamics Branch, Code 921, NASA/GSFC,Greenbelt, MD 20771.