The Milankovitch cycles pt I: a climate history written in the celestial sphere

Mark Osborne
Jun 19
14 min read

Updated: Jul 1

As sure as night follows day, we know the Earth rotates on its axis. Clocks spring forward, fall back around the equinox, and the solstice that marks the longest day of Summer and shortest in Winter, tells us the Earth orbits the sun.

Of course that was only made obvious when Copernicus (1543) [1], from celestial observation, placed the sun at the centre of the solar system. Before that, sunrise after the dark of night, fruitful summers and freezing winters were the hands of the gods!

Less obvious are the subtle peroidic changes in these cycles, known as the Milankovitch cycles. Their effects on Earth's climate are dramatic, ultimately dictacting the timing of ice ages and warmer interglacial periods. Luckily, these cycles operate on long geological time scales of many thousands of years.

This dive aims to sit somewhere between qualitative descriptions and the rigourous mechanics of the Milankovitch cycles. A bit of history and just enough maths to shine a light on the science behind our Earth's evolving climate.

From spinning tops to wobble boards

The view of Earth spinning on an axis was first conjectured by Greek astronomers of Pythogoros and Platonic camps, Philopaus and Heraclides [C4th BC], in a challenge to the ancient picture of the heavens spinning about a fixed Earth, as painted by Aristotle.

Pytheas shadows cast at the solstice & equinox

Obliquity, or axial tilt was first measured by Pytheas of Marseille (350 BC) [2], using the trigonometry of shadows cast by a gnomon (upright of a sundial) at his name place. At the Summer solstice noon, Pytheas quoted a gnomon height to shadow length ratio 209/600, giving the intermal ray-to-earth, projection angle

By alternate angles, this gives the latitude of the location with respect to the orbital plane, but not to the equator.

The same measurement at the equinox then gives the latitude from the equator directly where it intersects the orbital plane. The gnomon now casts a longer shadow, with a height to length ratio round 0.94 (~209/222), a projection angle, and hence latitiude of 44.3°.

The difference gives an obliquity of 24.1°, a value larger than today's tilt of 23.5°, but close to that 2000+ years past, given the obliquity is now known to have been in decline over the last 10000 years. A difference of < 0.6° from the shadow of stick, remarkable!

Axial Precession, the slow roll of Earth's axis around a fixed pole, like a kids spinning top, was first identified by Hipparchus of Rhodes (190 BC).

The plane of the equator necessarily wobbles in concert with the axis, such that the position of the equinox on Earth's orbit, shifts slowly over time. The constellations at the March equinox are seen to creep westward with this precession of the equinox.

Westward shift of the equinox due to axial precession. Apdapted from Dbachmann CC BY-SA 3.0 link

Hipparchus (129 BC) compared location of the star Spica in the northern night sky with records of Timocharis (293 BC), and concluded the stars made a longitudinal shift around 1° in a century, a full revolution in 36000 years.

Despite such understanding of the heavens, the geocentric view of Ptolemy's Algemist (150 AD) [3] prevailed for centuries, where day, night and seasons were dictated by motions of the celestial bodies about a central Earth (geo = Earth, helios = Sun in greek). The Algemist provides the only records of Hipparchus and development of the epicycle-deferent (sphere-on-sphere) model of celestial motions.

Translations of the Indian text Surya Siddhanta (800 AD) suggest refinement of the "trepidation" of the equinoxes with 54° shift-and-return over a 7200 year cycle, or axial precession of 54" pa (2X54X3600/7200) and period 24000 years [4].

Tables of Ulugh Beg (1347), refine star displacements to 1° in 70 years, an precession rate of 51.4" pa and period of 25 200 years. Translations of the Persion astronomers text indicate a first to measure the obliquity at 23.5° [5].

Copernicus (1543) [1] later showed the precession of the equinox to be variable about a mean 50.2" pa, close to modern measurements of an axial precssion period around 26000 years.

Around the same time Fracastoro (1538) [6] appears the first to acknowledge a north-to-south, latitudinal shift in the precession of the equinoxes from a change in obliquity.

