Kalman Filter Harmonic Bank for Vostok Ice Core Data Analysis and Climate Predictions

1. Introduction

Extensive climatic data for the past four ice ages and earlier, the period in which Homo sapiens evolved, is available in various scientific reports analyzing ice cores, commencing from the mid-1950s. There are 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. Currently, intensive ice core drilling is being conducted in other areas as well, so an even larger data set is anticipated. In previous papers [1, 2, 3, 4] the history and limitations of ice core drilling are described in detail. Our purpose in this paper is to employ Vostok ice core data for time and frequency related analyses to set a basis for improving prediction of climate variations. The data sets include derivations of relative temperature, carbon dioxide (CO₂), methane (CH₄), oxygen, dust and solar variation (insolation), during the past 420,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. References [5, 6] presented the first record for full glacial–interglacial period from an ice core drilled in the Russian Vostok station in Antarctica. However, it is important to establish criteria for controlling the quality of this new science that we claim could eventually lead to a Climate Prediction Engine (CPE), based on verified causes.

A 420,000 year record was constructed from the [5] study on a 3 km deep core of ice (Figure 1). Another source of similar ice core data is available from European EPICA drilling project [7] which lasted from 1998 until 2005. EPICA data are comparable with Vostok data. In this paper we focus on specific analysis related to only two of the Vostok data variables, namely relative temperature and CO₂ content. We also touch upon possible effect of the dust on various points along Vostok timeline. We will include other data at a later date. Data sets used are from [8] and, with corrections, [9]. The variation of atmospheric CO₂, temperature and dust are shown in Figure 1 together with our definition of four cycles (C1, C2, C3, C4) formed from the maxima of the variables. These cycles could also be defined starting from the variable minima. Before presenting very detailed time and frequency analysis of specific subset of Vostok data, we make some general observations related to climate fluctuations of CO₂, temperature and dust.

Variations in global climate measured by temperature change automatically involve the thermal energy content and heat capacity of the atmosphere. These highly variable systems can be contrasted with the energy content of conservative physical systems like planetary orbits, where the principle of least action defines the trajectories favored; conservative systems can be reliably modeled as an interaction between gravity and inertial forces, with the angular momentum and action held constant as a function of mass and the gravitational constant [10]. By contrast, the molecular states in the global ecosystem are dominated by transient energy flows, oscillating in a complex set of time scales as short as hourly to as long as geological epochs. Such variations in energy flows invite frequency analysis, characterized with Fourier transforms to provide underlying information.

A “snowball” Earth [11], more or less covered by ice and snow as the ice cores were fashioned clearly has diminished total energy in terms of molecular quantum states affecting molecular vibration, rotation, translation − as well as sustaining the more coherent circulatory motions in the atmosphere and the ocean, transferring thermal energy towards the poles. A warm Earth as at present has all these flows in higher quantum states, raising their entropy as we have quantified for atmospheric molecules elsewhere [12]. Variation in physical parameters such as temperature is obvious but complex in causes, with any longer-term trends overlying daily or seasonal trends as a function of latitude.

One state property that does not vary over time on Earth except geologically is atmospheric pressure, given that it is effectively based on the weight of the atmosphere, apart from variation for water, some 4% during the El Nino cycle [13]. The current IPCC consensus concludes that the most significant warming is from greenhouse gas content, such as that of CO₂, methane and nitrous oxide, but this is not necessarily the major controlling cause of increasing temperature, given how negotiable thermal energy is. Globally, we are conducting a large uncontrolled experiment to see if CO₂ is so important. Although originally thought that the CO₂ data from ice cores might be considered as proof of its causal role in global warming, but given that changes in CO₂ lag temperature changes its resultant release or absorption from solution in sea water is more likely the cause of correlations in the relationship.

An additional causative factor controlling the atmospheric pressure of CO₂ that seems to have been overlooked is varying acidity of waters, which can also respond to temperature. In warm periods with high oxygen levels, acidification by oxidation of reduced sulfur, nitrogen and carbon compounds is favored [14]. Alkaline conditions are favored anaerobically. The reduction of dimethyl sulfate, a major oxidant in sea water, can also lead to the evolution of reduced sulfur compounds, which can be converted to sulfuric acid in the atmosphere by ultraviolet oxidation. The alkaline pH of the ocean about 8.2 is effectively the same as that obtained in aerated water equilibrated with limestone (CaCO₃). Sources of alkaline carbonated salts in the ocean is basaltic volcanic rocks and soils, in contrast to more acidic granitic deposits that tend to acidify water.

