Chasing cyclicities

In this notebook, we demonstrate the use of the spectral and wavelet analysis features of Pyleoclim using a sea surface temperature reconstruction from Site ODP846 in the Eastern Equatorial Pacific. The goal is to explore the periodicities present in this 5 million years record and explore any coherency between insolation and the sea surface temperature record.

Data

The record is described in:

• Mix, A. C., J. Le, and N. J. Shackleton (1995), Benthic foraminiferal stable isotope stratigraphy from Site 846: 0–1.8 Ma, Proc. Ocean Drill. Program Sci. Results, 138, 839–847.

• Shackleton, N. J., Hall, M. A., & Pate, D. (1995). Pliocene stable isotope stratigraphy of ODP Site 846. Proc. Ocean Drill. Program Sci. Results, 138, 337-356.

• Lawrence, K. T., Liu, Z. H., & Herbert, T. D. (2006). Evolution of the eastern tropical Pacific through Plio-Pleistocne glaciation. Science, 312(5770), 79-83.

The data were aligned to the Benthic Stack of Lisiecki & Raymo (2005) (LR04) using the HMM-Match algorithm developed by Lin et al. (2014) as used by Khider et al, 2017. The latter is a probabilistic method that generates an ensemble of 1000 possible age models compatible with the chronostratigraphic constraints.

The dataset is stored in the Linked Paleo Data Format (LiPD, McKay and Emile-Geay (2016)) and available from the same GitHub repository as the notebook.

Data Exploration

Let’s import this file into Pyleoclim and make a informational dashboard about the dataset. The cell below will import the various packages needed throughout the notebook, including Pyleoclim.

import pyleoclim as pyleo
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import numpy as np
from pylipd import LiPD

# For insolation curves, let's use climlab
from climlab.solar.orbital import OrbitalTable
from climlab.solar.insolation import daily_insolation

pyleo.set_style('web') #set the figure style

import warnings
warnings.filterwarnings('ignore')

We can use the pyLiPD package to open and manipulate datasets published in the LiPD format.

D = LiPD()
D.load('./data/ODP846.Lawrence.2006.lpd')

Loading 1 LiPD files

100%|█████████████████████████████████████████████| 1/1 [00:00<00:00,  1.84it/s]

Loaded..

The next step is to isolate the time series of interest (in this case, the temperature record) and create a GeoSeries object in Pyleoclim. GeoSeries is itself a child of the Series object, which carries most of the functionalities of Pyleoclim. This means that all the methods applicable to Series are also applicable to GeoSeries (with some additional ones enabled by the additional metadata extracted from a LiPD file).

To extract individual timeseries using PyLiPD, you can use the method get_timeseries_essential, which will extract the time series information and associated metadata into a Pandas DataFrame. Each row will correspond to a timeseries.

df = D.get_timeseries_essentials()
df.head()

	dataSetName	archiveType	geo_meanLat	geo_meanLon	geo_meanElev	paleoData_variableName	paleoData_values	paleoData_units	paleoData_proxy	paleoData_proxyGeneral	time_variableName	time_values	time_units	depth_variableName	depth_values	depth_units