Aspsidal Precession, the wobble the orbital plane, was first identified by Hipparchus through observed anomalie in the Moon's apogee (farthest point from Earth) [3]. Ancient, Babylonian records were used to show it took 126007 days for the Luna eclipse to return to the same position in the night sky, within which time 4573 anomaly returns and 4612 orbits of Earth were noted.

Hipparchus recognised the difference arising from motion of the Moon's apsides (axes of orbit). The anomoly and orbital periods, 27.55 days (126007/4573) and 27.32, respectively, determine the precession period (T) from the "beat" frequency (ω) of the two cycles.

Modern measurements place the Moon's precession period at 8.85 years, while Earth's small, slow aspidal precession (relative its axial precession), was not identified until the precision of the 19th Century.

Primitive heliocentric views of the early Pythagorean astronomers were ultimately substantiated by Copernicus [1], who proposed that motions of the planets could only be explained by placing the sun at the centre of their orbits. Not until the following century, with Galileo's (1610) [7] invention of the astronomical telescope, was this view and motions of the five known planets Mercury, Venus, Mars, Jupter, Saturn confirmed [8].

At the same time Kepler published the first mathematical laws governing planetary motion (1609, 1619) [9], based on accurate sextant measurements of Mars' orbit by Tycho (Brahe) (1574) [10]. Newton (1687) followed up by identifying gravity as the force driving the motion [11].

Kepler's equal triangles for elliptical orbit of high eccentricity e = 0.75

Kepler's key finding? Triangles made from lines connecting the sun to the arc swept by a planet over a fixed time, are all equal regardless of distance between sun and planet. The planet must orbit faster at points closer to the sun than points farther from the sun, with a the crucial conclusion, that the orbit must be ellipitcal.

Eccentricity is defined by the semi-major and semi-minor axes of the elliptic orbit, governs the distance of closest approach to the sun, the perihelion and furthest, the aphelion, and determines the lengths of seasons. Kepler's value of 0.018 is today measured at 0.0167.

A calculated enlightenment

By C18th, accuracy of instruments and calculus accelerated discovery of new motions, with 3rd Astronomer Royal, Bradley [1747], first to identify a nutation of the obliquity; an 18 year, 9" wobble in the rotation axis from the pull of the moon [12]. His measurements allowed later derviations of an overall decline in axial tilt of around 0.45" pa [13]. Limits on the tilt came later.

Lagrange and Laplace (1766-1784) performed the first calculations of the six planetary orbits including "three body" interactions, between the planets themselves and the sun [14]. A planet's motion is then coupled to all others, forcing periodic changes in orbital inclination and the longitude of its node (pivot point), the eccentricity and longitude of the perihelion (point of closest approach).

The legacy? By treating eccentricity (and other periodic elements) as a vector changing in amplitude and direction over time, the Cartesian components (h and k) could be solved, each as trigonmetric series with amplitudes (α), frequencies (𝑔) and phases (β) defined by the strength of (gravitational) coupling between planets and its rate of change as planets orbit to and from each other.

Importantly changes in eccentricity were found to arise from changes in the semi-minor axes of the ellipitic orbit, while the semi-major axis remained invariant.

Eccentricity periodicity of Le Verrier vs Berger & table of highest amplitude components (Berger)

Uranus, discovered at his time by Herschel (1781) [15], was ultimately included in orbital calculations of Le Verrier (1843) that revealed (mass) corrected an eccentricity max and min of 0.06 and 0.003, respectively. Le Verrier (1848) went on to predict the existence Neptune, from irregularities observed in the orbit of Uranus [17].

Refined apsidal and axial precession rates of 11.66" pa (T = 360X3600/11.66 ~ 111150 years) and 50.1" pa (T = 25868 years) were reported, respectively by L.B. Francoeur (1812) and Delambre (1817) and Adhémar (1842) recognised the sourjon times combined for an overall "climate" precession period of 21000 years (360X3600/61.67) [18,19,20].

Calculations of the orbital eccentricity, inclination and longitude of the perihelion were subsequently corrected by Stockwell (1873), with the inclusion of all 8 planets in his celestial calcs [21].

