Nonlinear comparative analysis of Greenland and Antarctica ice cores data

The ice cores located in Greenland and Antarctica have been among the most extensively studied by the scientific community, as evidenced by various works (see Refs. 1–7 among others). These studies encompass the construction of age models, dating methods, ice core retrieval, isotopic stage recognition, and glacial/interglacial and stadial (cold)/interstadial (warm) periods. All these efforts aim to comprehend the climate change over the years using historical data that chronicle the ocean–atmosphere–cryosphere interaction. The comparison of paleoclimatic records to insolation variations (known as orbital tuning methods) is commonly applied to the entire ice cores, as long as the stratigraphy is preserved.^6,8

The records from Greenland ice cores have demonstrated the existence of repeated, prolonged, and abrupt climate changes in the Northern Hemisphere during the last ice age.^9,10 These changes represent one aspect of the global system capable of inducing significant alterations in climatic components such as ocean temperatures^11,12 and monsoon rainfall.¹³ Greenland records offer an archetypal insight into the abrupt climate variability during the last glacial cycle, characterized by rapid shifts between stadial and interstadial conditions, known as Dansgaard–Oeschger oscillations.¹⁴ 

In Antarctica, more than 70 lakes beneath the Antarctic ice sheet have been identified, with particular emphasis on Lake Vostok, a colossal subglacial lake covering an area of 14 000 m $2$⁠. The existence of Lake Vostok holds particular significance for microbiological investigations, aiming to uncover traces of ancient life forms on Earth.¹⁵ Furthermore, the Antarctic ice core from EPICA Dome C (EDC) provides exceptionally high-resolution records spanning changes in local climate, atmospheric composition, and greenhouse gas concentrations over the last nine terminations.^3,16–19 It is noteworthy that Antarctica holds the longest record (800 kyr, where 1000 years = 1 kyr).

Both Greenland and Antarctica provide considerable understanding of past climate and for future insight. Previous climate data analyses have focused on linear studies; though some nonlinear analyses have been carried out, this is rarely the case for paleoclimatic behavior. Such an approach may be more adequate given the intricate processes underlying climate change.^20–24 Due to the interaction among multiple elements shaping past climate, a nonlinear approach seems appropriate. Here, we have implemented a recently developed nonlinear sensitive method.²⁵ Typically, a starting point for the nonlinear analysis of time series involves phase space reconstruction to obtain an attractor topologically equivalent to the original one.²⁶ Additionally, estimators have been calculated to characterize the attractor’s dynamics; key estimators include Lyapunov exponents, attractor dimension, and Kolmogorov–Sinai entropy.^27–29 

A challenge in phase space reconstruction is the somewhat heuristic selection of a certain number of required parameters.³⁰ The success of phase space reconstruction also hinges on certain conditions rarely met by real data, such as stationarity, a substantial amount of data, low noise levels, and fixed data acquisition frequencies. Furthermore, the existence of a finite-dimensional attractor can only be assured for deterministic systems. This last restriction was pointed out since the pioneering work of Lorenz in 1968³¹ for climate studies; the question of whether climate has deterministic dynamics remains a matter of current research.

In this paper, we utilize the statistical index $J$⁠, which was previously defined in Ref. 25, to circumvent the aforementioned complications. This index was chosen because (1) it can be applied to time series from real data, without computationally demanding calculations, (2) it quantifies the degree of irregularity in time series without having to perform a phase space reconstruction, (3) it is sensitive to the presence of non-linear behaviors, (4) it detects signs of determinism with a well-established significance level without the necessity of comparisons with surrogate data, (5) it can be applied in a univariate fashion, and (6) it is robust to the presence of considerable amount of noise contamination in the data, and in our study, provides features for the comparison of Antarctica and Greenland regional behaviors. For the determination of J, we calculated the Hurst exponent and the detrended fluctuation analysis (DFA) exponent, which furthermore characterize the dynamics. Additionally, we also computed self-similarity by means of the series fractal dimension.

