Glacial Cycles, part 1

There’s really no doubt that astronomical cycles have influenced the growth and decay of ice on planet earth for the last 5 million years or so. The subject came up recently, and there seems to be a lot of confusion on the issue, so let’s take a closer look at the influence of astronomical factors on earth’s cryosphere.

First, how can we be sure that astronomical cycles have any influence at all? The answer is, that the growth and decay of ice shows cyclic influences with the same periods as astronomical cycles. The statistical significance is so high, and the synchronization so persistent, that frankly it’s pretty much impossible that there’s no relationship. And I think we can safely conclude that it’s not earth’s glacial changes that are causing the astronomical cycles.

For instance, here’s the data for delta-oxygen-18 from a stack of 57 ocean sediment cores, which is considered an excellent proxy for global ice volume, known as the “LR04 stack” (from Lisiecki, L.E., & Raymo, M.E. 2005. A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records. Paleoceanography, 20, 1003–1019, doi:10.1029/2004PA001071). Up is less ice (so generally warmer conditions), down is more (generally cooler), the distant past is on the left and the present is to the right:

Over the last 5 million years, the very-long-term trend is of a gradual increase in global ice. Let’s remove the very-long-term trend by subtracting a cubic polynomial fit, leaving this:

This emphasizes that over the last 5 million years, we’ve also seen an increase in the size of the fluctuations of global ice volume. The glacial cycles have been a lot bigger over the last million years, than they were prior to that.

If we Fourier analyze the detrended data we get this (sorry, the x-axis is mis-labeled “period” when it’s actually “frequency”):

I’ve indicated the important cylic (and pseudo-cyclic) peaks in the Fourier periodogram. The strongest is at frequency 0.024 cycles/kyr, period about 41,000 years. This is the same as the cycle of obliquity, the tilt angle of the earth’s spin axis relative to its orbit. We can see it plainly in a Fourier periodogram of earth’s obliquity over the last 5 million years:

There’s also a cluster of small peaks in the range 0.042 to 0.045 cycles/kyr (periods 22,000 to 24,000 years) and a small peak at 0.053 cycles/kyr (period 19,000 years) that are all coincident with periods in the changes of precession, the orientation of earth’s spin axis relative to the longitude of perihelion (closest approach to the sun) of earth’s orbit. We can see these peaks in the Fourier periodogram of the earth’s precession factor (we’ll define that precisely in the next installment) over the last 5 million years:

I also adorned the first periodogram with an arrow labelled “???”. This is a broadband response around frequency 0.01 cycle/kyr (period about 100,000 years), but its attribution is unclear.

Some may object that ice-core and sediment-core records are often “tuned” to orbital variations, i.e., the timing is adjusted to match orbital changes. However, the presence of orbital-cycle variations has (time and again) been confirmed in records which are not orbitally tuned (see e.g. Huybers, P. 2007. Glacial variability over the last two million years: an extended depth-derived age model, continuous obliquity pacing, and the Pleistocene progression. Quaternary Science Rev., 26, 37-55, doi:10.1016/j.quascirev.2006.07.013) — there’s just no doubt about it, these astronomical cycles are clearly imprinted on paleoclimate data. Both obliquity and precession have an undeniable effect on glacial changes.

We can identify the times at which these fluctuations are stronger and weaker with a wavelet analysis. Here’s the logarithm of the wavelet power (using the WWZ, or “weighted wavelet Z-transform”) as a function of time and frequency (NOTE: in this graph, unlike those preceding, time goes from right to left so the present day is at the far left):

There’s stronger response recently (especially in the last 800,000 years or so). In part it’s a purely statistical effect, due to the fact that we have more data points for more recent times. But there’s also a stronger response due to truly stronger periodic fluctuations recently, especially for the precession frequencies and the unidentified 100,000-year fluctuation. Only the 41,000-year fluctuation is ubiquitous in this data set, in fact before about 800,000 years ago it dominates glacial changes, which leads to that period being called the “41kyr world” (Raymo, M.E., and Nisancioglu, K. 2003. The 41 kyr world: Milankovitch’s other unsolved mystery. Paleoceanography, 18, 1011-1016, doi:10.1029/2002PA000791).

