Glacial cycles and orbital inclination

Richard A. Muller

January 1994

Lawrence Berkeley Laboratory Report LBL-35665

   Milankovitch1 attributed cycles in the earth's ice coverage to perturbations in the motion of the earth and the resulting changes of insolation (solar heating) in the northern hemisphere, which has most of the land mass. The strongest effects were expected to come from changes in the earth's obliquity (tilt of the earth's spin axis with respect to its orbit) and from a precessional term that accounts for the delay between summer solstice (when the pole faces the sun) and perihelion (when the earth is closest to the sun). Perturbations of the earth's motion come from gravitational effects of the planets and the moon, and can be calculated with precision back at least 10 million years.2,3

   The most compelling evidence for the theory is the presence of earth orbital frequencies in proxy climate variables such as CO2 and d18O.4,5 A careful spectral analysis by MacDonald6 shows that eccentricity, obliquity and precession frequencies match 8 of 9 statistically-significant d18O frequencies in an eastern Pacific core; the "only anomalous line" in the data with no theoretical counterpart is a d18O peak at 0.014 ± 0.001/kyr (period 71 ± 5 kyr).

   However the insolation mechanism of Milankovitch has difficulty accounting for the relative magnitudes of the cycles. The spectrum for the last 500 kyr is dominated by the 100 kyr period, usually attributed to oscillations of eccentricity e. But yearly average insolation is proportional to (1/2)e2, and the calculated 4% variations in e yield insolation changes of only 10-3. A linear calculation6 shows the solar insolation at a typical latitude (60o N) at the primary frequencies for (precession, obliquity, eccentricity) to be in the ratio (1, 0.2, 0.02), in sharp disagreement with the d18O data which show the climate variations in the ratios (1, 2.5, 11). The eccentricity term, rather than being 50 times weaker than precession, is 11 times larger, in disagreement by a factor of 550. Attempts to account for this discrepancy are summarized by Imbrie et al.5, who propose to solve the paradoxes with a "conceptual model" that depends on nonlinear behavior of ice sheets. They cite several additional features that a successful model should have, including the suppression of the 400 kyr cycle (strong in eccentricity calculations but absent in the climate data), and the variability of the amplitude of the 100 kyr cycle, which although dominant recently was weak in the period 1-3 Ma.5,7 . Broecker 8,9 and Winograd10 argue that new d18O measurements show warming associated with several glacial terminations begins before the start of above-average insolation, thus disproving Milankovitch's mechanism; however this conclusion is disputed by Imbrie et al.11. We do not try to resolve this controversy here, but instead offer a conceptual model that differs in a qualitative way from that of Milankovitch.

Orbital Inclination

   The inclination i of the earth's orbital plane has been assumed to have little effect on climate because of its negligible effect on insolation. (i should not be confused with the inclination of the earth's spin axis.) The variation of i ("nutation") as calculated by Berger2 is shown in figure 1; it is quasi-periodic with period 70 kyr (frequency 0.0144 /kyr). Also shown is Omega, the projection of the inclination vector on the reference plane, taken to be that of the earth's orbit in 1850. The increase of Omega is called the "advance of the line of nodes" in astronomical jargon. The parameters i and Omega are plotted in polar coordinates at the bottom.

   
Figure 1. Inclination and Omega (precession of inclination vector) for past 10 6 yr, based on Berger.2 Figure 2. Inclination and Omega using the invariable plane as a coordinate reference. The black circles show the inclination at the times of the abrupt glacial terminations; the uncertainties in the times are systematic and equal to or less than the diameters of the circles.

   Berger's choice of a reference coordinate based on the Earth's 1850 orbit is an arbitrary choice that doesn't affect insolation calculations. However the behavior of the orbit is clearer in coordinates based on the invariable plane of the solar system, perpendicular to the total angular momentum vector of the perturbing planets. We performed this transformation, and the results are shown in figure 2. The nutation is nearly sinusoidal with period 100 kyr rather than 70 kyr, the amplitude is smaller by two, and Omega now shows a virtually steady advance of one cycle per 70 kyr. The simplification that occurs can be understood by using polar coordinates, as shown in the bottoms of the figures. The 70 kyr precession of the vector (in the invariable plane), when measured in a coordinate system in which it is not centered (the frame of 1850), becomes mixed into the inclination to yield the apparent 70 kyr nutation.

Inclination rather than eccentricity

   We hypothesize that the 100 kyr variation in the orbital inclination is the dominant one seen in the climate data, rather than the 100 kyr period of eccentricity. This explains in a natural way the absence of a 400 kyr eccentricity cycle. The 70 kyr orbital precession (the advance of the line of nodes) matches the "anomalous" line of MacDonald.6