We have shown using equilibrium theory that a lowering in surface ocean pH of 0.01 units can lead to an increased CO₂ pressure in the atmosphere of 8–10 ppm by volume. So, a change of 100 ppm could occur locally if the ocean surface or water on land was acidified by 0.10 unit, say from pH 8.20 to 8.10 as observed recently. Much larger local changes could occur, from a catastrophic event such as major volcanic eruption. The annual oscillation in CO₂ pressure at 3500 m altitude on Mona Loa, Hawaii could be a result of the change in the high surface temperature of the nearby ocean between winter and summer, rather than imbalance between photosynthetic assimilation and respiratory evolution of CO₂, often claimed in general climate models. The CO₂ pressure in Hawaii peaks after winter in May when the sea temperature is about 23°C and reaches its minimum in late October when the sea is about 27°C. A peculiarity of calcite is that its solubility declines with higher temperature [15], so its precipitation removing CO₂ from the atmosphere could partly or fully explain the oscillation. By contrast, at the sampling site at Cape Grim in Tasmania, where the annual sea temperature has a mean temperature of 15°C with variation between 11 and 19, a range which does not lead necessarily to precipitation of calcite, there is only slight evidence of an oscillation. While the burning of fossil fuels would be the main current cause of pH variation in the ocean, any other acidifying processes in the atmosphere or on land such as from agriculture, could also be controlling influences.

For all of the above reasons regarding the complexity of causes of climate change, we consider it would be beneficial if the ice core data could be subjected to careful Fourier frequency analysis, yielding detailed data regarding long term mechanisms of variation in climate, also pointing to an appropriate dynamic modeling of the underlying processes. In this chapter we focus specifically on CO₂ Vostok data for C1 as defined in Figure 1, as an example of our CPE approach, extending this to CO₂ for C2, C3 and C4, or their combinations. Other data such as temperature, methane, dust, and any other available ice core data such as EPICA, can be analyzed similarly.

2. Methodology

The eminent statistician Fisher [16] was an early exponent of testing statistical significance by harmonic analysis. Some preliminary spectral analysis has been conducted on Vostok ice core data set as reported in [17, 18, 19]. In this paper, a detailed spectral analysis in R and Excel was applied to Barnola et al. data set [8, 9]. In such analysis, time series are decomposed into underlying sine and cosine functions to establish the most important frequencies. Various texts on Fast Fourier Transforms (FFT) were also sources for frequency determination. These approaches allow construction of periodograms quantifying the contributions of the individual frequencies to the time series regression. The methodology developed for bioinformatics [20] has general application for time series and was employed in this study. In this paper we perform our own time and frequency analysis using R and Excel. Trudinger et al. [21, 22] have also applied Kalman filter analysis to ice core data. This allowed a more rigorous analysis of CO₂ variability for the Law Dome ice cores of Antarctica over the recent 1000 years than the usual deconvolution method. These authors pointed out that the Kalman filter allows better calculation of uncertainties in the deduced sources. The uncertainties correspond to the selected range of frequencies. They claimed [22] that it allows investigation of statistical properties that are directly related to physical properties.

Our aim is to apply Kalman filter based methodology to define the logic for a Climate Prediction Engine (Figure 2) based on 4 data cycles, Figure 1. Cycle C1 is the most recent period, the current interglacial warming period nearing its maximum but still incomplete. The overall number of published data points for both variables (relative temperature and CO₂) is 363. The individual cycles are determined by locating maximum absolute values for relative temperature and CO₂. Individual cycles differ in the number of data points slightly. Also, because of a lag observed between the maxima for temperature and CO₂ data, the number of data points differs slightly in two data sets in each cycle. This is also confirmed in our time correlation analysis below. Table 1 summarizes the number of data points in each cycle. Note that each cycle is defined as top-to-top data values for both CO₂ and temperature, resulting in a different number of data points for CO₂ and temperature. However, the beginning of each interglacial where temperatures commences climbing might prove superior for some analyses, given that the controlling causes for reversal may be more consistent for these periods. Obviously, the current cycle is still evolving and new data may be added as required for further analysis. We can achieve that by appending the original Vostok data with additional current points, that might skew the previous natural data progression. Comparison of two data sets including current data may be very useful in estimating global near and far future effects of factors like temperature and CO₂ content − an ultimate goal in this research.