The article is organized as follows: In Sec. I, an introduction to previous studies by other researchers on Antarctica and Greenland climatology is provided, focusing primarily on linear analysis. Section II presents a detailed description of the data used in this study. In Sec. III, the calculation of the $J$ index is introduced, as well as an overview of usual complementary statistical analysis methods. Results are presented in Sec. IV. Section V is a summary and discussion of our work. Finally, Sec. VI provides concluding remarks.

In Fig. 1, the temperature fluctuation records from the Dome C ice core for Antarctica (derived from $δD$⁠), with a duration of 800 kyr BP, using the EDC3 age model¹⁶ are displayed. Also, the temperature fluctuations from the Greenland Ice Project, Summit Core (derived from $δ 18O$⁠), are shown, spanning a duration of 248 kyr BP with the SS09 age model. In the case of Greenland, conversion from $δ 18O$ to temperature was necessary, as a direct temperature record was unavailable, unlike in Antarctica. The conversion was performed using the formula: $T s=α+βδ+γ δ 2$⁠, where $δ$ represents the $δ 18O$ value of the ice. The employed coefficients were $α=−211.4 °$C, $β=−11.88 ° C % o$⁠, and $γ=−0.1925 ° C % o$⁠, sourced from Ref. 32. Both series were obtained from http://cdiac.ornl.gov/trends/co2/ice_core_co2.html and were interpolated every 100 years using a stochastic interpolation method developed by Ref. 33.

This section introduces the methods employed for the analysis of ice core data: (i) the $J$ index, sensitive to determinism nonlinearity and irregularity, (ii) the Hurst exponent computed to gauge “memory” and assess fractality, and (iii) DFA to identify correlations.

By construction, the $J$ index takes values between 0 and 1. The value 0 is achieved when all angles $δ$ are constant, while 1 is approached for completely random trajectories asymptotically as the number of data points increases because in this limit Fourier phases have a uniform distribution.²⁵ In the same article,²⁵ a curve of the J index calculated on noise signals is presented, enabling the identification of determinism signs in any type of time series based on its data count $N$ (see Fig. 4 in Ref. 25). Any value below the curve indicates the presence of some deterministic components with a significance level of $1%$⁠. Values above the curve cannot be distinguished from a noise signal.

The $J$ index can be applied to the entire dataset available in the time series and also to sections (windows) of these data. When applied to the entire series, only a global value is obtained. When calculated for data windows, the local evolution of the index along the series can be observed. In our calculations, we chose window sizes for Greenland and Antarctica, adequate for their comparison based on Hurst exponent determinations. Our choice of time delay is explained in the results in Sec. V.

The Hurst exponent exhibits long-term dependence and fractality of the time series.^34–36 Various methods have been employed for its calculation, the most well-known being rescaled range analysis $(R/S)$⁠, a non-parametric method introduced by Ref. 37, where $R$ is the global maximum range and $S$ is the standard deviation. The Hurst exponent varies between 0 and 1^38,39 and can indicate three types of behavior: (i) anti-persistent $(0<H<0.5)$⁠, (ii) persistent $(0.5<H<1)$⁠, and (iii) white noise $(H=0.5)$⁠. For further details, refer to  Appendix A 1.

Another technique that can be used to study the scaling properties of time series is detrended fluctuation analysis (DFA). The advantages of DFA are that it can detect long-term correlations for finite and non-stationary time series. It was first proposed by⁴⁰ while examining DNA nucleotide sequences. It eliminates trends by subtraction of polynomial fits at different data window resolutions $n$ into which the series is divided, thus quantifying fluctuations at different time scales and eliminating trends. The method has been successfully applied in the analysis of long-term climatic records.³⁹ It is widely used to relate scaling properties with fractal behavior in non-stationary time series, without being influenced by trends or non-stationarity in the series. The calculation of DFA is similar to (⁠ $R/S$⁠), except for the choice of the windows of size $n$ and the polynomial fit procedure, yielding a scaling relation $n α$⁠,⁴¹ where $α$ is equal to $H$ under certain conditions; for further details, refer to Appendix A 2.