The only plausible mechanism we have found that could relate glaciation to inclination is the accretion of interplanetary material. As the orbit of the earth changes, it passes through different parts of the sun's zodiacal ring, and encounters different regions of meteoroids and dust. This variability can, in turn, result in changes in the density of dust and aerosols in the upper atmosphere, as shown in Figure 3 below.

 
 Figure 3. The orbit of the hypothesized cloud of meteoroids and dust. When the inclination of the earth's orbit is low, it spends the entire year in this cloud. When it is high (as shown) it passes through only twice per year.

   Interplanetary dust accreting on the sun has previously been proposed as a driver of the ice ages.12,13 Clube14 discussed the possibility of accretion from a single large (and unknown) meteoroid stream affecting the Earth's climate, but he did not draw any connection to the periodicity of glacial cycles.

   The effects of high altitude dust and aerosols are known primarily from volcanic eruptions. Global cooling of 0.5-1 C was estimated from the eruption of Krakatoa, and measurable climate changes have been attributed to El Chichon and other recent eruptions that injected several megatons of material into the stratosphere.15,16 Large, explosive volcanic events occur typically once per century, so the average injection of volcanic material ¸ 100 kton/yr. Accretion could cause cooling (as volcanic eruptions indicate) or warming (if cometary particles inject water). Small particles (0.5 xm-diam.) take several years to reach the ground17; larger particles (10xm) take a few months; gases can reside for much longer.

   Measurements by Kyte and Wasson18 of iridium in oceanic sediments show that long-term global average flux from extraterrestrial materials for the period 35-70 Ma is 60-120 kiloton/yr, about the same as that from present-day volcanic eruptions. The density of meteoroids and dust the zodiacal ring could depend on the number of recent collisions among asteroids and comets. It would be interesting to see if the 100 kyr period could be resolved in the iridium data; however the shortest variation detectable by Kyte and Wasson (the average Nyquist period)¸ 600 kyr. We performed a spectral analysis of 299 iridium data measurements in the late Cretaceous samples of Alvarez et al.19, but found no peaks above expected noise fluctuations; however only 25% of the 8.8 ppt iridium in the sample is estimated to be cosmic in origin, and this severely limited our sensitivity. A definitive test requires sediment with low iridium background, low bioturbation, and frequent sampling. The best candidate for this may be in the Greenland ice cores.

Structure in the Zodiacal Ring

   For our model we must assume the existence of a solar ring with sufficient structure that accretion changes as the orbit of the Earth varies. A possible configuration is shown in figure 3. When i is low, the earth remains in the ring; when i is high, it passes through only twice per year. The Infrared Astronomical Satellite (IRAS) observed a band of emission near the ecliptic¸ 2 degrees wide, approximately the right thickness for our model, although IRAS is not sensitive to larger 100 xm to cm meteoroid-sized particles that could dominate. Jackson and Zook 23 showed that the earth and other planets act as "shepherds," stabilizing asteroidal dust into toroidal orbits; they suggested that these orbits could account for structure seen by IRAS. Dermott et al.24 verified that the IRAS data matches the calculated orbits. Dermott25 argued that narrow bands could extend into the region of the earth.

   The interplanetary dust has been thought to have a lifetime (limited by Poynting-Robertson drag) to be significantly less than the 100 kyr glacial period; see, for example, Grün.26 However Nishiizumi et al.27 measured the cosmic ray exposure age of dust particles to be between several times 105 years to several times 106 years. These long lifetimes may be due to planetary shepherding or to an unknown mechanism; however they are sufficiently long to allow dust to play a role in the observed 100 kyr cycles. Of course even if the dust had a short lifetime, the spatial structure of the orbits could endure much longer by continual replenishment from the asteroid belt.

Accretion measurements

   Since we are currently in a warm period then our model suggests that the present extraterrestrial flux should be small compared to the long term average. In a recent review, Rocchia et al.28 found 6 of 7 published determinations of recent influx measured¸ 10 kiloton/yr; the long-term flux was 6 to 12 times larger. They speculate that the recent low flux is a fluctuation, due to the absence of recent large impacts. Instead we suggest that the discrepancy is due to the fact that interglacials such as the present are periods of low extraterrestrial flux. Iridium in ice during the period 19 - 14 ka has been studied by LaViolette 29 found that total cosmic influx during that period were "one to two orders of magnitude higher" than present levels. However his absolute levels of present iridium levels disagree with those made by others [see Rocchia et al. 28] so we cannot consider this confirmation of our theory.

The 41 kyr obliquity period

   The second strongest line seen in the last million years is the 41 kyr period attributed by Milankovitch to changes in insolation associated with obliquity variations. We note that the accretion in the Northern Hemisphere depends on the tilt of the earth's axis with respect to its orbit, i.e. the obliquity, and so it is also possible that this period is driven at least partly by accretion rather than insolation. Measurements of iridium densities in the Greenland ice cores could also test these hypotheses.

Summary and Conclusions

   The 100 kyr period could come from oscillations in the orbital inclination, eccentricity, or both; however the lack of the 400 kyr period suggests that inclination dominates. Extraterrestrial accretion rather than solar insolation may be the forcing mechanism. The observed 70 kyr period may be due to the precession of the earth's orbital plane. The 41 kyr period could be due to insolation, as Milankovitch proposed, or from accretion in the Northern Hemisphere. Our model suggests that the recent measurements of low extraterrestrial fluxes may not be anomalous, but due to the earth being in a relatively debris-free region.

   We now summarize the strengths and weaknesses of the model:

Strengths:

Weaknesses:

   We note that weakness 2 is shared by most other models. The required dust and meteoroid streams, although not understood theoretically, are not ruled out experimentally. The greatest strength of the model is that it can be tested in the near future by looking at iridium levels in Greenland and Antarctic ice.