Time and frequency analysis of vostok ice core climate data

Received May 5, 2019 The periodicity of Vostok ice core climate temperature and gas concentration data indicate inherent long term past regularity of Earth’s climate, with a period of around 100,000 years, warming around 15,000 and cooling of around 85,000 years. At this point we are at the top of one of the warming periods. Vostok data cover around 430,000 years, ie 4 climate cycles (warming-cooling), of similar but not quite the same duration. In this paper we perform a detailed time and frequency analysis of these data for each of the cycles as well as their various combinations, including a full tested period of 430,000 years. Time correlation analysis allows for more accurate time lag estimate in each cycle already noted between temperature change and carbon dioxide content. We estimate these lags to lie between 10002500 years, longer than previously concluded. On the frequency side we perform Fast Fourier Analysis and identify full spectrum of harmonics for various cycles, and then perform energy analysis to identify which of the harmonics contributes the most. The idea is to reduce the computational load for further modeling and analysis using Kalman Filter based prediction method. Once the prediction model is defined (a follow up paper) data will be split into two segments, Learning and Testing, in preparation of a Machine Learning fine tunning methodology. We can use last three or last two or even just last cycle to learn on, and then the current on going cycle to test on. This will result in real time prediction of relevant climate data. Assuming causal time regularity, more of these cycles are employed in training, more the prediction error for the next cycle should be reduced. Hence it is critical to perform very detailed time and frequency analysis of Vostok data as a solid data base for the prediction model to follow. Keyword:


Introduction
Extensive climatic data on the past four ice ages and beyond is available from various studies, commencing from mid 1950s until now.Figure 1 shows various sites on Antarctica and Greenland where intensive ice core drilling has occurred since 1956, with several countries supporting more than a score of different drilling projects in the two areas.Currently, intensive ice core drilling is being conducted in other areas as well, so an even larger data set is anticipated.Reference [3] and [12] describe the history of ice core drilling in detail.Our purpose in this paper is to employ Vostok Station ice core data for a variety of time and frequency related analyses.The Vostok ice core data set includes derivations of relative temperature, carbon dioxide (CO 2 ), methane (CH 4 ), oxygen and solar variation (insolation), during the last 430,000 years.Because isotopic fractionation of oxygen-18 and deuterium in snowfall is temperature dependent and a strong spatial correlation exists between mean annual temperature and mean isotopic ratios it is possible to derive ice-core climate records [15].Reference [10] presented the first record to span and full glacial-interglacial from an ice core drilled in the Russian Vostok station in Antarctica.A 430,000 year record was constructed from the Petit et al. study on a 3 km deep core of ice (Fig. 2).Another source of similar ice core data is available from European EPIC drilling project (European Project for Ice Coring in Antarctica) which lasted from 1998 until 2005 (Fig. 3).In this paper we focus on specific analysis related to only two Vostok data variables, namely relative temperature and CO 2 content.Data used are as corrected in [4] and [5].The variation of atmospheric CO 2 and temperature is shown in Fig. 2.Although originally thought that the CO 2 data might be considered as proof of its causal role in global warming, it is now widely considered that CO 2 lags temperature change and its lower rate of solution in and release from sea water is more likely the cause of the relationship.
. Fig. 1 Ice Core drilling sites in Antarctica and Greenland.Note that two maps are not to scale, for Greenland is less than 1/6 th of the size of Antarctica.

Spectral Analysis Methodology
The famous statistician Fisher [7] was an early exponent of testing statistical significance in harmonic analysis.Some preliminary spectral analysis has been conducted on Vostok ice core data set in [21] but without major impact.Here, spectral analysis in R as described in [18] and [20] was applied in [4] and [5] for their data set.In such analyses, time series are decomposed into underlying sine and cosine functions into the most important different frequencies.Various texts on fast Fourier transforms were also sources for frequency determination.These approaches allow construction of periodograms quantifying the contributions of of the individual frequencies to the time series regression.The GeneCycle R Package [19] developed for bioinformatics has general application for time series and was employed in this study.Some recent time series analysis research deals with data classification as well as data series use in environment, production and financial research [24], [25], [26].In our paper we perform our own Vostok data analysis using R and Excel, in both time and frequency domains.Our aim is to use it for Machine Learning training and testing to follow.We analyzed all four Vostok cycles, C1, C2, C3 and C4 as indicated in Fig. 2 above.The cycles are determined by locating maximum values for relative temperature and CO 2 .Becuase of a lag determined between relative temperature and CO 2 content the number of data points is slightly different in two data sets.The overall number of data points for both variables is 363.Individual cycles differ in data points slightly.Table 1 summarises the number of data points in each cycle.Note that the number of data poiunts is not the same for each cycle.This will have important effect on our Machine Learning method and it has to be addressed properly.Also note that the cycle C1, the current cycle we live in, is still evolving and new data may be added if required.This data may be skewed by "global warming effects" which did not exist in previous cycles.Comparison of two data sets may be very useful in estimating global effects of temperature and CO2 content.In addition to the above, the uneven data distribution (non uniform data sampling times within each separate cycle) may pose some numerical issues as well for the analysis in the context of Machine Learning prediction methodology.This will be dealth with in our follow up paper.There are various methods to deal with this as described in [22] [23].One of the simpler methods, yet effective one, is to "fill in" missing data by some approximation method.In any case in this initial paper no additional data insertions are done for smplicity.Our aim is to gain more insight into Vostok data in time and frequency domains, and check the corresponding "energy" content in order to reduce number of significant harmonic components which will reduce complexity of our Machine Learning prediction model.