0	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	ukprime37	[0.821, 0.824, 0.828, 0.787, 0.777, 0.767, 0.7...	unitless	None	None	age	[5.228, 8.947, 11.966, 14.427, 16.502, 18.41, ...	kyr BP	depth	[0.16, 0.26, 0.36, 0.46, 0.56, 0.66, 0.76, 0.8...	m
1	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	c37 total	[2.37, 2.1, 1.87, 2.74, 3.75, 7.62, 7.86, 7.73...	nmol/kg	None	None	age	[5.228, 8.947, 11.966, 14.427, 16.502, 18.41, ...	kyr BP	depth	[0.16, 0.26, 0.36, 0.46, 0.56, 0.66, 0.76, 0.8...	m
2	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	u. peregrina d13c	[nan, nan, nan, nan, nan, nan, nan, nan, nan, ...	per mil PDB	None	None	None	None	None	depth cr	[0.12, 0.23, 0.33, 0.33, 0.43, 0.53, 0.63, 0.7...	rmcd
3	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	u. peregrina d13c	[nan, nan, nan, nan, nan, nan, nan, nan, nan, ...	per mil PDB	None	None	None	None	None	depth comp	[12.0, 23.0, 33.0, 33.0, 43.0, 53.0, 63.0, 73....	mcd
4	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	u. peregrina d13c	[nan, nan, nan, nan, nan, nan, nan, nan, nan, ...	per mil PDB	None	None	None	None	None	depth	[0.12, 0.23, 0.33, 0.33, 0.43, 0.53, 0.63, 0.7...	m

Let’s have a look at the variables present in the file:

df['paleoData_variableName'].unique()

array(['ukprime37', 'c37 total', 'u. peregrina d13c', 'temp prahl',
       'event', 'c. wuellerstorfi d13c', 'section',
       'c. wuellerstorfi d18o', 'site/hole', 'sample label',
       'temp muller', 'u. peregrina d18o', 'interval'], dtype=object)

Let’s pick the sea surface temperature record obtained from the Muller et al. (1998) calibration and create a GeoSeries object:

row = df[df['paleoData_variableName']=='temp muller'].reset_index() #select the proper row and re-index it to zero so we can select it easily

ts = pyleo.GeoSeries(time=row['time_values'][0], value = row['paleoData_values'][0], lat = row['geo_meanLat'][0],
                     lon = row['geo_meanLon'][0], elevation = row['geo_meanElev'][0], time_name= 'Age', time_unit = row['time_units'][0],
                     value_name = 'Sea Surface Temperature', value_unit = '$^\circ$C', label = row['dataSetName'][0], 
                     archiveType = row['archiveType'][0], control_archiveType=True, depth = row['depth_values'][0], 
                     depth_name = row['depth_variableName'][0], depth_unit = row['depth_units'],verbose = False)

Let’s have a look at the parameters passed to the GeoSeries object upon its creation:

• time: The values for the time axis. (Mandatory field)

• value: The values for the variable of interest, in this case sea surface temperature. (Mandatory field)

• lat: The latitude the archive was recovered from. (Mandatory field)

• lon: The longitude the archive was recovered from. (Mandatory field)

• elevation: The depth (or altitude) the archive was recovered from.

• time_name: THe name of the time axis. By convention, we use ‘Age’ for retrograde time and ‘Year’ for prograde.

• time_unit: The units associated with the time axis

• value_name: The name of the variable of interest

• value_unit: The units for the variable of interest

• label: The label for the timeseries. This is particularly useful if you plan to plot multiple timeseries to generate the legend automatically.

• archiveType: The type of geological archive the measurements can be made on.

• control_archiveType: If set to True, Pyleoclim will attempt to match your string to a standardized vocabulary of archiveType. In this notebook, we set it to True to enbale default plottinf style as the LiPD file was created before standardization. Note that this step is not necessary and is mostly useful for mapping.

• depth: the values for the depth in the archive where the measurements were made

• depth_name: The name for the depth in the archive

• depth_unit: The units of depth

• verbose: Here, we set it to False to prevent Pyleoclim from giving out too many warnings. If this is you first time using the package or if you are not familiar with all the parameters associated with the GeoSeries object, we recommend that you turn verbose to True.

Let’s create a dashboard, a functionality specific to GeoSeries, which will give us a plot of the time series, the location of the record and a spectrum along with some specific metadata.

ts.dashboard(plt_kwargs={'marker':None})

Performing spectral analysis on individual series: 100%|█| 200/200 [00:52<00:00,

(<Figure size 1100x800 with 4 Axes>,
 {'ts': <Axes: xlabel='Age [kyr BP]', ylabel='Sea Surface Temperature [$^\\circ$C]'>,
  'dts': <Axes: xlabel='Counts'>,
  'map': {'map': <GeoAxes: xlabel='lon', ylabel='lat'>},
  'spec': <Axes: xlabel='Period [kyr]', ylabel='PSD'>})

Sea surface temperatures have decreased throughout the last 5M years, with cycles superimposed on the long-term cooling trend. The periodogram indicates significant periodicities at orbital time scales, which merits further investigation.

Since the age model was obtained by alignment of the benthic \(\delta^{18}O\) to the LR04 curve, let’s also load it into a GeoSeries object for comparison. The data is stored in the chron tables, which we indicate by using the mode parameter:

df_c = D.get_timeseries_essentials(mode='chron')
df_c

	dataSetName	archiveType	geo_meanLat	geo_meanLon	geo_meanElev	chronData_variableName	chronData_values	chronData_units	time_variableName	time_values	time_units	depth_variableName	depth_values	depth_units
0	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	age	[3.645, 7.99, 11.18, 13.803, 15.886, 17.93, 19...	ky BP	age	[3.645, 7.99, 11.18, 13.803, 15.886, 17.93, 19...	ky BP	depth	[0.12, 0.23, 0.33, 0.43, 0.53, 0.63, 0.73, 0.8...	m
1	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	depth	[0.12, 0.23, 0.33, 0.43, 0.53, 0.63, 0.73, 0.8...	m	age	[3.645, 7.99, 11.18, 13.803, 15.886, 17.93, 19...	ky BP	depth	[0.12, 0.23, 0.33, 0.43, 0.53, 0.63, 0.73, 0.8...	m
2	ODP846.Lawrence.2006	Marine sediment	-3.1	-90.8	-3296.0	d18o	[3.38, 3.46, 3.765, 4.14, 4.47, 4.99, 4.99, 4....	permil	age	[3.645, 7.99, 11.18, 13.803, 15.886, 17.93, 19...	ky BP	depth	[0.12, 0.23, 0.33, 0.43, 0.53, 0.63, 0.73, 0.8...	m

Let’s grab the benthic \(\delta^{18}O\) data and store it into another GeoSeries:

row = df_c[df_c['chronData_variableName']=='d18o'].reset_index() #select the proper row and re-index it to zero so we can select it easily

ts_benthic = pyleo.GeoSeries(time=row['time_values'][0], value = row['chronData_values'][0], lat = row['geo_meanLat'][0],
                     lon = row['geo_meanLon'][0], elevation = row['geo_meanElev'][0], time_name= 'Age', time_unit = row['time_units'][0],
                     value_name = 'Benthic $\delta^{18}O$', value_unit = u'\u2030', label = row['dataSetName'][0], 
                     archiveType = row['archiveType'][0], control_archiveType=True, depth = row['depth_values'][0], 
                     depth_name = row['depth_variableName'][0], depth_unit = row['depth_units'],verbose = False)

And do a simple plot of the benthic record:

fig,ax = ts_benthic.plot(legend=False)
ax.invert_yaxis() #flip the axis by convention

The benthic \(\delta^{18}O\) record displays:

• a long-term cooling trend (\(\delta^{18}O\) gets more positive over time) characteristic of late Neogene and Quaternary and increasing ice volume

• shorter term shifts that represent the ice ages and that are used for the alignment of this benthic record with the benthic stack (LR04).

Speaking of the LR04 record, let’s load the benthic stack into a Series for comparison. The data is saved in csv format that can be loaded into a Pandas DataFrame.

df = pd.read_csv('./data/LR04.csv',skiprows=4)
df.head()

	Time (ka)	Benthic d18O (per mil)	Standard error (per mil)
0	0.0	3.23	0.03
1	1.0	3.23	0.04
2	2.0	3.18	0.03
3	3.0	3.29	0.03
4	4.0	3.30	0.03

And let’s create a Series from this DataFrame:

lr04 = pyleo.Series(time=df.iloc[:,0],value=df.iloc[:,1],value_name='Benthic $\delta^{18}O$',value_unit=u'\u2030',
                   time_name='Age',time_unit='ky BP',label='LR04')

Time axis values sorted in ascending order

And plot it:

fig,ax = lr04.plot(legend=False) 
ax.invert_yaxis() #invert the y-axis as per standard usage

Let’s plot the series together on a stackplot. This functionality is enabled through Pyleoclim’s MultipleSeries object, which can be created from a list of Series (and by extension LipdSeries) objects.

ms = pyleo.MultipleSeries([ts,ts_benthic,lr04])
fig,ax = ms.stackplot(figsize=[10,8],labels=['ODP846-SST','ODP846-Benthic','LR04'])
ax[1].invert_yaxis()
ax[2].invert_yaxis()

Spectral and Wavelet Analysis

Benthic record

Let’s quickly examine the spectral content of the benthic record. Because the age model was obtained by aligning the record to the LR04 stack, which is orbitally tuned, we would expect to recover the orbital periodicities in the record.

Pyleoclim is capable of running both spectral and wavelet analysis, which give slightly different information about the timeseries. Both allow to examine the frequency content of a time series but wavelet analysis offers the added opportunity to examine the sationarity of that signal over time. Pyleoclim provides easy access to several methods, allowing you to explore the trade-offs of each method as well as the robustness of the signals to the statistical method.

For spectral analysis, Pyleoclim offers a variety of methods (e.g., MTM, Welch, Lomb-Scargle…). Some are more appropriate to evenly-spaced datasets (e.g., Lomb-Scargle) than others. Let’s first plot the time difference between adjacent samples for the record.

Data exploration

dt=np.diff(ts_benthic.time) # calculate the difference between adjacent time
fig, ax = plt.subplots(figsize=[10,4])
sns.histplot(dt,kde=True,ax=ax) 
ax.set_xlabel('Age increments (ka)')

Text(0.5, 0, 'Age increments (ka)')

The data is unevenly-spaced and therefore, traditional methods such as MTM, may not be appropriate for this dataset. Let’s start with a Lomb-Scargle periodogram (Lomb (1976), Scargle (1982,1989)).

We should also decide on any pre-processing steps. In this case, the long-term trend is confounding and should be removed prior to analysis. Pyleoclim offers several detrending methods:

• linear

• constant (the mean is removed from the series)

• Filtering through a Savitzky-Golay filter and removal of the low-frequency component.

• Empirical mode decomposition (removal of the last mode, which corresponds to the lowest-frequency signal).

Of these, the last two seem to be the most appropriate for the data since the long-term cooling is not linear. Let’s run both and compare the results.

ts_benthic_emd = ts_benthic.detrend().standardize() #detrend using the EMD method, which is the default in Pyleoclim, and standardize
ts_benthic_emd.label = 'Detrended, EMD method' #label for plot
ts_benthic_sv = ts_benthic.detrend(method = 'savitzky-golay').standardize() #detrend using the Savitzky-Golay method and standardize
ts_benthic_sv.label = 'Detrended, Savitzky-Golay method' #label for plot
fig,ax = ts_benthic_emd.plot() #plot. 
ts_benthic_sv.plot(ax=ax)
ax.invert_yaxis()

The two detrending methods return very similar results, so we proceed with the default EMD method.

Analysis

The Lomb-Scargle periodogram can then be obtained with the following line of code. We use method chaining to create a workflow that first detrend the time series using the default EMD method, standardize, run spectral analysis (note that here we ask for the Lomb Scargle method using a frequency method called Lomg Scargle after the REDFIT alogrithm of Schulz and Mudelsee (2002) and 5 overlapping segments), and perform significance testing using an AR1 ensemble.

The Lomb-Scargle method implemented in Pyleoclim makes use of the Welch’s overlapped segment averaging (WOSA) technique (Welch, 1967, Schulz and Stattegger (1997)). The parameter n50 allows to set the number of segments with a 50% overlap. Feel free to change the settings and explore the effect of this parameter.

psd_benthic=ts_benthic.detrend().standardize().spectral(method='lomb_scargle', freq='lomb_scargle', settings={'n50':5,'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [04:44<00:0

psd_benthic.plot(xlim=[500,5]) #limit to 500k period

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [ky]', ylabel='PSD'>)

Unsurprisingly, the record displays periodicities around 40 kyr and 100 kyr. Why is this result not surprising? Remember that the age model was obtained by aligning the benthic record to the LR04 stack, which is itself tuned to those very Milankovitch periodicities. If anything, this proved that the Lomb-Scargle method can recover periodicities embedded in the record by construction.

Another way to estimate the spectral density is to use the Weighted Wavelet Z-Transform (WWZ) from Foster (1996) as implemented by Zhu et al. (2019) for Pyleoclim, which is also appropriate for unevenly-spaced data. Note that this method is wavelet-based; we, therefore, compute the power spectral density (PSD) from the scalogram.

This opens two possibilities: (1) Run the ‘wwz’ method through pyleoclim.Series.spectral and/or (2) Run the ‘wwz’ method for wavelet analysis and compute the PSD from the stored scalograms. Which one should be used depends on the goal at hand.

To see why, let us have a look at the parameters attached to the wwz method as called through the wavelet() and spectral() methods. One of the parameters is the decay constant c, which balances the time resolution and frequency resolution of the wavelet analysis. The smaller this constant is, the sharper the peaks. For spectral analysis, where the purpose is to find peaks, a smaller value for c is needed compared to wavelet analysis, where time-frequency localization is of the essence (Foster (1996), Wu et al., 2007). Therefore, if the purpose is to identify spectral peaks (periodicites) in the record, it may be useful to run the wwz method through spectral() (method 1). If the purpose is to obtain the spectral slope only, then method 2 is more computationally efficient.

Let’s compare what happens with the two methods. Let’s first compute the PSD from the spectral() mehod. This calculation may take a few hours to run.

psd_benthic_wwz = ts_benthic.detrend().standardize().spectral(method='wwz',settings={'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
Performing spectral analysis on individual series: 100%|█| 1000/1000 [58:45<00:0

psd_benthic_wwz.plot(xlim=[500,5]) #limit to 500k period

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [ky]', ylabel='PSD'>)

The PSD obtained from WWZ method agrees remarkably well with the PSD obtained using Lomb-Scargle. We will compare the two methods on the same plot after we obtain the PSD using the scalograms from wavelet(). To do so, you need to set the export_scal parameter to True.

Note that this computation takes some time and the resulting file is ~500 Mb

scal_benthic_wwz=ts_benthic.detrend().standardize().wavelet(method='wwz',settings={'standardize':False}).signif_test(number=1000, export_scal=True)

Performing wavelet analysis on individual series: 100%|█| 1000/1000 [53:52<00:00

scal_benthic_wwz.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'ODP846.Lawrence.2006 scalogram (WWZ) with 95% threshold'}, xlabel='Age [ky BP]', ylabel='Scale [kyrs]'>)

The scalogram reveals the non-stationary character of the periodicities seen in spectral analysis and offers additional information related to the well-known mid-Pleistocene transition from a “41 kyr world” to a “100 kyr world” somewhere around 0.8 Ma (Paillard (2001)). There is also a drop in power after 3 Ma, which could be a climate signal or due to core resolution/alignment. The LR04 stack is noisier since fewer records are part of the composite over that time period. In other words, the noise to signal ratio may prevent any climate-related cyclicities from being identified.

Next, let’s compute the PSD from the scalogram and plot it.

psd_benthic_wwz_fromscal = ts_benthic.detrend().standardize().spectral(method='wwz',scalogram=scal_benthic_wwz).signif_test(qs=[0.90,0.95,0.99], scalogram=scal_benthic_wwz)
psd_benthic_wwz_fromscal.plot()

Performing spectral analysis on individual series: 100%|█| 1000/1000 [00:05<00:0

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [ky]', ylabel='PSD'>)

Note the effect of c on the periodogram: the lower resolution creates a smoother PSD, with peaks less sharply defined.

Let’s compare the Lomb-Scargle and WWZ method. To this end, we recompute the Lomb-Scargle periodogram using the default frequency vector:

psd_benthic_default_freq=ts_benthic.detrend().standardize().spectral(method='lomb_scargle', settings={'n50':5,'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [04:39<00:0

psd_benthic_default_freq.plot(xlim=[500,5]) #limit to 500k period

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [ky]', ylabel='PSD'>)

We can plot the two spectra obtained from the wwz and Lomb-Scargle methods on the same figure:

fig,ax = psd_benthic_default_freq.plot(xlim=[500,5], signif_clr='slateblue') #limit to 500k period and change the color of the significance line. 
psd_benthic_wwz.plot(ax=ax, signif_clr='sandybrown',xlim=[500,5]) #change the color of the significance line. 
current_handles, current_labels = plt.gca().get_legend_handles_labels() #get the information about the legend
current_labels[0] = 'Lomb-Scargle' # change label of the first entry to describe the method rather than the record name
current_labels[-4] = 'WWZ' # change label of the second to last entry to describe the method rather than the record name
plt.legend(current_handles,current_labels) ##replace legend entries
ax.set_xlim([500,5]) # set the current limits  to the 500k periodicity

(500.0, 5.0)

This plot highlights the remakable agreement between the two methods.

Pyleoclim can also create a summary plot of the wavelet and spectral analysis:

fig,ax = ts_benthic.summary_plot(psd=psd_benthic_wwz,scalogram=scal_benthic_wwz, value_label = '$\delta^{18}O$',
                                title = 'ODP846 Benthic $\delta^{18}O$',
                                ts_plot_kwargs = {'legend':False}) #change the label so it fits on the plot
ax['ts'].invert_yaxis() #invert the y-axis for the timeseries plot (print ax to figure out the 'ts' label)

Comparison with Benthic LR04 stack

The loss of power before 3 million years is intriguing and merits further investigation. There are two possible (and non mutually exclusive) explanations: (1) The orbital signal has a weaker impact on climate, (2) the noise to signal ratio is preventing detection. The LR04 stack is composed of 57 globally distributed records. However, most of the high-resolution records cover the last 2 million years (Lisiecki and Raymo, 2005), making the alignment between 3-5 Ma more difficult.

Let’s compute the scalogram for LR04:

scal_LR04_wwz=lr04.detrend().standardize().wavelet(method='wwz',settings={'standardize':False}).signif_test(number=1000)

Performing wavelet analysis on individual series: 100%|█| 1000/1000 [1:00:34<00:

scal_LR04_wwz.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'LR04 scalogram (WWZ) with 95% threshold'}, xlabel='Age [ky BP]', ylabel='Scale [kyrs]'>)

The stack displays strong periodicitied in the 100 kyr and 40 kyr band up to about 3 Ma. Despite being oribitally tuned, the method cannot accurately capture the 40 kyr periodicities from the early part of the record. We would, therefore, not expect to see them in the ODP846 record either.

Planktonic record

Data exploration

Let’s repeat our process with the planktonic data. In this case, the variations in SST were not aligned to the LR04 stack and are therefore not orbitally tuned. However, the data is still unevenly-spaced:

dt=np.diff(ts.time)
fig, ax = plt.subplots(figsize=[10,4])
sns.histplot(dt,kde=True,ax=ax)
ax.set_xlabel('Age increments (ka)')

Text(0.5, 0, 'Age increments (ka)')

Again, let’s look at detrending methods for pre-processing:

ts_emd = ts.detrend().standardize()
ts_emd.label = 'Detrended through EMD method'
ts_sv = ts.detrend(method = 'savitzky-golay').standardize()
ts_sv.label = 'Detrended through Savitzky-Golay method'
fig,ax = ts_emd.plot()
ts_sv.plot(ax=ax)

<Axes: xlabel='Age [kyr BP]', ylabel='Sea Surface Temperature [$^\\circ$C]'>

Once again, the two are equivalent and we will proceed with the default EMD method.

Analysis

Let’s obtain the Lomb-Scargle periodogram:

psd=ts.detrend().standardize().spectral(method='lomb_scargle', freq='lomb_scargle', settings={'n50':5,'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [05:00<00:0

psd.plot(xlim=[500,5])

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [kyr]', ylabel='PSD'>)

The record displays significant periodicities around and 100 kyr and singificant peaks in the precessional and obliquity band (23-40 kyr).

Let’s examine the robustness of these findings with the WWZ method. As we have done with the benthic record, let’s first compute the PSD from the spectral method.

psd_wwz = ts.detrend().standardize().spectral(method='wwz',settings={'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [1:17:54<00

psd_wwz.plot(xlim=[500,5]) #limit to 500k period

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [kyr]', ylabel='PSD'>)

And let’s compute the scalogram to look at the intermittence of the periodicities:

Note that this calculation takes some time to complete.

scal_wwz=ts.detrend().standardize().wavelet(method='wwz',settings={'standardize':False}).signif_test(number=1000, export_scal=True)

Performing wavelet analysis on individual series: 100%|█| 1000/1000 [1:07:58<00:

scal_wwz.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'ODP846.Lawrence.2006 scalogram (WWZ) with 95% threshold'}, xlabel='Age [kyr BP]', ylabel='Scale [kyrs]'>)

The scalogram displays significant and intermittent periodicities in the 100 kyr band between 0 and 3 Ma (beyond the mid-Pleistocene transmission). There is also significant power between 20-40 kyr throughout the last 3 Ma.

Let’s compute the PSD from the scalogram:

psd_wwz_fromscal = ts.detrend().standardize().spectral(method='wwz',scalogram=scal_wwz).signif_test(number=1000, scalogram=scal_wwz)
psd_wwz_fromscal.plot(xlim=[500,5])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [00:00<00:0

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [kyr]', ylabel='PSD'>)

Notice the effect of the c parameter. The PSD is thoroughly smoothed out when obtained from the scalograms directly.

Let’s recompute the PSD with the Lomb-Scargle method and the same frequency vector as we used for the wwz method for direct comparison:

psd_default_freq=ts.detrend().standardize().spectral(method='lomb_scargle', settings={'n50':5,'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [04:57<00:0

psd_default_freq.plot(xlim=[500,5])

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [kyr]', ylabel='PSD'>)

The PSD obtained with the WWZ method shows the significant periodicites at 100 kyr, 40 kyr and 23 kyr. Let’s plot the two PSDs on the same figure for comparison.

fig,ax = psd_default_freq.plot(xlim=[500,5], signif_clr='slateblue') #limit to 500k period and change the color of the significance line. 
psd_wwz.plot(ax=ax, signif_clr='sandybrown',xlim=[500,5]) #change the color of the significance line. 
current_handles, current_labels = plt.gca().get_legend_handles_labels() #get the information about the legend
current_labels[0] = 'Lomb-Scargle' # change label of the first entry to describe the method rather than the record name
current_labels[-4] = 'WWZ' # change label of the second to last entry to describe the method rather than the record name
plt.legend(current_handles,current_labels) ##replace legend entries
ax.set_xlim([500,5]) # set the current limits  to the 500k periodicity

(500.0, 5.0)

Again, the two methods agree remarkably well.

Effect of Age Uncertainty

Age model uncertainty can results in shifted age-depth relationship that may affect the spectrum. If age models are present in the LiPD file, Pyleoclim can help you identify spurious peaks. You can use PyLiPD to grab the age models from the ensembleTable.

ensemble_df = D.get_ensemble_tables()
ensemble_df.head()

	datasetName	ensembleTable	ensembleVariableName	ensembleVariableValues	ensembleVariableUnits	ensembleDepthName	ensembleDepthValues	ensembleDepthUnits	notes
0	ODP846.Lawrence.2006	http://linked.earth/lipd/chron0model0ensemble0	age	[[4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0,...	kyr BP	depth	[0.12, 0.23, 0.33, 0.43, 0.53, 0.63, 0.73, 0.8...	m	None

You can use the function from_AgeEnsembleArray to create an EnsembleGeoSeries object:

sst_ens = pyleo.EnsembleSeries.from_AgeEnsembleArray(series=ts, age_array= ensemble_df.iloc[0]['ensembleVariableValues'],
                                                     age_depth = ensemble_df.iloc[0]['ensembleDepthValues'], 
                                                     value_depth = ts.depth, verbose=False)

EnsembleGeoSeries is a child of EnsembleSeries, therefore the methods of this class apply to the child. We can plot using the following code:

fig,ax = sst_ens.plot_traces()
ax.get_legend().remove()

The envelope plot is often more practical to look at ensembles. However, make sure that the plot only encompasses the common period represented in every member. Don’t worry, Pyleoclim has a function for this called common_time.

sst_ens.common_time().plot_envelope()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [kyr BP]', ylabel='Sea Surface Temperature [$^\\circ$C]'>)

Let’s run spectral analysis. For speed, we will only use the Lomb-Scargle method. For consistency, we will also work with the same frequency vector as for the original analysis, which we can directly obtain from the object:

freq = psd.frequency #obtain the frequency vector from the PSD object created after spectral analysis. 
psd_ens=sst_ens.spectral(settings={'freq':freq,'standardize':False})

Performing spectral analysis on individual series: 100%|█| 1000/1000 [04:50<00:0

We can now plot the resulting ensemble PSD:

fig,ax = psd_ens.plot_envelope()
ax.set_xlim([150,10])
ax.set_ylim([10,10000])
ax.get_legend().remove()

Most of the low periodicity noise observed originally is smoothed out in the median ensemble, which is not surprising. The record retains orbital periodicity (40-100kyr) that is expressed in most of the age models.

Relationship to Insolation

Let’s examine the relation between the orbital periodicities observed in the SST record Given that the signal only exhibits singificant power in the 0-3MA, let’s slice the record

ts_slice = ts.slice([0,3000])
ts_slice.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [kyr BP]', ylabel='Sea Surface Temperature [$^\\circ$C]'>)

ts_slice.detrend().standardize().plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [kyr BP]', ylabel='Sea Surface Temperature [$^\\circ$C]'>)

Hypothesis 1: The periodic signals are driven by insolation at 5S, where the record is located.

Let us use ClimLab to obtain the insolation at 5S.

# array with specified kyears
years = np.linspace(-3000, 0, 3001)

# subset of orbital parameters for specified time
orb = OrbitalTable.interp(kyear=years)

# insolation values for past 3 Myears at 5S at summer solstice (day 172)
S5 = daily_insolation(lat=-5, day=172, orb=orb)

# put in a Pyleoclim.Series object
sol = pyleo.Series(time=-years,value=S5,value_name='Insolation at 5 $^\circ$S',value_unit='$W/m^{2}$',
                   time_name='Age',time_unit='ky BP',label='Insolation')

sol.plot()

Time axis values sorted in ascending order

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [ky BP]', ylabel='Insolation at 5 $^\\circ$S [$W/m^{2}$]'>)

Since the record is evenly-spaced, we can use spectral methods appropriate for evenly-spaced time series such as MTM for spectral analysis and the continous wavelet transform. For the later, we use the method by Torrence and Compo (1998):

psd_sol=sol.standardize().spectral(method='mtm',settings={'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [00:46<00:0

psd_sol.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [ky]', ylabel='PSD'>)

scal_sol = sol.standardize().wavelet(method='cwt',settings={'standardize':False}).signif_test(number=1000)

Performing wavelet analysis on individual series: 100%|█| 1000/1000 [00:48<00:00

scal_sol.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Insolation scalogram (CWT) with 95% threshold'}, xlabel='Age [ky BP]', ylabel='Scale [kyrs]'>)

Bother spectral and wavelet analysis shows high, significant power in the 23-kyr band, which corresponds to the precessional cycle. Let’s run coherence analysis between the SST record and insolation at 5S:

coh = ts_slice.detrend().standardize().wavelet_coherence(sol.standardize(), method='wwz',settings={'standardize':False}).signif_test(number=1000)

Performing wavelet coherence on surrogate pairs: 100%|█| 1000/1000 [2:46:24<00:0

fig, ax = coh.dashboard()
ax['ts1'].set_ylim([-4,8]) #offset the lines
ax['ts2'].set_ylim([-6,3]) #offset the lines

(-6.0, 3.0)

XWT (cross-wavelet transform) informs about areas where there is high common power between the two series. The analysis reveals high common power in the precession band (~23kyr) but the phase angles are irregular. This is not surprising given the spectral analysis on the age ensembles, which shows large effects of age uncertainty at 20kyr scales (compared to 40-100kyr). Even if there was a regular behavior, the age uncertainty prevents us from capturing it in the analysis.

WTC (wavelet transform coherence) shows areas of common behavior between the two time series even if there is low power. The analysis reveals coherence in the 23kyr, 40kyr, 100kyr and 400kyr bands, consistent with orbital forcing of climate. The phase angles in the two upper bands are also regular and show and an in-phase behavior in the eccentrivity band (particularly ca 1Ma) and nearly in phase quadrature in the 400kyr band.

Let’s examine whether insolation at 65N has a greater impact on the SST record.

Hypothesis 2: The periodic signals are driven by insolation at 65N.

# insolation values for past 3 Myears at 65N at summer solstice (day 172)
S65 = daily_insolation(lat=65, day=172, orb=orb)

# put in a Pyleoclim.Series object
sol65 = pyleo.Series(time=-years,value=S65,value_name='Insolation at 65 $^\circ$N',value_unit='$W/m^{2}$',
                   time_name='Age',time_unit='ky BP',label='Insolation')

sol65.plot()

Time axis values sorted in ascending order

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Age [ky BP]', ylabel='Insolation at 65 $^\\circ$N [$W/m^{2}$]'>)

psd_sol65=sol65.standardize().spectral(method='mtm',settings={'standardize':False}).signif_test(number=1000,qs=[0.90,0.95,0.99])

Performing spectral analysis on individual series: 100%|█| 1000/1000 [00:45<00:0

psd_sol65.plot()

(<Figure size 1000x400 with 1 Axes>,
 <Axes: xlabel='Period [ky]', ylabel='PSD'>)

scal_sol65 = sol65.standardize().wavelet(method='cwt').signif_test(number=1000)

Performing wavelet analysis on individual series: 100%|█| 1000/1000 [00:48<00:00

scal_sol65.plot()

(<Figure size 1000x800 with 2 Axes>,
 <Axes: title={'center': 'Insolation scalogram (CWT) with 95% threshold'}, xlabel='Age [ky BP]', ylabel='Scale [kyrs]'>)

Spectral and wavelet analysis also shows high, significant power in the precession band. Unlike insolation at 5S, there is stronger power in the obliquity band.

coh65 = ts_slice.detrend().standardize().wavelet_coherence(sol65.standardize(), method='wwz',settings={'standardize':False}).signif_test(number=1000)

Performing wavelet coherence on surrogate pairs: 100%|█| 1000/1000 [3:01:55<00:0

fig, ax = coh65.dashboard()
ax['ts1'].set_ylim([-4,8]) #offset the lines
ax['ts2'].set_ylim([-6,3]) #offset the lines

(-6.0, 3.0)

The XWT analysis reveals high common power in the precession band (~23kyr) and interminently in the obliquity band (~40kyr) but the phase angles are irregular, which we again interpret as a consequence of age uncertainty

The WTC analysis reveals coherence in the 40kyr, 100kyr and 400kyr bands, consistent with orbital forcing of climate. The phase angles in the two upper bands are also regular. Interestingly, it shows an anti-phased behavior in the 100kyr band. This should not be interpreted as a 50kyr lag of sea surface response to the forcing since the eccentricity signal in insolation (the forcing) is not significant.

Conclusion

In this notebook, we demonstrate the use of the Pyleoclim package for spectral and wavelet analysis using a deep sea record from the Eastern Equatorial Pacific. We first confirmed the strong relationship between the benthic \(\delta^{18}O\) record and orbital cycle, which is note surprising since the age model for this record was obtained through alignment to the orbitally-tuned LR04 stack.

The independently obtained SST record also exhbitis significant (despite age uncertainty) power in the orbital band. It is interesting to note that the peaks are not well defined around the precession and obliquity cycles, but rather the periodogram shows several peaks within the 20-40kyr band in addition to the 100-kyr eccentricity cycle.

Cross-wavelet analysis shows high power in the precession band but age model uncertainties preclude any meaningful interpretation of the relative phasing.

References

Khider, D., Ahn, S., Lisiecki, L. E., Lawrence, C. E., & Kienast, M. (2017). The Role of Uncertainty in Estimating Lead/Lag Relationships in Marine Sedimentary Archives: A Case Study From the Tropical Pacific. Paleoceanography, 32(11), 1275-1290.

Lawrence, K. T., Liu, Z. H., & Herbert, T. D. (2006). Evolution of the eastern tropical Pacific through Plio-Pleistocne glaciation. Science, 312(5770), 79-83.

Lin, L., Khider, D., Lisiecki, L. E., & Lawrence, C. E. (2014). Probabilistic sequence alignment of stratigraphic records. Paleoceanography, 29(976-989), 976-989.

Lisiecki, L. E., & Raymo, M. E. (2005). A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanography, 20(PA1003).

Lomb, N. R. (1976). Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science 39, 447-462.

McKay, N. P., & Emile-Geay, J. (2016). Technical Note: The Linked Paleo Data framework – a common tongue for paleoclimatology. Climate of the Past, 12, 1093-1100.

Mix, A. C., Le, J., & Shackleton, N. J. (1995). Benthic foraminiferal stable isotope stratigraphy from Site 846: 0-1.8Ma. Proc. Ocean Drill. Program Sci. Results, 138, 839-847.

Müller, P. J., Kirst, G., Ruhland, G., von Storch, I., & Rosell-Melé, A. (1998). Calibration of the alkenone paleotemperature index Uk’37 based on core-tops from the eastern South Atlantic and the global ocean (60°N-60°S). Geochimica et cosmochimica acta, 62, 1757-1772.

Paillard, D. (2001). Glacial Cycles: Toward a new paradigm. Reviews of Geophysics, 39(3), 325-346.

Scargle, J. D. (1982). Studies in astronomical time series analysis. II. Statistical aspects of spectral analyis of unvenly spaced data. The Astrophysical Journal, 263(2), 835-853.

Scargle, J. D. (1989). Studies in astronomical time series analysis. III. Fourier transforms, aotocorrelation functions, and cross-correlation functions of unevenly-spaced data. . The Astrophysical Journal, 343(2), 874-887.

Shackleton, N. J., Hall, M. A., & Pate, D. (1995). Pliocene stable isotope stratigraphy of ODP Site 846. Proc. Ocean Drill. Program Sci. Results, 138, 337-356.

Schulz, M., & Mudelsee, M. (2002). REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. Computers and Geosciences, 28, 421-426.

Schulz, M., & Stattegger, K. (1997). SPECTRUM: spectral analysis of unevenly spaced time series. Computers and Geosciences, 23(9), 929-945.

Torrence, C., & Compo, G. P. (1998). A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 79, 61-78.

Welch, P. D. (1967). The use of Fast Fourier transform for the estimation of power spectra: A methid based on time averaging over short, modified periodograms. IEEE Transactions on Auio and Electroacoustics, 15(2), 70-73.

Wu, Z., Huang, N. E., Long, S. R., & Peng, C. K. (2007). On the trend, detrending, and variability of nonlinear and nonstationary time series. Proceeding of the National Academy of Sciences, 104, 14889–14894.

Zhu, F., Emile-Geay, J., McKay, N. P., Hakim, G. J., Khider, D., Ault, T. R., et al. (2019). Climate models can correctly simulate the continuum of global-average temperature variability. Proc Natl Acad Sci U S A, 116(8), 8728-8733.