Note that one such additional current data point may not make a large difference on the current C1, as far as harmonic analysis, which has long term time behavior built in, but the short term time effects may be more prominent. Our approach covers both short and long term effects and it can be used to perform sensitivity analysis using current climate data together with original Vostok data. The uneven data sampling times within each cycle) may pose some numerical issues as well for the analysis in the context of Kalman filter (KF) prediction methodology. There are various methods missing data problem and one of the simpler methods, yet effective one, is to enter missing data by some (linear or nonlinear) approximation method. We inject additional data using linear interpolation, also achieving Nyquist sampling requirements and in the process minimize the number of required harmonics used by the Kalman filter harmonic bank (KFHB). For C1 which originally had 80 data points for CO₂ we added additional data points for the total of 128 suitable for FFT but also corresponding to the length of 128,000 years so far in the cycle.

An aim is to gain more insight into Vostok data in time and frequency domains, and check the corresponding amplitude and energy content in order to reduce the number of significant harmonic components (Figure 2). Other authors used similar energy consideration but in our case, we propose the simple and effective mathematical model of a stochastic harmonic oscillator, based on which a KFHB can be built and used to further analyze Vostok data as well as predict near and far future data behavior. Note that besides C1 other cycles can be used to improve the precision by combining two or three cycles based on past data for C1, C2. C3. Amplitude and energy analysis performed here reduces a number of required individual Kalman Filter Harmonic Oscillators (KFHO) as individual blocks for the KFHB. It is important to note that using this approach rather than just combining amplitudes and harmonic components directly as described for example in [23], is to accommodate stochasticity of the data and the overall probabilistic nature of the prediction problem using KFs. This approach adds to the robustness of the method and provides better probabilistic accuracy both for short and long term periods, as well as allowing for prediction sensitivity analysis using various values, such as CO₂ levels as measured now and estimated for some future periods. Our methodology can be applied for both research as well as for policy making tools for the future climate related societal and technological decisions.

In this paper we perform time analysis of all cycles and C1 CO₂ frequency analysis (boldfaced in Table 1) to illustrate the KFHB approach (Section 7).

3. A climate prediction engine

The block diagram in Figure 2 summarizes our CPE. It seeks the ability to predict various quantities using both past data (from one or more cycles) as well as the current data available This also allows to run various sensitivity analysis by using a variety of future scenarios as far as CO₂ and temperature, for example. The CPE can be applied to any ice core data, Vostok, EPICA or any other and for any variable, CO₂, temperature, methane, O₂, or dust.

The CPE consists of three major blocks. First one is a data base of various ice core data, both original, conditioned and inserted as required to make the data more uniformly distributed across the time span of each cycle. The second block is data analysis, in time and frequency domains, including various correlation measures. The third block is the prediction engine which consists of a set of oscillators (KFHO) that produce a KFHB (Section 8). Various predictions as well as sensitivity analysis are performed in this block. Finally prediction parameters are fine-tuned by original ice core data vs. the prediction errors. This is primarily done to fine tune KFHO gains. Section 8 contains all the mathematical details of KFHB.

4. Average sampling time analysis

First, we considered the entire set of Vostok data regarding time sampling. This will yield some initial indication of important issues in dealing with the data and how to perform further Vostok data filling to minimize the effects of non-uniform data distribution. As Table 1 indicates the number of data points and the corresponding time periods for each cycle is quite different. This has to be taken into account when analysis is performed, as far as machine learning use of individual cycles for the benefit of estimating on going C1 data values. Table 2 summarizes approximate duration of each cycle based on our definition of four cycles.