Data cycles
In this Section we summarise the ice core data analyses.Figs.4-7 show data for relative temperature and CO 2 for each cycle, analyzed individualy as far as time correlations and also harmonic and energy content.

Duration and sampling times in cycles
As Table 1 indicates the number of data points for each cycle is quite different.This has to be taken into account when analysis is performed.For example some data will not be used in cross correlation calculations between temperature and CO 2 .The number of such data is small and will not affect the analysis very much.From Figs. 4-7 we can note general similarity between various cycles data sets.The difference is in the length of each cycle in thousands of years.Table 2 summarizes approximate duration of each cycle based on our analysis (Section 3.4).Note that CO 2 lags relative temperature in C4, C3, and C2 (Table 2).In C1 the data may be a little skewed due to number of recent data points where maximum values are not clearly identified in Vostok data (left most CO 2 data).Another feature to keep in mind is very different time differences when the data is obtained from ice core readings.In some cases difference between two data samples are only few years, or 400-500 years, but in some as much as 5,000 years.This adds to the imprecision of any analysis.For each cycle we calculated average time sampling values as summarized in Table 3.In our follow up paper this will be corrected.The differences in average sampling times obviously come from the number of data points collected and the duration of each individual cycle.To facilitate our Machine Learning approach all other possible data samples have been also created, i.e. any combination of available cycles not necesarily in time order as they transpired.
The idea is to add to the richness of all available data for training purposes of a Machine Learning approach.For example Table 4 indicates respectively average sampling times for Cycle C1 and C2 combined (noted as C12), Cycle 2 and C3 combined (C23), as well as Cycle 2, Cycle 3 and Cycle 4 combined (C234), plus total C1234.The total number of data points, duration and can be easily obtaind from Tables 1 and 2. Note from Table 4 that more cycles we add the more uniform data sampling times become.

