Computing the Sunlight That Started and Ended the Ice Ages

A Geometry Problem With Ice Sheets In It

Nothing in this article required looking at anything. That is the part worth sitting with for a moment. Ice sheets two miles thick advanced across North America and Europe and then collapsed, nine times in the last million years, and the timing of it can be attacked with a protractor and an ellipse. No thermometer. No core. Four numbers about where the Earth is and which way it is tipped, and out the other end comes a curve whose peaks fall within a couple of thousand years of when the ice went away.

The quantity is insolation: watts per square metre arriving at the top of the atmosphere, averaged over one rotation of the planet. It depends on latitude, on the time of year, and on three slow properties of the orbit. How elliptical the orbit is, which is eccentricity. How far the spin axis is tipped over, which is obliquity. And which season the Earth happens to be passing closest to the Sun in, which is precession. All three drift on timescales of tens of thousands of years, and Milankovitch's argument, made in the 1920s and 1930s and set out at length in 1941 [1], was that the sunlight reaching the high northern summer is what decides whether last winter's snow survives until autumn.

That last figure is the reason this article is not a victory lap. The most conspicuous rhythm in the ice record of the last 800,000 years is a sawtooth roughly 100,000 years long. The insolation curve we computed does not have it. Not weakly. It has seven hundredths of one percent of its power in that band, which is to say it is absent, and everybody who works on this has known so since Hays, Imbrie and Shackleton published the spectra of a deep-sea core in 1976 and found the orbital frequencies sitting in the sediment exactly where theory said, plus one big peak that had no business being there [2].

So this is a study with a clean result and an open ending. The geometry is exact, the arithmetic checks out to thirteen decimal places against closed forms, the orbital periods fall out of a Fourier transform without being asked for, the alignment with the terminations is real and is also partly an artefact of how those terminations were dated. We will take all of that in turn.

What the Formula Actually Says

Start with a flat patch of ground facing straight up. The power it collects is the solar constant, reduced by the inverse square of the current distance to the Sun, times the cosine of the angle between the Sun and the vertical. Average that over twenty-four hours, throwing away the hours when the Sun is below the horizon, and you have the daily-mean insolation:

Three pieces do all the work. The declination \(\delta\), which is how far north or south of the equator the Sun stands, given by \(\sin\delta = \sin\varepsilon \sin\lambda\) with \(\varepsilon\) the obliquity and \(\lambda\) the Sun's longitude around the year. The sunset hour angle \(H_0 = \arccos(-\tan\phi\tan\delta)\), which is half the length of the day in radians. And the distance \(r/a = (1-e^2)/(1 - e\cos(\lambda - \varpi))\), where \(e\) is the eccentricity and \(\varpi\) locates perihelion.

The arccosine is the elegant bit. When \(-\tan\phi\tan\delta\) exceeds 1 there is no solution, which is polar night, and clamping it to 1 gives \(H_0 = 0\) and \(Q = 0\) with no special case required. When it falls below \(-1\) the Sun never sets, and clamping gives \(H_0 = \pi\) and \(Q = S_0 (a/r)^2 \sin\phi\sin\delta\). Two of the more dramatic facts about the polar regions come out of one clipped arccosine.

Converting \(\lambda\) into a calendar date needs Kepler's equation, because the Earth does not move around the ellipse at a constant rate. We solve \(M = E - e\sin E\) by Newton iteration and measure days from the northward equinox. The headline number avoids the whole question by using the June solstice, defined as \(\lambda = 90^\circ\) exactly, which is an astronomical event and carries no calendar convention at all.

The solar constant enters as a simple multiplier, and we used 1361 W/m² from the TIM radiometer recalibration of Kopp and Lean [3]. Most of the older insolation literature uses 1365 or 1368; Laskar's own parameter file ships with 1368. Every absolute number in this article scales linearly with that choice, and every relative number, every percentage and every spectral share, is untouched by it.

The Arithmetic

Plain numbers, then. This section has no argument in it.

Today: eccentricity 0.0167024, obliquity 23.43929°, longitude of perihelion 102.9179°. Those are row zero of the La2004 file.