The calculation of the Hurst exponent enabled us to determine a window size for each climatic region, adequate for the comparison of the $J$ index between Greenland and Antarctica data. The value of the window size was chosen in order to guarantee an invariance of $H$ under step size changes for each of the climatic regions and, hence, a similar degree of stationarity for the calculation of $J$⁠. Additionally, $H$ indicates statistical persistence or anti-persistence in the data, i.e., memory content, for each of the regions. Furthermore, with the determination of $H$⁠, we were able to observe differences between the regions concerning fractal behavior. The results are presented in Tables I and II (see Appendix A 1).

From Tables I and II, we chose for Antarctica a window of 500, for the calculation of $H=0.666$ obtained from taking its average value over the step sizes 5, 10, 15, 20, and 25. With this window, changes in $H$ under modification of step size were the smallest, of the order of 10 $− 3$⁠. This result provides indications of an average persistent behavior. For Greenland, following the same procedure, we obtained an average value $H=0.227$ for a window of 200; here, changes in $H$ were of the order of 10 $− 4$⁠. In this case, the value of $H$ is indicative of an average anti-persistent behavior. Based on the above, we can say that for Antarctica, it is more predictable due to higher memory content. In contrast, for Greenland, there is little dependence on currently observed values with those previously recorded, making its prediction more challenging. Notice that despite both having low temperatures, their difference in Hurst exponents detect different statistical behaviors: Antarctica showing correlated fluctuations and Greenland anticorrelated fluctuations.

Additionally, we obtained the Hausdorff dimension $D$ using Equation (B1) (see  Appendix B) for both stations. In the case of Antarctica, $D=1.334$⁠, and for Greenland, $D=1.773$⁠. In order to ensure these fractal dimensions, we also calculated $D$ using other methods presented in  Appendix B. Tables III and IV show, rouding up, that the fractal dimension for Antarctica ranges between 1.3 and 1.4, while for Greenland, it ranges between 1.3 and 1.7. These values indicate a fractal behavior in both regions, showing greater roughness or complexity in Greenland and, consequently, as mentioned earlier, greater challenges in predicting its future behavior.

To gauge the effect of non-stationarity, we applied the DFA method in both regions with a third degree polynomial fit. The value of $α$ for Antarctica was $α=1.560$ (⁠ $H=0.666$⁠), which gives us very close values between both methods (0.106 difference). In the case of Greenland, the value was $α=1.20$ (⁠ $H=0.227$⁠) with 0.027 difference. DFA plots are shown in Fig. 6 (see  Appendix A 2). Both the Hurst exponent and $α$ yielded similar values in both regions, indicating persistence for Antarctica and anti-persistence for Greenland. This behavior is interesting because, as suggested by Ref. 42, such a series can hide long-term persistence of complex processes. This could be interpreted in terms of regulatory processes of climate variability that have changed and evolved over time. Hence, forecasting models have been developed for various climatic events or processes, always aiming to predict beyond a certain threshold (see, for example, Refs. 43–45).

First, the data were analyzed globally to identify signs of determinism. Subsequently, a window-based study (determined through the calculation of the Hurst exponent) was conducted to locate local changes in dynamics from the data.

To detect significant deviations between a finite time series with $N$ data points and a completely random series, a statistical comparison with data representing the null hypothesis of complete lack of deterministic features is necessary.²⁵ This is achieved by creating a random sequence of vectors $( ξ 1(f), ξ 2(f))$ with a length of $N$⁠, where $ξ i(f)$ are uniformly distributed random numbers between 0 and $2π$⁠. Then, the $J$ index is calculated for this random phase sequence of length $N$⁠. This procedure is repeated for 100 initial conditions chosen from a uniform distribution and the minimum value of $J$⁠; for this ensemble, $J ξ$ is recorded. If the value of $J$ for the series of interest is below $J ξ$⁠, then the series has some deterministic components with a significance value of 0.01; otherwise, it indicates that we cannot distinguish it from noise (for more details, see discussion in Ref. 25). For the series corresponding to Greenland, $N=2500$⁠, with a delay $τ=1$⁠, $J ξ G=0.9896$ (the subscript $G$ distinguishes the value of $J ξ$ in Greenland). In the case of Antarctica, $N=8020$⁠, and $J ξ A=0.9642$ (the subscript $A$ distinguishes the value of $J ξ$ in Antarctica) for a delay $τ=1$⁠.