By the 20th century, Pilgrim (1904) completed calculations on all secular pertubations (long term changes) in Earth's orbital properties, over extended periods (-1000000 to +40000 years), to reveal min & max of the obliquity of 22.1° and 24.5°, respectively [22].

With the current 23.5° obliquity near-midway and decreasing at a mean rate around 0.42" pa, the obliquity period is ~ 41000 years (2X(24.5-22.1)*3600/0.42).

From ice-age to nice-age & back again

The Milankovitch cycles: past and future orbital elements & resulting insolation Q. Data generated from Laskar (La2004)[27]

Milankovitch (1920) was the first to calculate changes in the insolation from a over 600K years before 1800, from the most accurate values of eccentricity, obliquity and precession at the time [23,24,25].

A 100 kiloyear period on the eccentricity revealed in the Milkankovitch cycles was was confirmed in a +/-5 million year projection by Berger (1976) and a final 405 kiloyear modulation was found superposed in calculations of Laskar (2010) extending back 50 million years [27].

The 100K cycle is seen in the mean rate of secular frequencies given by Berger (fig), while the gravitational pull of Venus and Jupiter on Earth's orbit, give rise to the 405K beat period from the difference in frequencies.

Finally in a full circle moment, location of the shortest distance from Earth to sun, the longitude of the perihelion, is defined with respect to the "first point of Aries", set by Hipparchus (130 BC) at the March equinox [28].

Heliocentric precessions showing first point in Aries set by Hipparchus & in Pisces today

In the Earth centred view of the time, the sidereal year started here at a geocentric longitude of 0°. As the sun moved counterclockwise, the north entered summer around 90°, reached the autumn equinox near 180° and the winter solstice about 270°, where the south approaches the height of summer when Earth is closest to the sun, the perihelion.

Today, precession of the equinox at 50.1" pa has moved this first point ~30° west from Aries to Pisces in over 2000 years (50.1*(2026+130)/3600).

And the heliocentric system now flips the first point 180° such that the perihelion resides around 270-180 ~ 90°.

The equinox lands around March 20 (19 to 21), while perihelion occurs around January 3 (2 to 5), a difference of 76 days or 75° (76X360/365). Perihelion is then 360-75 ~ 285° and 105° in geocentric and heliocentric frames, respectively. A mean anomaly is applied varying between +/- 1.915°, as a sine dependent on the number of Julian calendar days after 2000 (26X365.25 + 2451545) and currently -1.88° [29]. The final longitude of the perihelion is then 105-1.88 ~ 103°.

That's it, the orbital elements defining the Milankovitch cycles are covered and collated in Table 1 below!

Cycle	Today	Max	Min	Period
Slow 𝑒	0.0167	0.012	-0.012	405000
Fast 𝑒	0.0167	0.058	0.005	100000
Obliquity ϵ	23.4	24.5	22.1	41000
Perihelion (apsidal)	103	0	360	112000
Equinox (axial)	First Point in Pisces moves to Aquarias by 2600			26000
Climate cycle ϖ	Summer solstice aligned with perihelion in 1249			21000

Table 1

Shedding light on our sun driven climate

So how do these orbital elements and the periods over which they change lead to the climate cycles, from ice ages to warm interglacial periods, that we see printed in the geological history of Earth?

The ancients looked to the sun gods. From Helios (Greece) and Sol (Norse) in Europe, to Ra (Egypt) in Africa, from Mithra (Persia), Syrus (Hindu), Xihe (China), Amaterusa (Japan) across asia, to Kinich Ahau (Mayan) and Inta (Inca) in the Americas, the sun has been worshipped for bringing light, warmth, life, harvests.

The amount of warming the sun brings depends on the earth-sun distance, with Herschel (1832) the first to recognise that the heat received from the sun would diminish with a decreasing eccentricity of earth's orbit [30].

The heat is inversely dependent on mean earth-sun distance (strictly the square, Lambert (1760)[31]) over a revolution. As the minor axis increases and earth's orbit becomes more circular with decreasing eccentricity, so the mean distance increases and total heat decreases.