The June solstice falls 92.758 days after the northward equinox. The Earth is then 1.016265 semi-major axes from the Sun. The declination is 23.43929°. At 65°N the sunset hour angle is 2.7646 radians, so the day is 21 hours 7 minutes long. Multiply it out: 477.936747 W/m².

The December solstice falls 276.248 days after the equinox, at a distance of 0.983707. Northern summer, equinox to equinox, runs 186.406 days. Southern summer runs 178.836. The difference is 7.570 days, and it is not a rounding artefact; it is why northern autumn feels short.

Annual means, today: 415.596 W/m² at the equator, 172.349 at either pole, 340.298 averaged over the globe. That last one is the famous \(S_0/4\), which is 340.250, raised by a factor \(1/\sqrt{1-e^2} = 1.00013951\) because an eccentric orbit spends more time far away but the inverse square law more than makes up for it.

Over the past million years: eccentricity between 0.004155 and 0.057815. Obliquity between 22.0762° and 24.4546°. The precession index \(e\sin\tilde\omega\) between \(-0.05768\) and \(+0.05573\). The resulting 65°N solstice insolation between 430.528 and 558.777 W/m², mean 493.164, standard deviation 23.488.

Today's 477.94 sits low. Only 28.1% of the past million years were dimmer.

Checking It Against Things We Did Not Compute

A formula that returns plausible numbers is worth nothing. The question is whether it returns the right ones, and the only way to find out is to compute the same quantity twice by routes that share no code.

Four closed forms are available. Integrating \(Q\) over the whole sphere and the whole year has to give \(S_0 / (4\sqrt{1-e^2})\), and the \((a/r)^2\) in the formula cancels against the Kepler weight \(dt/d\lambda\) exactly, which is why eccentricity enters only through that square root. At the equator the hour angle is always \(\pi/2\), the integrand collapses to \(\cos\delta\), and the annual mean becomes \(2 S_0 E(\sin\varepsilon) / (\pi^2\sqrt{1-e^2})\) with \(E\) the complete elliptic integral of the second kind. At either pole the integrand is \(\pi\sin\delta\) for half the year and zero for the other half, giving \(S_0\sin\varepsilon/(\pi\sqrt{1-e^2})\). And at the equinoxes and solstices the general formula reduces to one term.

We evaluated the elliptic integral by the arithmetic-geometric mean, which shares not one line with the quadrature it is checking. Both routes gave 1.506684231543093.

That still leaves one gap. All four closed forms descend from the same daily-mean integral, so a mistake in that integral would corrupt the check and the thing being checked identically. The fix is a Monte Carlo estimate that never uses it: scatter 4,000,000 points uniformly over the sphere and uniformly in time, compute the solar zenith angle at each from scratch, average the flux. Uniform in time means uniform in mean anomaly, which is what mean anomaly is for. It came back at 339.968 W/m² against the analytic 340.297, a gap of 1.50 standard errors, which for a random estimator is not a disagreement at all.

Then, last, the reference implementation. Laskar's La2004 distribution ships the Fortran subroutine cwj that computes exactly this quantity, and it does not clip the hour angle: it branches explicitly on three cases. We transcribed it, ran both over a grid of 181 latitudes by 721 solar longitudes, and took the largest disagreement anywhere.

Nineteen rows, nineteen agreements. The only one carrying a real gap is the Monte Carlo, and its gap is smaller than its own noise. The 65°N row looks like a 0.43% discrepancy until you notice it is entirely the solar constant: put 1365 back in and it closes to 0.14%, which is as close as a value published to three significant figures can be approached.

Working Notes From the Club Table

One Parameter at a Time

The three orbital parameters do not contribute equally, and they do not contribute independently. Holding two at their million-year mean and sweeping the third across its full million-year range gives a clean ranking, and one surprise.

Obliquity is the easy one. Push it from 22.076° to 24.455° and the 65°N solstice goes from 446.38 to 484.35 W/m², a span of 37.97. A steeper tilt means a higher summer sun and a longer summer day at high latitude, and the effect is monotone and unexciting.

