An Ice model driven by NH insolation but adjusted for dust albedo does a pretty good job at reproducing the last 8 glacial cycles. Details are described below.
Simulated Ice Volume model in blue compared to measured LR04 Ice Volume data.
Prior to the Mid Pleistocene Transition (MPT) roughly one million years ago, glacial cycles followed a simple 41Ky obliquity cycle. This can be understood if Northern and Southern Hemisphere ice sheets waxed and wained independently, because the average of the NH summer and SH summer peak insolation follows obliquity. Raymo et al. [1] suggested that the East Antarctic Ice sheet melted back to bare land prior to 1M years ago. The ebb and flow of the southern ice sheets (Patagonia included) with SH precession simply offset that in the Arctic because the precession component is always out of synch between North-South summers. The average of both just leaves the in-synch obliquity component which is why it dominates the Global Ice volume data before MPT.
Hubers (2) also pointed out that the integrated summer insolation which he calls total “caloric summer” also follows obliquity. The reason for this is due to Kepler’s law. Conservation of angular momentum during the earth’s elliptical orbit ensures that each hemisphere’s precession enhanced maximum insolation has a shorter summer by several days. So when calculating the total energy received at 65N or the poles, the precession term simply cancels out. Both the above effects are summarised in figure 1., which also highlights the role played by eccentricity in modulating peak insolation at each pole.
Figure 1. Maximum and total solar insolation calculated at the poles during last 600,000 years. Both the total annual insolation and the N-S asymmetry clearly show the underlying effect of the 41,000 obliquity signal. STOT is equivalent to total calorific summer and (Smax NP – Smax SP) leaves just the obliquity signal which affects both hemispheres equally.
Raymo proposed that once Antarctica became permanently covered in ice after MPT the symmetry was broken, and NH insolation thereafter became the dominant driver of glacial cycles. So why didn’t ice ages simply follow the NH Peak 65N insolation? This is the unsolved mystery of ice ages as to why they apparently switched instead to a 100Ky “eccentricity” cycle. Eccentricity modulates the strength of the precession term and each interglacial does indeed coincide with eccentricity maxima. However, no-one has properly explained why other larger NH insolation peaks fail to have much effect, while far smaller peaks eventually succeed in melting back the ice sheets (figure 2).
Figure 2. The last 800,000 years of glacial cycles. Ice Volume (LR04) in grey, Eccentricity in black, EPICA (Antarctic) temp (red), dust(purple), CO2(yellow). The current interglacial is most similar to that 400,000 years ago which also coincided with low eccentricity. Click to expand
In 2006 Gerard Roe published a paper “In defence of Milankovitz” [3] where he showed a good qualitative agreement between the rate of change of the SPECMAP[4] Ice volume fit with Arctic summer melting, as measured by maximum insolation at 65N. The agreement though was too good, partly because the SPECMAP fit itself had used the 65N insolation for dating the stack. The LR04 stack[5] however, contains many more geographically spaced cores which are independently timed. When I checked the LRO4 derived DV/DT against N-pole summer insolation I found that there was indeed a correlation but the agreement was not nearly as good as Roe’s. However what Roe did demonstrate is that NH peak insolation is a major driver in the ablation and growth of the ice sheets, but cannot on its own explain the ice volume data. So what is the missing link needed to explain the dynamics of recent ice ages?
Ellis and Palmer[6] proposed that the missing link was the dust albedo effect. As glaciations deepen and CO2 levels fall below 200ppm. The arid conditions cause less snow which combined with CO2 starvation cause die back of boreal forests and increased dust deposits over the ice sheets lasting thousands of years. You can see this effect in figure 2. They propose that this decreases the ice surface albedo sufficiently such that the next NH peak summer is sufficient to rapidly melt back the ice sheets caused by a chain reaction ice albedo feedback. One difficulty with this explanation is that Greenland dust data is only available for the last 200,000 years and so the assumption needs to be made that Antarctic dust is also a proxy for NH glacial deposits. This seems a reasonable assumption at least for for the last glacial period (figure 3).
Figure 3: Comparison of Greenland GRIP data with EPICA Antarctic data.
Willeit and Ganopolski [7] have also recently highlighted the importance of dust albedo and snow ageing during the last glacial cycle.
“The surface energy and mass balance of ice sheets strongly depends on the amount of solar radiation absorbed at the surface, which is mainly controlled by the albedo of snow and ice.
Both the snow ageing effect and the effect of dust deposition on snow albedo play a fundamental role in reducing surface albedo, particularly in the ablation areas”
They used an “Earth system model” – CLIMER-2 and were able to reasonably well model the last glacial cycle. Yet the details still remain obscure because such complex climate models are black boxes. Surely there must be a simpler physical explanation for the slow growth and sudden collapse of ice ages.
Does dust enhanced albedo play the critical role in recent glacial cycles? To answer this question, I decided to try to develop a simple dust model based loosely around the original Imrie and Imrie Ice model [8]. First I tried a very simple equation.
where b is a nonlinearity constant between ice growth and decay, S is NH insolation, 
is the ice albedo which depends on dust levels.
I have tried two versions for so far.
a) 
and
b) 
The first assumes that a critical accumulation of dust over one glaciation is the prime driver once NH induced ablation starts, whereas the second assumes that a critical threshold 15,000 years ago is required. Scale factors are needed to match insolation in W/m with Ice Volume in Benthic 
. I used b = 0.03. The model is written in Python.
Results
First we look at the integral dust results a) as shown in Figure 4. The dust albedo model improves the agreement between Insolation and -DV/DT
Figure 4. Shown below are the model results in blue compared to -DV/DT data in red. The upper plot shows again -DV/DT in red compared to the 65N Insolation. Below this is the LRO4 ice volume data and dust data scaled for comparison.
The agreement of the albedo weighted NH insolation with -DV/DT is now much better. But the real test of the model though is to reproduce the Ice Volume results by integrating forwards in time based on the start value 800,000 years ago. The only input for this is the new NH insolation and the recorded EPICA dust data. Figure 5 shows the result.
Figure 5. A comparison of the Model A simulated Ice volume by integrating -DV/DT forward in time This is compared to the the actual LRO4 Ice volume data. Overlaid in red is -DV/DT proving that the integration is sound.
In general the agreement is remarkably good considering the time scale. However we see two problems caused by overcompensation for dust albedo in the early cycles. Perhaps Antarctic dust was not fully coupled to that in NH.
The results from model b) which includes a 15,000 year time lag in dust contributions are shown in figures 6 and 7.
Figure 6. Overview for the delayed dust albedo model. The comparison of -DV/DT for ice volume and Dust induced albedo is shown below.
Figure 7: Simulated Ice Volume results by integrating the delayed dust albedo model B forward in time. This is compared to the actual LR04 ice volume data in black/red. as in Figure 5.
The agreement is better, but we see potentially another problem in that the holocene is now not so well modelled due to latent effects of dust from LGM. However, the timings and basic shape of glacial cycles are reasonably well represented. So I am fairly hopeful that a solution to the mystery of ice ages can be found.
This is still a work in progress but you can download the current software and all relevant data here
References
- M. E. Raymo et al. Plio-Pleistocene Ice Volume, Antarctic Climate, and the Global ?18O Record. Science 28 Jul 2006: Vol. 313, Issue 5786, pp. 492-495
- Peter Huybers. Early Pleistocene Glacial Cycles and the Integrated Summer Insolation Forcing, Science 313, 508 (2006)
- Gerard Roe, In defense of Milankovitch, GEOPHY RES LETTS, VOL. 33, L2470
- Imbrie, J., Hays, J.D., Martinson, D.G., McIntyre, A., Mix, A.C., Morley, J.J., Pisias, N.G., Prell, W.L., and Shackleton, N.J., 1984. The orbital theory of Pleistocene climate: Support from a revised chronology of the marine ? 18O record. In Berger, A., et al., (eds), Milankovitch and Climate. Hingham, MA: D. Reidel, pp. 269–305.
- E. Lisiecki &Maureen E. Raymo, A Pliocene-Pleistocene stack of 57 globally distributed benthic D18O records, PALEOCEANOGRAPHY, VOL. 20, PA1003, 2005
- Ellis, Ralph and Michael Palmer (2016) “Modulation of Ice Ages via precession and dust-albedo feedbacks” Geoscience Frontiershttp://www.sciencedirect.com/science/article/pii/S1674987116300305
- Matteo Willeit and Andrey Ganopolski, The importance of snow albedo for ice sheet evolution over the last glacial cycle, Clim. Past, 14, 697–707, 2018
- John Imbrie and John Z. Imbrie, Modeling the Climatic Response to Orbital Variations, Science, Vol. 207, No. 4434 (Feb. 29, 1980), pp. 943-953
is the ice albedo which depends on dust levels.
. I used b = 0.03. The model is written in Python.