Detailed visual inspection of maximum values in Figures 2 and 3 indicates that CO₂ lags relative temperature in C4, C3, and C2. In C1, the data may be a little skewed due to the number of recent data points where maximum values are still not clearly identified in Vostok data (left-most CO₂ data, Figure 1), as C1 is still evolving. Better lag estimates can be obtained by cross correlation analysis as done in Section 4 below. Another feature is the varied time differences when the data is obtained from ice core readings. In some cases difference between two data samples are only 400–500 years, but in some as much as 5000 years. To start the analysis for each cycle we calculated average time sampling values as summarized in Table 3. In Section 5 this will be corrected to some extent by data filling and determining an ideal Nyquist and data sampling rate.

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 general approach all other combination of available cycles data are also considered. The idea is enrich with the varied of all available data for prediction purposes and also training purposes of machine learning approach. For example Table 4 indicates respectively average sampling times for C1 and C2 combined (C12), C2 and C3 combined (C23), as well as C2, C3 and C4 combined (C234), plus the total data set C1234. Note from Table 4 that more cycles we add the more uniform average data sampling times become, until it becomes equivalent for both CO₂ and temperature. Further refinement of average sampling time is given in Section 5.2 bellow.

5. Time correlation analysis

As a starting point in time cross correlation analysis we examine correlation coefficients − single numbers that can be used as simple measure of cross correlation intensity between two variables. There are several coefficients named after their inventors such as Spearman, Pearson and Kendal [24] and they all indicate certain statistical properties that can relate two data series. Table 5 summarizes standard cross coefficient between relative temperature and CO₂ 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.8 for the entire set. If we split the cycles into up and down sub cycles we obtain Table 6 which indicates cross correlation coefficients for up and down cycle parts. Overall these coefficients indicate bigger spread between up and down sub cycles, and are very sensitive to where the break between up and down parts is chosen.

In general, one of the coefficients (up or down) is considerably larger than the overall single cycle coefficient. This might indicate that the usefulness of the individual up and down cross correlation analysis may be limited of the current C1 cycle. In the context of machine learning methodology this points to putting less emphasis on the cycles from longer in the past compared to the ongoing C1 cycle. To complete this analysis, the down-up period (boldfaced) ratios in Table 6 indicate that the descending period is on average 6 to 8.6 times longer than the ascending period during interglacials, given that the cooling period is that much longer than the warming part. See also Figure 9 for C1 CO₂.

Time correlation analysis produces a variety of useful information about periodicity and correlation strength among data samples of a given quantity. In particular autocorrelations produce the measure of self correlation of a data series, and the cross correlations indicate how two different data sets are correlated. We used standard programs such as R and Excel to generate those. One property of crosscorrelations R_xy is very useful in analysis of relative temperature vs. CO₂ 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₂ content. Figure 3 (with numerical values around zero lag) illustrates short term calculations for cross correlation between the two variables for the total C1234 cycle. We can read the value of the delay τ_delay between temperature and CO₂ as approximately equal to two lag units. Since the calculation is done for the entire data set, from Table 4 we can read an average sampling time for C1234 as 1149 years, hence we can make a rough upper limit approximation of the size of the delay for the entire data set to be:

The calculation is approximate, primarily because the non-uniform distribution of the Vostok data. Some data is separated by only hundred years rather than thousands. This points to a need to make the data more uniform by inserting additional data. We address how to harmonize these data in Section 5. The total delay appears to be of order of 2000 or more years and not 100–200 years. That may be an important finding which can influence our thinking about the role of CO₂ increase caused by human actions. Figure 4 indicates entire data set auto and cross correlation. It is clear that the data exhibits some periodicity. To determine average time delays between relative temperature and CO₂ for individual cycles we examine Figures 5–8 which also have numerical values around zero lag. The first two diagrams in Figure 5 are autocorrelations and they also indicate certain periodicity within the each cycle but obviously not as well as the entire data set. The third diagram shows cross correlation between two variables. The time delay can be read from cross correlation and for C1 it is less than one lag period but maybe more than zero lag, due to the non uniformity of data. We can estimate it as less than half of one period lag. From Table 3 for both relative temperature and CO₂ the average data sampling times are 1703 and 1605 years putting the absolute delay at around 800–850 years or less. Similar approximate estimates can be done for other cycles. Reading from C2 and C3 cross correlations, Figures 7 and 8 and Table 3, by the same consideration the delay is less than half the cycle, which translates into 400–435 years on average, or less. For C4 (Figure 8 and Table 3) the delay appears to be of order of one cycle sampling time, a delay around 1800 years. To get a more precise approximations we would need more uniform data and finer resolution around the zero lag where the cross correlation is at its maximum. Note that on the individual cycle up and down parts this delay may differ from the average cycle level, also indicated by larger correlation coefficients spread in Table 6. We can identify no specific pattern in these coefficients regarding up and down sub cycles having larger or smaller coefficient. Overall, there is a significant CO₂ delay across total Vostok data C1234 compared to the relative temperature. For individual cycles a more detailed cross correlation analysis should be done, especially following data insertion. We made very rough estimates above. For practical reasons, the most important figure to keep in mind is the current cycle C1 delay which points to 800–850 years or less.