Eccentricity alone gives 46.98 W/m², running from 487.97 at \(e = 0.004155\) down to 440.99 at \(e = 0.057815\). That sweep holds perihelion at today's angle, which happens to put the Earth near aphelion in June, so a rounder orbit brings northern summer closer to the Sun. Turn perihelion halfway round the circle and the same sweep runs the other way.

Which brings the surprise, except it should not be a surprise. Sweeping the perihelion angle through a full circle at mean eccentricity moves the answer by 55.46 W/m². At the largest eccentricity of the last million years it moves it by 114.53. At zero eccentricity it moves it by exactly nothing, because a circle has no perihelion and the angle to it is undefined.

So the precession term is not a signal in its own right. It is a carrier at roughly 23,000 years whose amplitude is eccentricity, which is exactly why the standard index is written \(e\sin\tilde\omega\) with the two multiplied together. Over the actual million years, letting one parameter vary in time and freezing the others, the standard deviations come out at 11.24 W/m² for eccentricity alone, 8.93 for obliquity alone, and 19.63 for the perihelion angle at mean eccentricity, against 23.49 for the real curve with all three moving. Those three variances sum to 1.072 times the total, which they should not do if the terms were independent, and they are not.

A Million Years of June

Run the whole thing backwards, one value per thousand years, and the curve in the bottom panel of Figure 3 is what Milankovitch spent twenty years computing by hand.

The extremes are worth naming. The brightest June at 65°N in the past million years was 958,000 years ago, at 558.78 W/m², with eccentricity at 0.0545, obliquity at 23.78° and the precession index at its most positive. The dimmest was 231,000 years ago, at 430.53, with eccentricity nearly as high at 0.0465 but obliquity down at 22.09° and the precession index almost as negative as it gets. The whole swing is 128.25 W/m², or 26.0% of the mean.

That 128 W/m² is a large number by the standards of anything in climate, and it is also completely different in kind from the numbers it invites comparison with, because it is neither global nor annual. Averaged over the sphere and over the year, the orbital variation is very nearly nothing, since eccentricity enters the global annual mean only through \(1/\sqrt{1-e^2}\). Across the full million-year range of \(e\), from 0.004155 to 0.057815, that moves the global annual mean from 340.253 to 340.820 W/m². A swing of 0.567. Milankovitch forcing does not change how much sunlight the Earth receives in a year. It moves sunlight between latitudes and between seasons and leaves the annual total where it was, which is why a theory of the ice ages built on it has to be a theory about where and when, never about how much.

What the Fourier Transform Says, and What It Does Not

Take the periodogram of that curve and the orbital periods appear without being asked for.

The precession band from 17 to 26 kyr holds 81.26% of the variance. The obliquity band from 35 to 55 kyr holds 17.74%. Together with the eccentricity band the three account for 99.08% of everything, which is a decent sign that the band edges were chosen sensibly and that there is nothing else in the signal.

The individual lines are where the check bites. The four strongest peaks come in at 23.83, 22.24 and 18.89 kyr, which is the precession triplet, and 41.71 kyr, which is obliquity. A thousand-kyr record cannot place a peak near 41 kyr better than about 0.85 kyr, so we ran the same analysis on longer stretches of the same file. Over the most recent 5 Myr the obliquity line lands at 40.992 kyr and the precession line at 23.701 kyr. Over 51 Myr the eccentricity lines land at 404.770 ± 1.606, 94.974 ± 0.088, 124.090 ± 0.151 and 98.839 ± 0.096 kyr. Nobody tuned anything. These are the periods the literature quotes [4][6][8], and they drop out of a Fourier transform of somebody else's orbital integration.

Now the part that does not work. The 75 to 135 kyr band in the insolation curve holds 0.07% of the power, and its largest peak in that band sits at 125 kyr. Widen the window to 5 Myr and the 300 to 500 kyr band, where eccentricity keeps 40.57% of its own power, holds 0.06% of the insolation curve's. Eccentricity is loud in the orbit and almost silent in the sunlight, for the straightforward reason that it enters the daily insolation at first order only as the multiplier on precession, and multiplying a slow envelope onto a fast carrier puts power into sidebands around the carrier rather than at the envelope's own frequency.

