Exploring EPICA Dome C records with Pyleoclim

Authors¶

Julien Emile-Geay and Deborah Khider

Preamble¶

Goals¶

In this demo, we use Pyleoclim to explore the mystery of Ice Ages as documented from the EPICA Dome C ice core.

Technical skills involved:

• loading data using pandas

• applying and interpreting spectral analysis in Pyleoclim

• applying and interpreting singular spectrum analysis in Pyleoclim

• applying and interpreting Wavelet Transform Coherency analysis

Data¶

The Deuterium data $\delta D$ was published by Jouzel et al. (2007). It is a proxy for the temperature of snow formation above the site, taken to be representative of Antarctic temperature over such timescales. The demo also makes use of $CO_2$ measurements came from Luthi et al (2008), corrected by Bereiter et al. (2015).

For a video version, check out this link, published in Dec 2021 under version 0.7.4.

Reading time¶

15min

Keywords¶

Pyleoclim, EPICA Dome C, Spectral Analysis, Singular Spectrum Analysis, Wavelet Transform Coherency

Let us first load necessary packages:

Loading and visualizing $\delta D$ series¶

Let’s load the data into a pandas DataFrame:

Let’s have a look at the DataFrame and see if there are any missing values:

There seems to be missing values in the Deuterium data. This is not an issue to work with Pyleoclim since the first step is to create an object called Series (more info here), which will remove any missing values and sort the values in ascending order with default parameters:

We can plot this with a single line of code, which greatly simplifies life:

Notice how we put in quite a bit of metadata (information about the data) when we created the Series object, and that this information is now being used by the code to properly label axes for you. Structuring and labeling your data is a little more work upfront but it pays big dividends down the line. Let’s see that now.

First, take a look at the distribution of time increments:

The data are not evenly spaced, which is a challenge for spectral analysis because most timeseries analysis methods (e.g. Fourier analysis) make the implicit assumption that data are evenly spaced, so applying those tools to unevenly-spaced data will result in methods not behaving as they are supposed to.

Spectral Analysis¶

Thankfully, Pyleoclim can deal with that challenge (unevenly-spaced data) very easily. Two of the spectal methods in the package can handle an uneven time axis:

• Lomb-Scargle periodogram

• Weighted-Wavelet Z-transform

The former is fast, but far from optimal when interested in looking at the slope of the spectrum. WWZ is usually better, but can be very slow for large datasets.

Lomb-Scargle Periodogram¶

We start with the Lomb-Scargle periodogram, which may be accessed via the spectral method associated with any Series object:

Notice how we chained together two commands to do this: we first standardized the data (not strictly necessary, but useful later) and computed the spectrum. The output (called a periodogram) is stored in a variable (we’ll see this is a good idea why in a minute). Then we plotted it. Notice also how the plot contains all the labels you would want, and the frequency axis is labeled in terms of periods, which are much more intuitive than their inverse (periods).

We could also have chained all three commands to directly create the plot. This can be useful for data exploration; but in this case, the periodogram is not saved as a variable and is therefore not available for later use.

Which of these peaks are significant? To do so, the method signif_test can be applied to the PSD object psd_ls, the structure that stored the result from the LS method:

We see that the orbital periods poke above the 95% threshold, but not by much. Certainly, none of the high-frequency wiggles are. It would make more sense to interpolate the data at 0.5ky resolution and redo this. In fact, you can do this all in a single line of code:

The result is still a bit noisy, but at least we’ve cut all that high-frequency junk, and because the timeseries are shorter, the computations are much faster. We can see peaks poking above the AR(1) significance level at the classical orbital periodicities: 21 ky (precession), 41 ky (obliquity) and 100 ky (eccentricity), with a number of higher frequency peaks breaching that limit as well.

Now, let’s see what the Weighted Wavelet Z-transform (WWZ) makes of the same data:

Notice how we saved the output of the calculation into a structure called psd_wwz before passing that to the plot function; this is to save time, as these calculation can get lengthy. Indeed:

The same orbital periods around 20, 40 and 100 kyr pop out, with a few high-frequency bumps grazing the 95% threshold. This may disappear if we ran more AR(1) benchmarks (say, 1000), but running 50 already took some time. You’d certainly want to invest that time for paper/talk, but in this demo we have other fish to fry.

Let us note in passing that we could plot both of these spectra on the same plot, and add Multitaper Method (MTM) for good measure:

We see that MTM’s bias variance tradeoff (with the default time-bandwith product used here) results in very attenuated, broad peaks at orbital periods. WWZ (again with default parameters) results in sharper peaks, with a secondary one around 80ky. Lomb-Scargle is in between those two extremes, and markedly noisier at high frequencies (no matter what is done to its parameters - try it)

Singular Spectrum Analysis¶

If we really wanted to focus on orbital periods, we might want to clean the record still further. One way to do this is to apply Singular Spectrum Analysis (SSA) to a standardized, interpolated version of the series (SSA will not work on unevenly spaced data).

The output SSARes object, dDssa, contains a bunch of usable information, with some associated methods. The first thing that one does with a decomposition like this is to look at the eigenvalue spectrum via a “scree plot”:

Clearly, the first 8 eigenvalues are a cut above the rest. We should be able to reconstruct most behavior of interest using just those 8. Indeed, if we sum their variance we get to over 99%:

Let’s extract those first 8 modes using RCmat, the matrix containing the reconstructed SSA modes, and plot it alongside the original, interpolated data:

Let’s compute its WWZ spectrum and plot it along the old one:

We see that SSA was able to capture the most essential swings of the original series, while filtering out the high frequency noise. In this respect, it acted very much like a lowpass filter, which we could also have applied like so:

Wavelet Analysis¶

