The Azimuth Project
Milankovitch cycle



The Milankovitch cycles are periodic or quasiperiodic changes in the parameters of the Earth’s orbit and tilt, which in turn affect the climate. The three major types of Milankovitch cycle are:

(changes greatly exaggerated)

These changes do not effect the overall annual amount of solar radiation hitting the Earth, but they affect the strength of the seasons in different ways at different latitudes. It is widely believed that they are partially responsible for the glacial cycles. However, the details of how this occurs are complex and poorly understood, as the following graph makes clear. Some form of amplification, e.g. stochastic resonance, may be involved.


A brief history of Milankovitch cycles and their effect on glaciation is given in:

Let us quote a bit:

Joseph Adhémar (1842) formulated the conjecture that the precession of Earth’s orbit caused ice ages, because of its effects on the seasonal distribution of incoming solar radiation (insolation). He realised that one season had to be critical given that the changes in winter and summer insolation caused by precession cancel each other (Herschell, 1830). James Croll (1875) further developed the theory and appreciated the importance of the three Earth orbital elements determining insolation: eccentricity, obliquity (angle between the equator and the ecliptic) and the longitude of the perihelion.

Milutin Milankovitch (1920, 1941) is quoted as the first one to have established a self-contained mathematical theory linking orbital elements to insolation, and insolation to climate changes. On a suggestion of V. Köppen, he calculated the secular evolution of the summer radiation of the northern latitude (published in Köppen and Wegener, 1925), following the hypothesis that glacial cycles are driven by the effects of changes in summer insolation on snow ablation rate during this season (the conjecture had in fact been earlier expressed by Murphy, 1876, and it contradicted Croll’s views). Summer insolation is largest when obliquity is large and/or when summer in the northern hemisphere corresponds to the time of perihelion passing. Milankovitch later substantiated the hypothesis by explicitly calculating the effects of radiation changes on the position of the snow line.

By the mid-seventies the Earth’s orbital parameters were known over the last million years to a good accuracy thanks to the work of several generations of astronomers who based themselves on foundations lain by Laplace and Le Verrier. A decisive step was made by Berger (1978), who expressed in an analytical form the Fourier decomposition of the Earth’s orbital parameters relevant for the astronomical theory of palaeoclimates. This work constitutes the first demonstration that the spectrum of climatic precession is dominated by periods of 19, 22 and 24 kyr, that of obliquity, by a period 41 kyr, and eccentricity has periods of 400 kyr, 125 kyr and 96 kyr. The most accurate solution to date is the La04 (Laskar et al., 2004).


There are three major forms of Milankovitch cycle:

Eccentricity: The Earth’s orbit is an ellipse, and the eccentricity of this ellipse says how far it is from being circular. But the Earth’s orbit slowly changes shape: it varies from being nearly circular (eccentricity of 0.005) to being more strongly elliptical (eccentricity of 0.058), with a mean eccentricity of 0.028. There are several periodic components to these variations. The strongest occurs with a period of 413,000 years, and changes the eccentricity by ±0.012. Two other components have periods of 95,000 and 125,000 years.

The eccentricity affects the percentage difference in incoming solar radiation between the perihelion, when the Earth is closest to the Sun, and aphelion, when it is farthest from the Sun. This works as follows. The percentage difference between the Earth’s distance from the Sun at perihelion and aphelion is twice the eccentricity, and the percentage change in incoming solar radiation is about twice that. The first fact follows from the definition of eccentricity, while the second follows from differentiating the inverse-square relationship between brightness and distance.

Right now the eccentricity is 0.0167, or 1.67%. Thus, the distance from the Earth to Sun varies 3.34% over the course of a year. This in turn gives an annual variation in incoming solar radiation of about 6.68%. Note that this is not the cause of the seasons, which arise due to the Earth’s tilt, and occur at different times in the northern and southern hemispheres.

Obliquity: The angle of the Earth’s axial tilt with respect to the plane of its orbit, called the obliquity, varies between 22.1° and 24.5° in a roughly periodic way, with a period of 41,000 years. When the obliquity is high, the strength of seasonal variations is stronger.