The ice record does not agree. The dominant rhythm since about 800 kyr ago is a 100 kyr sawtooth [7][9], and the interglacials that punctuate it have been catalogued and compared one by one [18]. Whatever produces it, it is not a 100 kyr line in the forcing, because there essentially is not one. Explanations on offer include a threshold mechanism that skips insolation peaks until the ice sheet is old and large enough [10], every second or third obliquity cycle being selected [11], ice-sheet hysteresis with isostatic rebound [12], and a simple rule based on how long it has been since the last interglacial [13]. None of those is something our calculation could produce, because none of them lives in the sunlight. They live in the ice.

The Strongest Objection: the Clock Was Set By the Thing We Are Testing

Here is the case against this section's own headline, made before the headline is.

We report that nine glacial terminations sit an average of 1.78 kyr from the nearest maximum of our insolation curve, against 5.52 kyr for dates thrown at random, and that no set of nine random dates in 20,000 draws matched it. Taken at face value that is a striking confirmation of Milankovitch. Taken with any care at all, it is close to circular.

The LR04 stack is an orbitally tuned age model [7]. Its authors did not date those cores independently and then notice they lined up with insolation. They aligned the benthic oxygen isotope record to a simple ice model driven by 21 June insolation at 65°N, which is the quantity in our bottom panel, and read the ages off. So when we measure the distance from a termination to an insolation maximum we are partly measuring the tuning procedure. The residual of 1.78 kyr is, to a significant extent, a statement about how tightly Lisiecki and Raymo's algorithm was allowed to pull.

How much is left once the circularity is subtracted? Some, and we cannot say how much from inside this study. Three things push back. Terminations II and V and the more recent ones have been dated independently by uranium-thorium on speleothems, which owes nothing to orbital tuning, and those ages agree with the tuned ones to within the tuning residual [14]. The tuning target is a single insolation curve at one latitude on one day, so it is mildly informative that the caloric summer half year, a different curve nothing was tuned to, still gives a mean offset of 2.89 kyr. And the alignment of the frequencies, which is the Hays, Imbrie and Shackleton result [2], was established on records with independent magnetostratigraphic age control and does not depend on tuning at all.

None of that rescues the specific number 1.78. We report it, and we report that it cannot be read as an independent test, and a reader who wants an independent test should go to the radiometrically dated speleothem record rather than to us.

There is a second objection, smaller and cleaner. Our curve has 46 local maxima over the past million years, spaced 21.5 kyr apart on average. There were nine terminations. Even if every termination sits on a maximum, 37 maxima produced nothing whatsoever, so the curve cannot be a sufficient cause of anything. At best it is a trigger that something else has to arm.

Which Summer? Three Answers, Three Cultures

Everything above takes "summer insolation at 65°N" to mean the daily-mean value on the June solstice. That is one choice. It is not Milankovitch's.

He used the caloric summer half year: the 182.62 days of the year that receive the most sunlight, averaged. Huybers later argued for the total energy received on every day whose insolation clears a melting threshold, on the grounds that what melts ice is joules, and that a longer summer is worth as much as a brighter one [15]. We computed all three on the same orbital elements.

The numbers are not close. Correlated against the orbital parameters over the million years, the solstice value has \(r = 0.915\) with the precession index and \(0.395\) with obliquity. The integrated summer energy has \(r = 0.393\) with precession and \(0.916\) with obliquity. The two quantities have almost exactly swapped. The caloric half year lands in between, at \(0.698\) and \(0.712\).

The reason is a cancellation that Huybers pointed out. When precession brings perihelion into northern summer the days get brighter, and by Kepler's second law that summer also gets shorter. Brighter times shorter is nearly constant, so anything that integrates over the season loses most of the precession signal, while a single-day snapshot keeps all of it. Obliquity has no such cancellation: a bigger tilt makes the summer brighter without shortening it.

This matters for more than bookkeeping. Before about 1 Myr ago the ice record is dominated by a 41 kyr rhythm, not a 100 kyr one, and the obliquity-heavy definitions of summer are much better at explaining that [15][16]. Whether the early Pleistocene looks like an obliquity world or a precession world depends on which integral you take, and both integrals are of the same sunlight.

What We Computed and What We Did Not

We computed watts per square metre arriving at the top of a transparent atmosphere. That is all.