Note at the end of this Section that the way various cycles are defined (Figure 1) also affects the analysis. In a follow up work we aim at repeating this analysis by defining cycles from minim to minimum CO2 values. Choice of maximum or minimum values may assist in determining what triggers cycle reversals. We have early indications that the dust might have a significant role in this reversal.

6. Cycle C1 CO₂ frequency analysis

As an illustration of frequency analysis we examine Cycle C1 CO₂ data in details. Similar approach can be used for other data, such as temperature, methane, oxygen and dust. The first concern is to make the C1 CO₂ data more homogeneous hence the sampling time can be more precise and useful.

6.1 Data filling

The original Cycle 1 Vostok data contains a total of 80 CO₂ data points are very unevenly distributed in time. This poses issues with FFT or DFT harmonic analysis. Hence we corrected the situation somewhat by inserting certain number of additional “data points” using linear approximation between the original data points with the biggest time difference between the neighboring points, hence generating a total of 128 original and inserted data points. In particular, the very first point we added to the original Vostok data is the current value of CO₂ which we stated at 350 [25] for the year 1990, as an example. Current levels are estimated at over 400 ppm. Our CPE model can be used for a variety of sensitivity analysis by playing “what if” and processing the data through CPE.

Next data is the first original Vostok data, 2362 years in the past. This will correspond to 64 FFT harmonics described next. Note that the time period of C1 CO₂ data is 128,399 years (Table 2), time difference between the oldest C1 data and the current time, Figure 9 above.

7. Kalman filter harmonic bank

In an illustrative example in this paper we focus on C1 CO₂ harmonic analysis and the corresponding KFHO for the strongest, amplitude wise, first harmonic H₁ in Table 7. The same analysis can be repeated for any harmonic as we elaborate in Sections bellow. We proceed with the development of KFHB which models a series of harmonic components obtained by the energy analysis, which approximates to a reasonable degree (say 90% + or higher cumulative amplitude percentage) the original Vostok data. First step is to define a general Markov model [28] for a generic harmonic oscillator in discrete time described in the following Section. Other authors also used the Kalman filter approach to analyze ice core data, but with a different emphasis, [21, 22] as mentioned in Introduction.

7.2 Kalman filter harmonic oscillator

To facilitate the next step, we rewrite (8) as:

x1t+1x2t+1=−a0−110x1tx2t+w1tw2tE10

where x2(t+1) is just an auxiliary notation for x1tand it is not x2tin (3). Then the standard KF equations in the above case produce [2]:

Prediction Step:

x̂1t+1/t=−a0x̂1t/t−x̂2t/tE11

x̂2t+1/t=x̂1t/tE12

Correction Step:

x̂1t+1/t+1=x̂1t+1/t+K11t+1y∼0t+1E13

x̂2t+1/t+1=x̂2t+1/t+K21t+1y∼0t+1E14

In (13), (14) above, y∼0t=y0t−x̂1t/t−1is Innovation Sequence and filter gains are:

K11t=P11t/t−1/[P11t/t−1+R]E15

K21t=P12t/t−1/P11t/t−1+RE16

