Milankovitch cycles and glaciations

Milanković recognized almost 100 years ago that the quasi-periodic variations of earths orbital parameters influence its climate and cause ice ages and interglacials.
Calders and Roe showed how northern hemisphere summer insolation corresponds to the rate of change of ice volume and Huybers showed that deglaciations happened during times of high precession and especially obliquity, but in the late Pleistocene some obliquity cycles are skipped, resulting in 80 and 120 thousand year periods.

Earth obliquity and precession anomalies 500 thousand years in the past and future, based on Berger 1978 and 1991.

When high obliquity and precession coincide, the high northern latitudes get >40 W/m² more peak summer insolation than on average and the ice sheets shrink or disintegrate, while below average obliquity lets the ice sheets grow and global temperatures fall.

The orbital forcing function is a linear combination of obliquity and climatic precession.
The red peaks are deglaciations. The overall correlation of the two graphs is 65%.

Model input is obliquity and precession and model output is the inverted δ¹⁸O record, with zero mean during the Pleistocene, from Lisiecki and Raymo 2004 and Huybers 2007.
Lisiecki and Raymo use orbital tuning to constrain the age of the benthic records, while Huybers explicitly avoids this, consequently the two datasets are occasionally completely out of phase, but generally in good agreement, especially in the late Pleistocene.
As fitness function we take the product of the sum of squared errors (SSE) between the model and the two reference records from 2580 thousand years before present, with 1000 year timesteps.

Backcast of the model for the last million years, and reference data.
The early Pleistocene is also used for optimization, but not shown.

Forecast of the model for the next million years.
The error bounds should not be taken too seriously, as it's just the statistical error assuming the model is 100% correct, which it is obviously not.

For a better error estimate, one should compare several models:

Back- and forecast of 7 different models.
The models suggest there will be gradual cooling for ~60 thousand years, during the next two obliquity minima, with a potential partial melt in ~30 thousand years, while the subsequent obliquity maximum will probably get us out of the next ice age.

Download the consolidated data, including orbital parameters, insolation calculations, reference data and model output:
Milankovitch.xlsx