My name is Blake S. Pollard. I am a physics graduate student working under Professor Baez at the University of California, Riverside. I studied Applied Physics as an undergraduate at Columbia University. As an undergraduate my research was more on the environmental side; working as a researcher at the Water Center, a part of the Earth Institute at Columbia University, I developed methods using time-series satellite data to keep track of irrigated agriculture over northwestern India for the past decade. I am passionate about physics, but have the desire to apply my skills in more terrestrial settings. That is why I decided to come to UC Riverside and work with Professor Baez on some potentially more practical cross-disciplinary problems. Before starting work on my PhD I spent a year surfing in Hawaii, where I also worked in experimental particle physics at the University of Hawaii at Manoa. My current interests (besides passing my classes) lie in exploring potential applications of the analogy between information and entropy, as well as in understanding parralells between, statistical, stochastic, and quantum mechanics.

Glacial cycles are one essential feature of Earth's climate dynamics over timescales on the order of 100's of kiloyears (kyr). It is often accepted as common knowledge that these glacial cycles are in some way forced by variations in the Earth's orbit. In particular many have argued that the approximate 100 kyr period of glacial cycles corresponds to variations in the Earth's eccentricity. As we saw in Professor Baez's earlier posts, while the variation of eccentricity does affect the total insolation arriving to Earth, this variation is small. Thus many have proposed the existence of a nonlinear mechanism by which such small variations become amplified enough to drive the glacial cycles. Others have propsed that eccentricity is not primarily responsible for the 100 kyr period of the glacial cycles.  

Here is a brief summary of some time series analysis I performed in order to better understand the relationship between the Earth's Ice Ages and the Milankovich cycles.

The orbital data is due to Berger et. al. (1991).

The temperature proxy is based on changes in deuterium concentrations from the EPICA Antarctic ice core dating back over 800 kyr. Both are available via ftp from: ftp://ftp.ncdc.noaa.gov/pub/data

Jouzel, J., et al. 2007.
EPICA Dome C Ice Core 800KYr Deuterium Data and Temperature Estimates.
IGBP PAGES/World Data Center for Paleoclimatology
Data Contribution Series # 2007-091.
NOAA/NCDC Paleoclimatology Program, Boulder CO, USA.

The orbital data also includes an estimate of the solar insolation derived from the orbital parameters, which is plotted below against the temperature proxy from the ice core.

{1InsolTempPlot.png}

I'm going to focus on the data from the orbital parameters themselves, plots of which are shown below, obliquity (tilt), precession (direction tilted axis is pointing), and eccentricity (deviation from circular).

{2OBL} {3PREC}

{4ECC}

Muller and MacDonald have argued (paper, and their very relevant book) that another astronomical parameter, the orbital plane  inclination, the angle between the plane Earth's orbit and the 'invariant plane' of the solar system. This invariant plane of the solar system depends on the angular momenta of the planets, but roughly coincides with the plane of Jupiter's orbit from what I understand. Here is a plot of the orbital plane inclination for the past 800 kyr, the data is from Professor Muller's website.

{5Oplane}

One can see from these plots, or from some spectral analysis, that the main periodicities of the orbital parameters are: Obliquity ~ 42 kyr, Precession ~ 21 kyr, Eccentricity ~100 kyr orbital plane~ 100 kyr .

Of course the curves clearly are not simple sine waves with those frequencies. Fourier transforms give information regarding the relative power of different frequencies occurring in a time series, but there is no information left regarding the time dependence of the principle frequencies as the time dependence is integrated out in the Fourier Transform.

The Gabor Transform is a generalization of the Fourier Transform, sometimes referred to as the 'windowed' Fourier Transform. For the Fourier Transform, $latex F(w) = \dfrac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(t) e^{-iwt}$, one may think of $latex e^{-iwt} $, the 'kernel function', as the guy acting as your basis element in both spaces. For the Gabor Transform instead of $latex e^{-iwt}$ one defines a family of functions, $latex g_{(b,\omega)}(t) = e^{i\omega(t-b)}g(t-b) $ where $latex g(t-b)\epsilon L^{2}(\mathbb{R})$, is called the window function. Typical windows are square windows, triangular (Bartlett), and the most common is the Gaussian, $latex g(t)= e^{-kt^2} $, which is used in the analysis below. The Gabor Transform of a function f(t) is then given by $latex G_{f}(b,w) = \int f(t) \overline{g(t-b)} e^{-iw(t-b)}dt$Note the output of a Gabor Transform, like the Fourier Transform, is a complex number. The modulus of this number indicates the strength of a particular frequency in the signal, while the phase carries information about the... well, phase.