The corresponding Prediction Step and Correction Step variances and covariances of the estimation error x∼1t/t−1=x1t−x̂1t/t−1, and x∼1t/t=x1t−x̂1t/t, and similarly for the state x2t, are:

Prediction Step:

p11(t+1/t)=a02p11t/t+2a0p12t/t+p22t/t+Q11.E17

p12(t+1/t)=−a0p11t/t−p12t/t+Q12E18

p22t+1/t=p11t/t+Q22E19

Correction Step:

p11t+1/t+1=[1-K11t]p11t+1/tE20

p12t+1/t+1=[1-K11t]p12t+1/tE21

p22t+1/t+1=−K21tp12t+1/t+p22t+1/tE22

The initial conditions for the above equations are:

p111/0,p121/0,andp221/0E23

and they are determined by the initial Kalman Filter design. One way to determine them is to use matrix Qvalues:

p111/0=Q11,p121/0=Q12,andp221/0=Q22E24

Here the values of Q11and Q22are assumed to be of same the order because they represent uncertainty in modeling x1tand x2(t), and they are just one step apart values of the same state. Simple correlation analysis of x1tand x2(t) indicates that Q12is of order of −a0Q11/4. Figure 10 bellow shows a block diagram of a single ω0KFHO. Here we have the total state vector corresponding to specific harmonic ω0as:

x̂0t/t−1=x̂1t/t−1x̂2t/t−1TE25

7.3 Kalman filter harmonic bank

Once we define single harmonic KF as in Figure 10 we can proceed and construct a KFHB as an assemblage of a number of individual harmonic filters in parallel with the combine outputs to form the original signal (data). We assume a set of harmonics ω=ω1ω2ω3…ωNand for each of ωi, i=1,2,…,Nwe define a separate KFHO as described above. Note that harmonics are related to each other via:

ωi=iω1,i=1,2,…,NE26

and the total signal output (such as Votok data) is the sum of individual harmonic ωi=2πfioutputs:

yt=∑1NyitE27

Figure 11 indicates this arrangement with KFHO where we have:

x̂t/t−1=∑1Nx̂it/t−1E28

representing the total KFHB predicted state estimate. These can be used for short and long term prediction purposes for CO2, temperature or other variables of interest. We can similarly define a set of corrected state estimates which are more precise than predicted with the real data available.

x̂t/t=∑1Nx̂it/tE29

General idea here as compared to a simple sum of harmonic cosine signals (following inverse Fourier Transform) as in (27) with:

yt=∑1Nyit=∑1NAicosωit+θiE30

is to accommodate stochasticity of the underlying Vostok measurement data as well as a simple linear structure of Kalman Filters Harmonic Oscillator, and its ability to make predictions for the signal future values in the probabilistic (more realistic) environment.

8. Vostok cycle C1 CO₂ Kalman filter harmonic oscillators

8.1 First harmonic H₁

The harmonic analysis of Cycle C1 Vostok data is presented in Section 5. In this Section we proceed and build a specific first harmonic H₁ KFHO, as an example of CPE set up, with the initial conditions:

x̂10/0=A1cosθ1,x̂20/0=−2πf1A1sinθ1E31

were A₁ and θ₁ are first harmonic amplitude and phase. In Table 10 we show H₁ KFHO gains from Eqs. (15), (16) for the first seven time samples. The boldfaced values indicate the moment (after 5 time samples) when the gains become constant. These constant gains can be used in the KFHO design for its simplicity. Table 11 summarizes various filter parameters. The results for predicted and corrected state estimates x̂1and x̂2for 128 time samples are shown in Figure 12. The reference harmonic H₁ data is calculated using standard cosine function and it is the “measurement” as in (9) (Table 12).

Optimal	Constant Gain
k1(+)	k2(+)
N/A	N/A
0.96149	0.48075
0.93616	0.29119
0.93633	0.29926
0.93578	0.30137
0.9357	0.30151
0.9357	0.30151
0.9357	0.30151

Table 10.

Filter gains.

Boldfaced entries indicate constant Kalman Filter gains.

Q11	100
Q22	100
Q12	49.759
R	25
a	−1.99759

Table 11.

Parameters.