Figure 2 shows the results of $J$ index for the Antarctica and Greenland series for different values of $τ$⁠. Notice that $J$ is not statistically dependent on the delay value $τ$⁠. The red line indicates the value of $J ξ A$ for Antarctica. The Greenland value $J ξ G$ is not shown because it basically overlaps. Dashed lines show the average of $J$ for each series: average $J$ for Greenland is 0.9683 with a standard deviation of 0.0173, while for Antarctica,the $J$ average is 0.8472 with a standard deviation of 0.01. In both cases, the value of $J ξ$ falls above the $J$ values within the standard deviation. In the case of Antarctica, the gap is much larger.

For the dataset shown in Fig. 1, the $J$ index was calculated in a univariate manner employing segments of 500 records (Antarctica) and 200 records (Greenland), with overlapping windows every 100 data points over the entire time series sequence and lag $τ=1$⁠. In Ref. 25, it is shown that the lag value is not a relevant parameter in the estimation of the $J$ index. Concerning Antarctica, a decrease in $J$ can be observed around $∼130$ kyr and 400 kyr BP, as shown in Fig. 3. For Greenland, only the drop around $∼130$ kyr BP is noticeable, as this region covers a shorter record (see Fig. 4).

Taking into consideration the values of $J$⁠, and identifying declining peaks around $∼130$ kyr BP, we decided to compare the dynamics of both regions at this point. It was possible to observe concordance in the nonlinear dynamic behaviors of both hemispheres (see Fig. 5).

Our results on the Hurst exponent show that, despite both Antarctica and Greenland having cold climates, their value of $H$ differs. This type of discrepancy is in accordance with the findings of Refs. 39 and 46–48 on the determination of the Hurst exponent for various precipitation and temperature series in different regions of the Earth. For example, the precipitation series in a tropical climate (Tamil Nadu, India; Pastaza, Ecuador; Regions of Ethiopia; Australia; Venezuela; Brazil; and Queensland, Australia) exhibited $H<0.5$⁠, except for Venezuela where $H=0.7$⁠. In the case of dry climates (Canary Islands; Zacatecas, Mexico; regions of Israel and Saudi Arabia) (see Refs. 46 and 49–51), similar results were obtained (i.e., some stations had $H<0.5$ while others had $H>0.5$⁠). In the work cited above in this subsection, it is concluded that as precipitation decreases, the value of $H$ increases. Our values of $H$ coincide with this trend, taking into account the findings of Ref. 52 that Antarctica (⁠ $H=0.666$⁠) is an area where the snow accumulation rate is lower compared to Greenland (⁠ $H=0.227$⁠). Moreover, the value of $H$ for Antarctica is close to what is⁴⁷ reported for cold climates (⁠ $H=0.61$⁠) using precipitation time series. As mentioned by Ref. 39, $H$ can explain the spatiotemporal behavior of variables involved in a time series, and its values are correlated with the climatic conditions of the area under study.

On the other hand, from the fractal dimension calculation, both series exhibit self-similar behavior. Given the data availability, in the case of Greenland, this behavior is observed over 240 kyr, and for Antarctica, it extends to 800 kyr. Based on the values of $D$⁠, Greenland shows greater complexity in its climatic variability. This could be attributed to the fact that climate change events were faster in the Northern Hemisphere,⁵³ where warm events known as Dansgaard–Oeschger and cold events known as Heinrich events are located. This observation, although rooted in the past climate, is confirmed in a special IPCC report,⁵⁴ and a NASA study in collaboration with the National Snow and Ice Data Center (NSIDC). The Arctic sea ice reaches its minimum extent in September, decreasing at a rate of 12.6% per decade since 1979. The NSIDC study reported that the reduction in sea ice extent in this region during the first half of August 2023 was faster than the average. Meanwhile, the sea ice in Antarctica has experienced few net changes since 1979.⁵⁵ Hence, the Northern Hemisphere requires greater attention given its significant contribution to climate change. It also reinforces our results regarding the greater complexity found in the region, due to more significant wear and tear on its sea ice.