Time correlation analysis
Time correlation analysis produces a variety of useful information about periodicity and correlation strenght among data samples of a given quantity.In particular autocorrelations produce the measure of self correlation of a data series.Standard definition of the autocorrelation for a data time series x(t), t = 1,2,...,N, such as relative temperature or CO 2 content, is given as: where m stands for the lag (delay), and m=1,2,...,N-1.Standard definition of crosscorrelation for two discrete time data series x(t) and y(t), such as relative temperature and CO 2 concentration in our case, is given as: R xy (m) = ∑ t x(t)y(t+m), m= -N+1,...,-2,-1,0,1,2,...,N-1 Here m stands fot the lag (delay) and t =1,2,...,N.The above formulas are implemented in various ways depending on the software tool used, such as R, Matlab or Python.In our work reported in this paper we used mostly R for coding and analysis.One property of crosscorrelations is very useful in analysis of relative temperature vs CO 2 content and that is the estimate of the time delay between the two.In general one needs to locate the maximum value point for R xy and locate the corresponding argument, time lag in our case, between relative temperature and CO 2 content i.e.: Fig. 8 illustrates short term calculations for autocorrelations for relative temperature and CO 2 (first two diagrams) as well as the corresponding cross correlation between the two variables (third diagram) for the entire Vostok data data set.We can read the value of τ delay between temperature and CO 2 as approximately equal to 2 lag units.Since the calculation is done for the entire Vostok data set, from Table 4 we can read an average sampling time for C1234 as 1,149 years, hence we can make an upper limit approximation of the size of the delay to be: τ delay < 2 times 1,149 years = 2,298 years ( One needs to note here that the calculation is not very precize primarily due to very non uniform dsitribution of Vostok data.We will address how to harmonize these data in our subsequent work.The point is that the delay is of order of 2,000 years and not of 100 or 200 years.That is an important finding which can influence our thinking about global climate change caused by natural processes that happened prior to CO 2 increase caused by human actions.Instead of calculating cross correlation point by point, one can also calculate correlation coefficient, which are single numbers and they can be used as a simple measures of cross correlation intensity between two variables.There are several coefficients named after their inventors such as Spearman, Pearson and Kendal [6] and they all indicate certain statistical properties which tie two data series together.Table 5 summarizes standard cross coefficient (x) between relative temperature and CO 2 for individual cycles as well as for the entire Vostok data set.The intensity of the cross correlation is quite high, on average more than 0.82 for the entire set.Continuing with our analysis, Fig. 9 indicates entire (long) term auto and cross correlation.It is clearly visible that we are dealing with a periodic process with non uniform periods.This is also confirmed in harmonic analysis in Section 5.In order to determine average time delays between relative temperature and CO 2 for individual we look at Figs. 10-13.The first two diagrams (in Fig. 10 part a) are autocorrelations and they also indicate certain periodicity within the each cycle but obviously not as much on the entire data set.Again this will be confirmed with our harmonic analysis.The third (fourth in Fig. 10) diagram shows short term and long term cross correlation, and they indicates periodicity between two varaibles.The time delay can be read from short term cross correlation.For Cycle 1 it is less than one lag but more than zero lag, we can estimate it at a bit less than half of one lag.From Table 3 for both relative temperature and CO 2 the average data sampling times are 1,703 and 1,605 years which puts the delay at around 800 -850 years.Similar calculations can be done for other cycles as well.To get a more precise approximation we would need more data and finer resolution around the zero lag where the cross correlation is maximum.It is also worth noting that that the maximum values of autocorrelation and cross correlation indicate average energy in the data series itself or between the two series.Next to each harmonic there is an indication (in parenthesis) of harmonic.For example 1 is the basic one, 2 indicates half of it (in years), and so on.For example if the 1st harmonic is 75,000 years, the second will be half of that, i.e. 37,500 years, 10th harmonic will be equal to 7,500 years, and so on.Considering frequency values, 75,000 years corresponds to 1/75,000 = 0.000013333 = 1.3333x10 -5 Hz.
For the entire Vostok data set some of the harmonics, per our anlysis in Section 5, are of order of 375,000 years which corresponds to 1/375,000 = 0.0000026666 = 2.6666x10 -6 Hz.Typically, first several harmonics carry most of the energy, as given in Section 4 on energy considerations.

Harmonic Analysis
The key item in our Machine Learning Vostok ice core data preprocessing summarized in this paper lies in both time correlation analysis (Section 3.3) as well harmonic analysis in this Section.The methodology used is based on calculating FFT (Fast Fourier Transform) using an R language tool and then recalculating to real frequencies.The "real" frequencies (harmonics) are approximations due to a non uniform data points distributions.There are methods described in literature to deal with certain non uniform data sets, but we did not pursue that in this paper.Our approach was to use an average data sampling due to a wide sampling time variability in Vostok data set.Our follow up work addresses nonuniformity of data with filling for missing data.Figs.14-18 indicate FFT diagrams for both relative temperature as well as CO 2 .The frequencies on horizhontal axes are normalized FFT frequencies which need to recalculated to correspond to real frequencies assuming correctness of average sampling times.We also included data sets for individual cycles to have a feel for their variability within each cycle.Note that vertical lines indicate Fourier coefficients which indicate energy level of a specific harmonics.It is obvious from all FFT diagrams that the energy is concentrated in very small frequency range which has several harmonics carrying majority of the energy content of a specific sdata set.Energy analysis is presented in Section 4 with an indication for each cycle or combination of cycles as to the energy distribution across various harmonics.This distribution is a base for making a reasonable approximation in our Machine Learning approach and it will be used to train and learn from past data to predict the future ones.In this case "future" refers to both actual past data which will be used to test how the "training" is sucessfull as well as real future of new data sets to be predicted by our methodology.In order to calculate real frequencies from the FFT ones indicated in Figs. 14 -18 (obtained from data in Figures 4-7) we perform the following simple calculation: f R = f FFT times T s (6) where f R is a real frequency in Hertz, f FFT is normalized FFT frequency and T s is an average sampling time, different for each different cycle scenario.Normalized FFT frequencies go from zero up to 0.5, and FFT algorithm all produce number of normalized frequencies which correspond to half of the number of data points available for analysis.The time periods which correspond to each individual harmonics are calculated by simple inverse operation, i.e.T R = 1/f R (7) FFT assumes that number of data points is always a power of 2, such as 2 6 = 64 data points, or 2 8 = 256.If data set does not have precise power of 2 data points, then certain number of zeroes are appended to the end of the series to make up to power of 2. In our case, since the total number of data points is 363, and the next higher power of 2 above 363 is 2 9 = 512, a total of 149 zeros are appended to the end of both relative temperature and CO 2 .This is all taken care of by R statistical software tool.As an example Table 6 shows the first five most energy containing harmonics for temperature and CO 2 for the entire Vostok data set as well as for four individual cycles.Due to very small frequency values, Table 6 shows the corresponding time duration of one period of each harmonics per formula 7 (harmonic wave length).This is also more natural to visualize and associate to the number of (thousands of years, as well as to the each cycle duration.The numbers are rounded for simplicity.Several observations can be made from Table 6.First, it is not always the case the first harmonic is the most "energetic".For example for the entire Vostok data set cycle (428,000 years approximately) for relative temperature (Table 6), the most energy is carried in the 5th harmonic followed by the 3rd, whereas for CO 2 it is vice versa, and so on.More is described in Section 4 on energy considerations.Second oof all, in most cases both temperature and CO 2 carry very similar, few times equivalent harmonics, but not necessarily with the same enrgy.We also note from the Figs.14 -18 that longer the cycle or cycle combinations are, more harmonics are present which carry considerable energy in the total signal.Another observation is that, for example, for Cycle 4 temperature shows only two very significant harmonics and CO 2 only one.That could due to lack of data points in Vostok set (Table 1), as there are only around 50 data points collected, whereas for other cycles it is 100 or more, for Cycle 2 and Cycle 3, and around 75 for Cycle 1.It is worth nothing that Cycle 1 is still ongoing and more data in the future will be collected.Plus it is not clear if we are on top of Cycle 1, it will take many years to test this.One of important items for future research is to try to identify known climate related events, including Milankovitch Cycles of certain duration which may coincide with the harmonics in our analysis.

Preliminary Time_Frequency Analysis
To conclude this Section we also show very preliminary results on simultaneous time-frequency analysis of Vostok ice core data using Short Term Fourier Transform (STFT).The idea is to se the results for short and long term predictions for both relative temperature and CO 2 via so called Uncertainty Principle which relates resolution in time (short term) vs frequency (long term).One version of the Principle [9] reads as: where loosely speaking Δt and Δf represent essential time and frequency range around a certain moment in the data series time.There are other versions of the Principle as well.This principle resembles classis Heisenberg Uncertainty Principle in physics and comes from very similar mathematical analysis [9].

Energy Consideration
This section summarizes Vostok ice core data energy analysis.The idea is based on choosing only the highest energy harmonics to simplify future Machine Learning algorithm development.For this we set, as an example, a limit that any number of harmonics need to add up to more than 90% of the total data set energy for each cycle or any of their combinations.We expect that this will result in a considerable reduction of required harmonic components for further inclusion in any training and testing approaches.The FFT algorithm in R produces a set of Fourier transform coeffficients which indicate energy level of each harmonic.They correspond to Table 6 harmonics for each cycle as as the entire Vostok data set.The values are higher for CO 2 compared to relative temperature due to larger values of original data points.As mentioned earlier, these values are Fourier Transform coefficients.The total energy in the data set can be me3asured using harmonic components or using time data samples.Per Plancherel Theorem (or Parsseval) both approaches produce the same energy content, and for FFT and DFT (Discrete Fourier Transform) this is stated as in: , n and k = 0,1,...,N-1 (9) where X(k) is the DFT of x(t), both data sets of of length N. The formula applies to both relative temperature as well as CO 2 as they are examples of data sequencies.The Theorem above simply states that there is a oneto-one correspondence between original data set x(t) ND ITS Fourier Transform version in frequency domain.Hence total signal energy in time is equivalent to total ener gy in frequency domain.Because of that we can freely look for energy distribution in frequency domain and use it to choose the most "energetic" harmonics in temperature or CO 2 data set.As the number of these components will be just few of all present but carrying most of the energy (typically 5 components carry more than 90% of the total energy) per Table 7 bellow, this will simplify choice in Machine Learning approach of an appropriate method to utilize these harmonics.We are working to develop Kalman Filter like methodology to employ these results in the energy domain.

Results summary
Various results of the analysis in this paper are summarized in Table 8 for ease of getting a general and complete view of the results.As stated earlier, all possible combinations of four cycle data were analyzed with the idea to enrich the data set for Machine Learning work in progress.Table 8 shows results for four cycles plus few more cycle combinations, as examples.Laft half of the table is devoted to relative temperature and the nrigh halh to CO 2 .C1, C2 etc indicate which cycle the data refers to.On the upper left there is an indication of the length of the cycle period in thousand of years, and also the number of data points for each case, together with the information as to how many "energy" harmonics were chosen.For example, for C1 temperature portion this is 5 harmonics out of 37, which is 13.5%.This corresponds to 93% of the total energy.There is also an indication of average sample time, for C1 temperature it is 1,703 years.Next three sub columns to the right indicate presence of specific harmonic numbers (1, ½ or 1/3, etc.), plus normalized FFT harmonics in time, and finally real harmonics estimate in 1000s of years.This is then repeated for all other cases of temperature (left side of the Table ) as well as for CO 2 on the right side.This way one we can quickly compare temperature data with that of CO 2 as far as how many harmonics each requires, what is the enrgy percentage of the total energy, as well as average time sampling in years, as well as the cycle duration.The table can be expanded to incorporate all other possible cycle combinations for the sake of enriching Machine Learning data sset for training, testing and predicting purposes.We expect more data to be beneficial in order to reduce prediction errors.As the errors get calculated, they can be added to this table for completeness and quick analysis.

Conclusion, Future work and Machine Learning
In this paper we analyzed Vostok ice core data using (i) time correlations, (ii) harmonic analysis, as well as (iii) energy consideration.The general approach split the Vostok data set into 4 smaller sets, as per climate periodiciy indicated in the set.The general outcome is a choice of set of high energy harmonics for all cycles and any of their combinations for Kalman Filter based predictor for Machine Learning prediction purposes.The minimal choice of harmonics will allow us to devise a reasonably simple Machine Learning algorithm for training, testing and prediction purposes.For example long data set in C234 can be used as a Training Data Set, so can C23 (shorter though), in order to "predict" C1, calculate prediction error E1,234 (Error in predicting C1 given C2, C3 and C4 cycles) or E1,23 (preedicting C1 given onbly C2 and C3 data).This applies to both temperature and CO 2 .We will repeat the above analysis for other components in Vostok data set, such as methane, oxygen and insolation.), it is conceivable we will be able to reduce the erros even further, all for the benefit of predicting NOW and in the FUTURE.If a precise numerical correlations are found between Vostok and Milankovitch data sets [14], prediction using our Machine Learning approach can be combined with Milankovitch future cycles as they can be correctly predicted.When chosing harmonics set, various energy targets, such as at 90%, 93%, 95% or 99%, at the expense of incorporating more harmonic components.Another by-product of the temperature data analysis is to estimate the likely water vapour content of atmospheres of different eras.Given that we have proposed positive forcing from irrigation water [11] in addition to other primary sources of warming such as the Milankovitch astronomical cycles.This may prove a more reliable means of correlation using the link already established between water vapour responsible for more than 80% of the heating of air [11].

Fig. 19 Fig. 21 Fig. 22 Fig. 23
Fig. 19 Cycles 1 and 2, temperature and CO 2 harmonic components, FFT normalized Figs. 24 -27 (normalized frequency vs time) show results for STFT for Cycles 1 through 4, for temperature (first figure), CO 2 (second figure) and cross spectrum between the two (third figure).Darker areas indicate higher harmonic energy.Note dark areas in cross correlation are concentrated around certain time and indicate CO2 delay with respect to temperature (for example longer for C1 compared to C4).Note also more energy (dark area) on the right of second diagram in Figure 24 (CO2) comapred to the first diagram Figure 24 (temperature).Compare this with Figure 4 earlier indicating lots of CO2 data "activity" vs much less for temperature towards the end of the C1 record.These kinds of considerations can be used to further fine tune our algorithm for short and long term training and prediction.As stated in Abstract we are working on Kalman Filter Predictor as a base of Machine Learning prediction methodology which we will report in a follow up paper.

Table 1
Number of data points for each cycle

Table 2
Cycles approximate duration in years

Table 3
Average data sampling time in years for each cycle

Table 4
C12, C23, C234 and total C1234 average data sampling times in years

Table 5
Cross correlation coefficient for individual and entire cycle

Table 7
First 5 "energy" harmonics for individual and entire Vostok cycles

Table 8 .
Summary of results PEN Vol. 7, No. 2, August 2019, pp.907-923 Similarly for European EPICA data set as well as set of cycles indicated in Milankovich theory.Hence, once we predict C1 using C123 or C23, we obtain errors E1,234 and E1,23.Intuitively we can expect that E1,23 > E1,234, i.e. training based on larger data set ideally would produce smaller test and prediction errors.With this in mind we can use C1234 to predict time now, or in the near or far away climate future of the data set.Using other C combinations not necesarilly next to each otrher in time (such as C13, C14, C24, C123, C134