x1(0)	x2(0)
Initial	Initial
231.0068	0
26.42065	−1.20707
13.48622	−0.81265
4.1098	1.643843
6.050971	−0.45685
0.754028	−1.87793
4.844518	−0.7715
2.016297	1.070199

Table 12.

Initial values.

8.2 Other harmonics (H₂, H₄, H₃, H₇, H₁₀, H₅)

The other harmonic calculations are done likewise, with specific amplitude and phase used and we will not give the obvious details here. The initial conditions are calculated in the equivalent way as given in (31) using different values for amplitude, frequency and the phase angle, using Table 7 harmonic values. Once calculated, initial conditions drive the KFHO for each harmonic, with the various harmonics KFHO parameters similar to the Tables 12–11 above for H₁ and the other harmonic signals like in Figure 12 but with different frequencies and phases. Next all of the predicted and corrected x₁ and x₂ estimates are combined per KFHB output in Eqs. (28) and (29). Resulting values are compared to the known original and inserted Vostok data and appropriate MAPE are generated to check the performance of the method. Table 9 summarizes resulting MAPE for all 7 harmonics. Note that the results are based on constant filter gains in KFHB similar to the values in Table 10 which are equivalent to the optimal time varying gains to at least 2nd decimal point. Parameters in Table 11 were not optimized in any way and they can be tuned further based on errors obtained. This may be done so Q values would be adjusted per specific harmonic amplitude. We made a simple assumption that the standard model deviations and variances are of order of 100 for both states and all harmonics, and cross variance of order of −a0Q11/2, as noted earlier. Table also indicates what we expected, namely that the corrected errors are smaller than the predicted ones.

Note that the values for y are just specific harmonic amplitudes calculated based on maximum amplitude and phase angle, per Eq. (20). The prediction errors correspond to one sampling interval, in our case 1011 years. This is large enough for long terms calculations and prediction. If we chose to have 2 or more sampling intervals the prediction errors will obviously increase.

As far as corrected errors, they correspond to the situation where we have a specific CO₂ value to use for KFHB correction. For example if we chose to consider a current CO₂ value we can append it to the beginning of the Vostok data, and we actually did that by adding 350 PPM value for CO₂ ‘now’ in front of the ‘newest’ Vostok data, more than 2000 years old. The CPE as we envisioned it allows for all sorts of scenarios, sensitivity analysis, ‘what if’ scenarios for all the variables, CO₂, temperature, methane, and so on, EPICA and other ice core measurements, both for short term (100+ years) as well as long term (1000 years+). To give a visual impression of Vostok data approximated by first seven harmonics, see Figure 13 above. Approximation by 7 strongest amplitude harmonics produces very smoothed Vostok data (in blue) without capturing very abrupt changes on a smaller time scale, as in original Vostok data (in red). This can be improved by addition of more harmonics. The advantage which our CPE offers is its flexibility to treat both natural induced values of CO₂ and other variables as well as modern times human produced effects attributed to global warming.

9. Further considerations

Other cycle contributions can be examined as a base for KFHB based machine learning training and testing. This is important due to different duration of various cycles and the uncertainty about the duration of the ongoing Cycle 1. Then older cycles may carry less useful data for the predictions about the current cycle, especially in a view of current global warming effects due to human activities not present in the climate past. On the other hand the old data carries useful short and long term information which can be used judiciously in the KFHB fine tuning, in particular due to its inherent periodicity. Hence multicycle analysis can be very useful for machine learning implementation of our CPE based on KFHB idea. The minimal choice of harmonics will allow us to devise a reasonably simple machine learning algorithm for training, testing and prediction purposes based on KFHB. For example the longer data sets in C234 can be used as a training data set, as can shorter C23, in order to predict the completion of C1, calculating prediction error E1,234 (Error in predicting C1 given C2, C3 and C4 cycles) or E1,23 (predicting C1 given only C2 and C3 data). This applies to both temperature and CO₂. This can be repeated for other components in Vostok data set, such as methane, oxygen and insolation. Similarly for the European EPICA data set as well as set of cycles indicated in Milankovich theory [31]. 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. This has to be confirmed by the further analysis, in particular due to different lengths of the various cycles.