There is no atmosphere in this model, no cloud, no albedo and therefore no ice-albedo feedback, which is the largest amplifier in the real system and the reason a 20 W/m² change in summer sunlight can move an ice sheet at all. There is no carbon dioxide, and the glacial cycles are accompanied by a swing of about a third in atmospheric CO₂ that our calculation knows nothing about. There is no ocean, no heat transport and no thermal inertia, so there is no lag between a forcing and a response, when a real ice sheet takes millennia to build and millennia to fall apart [17]. There is no ice sheet at all, so no isostatic rebound and no elevation feedback [12]. The solar constant is frozen at today's value for a million years.

What the calculation does establish is the shape of the forcing, and that shape is not negotiable. Summer sunlight at high northern latitudes really did swing by 128 W/m². It really did do so mostly at 23 and 41 thousand years. There really is almost nothing at 100 thousand years in it. Any theory of the ice ages has to start from those three facts, and the first two are the reason Milankovitch was right and the third is the reason the argument has run for a hundred years and is not over.

Where a different choice would change the answer

The definition of summer. Of everything in this study this moves the conclusion most, and section 10 is about nothing else. Solstice insolation says precession runs the show, 81% to 18%. Integrated summer energy says obliquity does, 83% to 14%. We led with the solstice because it is the conventional Milankovitch index and because it is the quantity LR04 was tuned to, and a reader who thinks melting is about joules rather than peak brightness should read Figure 5A right to left.

The latitude. 65°N is a convention, chosen because it is roughly where the great northern ice sheets had their southern margins. Nothing in the arithmetic picks it. We reran the whole spectral decomposition at five other latitudes: at 45°N the precession share climbs to 94.80% and obliquity falls to 4.57%, at 30°N it is 98.31% against 1.15%, and at 85°N it is 75.79% against 23.06%. Obliquity's leverage grows towards the pole and precession's does not, so the choice of latitude is itself an argument about which orbital parameter matters. One curiosity fell out of that sweep: above the Arctic Circle every latitude gives exactly the same band shares, because under the midnight sun the formula collapses to \(S_0 (a/r)^2 \sin\phi \sin\delta\) and \(\sin\phi\) is then a constant multiplier that cancels out of every ratio.

The orbital solution. We used La2004 [8]. La2010 revises the planetary ephemeris and changes eccentricity noticeably before about 50 Myr ago [16]. Over one million years the two agree to far better than anything here depends on, and swapping them would change no number in this article at the precision we quote.

The solar constant. 1361 against Berger's 1365 or Laskar's 1368 shifts every absolute watt figure by up to half a percent and shifts no percentage, no correlation and no spectral share at all. The 128.25 W/m² range would read 128.63 at 1365.

The obliquity reference frame. La2004 tabulates obliquity relative to the fixed ecliptic of J2000 rather than the ecliptic of date. The two differ by a small amount that grows with time, and using the other convention would move the absolute insolation values slightly while leaving the spectral structure alone.

Reproducing this

Two commands, and the first one is a download. You need Python 3.12 and numpy, nothing else.

curl -o analysis/data/INSOLN.LA2004.BTL.ASC \
  http://vo.imcce.fr/insola/earth/online/earth/La2004/INSOLN.LA2004.BTL.ASC
python milankovitch-insolation.py > milankovitch-insolation-output.txt

Expect about 21 seconds. Ours took 20.9 s on numpy 2.4.2 and Python 3.12.3. The orbital file is 4,590,090 bytes with SHA-256 beginning 3f13b9f8e69085ba, and the script prints the full digest so a reader can confirm they have the same data we did; it was retrieved on 14 September 2026. The termination ages are read from the LR04 age-model table at lorraine-lisiecki.com/LR04_MISboundaries.txt, cached in the same directory, and are not typed into the source. The seed 20260320 governs the Monte Carlo integration in section 4 and the random-date null in section 9, and nothing else; everything else in the file is deterministic arithmetic. Numbers differing from ours by more than the printed standard errors mean something is wrong, and we would like to hear about it. The interactive model runs the same insolation formula in your browser, and its default settings reproduce the headline numbers above.

References