The inverse square dependence on earth-sun distance gives a current 7% difference in heating between perihelion (shortest) and aphelion (longest), but as high as 25% and low as 2%, as the eccentricity cycles between max and min limits over 100K year timescales.

Insolation, the total solar energy received at the top of Earth's atmosphere per unit area and time, was first determined by Pouillet (1838)[32] from the temperature change of a fixed container of water exposed to the sun (the pyrheliometer) with extrapolation to the top of the atmosphere assuming exponential attenuation (as per Bougeur(1729)[33]).

Measured then at 1.7366 calories/cm2/min ~ 1211 W/m2 (4.184 J per cal), the solar constant is now measured by satellite at a mean of 1361 W/m2 and varying by only 0.1% over the 11 year sunspot cycle.

The solar constant then allows calculation of insolation at any point across the Earth and its variation across seasons (an orbit of the sun) and long term, through the 21 kiloyear precessions of the equinox, 41 kiloyear cycles of obliquity and 100+ kiloyear changes in eccentricty.

Simple geometry gets you so far, but the projection of a 2D disc of illumination on a 3D spherical earth, ultimately requires some complex derivation.

In a toy model (see fig), the curvature of the Earth can be thought as "spreading" the solar consant, S0, over an increasing projected area with increasing latitude.

For a latitude ϕ and obliquity ϵ, the irradiance is reduced by cos(ϕ +/- ϵ) at aphelion/perihelion.

The total insolation then depends on the length of day, the fraction of a 24hr (2π) rotation spent exposed to the sun.

Tilt of the rotation axis makes this half a rotation +/- an arc length east and west that depends on tilt ϵ and increases from 0 at a latitude matching the tilt (12hr day) and π at latitudes above ϕ = π/2 - ϵ (24hr sunlight) and below ϕ = -π/2 + ϵ (no sunrise).

Between aphelion and perhelion, the arc lengths are further modulated by the ecliptic longtiude, θ, locating the time of year on Earth's orbit. Approaching the equinoxes (θ = 0, π), where day-night lengths are equal, the arcs reduce to 0 as sinθ, leaving the half arc and 12hr day.

Finally the insolation is ultimately modulated by the inverse square of the Earth-sun distance RE with resepct to the orbital mean R0 to which the solar constant applies.

In practice the length of day is determined from integration of the solar zenith angle, Θ, between hour angles, +/- h0, representing sunrise and sunset.

The zenith angle is equivalent to that between a perpendicular to the Earth surface and the direction of suns rays and derived from the spherical law of cosines.

Since the perpendicular depends on the latitude ϕ, and the direction of the sun depends on its declination δ and hour angle, h, in the day, so these are found in the cosine expression.

The mismatch between orbital and axial precessions adds a phase, between "first point in Aries" (θ = 0) , at the spring equinox and perihelion, as defined by the longitude of the latter, ϖ.

The limits on the hour angles depend, in the same way as the toy model, on latitude, ϕ, and obliquity, ϵ, with modulation by the ecliptic longitude, θ. So again, at full tilt toward the sun, the high arctic never sleeps (maximum insolation) and the antarctic never sees sunrise (zero insolation).

Integration between limits gives the insolation

Insolation across Earth latitude and a year orbit

For current orbital elements, 𝑒=0.0167, ϵ=23.4°, ϖ=103°, the contour map of the daily insolation (W/m2) over the year (0 < θ < 360°), from the spring equinox (θ=0), shows the expected pattern of solar heat with only a small north-south asymmetry in peak solar due to eccentricity; the southern summer being at perihelion, north at aphelion.

Given the eccentricity, obliquity and precession change with time, so the insolation will vary accordingly.

The Milankovitch Cycles revealed

Milankovitch identifed that the large land mass in the northern hemisphere, would be more responsive to changes insolation, through a lower heat capacity and faster land-air heat exchange than the ocean that fills the south.

Then at a latitude, ϕ = 65°N and summer solstice, θ = 90°, the hour angle aproaches h0 ~ π, and for a small eccentricity, 𝑒 << 1 and obliquity, sinϵ ~ ϵ, the insolation reduces to