We can even use the wavelet analysis to estimate changes over time in the strength of the astronomical cycles. We can do this for obliquity and precession themselves, to gauge the changing strength of the driving influence, and to the LR04 stack to gauge the strength of the glacial response. Here, for example, is normalized amplitude of the 41kyr obliquity changes compared to normalized amplitude of the 41kyr response in the LR04 stack (time goes from left to right, with the present at the far right):

There’s a lot of commonality to the changes, not only on long (million-year) time scales but even on very brief (100,000-year) time scales as well (interesting to refer to 100,000-year time scales as “brief”). However, there’s less match between the amplitude changes in the most recent million years or so, yet another indicator that the response of glacial changes to astronomical cycles is itself changing over time.

So without doubt astronomical cycles (specifically, obliquity and precession) have profoundly affected glacial changes — essentially, Milankovitch’s overall thesis has been proved correct. The question remains, how do they do so? And why does the 41kyr obliquity cycle dominate most of the time? And wherefore has the glacial response changed over time? And what is the cause of that 100,000-year pseudoperiodic response over the last 800,000 years or so? Stay tuned — we’ll look into those questions in our next installment.

24 responses to “Glacial Cycles, part 1

  1. Your third graph x-axis is incorrectly labeled as Period instead of Frequency.

    [Response: Oops.]

  2. Can you say something about the source of paleo orbital data? How do we know what the obliquity, precession, etc, was at some given time in the past?
    I know we can detect current 41kyr nutation in astronomical observations. But that doesn’t apply to, say, 3Mya.

    [Response: The paleo orbit data is based on computing the gravitational dynamics of the solar system (essentially, solving the N-body problem numerically). I used the data from Berger & Loutre]

  3. You mentioned the tuning of these stacks, but can you elaborate on why you think it is okay to ignore it? This procedure is sometimes derided as “wiggle matching,” and I have literally seen it done in a way that the obliquity cycle is overlain with delta-18-O and peaks are moved to match the obliquity peaks. This wiggle matching produces a “model” of sedimentation rate, but sometimes when you see the implied sedimentation rate it is jumpy, which doesn’t seem physical. That said, maybe the tuning is a good way to transform from depth to time, but then to do time series analysis on the result and claim anything about the frequencies associated with the tuning target seems misguided to say the least. Or maybe I’m missing something?

    [Response: Because the match between sediment core fluctuations and the obliquity/precession cycles has been well established using cores that are *not* orbitally tuned. And, orbital tuning won’t create exactly the *combination* of frequencies (as is clearest with precession cycles) that match the astronomical cycles. Nor can it create the match between changes in the *amplitude* of astronomical cycles and glacial response.

    In my opinion, the LR04 stack is one of the most complete records of glacial changes over the last 5 million years, so it’s a good choice for this exposition.]

    • I just went back and re-read Lisiecki & Raymo 2005 to refresh my memory of their method. I think they did a great job in combining the many sediment cores into a single stack. I’m still just a bit concerned with interpreting time series analysis of the stack resulting from their tuning process, at least around the frequencies that are present in the tuning target. I agree with your first point: there’s no doubt that these astronomical variations influence the global ice volume, and therefore there’s reason to believe that tuning to the insolation is justified. I also think your third point is a good one, the matching of amplitude is useful evidence that the signals are connected. I disagree with your second point though. The orbital tuning will absolutely produce spectra with the combination of frequencies that match the astronomical cycles. If you tune only to obliquity and get a lot of power at the precession frequencies, that is evidence that precession is really a significant signal, but if you tune to the total insolation at 21 June 65N, then there is no doubt that the result will have obliquity and precession signals. You can test this by tuning noise to a composite of oscillating signals; the tuning will produce spectral peaks at the same frequencies where the composite has peaks.

      [Response: Upon reflection I conclude that you’re right, it will produce the same mixture of frequencies as the signal to which it’s tuned.

      But of course it won’t reproduce the amplitude modulations, and yes the astronomical cycles are confirmed in non-orbitally-tuned records. We should also bear in mind that there are limits to how far one can move the wiggles to do the matching, beyond which there’s evidence of problems, but as far as I know no such evidence exists. The conclusion is, as you suggest, that the orbital tuning is justified.]

  4. While I agree that orbital tuning is both warranted and useful, for purely illustrative purposed like this blog post using a non-tuned dataset will be much more powerful. It’s all very well saying that tuning makes no difference – far easier to exclude the problem and cut off the bleating before it even starts.

    I suppose you just went for the longest high quality record, rather than picking a shorter, less certain record.

    Sadly, the amplitude discussion completely lost me, but that’s nobody’s fault but mine.

    [Response: Can anyone suggest an alternative of comparable (or close to) quality and duration, which is publicly available and not orbitally tuned?

    As for the bleating, that will continue regardless.]

  5. Tamino, I agree with you about the quality of Lisiecki and Raymo 2005. I also agree that the fact that an effective stack can be established by orbital tuning establishes an argument from coherence that the pattern is found in the original data rather than simply imposed on it. Nevertheless, it only provides weak support for that conclusion. This is particularly true in that much of the data in the stack is not continuos.

    Further, I do not think your point about the relative strengths of the signals is well made. The stacks in Lisiecki and Raymo were tuned to a model of 21 June insolation at 65 degrees North. Obliquity is the major contributor to variation in that insolation. In tuning the stacks to that model, therefore, the effect would be similar to tuning to a model of obliquity, and only using a model of precession to break ties. That would account for the relative strengths of the obliquity and precessional signals in the data.

    Further, the fact that the stack was orbitally tuned detracts from the otherwise excellent post for purposes of exposition.

    Is it possible, therefore, for you to show a similar analysis of paleo data that has not been orbitally tuned – perhaps one of the long ice cores from Vostock?

    [Response: I’m afraid you’re mistaken. First, orbital tuning can increase the response by matching the signal times, but there’s a limit to how far it can go — it simply can’t inflate the physical amplitude of the response, which dominates the period analysis once the tuning is done, and it certainly can’t impose temporal variations in the physical amplitude of the response.

    But mainly, obliquity is *not* the major contributor to variation in Jun21 insolation. Precession is overwhelmingly dominant, while for most of the Plio-Pleistocene it’s obliquity which dominates glacial changes. This is a clear result from my own analysis, and is testified by Huybers (2006, Science, 313, 508-511), who says:


    For example, average insolation on the 21st day of June at 65-N has 80% of its variance at the precession periods (1/21 ky +/- 1/100 ky). The caloric summer half-year at 65-N, defined as the energy received during the half of the year with the greatest insolation intensity (4), also has more than half its variance in the precession bands. But a major problem exists for the standard orbital hypothesis of glaciation: Late Pliocene and early Pleistocene glacial cycles occur at intervals of 40 ky (8–11), matching the obliquity period, but have negligible 20-ky variability.


    In fact that’s one of the main motivations for Huybers to propose his “integrated summer insolation forcing” model.

    But what the heck, I’ll add Vostok to the mix.]

  6. As an alternative to Vostock data, and if you want a record of similar length, several of the sediment records used in Lisiecki and Raymo extend for the full 5 million years. In particular, ODP 659, 846, 982 and 1143 may be suitable for the task.

  7. David B. Benson

    Thanks. First time to see the Pliocene varations for me.

  8. I just wanted to say that I’m impressed (and pleased!) with the level of thoughtfulness and knowledge shown in reader comments. Lucky me! But be prepared, that in the following posts I don’t want to leave those who are not so well-informed in the dust.

    Even so: compare to the comment threads at WUWT (particularly about Marin Hertzberg’s post).

  9. Tamino, you said:

    “I also adorned the first periodogram with an arrow labelled “???”. This is a broadband response around frequency 0.01 cycle/kyr (period about 100,000 years), but its attribution is unclear”

    Well, this 100 000 year cycle is the ECCENTRICITY CYCLE of the Earth Orbit around the Sun: The orbit oscillates between a more elliptical and a more circular orbit every (approximately) 100 000 years.

    A site that show the astronomical cycles, solar insolation, CO2, CH4 and temperature record is the Vostok Viewer:

    Just activate the option “eccentricity” and a blue line showing Earth eccentricty vs. time will be shown.

    [Response: The attribution to eccentricity is possible, but in doubt, as I’ll discuss in an upcoming post.]

  10. From Peru: You’ve complimented me by sharing my project (the vostok viewer) on this blog, but I should add that that project began here a couple years ago when I asked Tamino to define climatic precession. I created an illustration to help me understand his answer. Later Tamino showed me where to find Laskar’s orbital data. So, my project exists because I follow (and have a lot to learn from) this blog. Please consider my project an illustration; it is not an analysis, and not of the caliber we get here. And thank you, Tamino.

  11. On a slight tangent, there was an excellent documentary called “Men of Rock” presented by the geologist Prof Ian Stewart, which profiled the work of Scottish scientist James Croll who originally came up with the theory of climatic change from orbital forcing. Croll was actually a janitor at Edinburgh University and was self-taught in celestial mechanics and very interested in Louis Agassiz theory of glacial periods. It was actually Croll that first proposed the relationship between glacial cycles and orbital insolation. Unfortunately, he came up against vociferous opposition from the scientific establishment at the time and it wasn’t until much later that Milankovitch picked up Croll’s baton and developed his ideas further.

    Milankovitch-Croll cycles anyone?

    • Hmmm….this paper argues otherwise:

      Click to access 3fleming_croll.pdf

      about the reception of Croll’s work.

      Note also that it also ran into opposition because one of its predictions did not fit the facts: According to Croll, only one hemisphere should show glaciation, not both (there also was a timing issue incorrectly predicted, IIRC).

    • IIRC, both Arrhenius and Ekholm cited Croll’s work–I know Ekholm did, and believe his friend Arrhenius did as well.

      For info on the rather eventful life and career of Nils Ekholm, see:

    • Steve,

      You beat me to it. Perhaps it should be Croll-Milankovich (Milanković)? Croll’s ‘Climate and Time, in Their Geological Relations’ (1875) explained how orbital variations could cause cycles of ice ages. Milanković improved on Croll but Croll should be credited with the idea as well IMO.

  12. And I think we can safely conclude that it’s not earth’s glacial changes that are causing the astronomical cycles.

    Er, don’t be too sure of that. Remember Jim Goodridge’s post at WUWT, where he accidentally proved that global warming causes sunspots?

  13. On a more serious note, this is nice post, and I agree that the comments are a pleasure to read. I’m very much looking forward to the next installment.

  14. Tamino, can you help me understand my error.

    With regard to orbital tuning, it seems logical to me that if I took a large number of short random sequences and stacked them to form a long sequence taking care that the largest peak in any sequence always matched a 41k year cycle, a fourier analysis would pick out the 41k year cycle very strongly. This would be because the large peaks all align on the 41k period, while the many other random variations would not align, and hence would cancel each other out. Even though the largest peak in any given series would not by much larger than other positive fluctuations, it would be significantly larger in the stack because it would not have been matched with any negative fluctuations. It also seems logical that the shorter the random sequences compared to the 41k period, the stronger this artifact would be.

    In believing this, I think I am assuming that power in the fourier analysis graph is a function both of the reliability of the signal (is the peak always present at that period) as well as the relative magnitude of the signal. Of course, I cannot perform fourier analyses, so my assumptions about “what is logical” may be entirely in error. Are there any specific errors in my analysis in the paragraph above, or in my interpretation of “power”in fourier analysis?

    If I am correct about that, then the more the 41k period is used to determine relative position of sequences within the chronology, or the duration of individual sequences; the weaker will be a strong 41k signal from fourier analysis as evidence that the 41k signal was not an artifact. (I am aware that my model of “tuning the stack” is a travesty of what was actually done by Lisiecky and Raymo, which is why I agreed that establishing an effective stack does support the idea that the 41k signal in their stack is a natural signal from considerations of coherence. I am trying to understand why my point that the orbital tuning weakens the evidence is false.)

    [Response: I think we’re in less disagreement than you suspect. Orbital tuning *does* weaken the evidence for the reality of an orbital cycle.

    But, although it will increase the power in a power spectrum, we obviously agree that it won’t increase the physical amplitude of a signal, or the variance of the data even if they’re nothing but random noise. And that places a limit on how much the tuning can inflate the periodic response, especially in its estimated physical amplitude. In a sense, there’s more “evidence” from an orbitally tuned stack in the amplitude spectrum from Fourier analysis, than there is in the power spectrum (only rarely does anyone plot an amplitude spectrum, although I’m a big fan of them).

    Of course the real “proof” is in the analysis of data which are *not* orbitally tuned — and that evidence is abundant.]

  15. Tamino, thanks for your responce to my comments in this post and in part 1b.