Lastly, we play with wavelet analysis, which may be used to “unfold” a spectrum and look at its evolution over time.

We can now use the summary_plot functionality, which stacks together the timeseries itself, its scalogram, and the power spectral density (PSD) obtained from summing the wavelet coefficients over time. By default, the function will calculate the periodogram and scalogram for you. But since we have already done so, we can pass the results of these analyses directly to the plotting function:

The scalogram reveals how spectral power (technically, wavelet power) changes over time. Integrating this wavelet density over time yields something close to the spectral density (PSD), which this summary plot displays to the right of the scalogram. Integrating over scale recovers the original signal, which is plotted on top. However, note that the PSD and scalogram were obtained from different methods here, so they are not exactly comparable as indicated by the labels. For figures intended for publication, make sure that your labels reflect the methods use for transparency. If space of the essence, do so in the caption.

There is another difference between this PSD and scalogram: the PSD came with an estimate of significance ; the scalogram did not. However, this is easy to add in; we just apply the signif_test() method to the scalogram, and pass its result to summary_plot():

The white lines delineate regions of the scalogram that are significant against an AR(1) benchmark, so encircle “islands” of notable power. We see that the 100kyr periodicity is particularly pronounced around 300-600 kyr BP, while the 40 and 20kyr periodicities are more pronounced in the later portions (since 400 ky BP). This may be because of the compaction of the ice at depth, which you wouldn’t know from analyzing just this dataset. Paleoclimate timeseries must always be interpreted with those possible artifacts in mind.

Now, notice that we specified the method ‘ar1asym’ for the asssessment of significance; this tells the code to use an asymptotic approximation to the distribution of the AR(1) benchmarks, making it lightning-fast. We could also have used direct simulation (the default method, because it applies to other wavelet methods like wwz), and indeed this would yield a very similar result:

Assessing Memory¶

One may also assess the scaling properties of a timeseries with Pyleoclim. Scaling exponents are a very useful way to characterize timeseries, and are related to the concept of “memory”. One method of doing so, is to fit a regression line to the PSD and estimate its slope, indicative of fundamental scaling laws in the continuum of climate variability. Pyleoclim makes this easy:

Note that this slope is related to the Hurst exponent: $H = \frac{\beta-1}{2}$

This value is close to 0.5, the expected value for a Brownian motion type of stochastic process.

Temperature and $CO_2$¶

Now let us load the $CO_2$ composite from this and other neighboring sites around Antarctica:

We see very similar Ice Ages as in the deuterium data and of course a precipitous rise since the Industrial Revolution. To plot the two series side by side, we can put them into a MultipleSeries object and use the default plot function:

By default, the MultipleSeries class assumes commensurate units, which is not really the case here. Fear not, we can just standardize the series:

We could also plot them as a stack, like so:

Now, one might be interested in lead/lag relationships between those two series. Before that, a brief primer: the temperature proxy $\delta D$ is measured on the ice, whereas $CO_2$ is measured in the air trapped in the ice. Because bubbles close only progressively as the firn gets compacted, the air can be younger than the surrounding ice by several hundred years. The ice core community has worked diligently on this for decades and have made very clever adjustments to correct for this effect, but that is to be kept in mind when comparing those two data streams.

With that in mind, let us apply wavelet transform coherency to identify phase relationships between the two series at various scales. The wavelet_coherence method is attached to the Series object in Pyleoclim and takes as input another Series object:

Arrows pointing to the right indicate a phase angle close to zero at all scales longer than precession. We can quantify that by computing the average phase angle in a given range of scales, like so:

This means that on orbital timescales, the two series are essentially in phase ; there is no consistent lead or lag between the two series. This is remarkable given the dating challenges mentioned earlier, and is widely interpreted to mean that on such timescales, atmospheric $CO_2$ is a part of a positive feedback loop amplifying orbitally-triggered changes in temperature. However, the variance of this angle is fairly high, and by this test it does not register as a very consistent signal. Lastly, note that in the modern era, atmospheric $CO_2$ is of course a forcing (the largest climate forcing, currently), acting on much, much faster timescales for which the climate system is still catching up. You would not know it from this analysis, but it’s important to state out loud, given that climate deniers have an annoying tendency to cherry-pick the paleoclimate record in support of their feeble arguments.

One might be tempted to interpret phase angles at shorter scales. Let us first see if they are significant against AR(1) benchmarks, As before this can be done in one function call, though some patience is required to obtain the result.

These phase differences are indeed significant according to this test, but mostly inconsistent (at a given scale, the arrows veer one way or the other). This would not be a good candidates for identifying a robust lead or lag. One would find the same result with phase_stats().

Another consideration is that coherency is like the correlation coefficient in wavelet space: it tells you how similar two series are for a given time and scale, yet it says nothing about what fraction of the total variability is shared. This is better measured by the cross-wavelet transform (XWT), which highlights areas of high common power. Both of those, along with the original series, can be visualized with one swift function call:

Here we see that the orbital bands are the only ones to show up consistently throughout the record, but precession and obliquity fade before 200 ky BP in their common power, and the XWT band near 120 ka drifts to 80ka back in time. This means that, out of the large swathes of WTC one would consider “coherent”, only those areas highlighted in XWT in green-yellow are likely important. That’s where we stop the banter, though: time to go ask a glaciologist about what is going on in those datasets. I know a few good ones.

This concludes our mini tour of paleoclimate timeseries wizardry with Pyleoclim. Be sure to check our documentation for more, and remember that magic tricks are not science until backed up by knowledge.

Finally, please watch the PaleoBooks Gallery for more examples. Feel free to get in touch on Discourse if you want to suggest (or make!) other demos/tutorials.