where the orbital elements are strictly functions of time as they vary over their respective periods and the dependency of insolation on the key components 𝑒(t), sinϖ(t) and ϵ(t), as seen in Milankovitch-type analyses (fig above), is made clear.

Correlations between insolaiton, orbtial elements and proxy temperature records

So what? How's this help understand long term changes in climate of the past and climate change of the future.

Ice cores provide a best reference for Earth temperatures, extending over 400K years past (Vostock [34]).

Analysis of the last two ice ages, ending with interglacial warming, around -20K year, shows generally good correlation between solar heating of north (Q at 65°N and the temperatue anomly (ΔT from the 1850-1900 mean).

Associations are observed between the short initial cycle of rapid warming, cooling and warming (25K years) and the large peak-to-trough amplitudes in the insolation that follows precession.

Warming is then correlated with shorter, but hotter summers where the solstice aligns with perihelion, while cooling occurs during longer, but cooler summers, as precession brings the solstice in alignment with aphelion.

Duration of the ice-age (100K year) and a slower and deeper decline in temperatures trends with the reduction in eccentricity, with low warming-cooling cycles (40K years) matching obliquity.

Temperature swings appear larger and faster at high eccentricity (-225K years of the 405K cycle) and the current period relative stability and slow cooling is expected to persist while eccentricity decreases to its minimum in around 20K years.

Currently the north tilts toward the sun at aphelion, for long, warm summers, and away at perihelion, for short, cool winters. In 13000 years of the 21K year precession cycle, the north will tilt away from the sun at aphelion, for long, freezing winters and at perihelion, short, intense summers.

The mechanism by which the climate tips into an ice age is not well understood. Early theories of Adhémar (1842), based on precession, and Croll (1864), adding effects of eccentricity on the changing lengths of seasons, landed on long cold winters at aphelion as the trigger for glaciation [20, 35]. Observations of "arctic cricle" snow lines by Buch (1800), Forbes (1853) and Murphy (1869), drew conclusions the advance of the snow line followed from longer, cooler summers [36].

Milankovitch showed that, where precession brings the summer solstice in alignment with aphelion during periods of peak in eccentricity and falling obliquity, the slow orbit far from the sun, and lower angle of solar projection combine to make summers sufficiently cool and long to prevent polar ice melt.

Velocity & angular momentum relationship on the elliptical orbit

Time spent at aphelion is inversely dependent on speed at the turn point. Then from conversation of angular momentum at aphelion and perihelion, the ratio of time spent at these extrema is simply the ratio of their distances (Table 2).

Thus at maximum eccentricity, the summer at aphelion lasts 158/141 ~ 1.12X longer than winter; an extra 11 days (0.12X365/2.12) compared to just a 5 day (152/147 ~ 1.03X, 0.03*365/2.03) difference today, and less than 2 days at minimum eccentricity.

Eccentricity Limit	𝑒	Semi-major a Mkm	Perihelion rp Mkm	Aphelion ra Mkm
Max	0.057	149.597	141.070	158.124
Now	0.0167	149.597	147.099	152.095
Min	0.005	149.597	148.849	150.345

Table 2

A minimum insolation of around 430 Wm-2 (-120K & -230K yrs) is reached when obliquity dips to 22.1° and a maximum near 580 Wm-2 when precession aligns the summer solstice with perihelion (ϖ = 90) and obliquity is maxed out.

The matrix of solar irradiance and the "normalised" product of irradiance x summer to winter lengths ratio (Table 3), for the eccentricity, obliquity and precession extremes (summer at aphelion ϖ = 270° & perihelion ϖ = 90°), highlights the key role that obliquity plays in the large swings in solar heating that control long term climate changes.

Changes in insolation due to eccentricity and precession of summer between aphelion and perihelion are normalised by changes in season length, leaving the extent of Earth's tilt on its axis, the obliquity as the standout element. The central conclusion of Milankovitch!

ϖ	𝑒	ϵ 24.5°	ϵ 22.1°	summer: winter	ϵ 24.5°	ϵ 22.1°
270	0.05	472	426	1.11	521	470
270	0.005	519	468	1.01	524	473
90	0.05	577	£520	0.90	522	471
90	0.005	529	478	0.99	524	473