While various studies using detrended fluctuation analysis have been conducted on temperature variations in different parts of the world, such as the work by Ref. 56 in Spain, similar studies have not been reported in the polar regions.

We have shown that values of $H$ and DFA not only helped find adequate windows for $J$ calculations, but also provided a clearer picture of the climatic behavior of these regions. They emphasize the need to consider regional processes in both areas, rather than assuming that two cold-climate regions will share the same characteristics. For example, our study provides information not previously reported about the temperature in Antarctica and Greenland such as fractality. We have also contributed knowledge about data memory features and predictability.

Regarding the results obtained from $J$ analysis, this index is capable of directly revealing important aspects of the system’s dynamics without the need for additional studies. For both temperature series, a significant drop in $J$ around $∼130$ kyr BP was encountered. This “simultaneity” in geological timing coincides with a substantial warming during Marine Isotope Stage 5e (MIS 5e), drawing attention to the general events during MIS 5 (⁠ $∼74$–130 kyr BP). According to Ref. 57, this period marks the Pelukian transgression, a rise in sea levels. Luminescence data of Pelukian deposits place the high sea level (⁠ $∼3.5$ m) at approximately 123 500 years ago⁵⁸ (which for paleoclimate issues is not very far from the value of $∼130$ kyr BP that we obtained), correlating with Marine Isotope Stage 5e (MIS 5e). Globally, ocean surface temperatures during the peak of this last interglacial (substage 5e) averaged 1–3  $°$C higher than today.⁵⁹ Particularly, these high sea levels are expected to be observed in MIS 5, 7, 9, and 11. Indeed, during MIS 11 (⁠ $∼400$ kyr BP), another marine transgression called Anvilian occurred. It is during MIS 11 that the highest sea levels are reported, implying the disappearance of the Greenland ice sheet,⁶⁰ the entire West Antarctic ice sheet, and part of the East Antarctic ice sheet.^57,61 In our case, due to the length of the Antarctica time series, we were also able to observe a pronounced change in the value of $J$ during MIS 11, approximately $∼400$ kyr BP.

However, according to Ref. 57, it is not clear how a longer but not warmer interglacial, as suggested by MIS 11 indicators in the North Atlantic, far from Greenland, would have caused a notable or even complete loss of the Greenland ice sheet. Koop et al.⁶² analyzed a global compilation of MIS 5e sea-level indicators using a Bayesian statistical approach, arguing that it was 95% probable that both the Antarctic and Greenland ice sheets lost at least 2.5 m of sea level during MIS 5e (though not necessarily simultaneously). This aligns well with the pronounced drop in $J$ observed around $∼130$ kyr BP for both series. The main heat forcings during MIS 5e were seasonal and latitudinal changes in Milankovitch solar insolation at the top of the atmosphere,⁶³ resulting in anomalously high summer insolation in the Northern Hemisphere during the first half of MIS 5e.^64,65 This led to the loss of sea ice and snow, triggering positive albedo feedback that reduced sunlight reflection.⁶⁴ 

According to Ref. 57, Greenland’s paleoclimate data indicate very large warming, but much of that warming was likely due to the loss of the ice sheet itself.⁶⁴ Data from across the North Atlantic for MIS 11 and other interglacials do not show significantly higher temperatures than during MIS 5e, allowing for the possibility of maintaining MIS 5e levels for a more extended period and, hence, producing ice sheet loss. This implies that additional warming within the error bounds of other records, based on the assessment that MIS 5e, was long enough for much of the ice sheet’s response to be completed, requiring additional heat to cause further retreat.