Besides our aim at producing an effective CPE methodology, the harmonic analysis spurred a variety of related thinking and ideas which we also summarize in this paper. Some further ones follow. Climate change on the time scales of the ice cores has been considered as consistent with the IPCC’s hypotheses [11], focusing on permanent greenhouse gases, particularly CO_2, methane and nitrous oxide. This has included the role of increasing water vapor but viewed only as a secondary amplifying factor. The effect of temperature on water vapor pressure is shown by the Clausius-Clapeyron equation, dictating an exponential increase in vapor pressure as temperature rises. But the role of water vapor in modern global warming is only considered in GCMs as a derivative of primary warming by permanent greenhouse gases. This may be in error, as modern irrigation is now adding an extra 4% of water to land surfaces from 1960. It is possible to estimate water vapor content of atmospheres of different eras from temperature data. We have also hypothesized a positive forcing from irrigation water [32, 33] in addition to other primary sources of warming such as the Milankovich astronomical cycles. This may prove a more reliable means of correlation using the link already established between water vapor responsible for more than 80% of the air heating.

One feature of some of the ice core analyses is the irregular but rising increasing levels of dust as the planet became colder (Figure 1), possibly absorbing rather than scattering solar radiation. From frequency analysis, a marked impulse effect can be recognized, given that the peak of the dust samples in ice cores clearly coincides with the commencement of the interglacials, suggesting a role in initiating the warming process. Dust-climate feedbacks have recently been highlighted as having a role in the final stages of the glacial cooling process [34]; but we have surmised, the opposite that dust on the surface of snow and ice has decreased the albedo, capturing more radiant heat from the Sun. If true, a process in which a decreasing albedo effect from widespread dust in very cold, dry conditions (similar to the Antarctic now) initiates warming that is accentuated by increasing water vapor can be proposed. The rapid upward change of temperature in the interglacial would be directly consistent with the exponential increase in water vapor with temperature. In that case, the release of CO₂ could properly be seen as an effect of Henry’s Law [35], rather than cause.

Consequently, the climate sensitivity for CO₂ may be overestimated. It would be of interest to know the scale of wind velocities in the glacials, since the inertial force of wind motion is opposed to gravity in the case of dust suspension and its carriage to high albedo surfaces. When dust particles are dry and disaggregated, current dust storms carrying particles above are most prominent. Furthermore, we have proposed [33, 36] that the atmosphere can store an order of magnitude more thermal energy because of the additional degree of freedom of motion in vertical motion. That is the motion involved in circulating air in anticyclones-cyclones. In a tropical cyclone, condensing water vapor is considered as providing energy to drive the cyclonic motion, using heat derived from the surface of the ocean. This follows from our finding [12] that configurationally entropy (the inverse of free energy) is a logarithmic function of the physical action, a scalar property related to angular momentum but including the dimensionless angular motion [10]. If confirmed this could mean that weather extremes such as very hot days with greater fire risk could be caused by collisions between anticyclones with extra thermal energy released as heat as the laminar flow of air becomes turbulent.

10. Conclusion

In this paper we have analyzed Vostok ice core data using (i) time correlations, (ii) harmonic analysis, as well as (iii) amplitude and energy consideration, and proposed a (iv) general prediction approach using KFHB methodology. In particular we focused on Cycle 1 of CO₂ data in frequency domain as a representative example. The general approach is to split Vostok data set into 4 smaller sets, as per climate periodicity indicated in the set. The outcome is a choice of set of high energy harmonics for all cycles and any of their combinations for designing CPE based on KFHB which is a linear combination of several individual KFHOs, for effective data prediction purposes. This can be incorporated into a practical machine learning methodology for training and testing, as well as data prediction using collected Vostok or other available climate data sets. We believe that our CPE based approach offers advantages in its simplicity and for short as well as long term prediction abilities via KFHB approach which can produce optimal results in the context of stochastic data set. Several issues remain to be solved, in particular (i) uneven cycle lengths as well as better approach to (ii) non uniformity of ice core data. We are addressing both in our ongoing work. Our analysis also seeks to find evidence regarding causes and results in climate science, an essential requirement for more certainty in weather and climate predictions. We consider that current GCMs have considerable uncertainties in such predictions.