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

RSS- & UAH-AMSU prediction

Welcome to my new blog!
Updates will be infrequent and mainly concerned with results of my hobby of (automated) time series analysis.

Let's start with the global temperature anomaly, as calculated from satellite measurements since 1979.

RSS- & UAH AMSU data with zero mean, Excels linear regression and 4th order polynomial fit,

as well as a simple model of linear trend + single cosine. The asterisks indicate El Niño events.

Most of the medium term variance of ~0.5 °C can be attributed to ENSO and the big volcanic eruptions, like El Chichón 1982 and Mount Pinatubo 1991. Curiously the cooling caused by Mount Pinatubo fell into a period which should have had a La Niña episode, but didn't actually, so the pattern is probably more regular than it would have been without the volcano.

The overall correlation between between the RSS and UAH series is high at 0.957,

but over the last 12 months, UAH has been on average 0.05 °C warmer than RSS.

Difference between the two datasets.

The autocorrelation of the series is significant and confirms the ENSO signal with a period of ~3.7 years.

Auto- and cross-correlation of the detrended data.

Now for the interesting part of using genetic algorithms to automatically build models and use them for forecasting. The models (neural net like recursive algorithms) use only the timeseries as input and try to predict one step ahead. During optimization, the input is mostly clamped to the data but sometimes replaced by the model output for a short period. The model fitness is usually the mean squared deviation from the data, but could take into account other desirable properties like the first derivative.
The candidate models have a high potential number of freedoms, and are thus prone to overfitting (modelling the noise), so we're not really interested in the global minimum, but the genetic optimization usually converges on a local minimum with a low effective number of freedoms, which models an essential property of the data and is thus useful for forecasting.


The final model is then run many times with added noise and the results are averaged, to get a representative output and estimate the error.
Here are the results of 3 optimization runs:

The models almost always exploit the ENSO regularity, which the manual analysis indicated, but model c also picks up some higher frequency details, and consequently has the lowest error bounds and highest R²=0.84, though this measure is not too meaningful, as it depends on the forecasting horizon and is calculated on the same data used during fitting.
To verify a model, its performance should be evaluated on unseen data, but this timeseries was so short, that I decided to use all existing data during optimization and the verification must be done on future data.
Model C is probably overfit, and model A and B show that the error bounds in the far future are most likely underestimated, as small timing differences lead to a big phase shift. Getting the timing right is often more difficult than the magnitude and direction.