Right now the obliquity is 23.44°, roughly halfway between its extreme values. It is decreasing, and will reach its minimum value around the year 10,000 CE.

Precession: The slow turning in the direction of the Earth’s axis of rotation relative to the fixed stars, called precession, has a period of roughly 26,000 years. As precession occurs, the seasons drift in and out of phase with the perihelion and aphelion of the Earth’s orbit.

Right now the perihelion occurs during the southern hemisphere’s summer, while the aphelion is reached during the southern winter. This tends to make the southern hemisphere seasons more extreme than the northern hemisphere seasons.

The gradual precession of the Earth is not due to the same physical mechanism as the wobbling of the top. That wobbling does occur, but it has a period of only 427 days. The 26,000-year precession is due to tidal interactions between the Earth, Sun and Moon. For details, see:

Glacial cycles

The relation between the Milankovitch cycles and glacial cycles is complex, fascinating, but also poorly understood and controversial. One can see why here:

The above graph was taken from Wikimedia Commons. The Milankovitch cycles were computed using Fortran programs available here:

At this website you can create your own table of the eccentricity, obliquity and longitude of perhelion for selected years by entering the minimum year, maximum year, and yearly increment.

Stochastic resonance

There have been interesting attempts to argue that stochastic resonance is important in amplifying the small effects of Milankovitch cycles to create glacial cycles. Work along these lines began with Benzi, Patera, Sutera, and Vulpiani:

For a review of this line of work see:

Also see:

For a fascinating but highly mathematical treatment using energy balance models see:

For more, see Stochastic resonance.

The effect of changes in eccentricity

Changes in obliquity clearly have very little effect on the global, annual average insolation, since the Earth is nearly spherical. It is more tricky to determine how changes in the eccentricity of the Earth’s orbit affect this quantity. Let us work this out in detail, following a calculation presented by Greg Egan on the Azimuth blog. While the result is surely not new, Egan’s approach makes nice use of the fact that both gravity and solar radiation obey an inverse-square law!

Here is his calculation:

The angular velocity dθdt=Jmr 2\frac{d \theta}{d t} = \frac{J}{m r^2}, where JJ is the constant orbital angular momentum of the planet and mm is its mass, so if the radiant energy delivered per unit time to the planet is dUdt=Cr 2\frac{d U}{d t} = \frac{C}{r^2} for some constant CC, the energy delivered per unit of angular progress around the orbit is

dUdθ=Cr 2dtdθ=CmJ.\frac{d U}{d \theta} = \frac{C}{r^2} \frac{d t}{d \theta} = \frac{C m}{J}.

So, the total energy delivered in one period will be U=2πCmJU=\frac{2\pi C m}{J}.

How can we relate the orbital angular momentum JJ to the shape of the orbit? If you equate the total energy of the planet, kinetic 12mv 2\frac{1}{2}m v^2 plus potential GMmr -\frac{G M m}{r}, at its aphelion r 1r_1 and perihelion r 2r_2, and use JJ to get the velocity in the kinetic energy term from its distance, v=Jmrv=\frac{J}{m r}, when we solve for JJ we get:

J=m2GMr 1r 2r 1+r 2=mbGMaJ = m \sqrt{\frac{2 G M r_1 r_2}{r_1+r_2}} = m b \sqrt{\frac{G M}{a}}

where a=12(r 1+r 2)a=\frac{1}{2} (r_1+r_2) is the semi-major axis of the orbit and b=r 1r 2b=\sqrt{r_1 r_2} is the semi-minor axis. But we can also relate JJ to the period of the orbit, TT, by integrating the rate at which orbital area is swept out by the planet, 12r 2dθdt=J2m\frac{1}{2} r^2 \frac{d \theta}{d t} = \frac{J}{2 m} over one orbit. Since the area of an ellipse is πab\pi a b, this gives us:

J=2πabmT.J = \frac{2 \pi a b m}{T}.

Equating these two expressions for JJ shows that the period is:

T=2πa 3GMT = 2 \pi \sqrt{\frac{a^3}{G M}}

So the period depends only on the semi-major axis; for a fixed value of aa, it’s independent of the eccentricity.

If we agree to hold the orbital period TT, and hence the semi-major axis aa, constant and only vary the eccentricity of the orbit, we have:

U=2πCmJ=2πCbaGMU=\frac{2\pi C m}{J} = \frac{2\pi C}{b} \sqrt{\frac{a}{G M}}

Expressing the semi-minor axis in terms of the semi-major axis and the eccentricity, b 2=a 2(1e 2)b^2 = a^2 (1-e^2), we get:

U=2πCGMa(1e 2)U=\frac{2\pi C}{\sqrt{G M a (1-e^2)}}

So to second order in ee, we have:

U=πCGMa(2+e 2)U = \frac{\pi C}{\sqrt{G M a}} (2+e^2)

The expressions simplify if we consider average rate of energy delivery over an orbit, which makes all the grungy constants related to gravitational dynamics go away:

UT=Ca 21e 2\frac{U}{T} = \frac{C}{a^2 \sqrt{1-e^2}}

or to second order in ee:

UT=Ca 2(1+12e 2)\frac{U}{T} = \frac{C}{a^2} (1+\frac{1}{2} e^2)

We can now work out how much the actual changes in the Earth’s orbit affect the amount of solar radiation it gets! In what follows we assume that the semi-major axis is indeed held constant while eccentricity changes; from the work of Laskar (see below) this is approximately correct. According to Wikipedia:

The shape of the Earth’s orbit varies in time between nearly circular (low eccentricity of 0.005) and mildly elliptical (high eccentricity of 0.058) with the mean eccentricity of 0.028.

The total energy the Earth gets each year from solar radiation is proportional to

11e 2\frac{1}{\sqrt{1-e^2}}

where ee is the eccentricity. When the eccentricity is at its lowest value, e=0.005e = 0.005, we get

11e 2=1.0000125\frac{1}{\sqrt{1-e^2}} = 1.0000125

When the eccentricity is at its highest value, e=0.058e = 0.058, we get

11e 2=1.00168626\frac{1}{\sqrt{1-e^2}} = 1.00168626

So, the change is about

1.00168626/1.0000125=1.001673731.00168626/1.0000125 = 1.00167373

In other words, a change of merely 0.167%.

That’s very small And the effect on the Earth’s temperature would naively be even less!

Naively, we can treat the Earth as a greybody. Since the temperature of a greybody is proportional to the fourth root of the power it receives, a 0.167% change in solar energy received per year corresponds to a percentage change in temperature roughly one fourth as big. That’s a 0.042% change in temperature. If we imagine starting with an Earth like ours, with an average temperature of roughly 290 kelvin, that’s a change of just 0.12 kelvin!

The upshot seems to be this: in a naive model without any amplifying effects, changes in the eccentricity of the Earth’s orbit would cause temperature changes of just 0.12 °C!

This is much less than the roughly 5 °C change we see between glacial and interglacial periods. So, if changes in eccentricity are important in glacial cycles, we have some explaining to do. Possible explanations include season-dependent phenomena and climate feedback effects.


News articles and overviews

This news item surveys 17 papers in an issue of Paleoceanography entitled ‘Milankovitch climate cycles through geologic time’:

This news item provides a quick summary of John Imbrie’s work on the SPECMAP project, which uses spectral analysis to look for effects of the Milankovitch cycles on different aspects of our climate:

Abstract: Milankovic climate oscillations help define climate sensitivity and assess potential human-made climate effects. We conclude that Earth in the warmest interglacial periods was less than 1°C warmer than in the Holocene and that goals of limiting human-made warming to 2°C and CO2 to 450 ppm are prescriptions for disaster. Polar warmth in prior interglacials and the Pliocene does not imply that a significant cushion remains between today’s climate and dangerous warming, rather that Earth today is poised to experience strong amplifying polar feedbacks in response to moderate additional warming.

Deglaciation, disintegration of ice sheets, is nonlinear, spurred by amplifying feedbacks. If warming reaches a level that forces deglaciation, the rate of sea level rise will depend on the doubling time for ice sheet mass loss. Gravity satellite data, although too brief to be conclusive, are consistent with a doubling time of 10 years or less, implying the possibility of multi-meter sea level rise this century. The emerging shift to accelerating ice sheet mass loss supports our conclusion that Earth’s temperature has returned to at least the Holocene maximum. Rapid reduction of fossil fuel emissions is required for humanity to succeed in preserving a planet resembling the one on which civilization developed.