For example the modulus of the Gabor transform of $latex f(t)=\sin(\dfrac{2\pi t}{100})$ is shown below. For these I used the package Rwave, orginally written in S by Rene Carmona and Bruno Torresani; R port by Brandon Whitcher.

{6CGTsint}

You can see that the line centered at a frequency of .01 corresponds to the function's period of 100 time units.

A Fourier Transform would do okay for such a function, but consider now a sine wave whose frequency increases linearly. As you can see below the Gabor Transform of such a function indicates the linear increase of frequency with time.

{7CGTsintsqr}

The window parameter in both of the above Gabor Transforms is 100 time units. Adjusting this parameter effects the vertical blurriness of the Gabor Transform. For example here is the same plot as a above, but with window parameters of 300, 200, 100, and 50 time units:

{8window300} {9window200}
{10window100} {11window50}

You can see a you make the window decreases the line gets sharper, but only to a point. When the window becomes approximately smaller than a given period of the signal the line starts to blur again.

Now I'll say a few things to warm you up to the temperature data (plot below). First you notice spikes occurring about every 100 kyr. You can also see that the last 5 of these spikes appear to be bigger/more dramatic than the ones occurring before 500 kyr ago.  Roughly speaking each of these spikes corresponds to rapid warming of the Earth, after which occurs slightly less rapid cooling, and then a slow decrease in temperature until the next spike occurs. These are the Earth's glacial cycles. At the bottom of the curve where the temperature is about $latex 4^{o} C$ cooler than the mean of this curve glaciers are forming and extending down across the northern hemisphere. The relatively warm periods on the top of the spikes, about $latex 10^{o} C$ hotter than the glacial periods are called the interglacials. You can see that we are currently in the middle of an interglacial, so the Earth is relatively warm compared to rest of the glacial cycles.

The following plot shows Earth's temperature, calculated from the EPICA ice core deuterium concentrations, over the past 800 kyr.

{12TEMP}

Now we'll take a look at the spectrum of the windowed Fourier Transform, or the Gabor Transform, of this data. The window 'size' for these plots is 300 kyr.

{13GTMtemp}

One can see a few interesting features in this plot, zooming in a bit.

{14GTMtempzoom}

We see one line at a frequency of about .024, with a sampling rate of 1 kyr, corresponds to a period of about 42 kyr, close to the period of obliquity. We also see a few things going on around a frequency of .01, corresponding to a 100 kyr period.

The band at .024 appears to be relatively horizontal, indicating an approximately constant frequency. Around the 100 kyr periods there is more going on. At a slightly higher frequency ~ .015, there appears to be a band of slowly increasing frequency. Also around the .01 it's hard to say what is really going on. It is possible that we see a combination of two frequency elements, one increasing, one decreasing, but almost symmetric. This may just be an artifact of the Gabor Transform or the window and frequency parameters.

The window size for the plots below is slightly smaller, about 250 kyr.
If we put the temperature and obliquity Gabor Transforms side by side:

{15cgtOBLtemp}

It's clear the lines at .024 line up pretty well.

Doing the same with eccentricity:

{16cgtECCtemp}

Eccentricity does not line up well with temperature in this exercise though both have bright bands above and below .01 . 

Now for temperature and orbital inclination:

{17cgtOPlanetemp}

One sees that the frequencies line up better for this than for eccentricity, but one has to keep in mind that there is a nonlinear transformation performed on the 'raw' orbital plane data to project this down into the 'invariant plane' of the solar system. While this is physically motivated, it surely nudges the spectrum.

The temperature data clearly has a component with a period of approximately 42 kyrs, matching well with obliquity. If you tilt your head a bit you can also see an indication of a fainter response at a frequency a bit above .04, corresponding roughly to period just below 25 kyrs, close to that of precession.

As far as the 100 kyr period goes, which is the periodicity of the glacial cycles, this analysis confirms much of what is known, namely that we can't say for sure. Eccentricity seems to line up well with a periodicity of approximately 100 kyrs, but at closer inspection there seems to be some discrepancies if you try to understand the glacial cycles as being forced by variations in eccentricity. The orbital plane inclination has a more similar Gabor transform modulus than does eccentricity. A good next step would be to look the relative phases of the orbital parameters versus the temperature, but that's all for now.