Surprisingly, these anomalous warming events, the main cause of melting major ice sheets, are reflected in the $J$ index. This demonstrates that $J$ is capable of capturing important, relevant dynamics of the study system; moreover, without information loss through filtering or embedding methods. More importantly, $J$ evidences deterministic dynamics for these two regions during the study period covered by the series. This encourages work on forecasting methods, at least for Greenland and Antarctica. We have provided arguments geared to solving Lorenz’s³¹ conundrum about climate: “It may or may not be deterministic.”

The study of Earth’s dynamics, characterized largely by nonlinear behavior, has become an important research area in nonlinear physics and complex systems.^66–68 It is essential to implement mathematical tools that can contribute to our better understanding of these dynamic systems. These nonlinear regime changes can be viewed as alterations in the predictability, regularity, and complexity of signals.²⁴ 

In this work, we successfully identified significant climatological events in both Greenland and Antarctica by examining nonlinear features using the $J$ index. By focusing on the nonlinear dynamics revealed by $J$⁠, we obtained important information, for instance, we found that both series have deterministic features, implying predictability. With the assistance of the Hurst exponent and the fractal dimension of the series, we found that prediction in Antarctica should be easier to accomplish due to its persistence, compared to Greenland with its anti-persistence. To better understand dynamics in the poles, it is necessary to carry out studies involving regional characteristics. Despite that both regions have cold climates, they exhibit different fractality and, consequently, different levels of complexity for their study, while remaining deterministic. This result appears to be attributed to the distinct effects that climate individually exerts on each region. Last, we must reconsider the forecasting methods for studying climate variability in both ice sheets. The encountered determinism calls for the development of more deterministic methods, possibly focusing more on the Northern Hemisphere. This is due to the pronounced degradation of the ice sheet in Greenland, contrasting with the seemingly stable behavior of the ice sheet in Antarctica.

This work was supported by the UNAM Postdoctoral Program (POSDOC), through a scholarship awarded to Dr. Berenice Rojo-Garibaldi. Alberto Isaac Aguilar thanks Conahcyt for a doctoral scholarship. The authors thank Professor Markus Müller for sharing his work on the proposal, development, and calculation of the J method, prior to its publication.

We also acknowledge support from CONACHCYT-Mexico, Grant CF-23377.

The authors have no conflicts to disclose.

Berenice Rojo-Garibaldi: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Alberto Isaac Aguilar-Hernández: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Gustavo Martínez–Mekler: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Anti-persistent series exhibit “mean-reversion” characteristics. If a time series value was high in the previous period, it is likely to decrease in the next period toward the mean value.^38,69,70 These series are also referred to as short-memory series; however, the term “anti-persistent” may obscure the fact that a series with $H<0.5$ has the characteristic of persisting around a central value.⁷⁰ While the term “short memory” may obscure the fact that the series behavior may be following a trend or memory over a much longer period or, more importantly, some form of complex memory.^42,71

In Tables I and II, we present the calculation of $H$ for different window values $n$⁠.

Repeating the calculation over segments of different size $n$⁠, if a scaling relationship such as $F∝ n α$ is obtained , $α$ is called a scaling exponent, a self-affinity parameter representing the long-term power-law correlation properties of the signals. This can be determined by $log⁡F(n)=αln⁡(n)+C$⁠. Depending on the value of alpha, the following holds:

1. $0<α<0.5$⁠. The process has memory and exhibits anti-correlation [it can be modeled by a fractal Gaussian noise (fGn) with $H=α$].

2. $0.5<α<1$⁠. The process has memory and exhibits positive correlations [it can be modeled by fractal Gaussian noise (fGn) with $H=α$].

3. $α=0.5$⁠. The process is indistinguishable from a memory-less random process (it can be modeled by $H=α$⁠)

4. $1<α<2$⁠. The process is non-stationary [it can be modeled as a fractal Brownian motion (fBm) with $H=α−1$].

We add the results obtained by this method with the aim of having a measure better suited for data with a certain degree of nonstationarity complementary to the usual scaling methods for estimating power-law correlation exponents in signals such as fractality and memory determinations, primarily obtained through the Hurst exponent.

In Fig. 6, the results obtained with DFA are shown for both regions.