Technical papers

Hays, Imbrie and Shackleton 1976

This was a major analysis that applied Fourier analysis to climate proxies. From several hundred sediment cores studied stratigraphically by the CLIMAP project, they selected two (RC11-120 and E49-18) from the Indian ocean and measured:

1) δ 18O\delta^{18}O, the oxygen isotopic composition of planktonic foraminifera,

2) Ts, an estimate of summer sea-surface temperatures at the core site, derived from a statistical analysis of radiolarian assemblages,

3) percentage of Cycladophora davisiana, the relative abundance of a radiolarian species not used in the estimation of Ts.

Identical samples were analyzed for the three variables at 10-cm intervals through each core.

Here are their findings for RC11-120:

and for E49-18:

After some processing described in the paper, here are the power spectra for these data:

As you can see, there seem to be peaks at frequences of 100 ka, 42 ka, 23 ka, and 19 ka.

Berger 1989

This piece of work reviews evidence for the Milankovitch explanation of glacial cycles, at least up to 1989:

2 (1989), 1-14.

Abstract: Spectral analysis of geological records show periodicities corresponding to those calculated for the eccentricity (400 and 100 ka), the obliquity (41 ka) and the climatic precession (23 and 19 ka). It is precisely the geological observations of this bi-partition of the precessional peak, confirmed to be real in astronomical computations, which was one of the most impressive of all tests for the Milankovitch theory. Concerning the question of whether or not the observed cycles account for most of the climatic variability having periods in the range predicted by the astronomical theory, substantial evidence (from cross-spectral analysis, coherency analysis, and modelling) is provided that, at least near frequencies of variation in obliquity and precession, a considerable fraction of the climate variance is driven in some way by insolation changes accompanying changes in the earth’s orbit. The variance components centered near a 100 ka cycle dominate most ice volume records and seem in phase with the eccentricity cycle, although the exceptional strength of this cycle needs a non-linear amplification by the glacial ice sheets themselves and associated feedbacks.

As the insolation spectra change with latitudes and with the type of parameters considered, the diversity of the spectra of different climatic proxy data recorded in different places of the world over different periods is used to better understand how the climate system responds to the insolation forcing. This study of the physical mechanisms involved is also achieved through the analysis of the log-log shapè of the geological records and through the comparison, in frequency domain, between simulated climatic time series and proxy data.

The evidence, both in the frequency and in the time domain, that orbital influences are felt by the climate system, implies that the astronomical theory may provide a clock with which to date Quaternary sediments with a precision several times greater than is now possible.

Berger cites the aforementioned paper of Hays, Imbrie and Shackleton as important in finding empirical evidence for long-term climate cycles. He writes:

By the 1970s, the grounds upon which the first strong test of the Milankovitch theory was going to be based, were settled. Judicious use of radiometric dating and other techniques gradually clarified the details of the time scale (Broecker et al., 1968); better instrumental methods came on the scene for using oxygen isotope data as ice age relics (Shackleton and Opdyke, 1973); ecological methods of core interpretation were perfected (Imbrie and Kipp, 1971); global climates in the past were reconstructed (CLIMAP, 1976); and astronomical calculations were checked and refined (Vernekar, 1972; Berger, 1976b). Unfortunately it is not possible in this paper to give credit to all scientists responsible for the new picture of ice ages and their cause (refer to Imbrie and Imbrie, 1979; Berger, 1980a, 1988; Imbrie, 1982). Scientists were thus ready to show that the climatic features in the geological record go through precisely the same rhythms as do the orbital parameters of the earth, with sufficiently close links in time (in phase) to confirm a relationship of cause and effect.

In August 1975, at the climate conference of the World Meteorological Organization in Norwich, the spectral characteristics of two time series independently built – one in geology and one in astronomy – were compared, for the first time in the history of climate investigation. The joint efforts of Hays, Imbrie, and Shackleton (published in full in Science, December 1976) demonstrated that geologic spectra contain substantial variance components at many frequencies similar to those obtained by Berger from astronomical computations (published in Nature, January 1977).

The principal feature of the variance spectra obtained from a composite deep-sea core from the Indian Ocean (RCl1-120 and E49-18), spanning the past 468 ka, was a characteristic red noise signal over periods ranging from about 100 ka to some thousand years. Superimposed on this signal were several peaks representing significant responses at orbital frequency (Fig. 1). On both isotopic (ice volume) and temperature radiolarian spectra, peaks appeared for cycles roughly 100 ka, 42 ka, 23 ka, and 19 ka long.

The obliquity signal with a period of about 42 ka was gratifyingly consistent and detectable for the past 300 ka. The related climatic variations lagged about 9 ka behind the rhythm of the changing tilt of the earth. This consistency was less, however, in the older part of the record, because of a presumed less accurate time scale. This lead Hays et al. (1976) to adjust the time scale slightly (the ‘tune-up’ being well within the uncertainties of the dating) to obtain the same phase relationship going back to 425 ka BP.

A 23 ka cycle, corresponding with the precession of the equinoxes, also appeared strongly in the geological record. Although the uncertainties about dating produced relatively greater uncertainties in the phase relationships for this higher frequency, and the earth’s orbit was almost circular between 350 ka and 450 ka, the tune-up adjustments still showed that the climatic factors change closely in step with the wobble over the whole period. The most delicate and impressive of all tests for the Milankovitch theory, however, came with a closer examination of the astronomical solution. The precession of the equinoxes was split into two frequency components with periods of 23 ka and 19 ka. This was calculated both by Berger (1977) from his theoretical investigation of the astronomical theory and observed by Imbrie in the Hays and Shackleton records (1976).

Surprisingly, however, the dominant cycle in Hays et al.’s data had a consistent period of 100 ka, with the coldest periods of ice age activity coinciding dramatically with the periods of near-circularity of the earth’s orbit. This is exactly the opposite of what Milankovitch claimed, because in his theory eccentricity could only modulate the size of the precession effect, but the conclusion is compatible as far as the total energy effect is concerned (lower eccentricity, lower total energy received from the sun).

As the 100 ka year cycle is extremely weak in the insolation data set, with obliquity and precession sharing almost entirely the total variance, one of the mysteries which remained when trying to relate paleoclimates to astronomical insolations was to find how this cycle could be enhanced to become the dominant cycle observed in geological records.

Laskar et al

Laskar seems to be an authority on celestial mechanics. Relevant papers of his include:

  1. J. Laskar, The chaotic motion of the solar system. A numerical estimate of the chaotic zones, Icarus 88 (1990), 266-291.

  2. J. Laskar, F. Joutel and F. Boudin, Orbital, precessional and insolation quantities for the Earth from -20 Myr to + 10Myr, Astronomy and Astrophysics 270 (1993), 522-533.

  3. J. Laskar, M. Gastineau, F. Joutel, P. Robutel, B. Levrard, and A.C.M. Correia, A long term numerical solution for the insolation quantities of Earth. Astronomy and Astrophysics 428 (2004), 261-285.

In (2) section 4.1, they state that for the purpose of computing insolation, they neglect the secular variations of the semi major axis. There is a diagram for the variation of the semi-major axis in (2) (fig. 11) which seems to justify this. In (3) section 7, they state that an early version of this fact was known to Laplace and Lagrange!

In these papers they describe a rather complete computation of orbit elements, obliquity and precession for several tens of millions of years. There is some uncertainties, for instance due to the tides on Earth. In (2) 7.2 it is taken into account that the magnitude of the tides might change during an ice age. This is further discussed in [3] 4.2, but I’m not entirely convinced that they really really know by how much. They do seem pretty convinced that their calculations are correct from -20 My to present. If you try to extend the calculations to earlier times, you will eventually be up against chaos. In (3) section 11 they state that it is hopeless to try to get to the Mesozoic, before - 65My.

The best place to start reading these papers is probably (3), especially because of the interesting diagrams and tables.

Bennett 1990

Here is a curious paper on Milankovitch cycles and evolutionary biology:

Abstract: The Quaternary ice ages were paced by astronomical cycles with periodicities of 20-100 ky (Milankovitch cycles). These cycles have been present throughout earth history. The Quaternary fossil record, marine and terrestrial, near to and remote from centers of glaciation, shows that communities of plants and animals are temporary, lasting only a few thousand years at the most. Response of populations to the climatic changes of Quaternary Milankovitch cycles can be taken as typical of the way populations have behaved throughout earth history. Milankovitch cycles thus force an instability of climate and other aspects of the biotic and abiotic environment on time scales much less than typical species durations (1-30 m.y.). Any microevolutionary change that accumulates on a time scale of thousands of years is likely to be lost as communities are reorganized following climatic changes. A four-tier hierarchy of time scales for evolutionary processes can be constructed as follows: ecological time (thousands of years), Milankovitch cycles (20-100 k.y.), geological time (millions of years), mass extinctions (approximately 26 m.y.). “Ecological time” and “geological time” are defined temporally as the intervals between events of the second and fourth tiers, respectively. Gould’s (1985) “paradox of the first tier” can be resolved, at least in part, through the undoing of Darwinian natural selection at the first tier by Milankovitch cycles at the second tier.

Gallée et al 1991-1992

These papers present a model that simulates the last glacial cycle starting from insolation data over the last 122 kiloyears. Ice albedo feedback plays a crucial role in amplifying the Milankovich cycles. The importance of many other effects is studied, with many interesting conclusions. The effect of changing CO2 concentrations is neglected except in Section 6 of the second paper, and even here the feedback cycle involving temperature and CO2 concentration is not discussed.

There are many interesting figures. Here is one of the most basic: a graph of June insolation as a function of latitude for the last 130,000 years:

Huybers 2009

In this paper, Huybers proposes the following: The eccentricity period does not matter, that’s a red herring. Throw it away, don’t think about it. It’s all about tilt (obliquity). The glaciation 40ky period indeed corresponds to the dominant obliquity period. Sometimes, this period gets doubled or tripled, giving periods of 80 ky or 120 ky. These two cases he subsumes under the common name 100\sim100 ky.

The problems he wants to solve is: why does doubling and tripling occur, and why do we get sequences of single periods (as in early Pleistocene) , and sequences of repeated periods (as in late Pleistocene and Holocene). That is, following a single period, we will most likely get a single period, and following a 100\sim 100ky period, we will most likely get a 100\sim100 ky period.

As a preliminary and as supporting evidence he makes one strange observation about the δ 18O\delta^{18}O data. The curve γ\gamma describing this, which presumably also measures the total amount of ice on the planet, has the following strange property: The derivative γ\gamma\prime is negatively correlated to γ\gamma, with a time lag of 9 ky. The long time lag is curious, and he does suggest several possible physical explanations for this correlation, without chosing his favourite. It seems (though it’s not stated absolutely clearly in the article), that this correlation is mainly valid in the ablation phase, that is during de-glaciation.

An interesting fact is that the correlation persists through the change from 40ky to 100\sim100 ky periods. This seems to imply that the physics of ablation is the same in both cases. I think that Huybers takes this as evidence that the forcing is also the same. Since eccentricity can’t be blamed when the period is 40ky, eccentricity can’t be guilty even when the period is 100ky.

Based on this observation, Huybers writes down a time dependent difference equation, whose time dependence codifies a periodic forcing, with period 40 ky. He states that the difference equation has chaotic behaviour, but also that there exists a periodic orbit with period 40ky. Close to this periodic orbit, there are solutions with periods 80ky. These solutions alternate between 40ky of mild glaciation and 40 ky of severe glaciation. There are also cycles of period 120 ky (I don’t know if they are also supposed to be close to the basic periodic orbit) , passing through 40ky of interglacial, 40 ky of mild glacial and finally 40 ky of full glacial. So it seems that we can produce sequences of cycles close to these periodic orbits which reproduce the geological record. Unfortunately I did not entirely understand the mathematics of this part of the article.

The system is chaotic, so we should expect random shifts between the various periodic orbits. There is no particular external reason for the shift from 40ky to 100ky.

Miscellaneous online references


Let’s list books chronologically and comment on them. Here are some:

category: climate