We analyze the temperature time series of the EPICA Dome C ice cores in Antarctica and of the Greenland project, Summit, with durations of 800 000 and 248 000 years, respectively, with a recent mathematical tool defined through the Fourier phases of the series, known as the $J$-index. This data driven index can differentiate between purely random dynamics and dynamics with a deterministic component. It is sensitive to nonlinear components and robust to the presence of noise. Our $J$-index data analysis shows that both Greenland and Antarctica climatic fluctuations possess deterministic traits and suggests the presence of an underlying nonlinear dynamics. Furthermore, in both regions, it reveals the simultaneous occurrence of an important global event known as the “Pelukian transgression.” For Antarctica, it also detects the marine isotopic stage 11. Additionally, our calculation of the time series Hurst exponents and our detrended fluctuation analysis show the presence of long-range persistent correlations for Antarctica and anti-persistent correlations for Greenland. For the latter case, our fractal dimension determinations are indicative of a more complex climatic dynamics in Greenland with respect to Antarctica. Our results are encouraging for further development of climate variability deterministic models for these regions.

The study of climate variability is of the highest priority. Here, we present a comparative paleoclimatic analysis of temperature variations between Antarctica and Greenland. To identify whether the underlying dynamics of these measurements are random or deterministic, and to what extent a nonlinear analysis is required, we calculate a recently developed data driven irregularity index “J,” defined in terms of Fourier phases. This index can differentiate a purely random dynamics from one with deterministic features. It is sensitive to nonlinearity and robust to the presence of noise. In our study, we show that both Greenland and Antarctica climatic fluctuations uncover a noticeable deterministic component and nonlinear traits. Together with Hurst exponent determinations, detrended fluctuation and fractal dimension analyses, we conclude that Greenland presents a more complex evolution, which should be expected given that geologically it has been subject to stronger climatological perturbations than Antarctica. Greenland presents higher fractal dimension, anticorrelated statistical scaling, reduced memory, and higher $J$ irregularity index. In the case of Antarctica, correlation is positive. The $J$ index also evidences a simultaneous transgression climatic event and a global warning instance. The perception of determinism is an incentive for the development of appropriate modeling, which will most likely be more elaborate for Greenland, given its higher content of nonlinearity, inhomogeneity, and irregularity.

## I. INTRODUCTION

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.^{13} 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.^{14}

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.^{15} 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.^{25} 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.^{26} 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.^{30} 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^{31} 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.

## II. DATA DESCRIPTION

In Fig. 1, the temperature fluctuation records from the Dome C ice core for Antarctica (derived from $\delta D$), with a duration of 800 kyr BP, using the EDC3 age model^{16} are displayed. Also, the temperature fluctuations from the Greenland Ice Project, Summit Core (derived from $ \delta 18O$), are shown, spanning a duration of 248 kyr BP with the SS09 age model. In the case of Greenland, conversion from $ \delta 18O$ to temperature was necessary, as a direct temperature record was unavailable, unlike in Antarctica. The conversion was performed using the formula: $ T s=\alpha +\beta \delta +\gamma \delta 2$, where $\delta $ represents the $ \delta 18O$ value of the ice. The employed coefficients were $\alpha =\u2212211.4 \xb0$C, $\beta =\u221211.88 \xb0 C % o$, and $\gamma =\u22120.1925 \xb0 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.

## III. METHODS

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.

### A. J index analysis

^{25}The main idea of the method is, given that the nonlinear properties of a system are manifested in the Fourier phase space, such data driven trajectories may unravel deterministic and nonlinear traits related to underlying dynamics. The methodology used in for the calculation of the J index in this study is the univariate case introduced in Ref. 25. From the climate series $ S 1$, a new series $ S 2= S 1(t+\tau )$ is generated using a lag $\tau $. Subsequently, the Fourier transform is applied to $ S 1$ and $ S 2$ storing their Fourier phases in $ \Theta 1( f k)$ and $ \Theta 2( f k)$ series, ranging from the lowest $ f 1$ to the Nyquist frequency $ f N$. Subsequently, the sequence of points ${ P 1( \Theta 1( f 1), \Theta 2( f 1)),\u2026, P N( \Theta 1( f N), \Theta 2( f N))}$ is created. Since the Fourier phases are angles between $0$ and $2\pi $, points $ P l$ lie on the surface of a torus in three dimensions. To quantify the changes in direction along the trajectory generated by the succession of points $ P l$, from the geodesic lines from point $ P l \u2212 1$ to point $ P l$ and from point $ P l$ to point $ P l + 1$, we define trajectory vectors, for $l=2,\u2026,N\u22121$. Then, the angle $\delta ( f l)$ is measured, which quantifies the change of direction between two consecutive path vectors. This is done on a covering space of the torus.

^{25}Using the sequence of angles $\delta ( f l)$, the $J$ index is computed using the following equation:

By construction, the $J$ index takes values between 0 and 1. The value 0 is achieved when all angles $\delta $ 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.^{25} In the same article,^{25} 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.

### B. Hurst exponent

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.

### C. Detrended fluctuation analysis (DFA)

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^{40} 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.^{39} 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 \alpha $,^{41} where $\alpha $ is equal to $H$ under certain conditions; for further details, refer to Appendix A 2.

## IV. RESULTS

### A. Hurst exponent calculation with R/S method

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 $ \u2212 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 $ \u2212 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.

### B. DFA for Antarctica and Greenland

To gauge the effect of non-stationarity, we applied the DFA method in both regions with a third degree polynomial fit. The value of $\alpha $ for Antarctica was $\alpha =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 $\alpha =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 $\alpha $ 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).

### C. J analysis for Antarctica and Greenland

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.

#### 1. Signs of determinism

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.^{25} This is achieved by creating a random sequence of vectors $( \xi 1(f), \xi 2(f))$ with a length of $N$, where $ \xi i(f)$ are uniformly distributed random numbers between 0 and $2\pi $. 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 \xi $ is recorded. If the value of $J$ for the series of interest is below $ J \xi $, 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 $\tau =1$, $ J \xi G=0.9896$ (the subscript $G$ distinguishes the value of $ J \xi $ in Greenland). In the case of Antarctica, $N=8020$, and $ J \xi A=0.9642$ (the subscript $A$ distinguishes the value of $ J \xi $ in Antarctica) for a delay $\tau =1$.

Figure 2 shows the results of $J$ index for the Antarctica and Greenland series for different values of $\tau $. Notice that $J$ is not statistically dependent on the delay value $\tau $. The red line indicates the value of $ J \xi A$ for Antarctica. The Greenland value $ J \xi 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 \xi $ falls above the $J$ values within the standard deviation. In the case of Antarctica, the gap is much larger.

#### 2. Windowed analysis

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 $\tau =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 $\u223c130$ kyr and 400 kyr BP, as shown in Fig. 3. For Greenland, only the drop around $\u223c130$ 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 $\u223c130$ 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).

## V. DISCUSSION

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^{47} 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,^{53} 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,^{54} 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.^{55} 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 $\u223c130$ 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 ( $\u223c74$–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 ( $\u223c3.5$ m) at approximately 123 500 years ago^{58} (which for paleoclimate issues is not very far from the value of $\u223c130$ 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 $ \xb0$C higher than today.^{59} Particularly, these high sea levels are expected to be observed in MIS 5, 7, 9, and 11. Indeed, during MIS 11 ( $\u223c400$ 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,^{60} 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 $\u223c400$ 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.*^{62} 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 $\u223c130$ 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,^{63} 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.^{64}

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.^{64} 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^{31} conundrum about climate: “It may or may not be deterministic.”

## VI. CONCLUSION

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.^{24}

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: RESCALED RANGE ANALYSIS (R/S)

#### 1. Hurst analysis

^{38}

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.^{70} 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}

Hurst exponent H
. | Window size n
. | Window step-increase . |
---|---|---|

0.7180 | 200 | 5 |

0.7214 | 200 | 10 |

0.7068 | 200 | 15 |

0.6532 | 200 | 20 |

0.6745 | 200 | 25 |

0.7432 | 300 | 5 |

0.7316 | 300 | 10 |

0.6781 | 300 | 15 |

0.6771 | 300 | 20 |

0.6835 | 300 | 25 |

0.6658 | 400 | 5 |

0.6546 | 400 | 10 |

0.6882 | 400 | 15 |

0.6473 | 400 | 20 |

0.7025 | 400 | 25 |

0.6680 | 500 | 5 |

0.6670 | 500 | 10 |

0.6656 | 500 | 15 |

0.6639 | 500 | 20 |

0.6645 | 500 | 25 |

Hurst exponent H
. | Window size n
. | Window step-increase . |
---|---|---|

0.7180 | 200 | 5 |

0.7214 | 200 | 10 |

0.7068 | 200 | 15 |

0.6532 | 200 | 20 |

0.6745 | 200 | 25 |

0.7432 | 300 | 5 |

0.7316 | 300 | 10 |

0.6781 | 300 | 15 |

0.6771 | 300 | 20 |

0.6835 | 300 | 25 |

0.6658 | 400 | 5 |

0.6546 | 400 | 10 |

0.6882 | 400 | 15 |

0.6473 | 400 | 20 |

0.7025 | 400 | 25 |

0.6680 | 500 | 5 |

0.6670 | 500 | 10 |

0.6656 | 500 | 15 |

0.6639 | 500 | 20 |

0.6645 | 500 | 25 |

Hurst exponent H
. | Window size n
. | Window step-increase . |
---|---|---|

0.2267 | 200 | 5 |

0.2277 | 200 | 10 |

0.2276 | 200 | 15 |

0.2290 | 200 | 20 |

0.2266 | 200 | 25 |

0.2320 | 300 | 5 |

0.3097 | 300 | 10 |

0.3525 | 300 | 15 |

0.2810 | 300 | 20 |

0.2830 | 300 | 25 |

0.2315 | 400 | 5 |

0.2213 | 400 | 10 |

0.2111 | 400 | 15 |

0.2831 | 400 | 20 |

0.2864 | 400 | 25 |

0.2726 | 500 | 5 |

0.2608 | 500 | 10 |

0.2632 | 500 | 15 |

0.2480 | 500 | 20 |

0.2595 | 500 | 25 |

Hurst exponent H
. | Window size n
. | Window step-increase . |
---|---|---|

0.2267 | 200 | 5 |

0.2277 | 200 | 10 |

0.2276 | 200 | 15 |

0.2290 | 200 | 20 |

0.2266 | 200 | 25 |

0.2320 | 300 | 5 |

0.3097 | 300 | 10 |

0.3525 | 300 | 15 |

0.2810 | 300 | 20 |

0.2830 | 300 | 25 |

0.2315 | 400 | 5 |

0.2213 | 400 | 10 |

0.2111 | 400 | 15 |

0.2831 | 400 | 20 |

0.2864 | 400 | 25 |

0.2726 | 500 | 5 |

0.2608 | 500 | 10 |

0.2632 | 500 | 15 |

0.2480 | 500 | 20 |

0.2595 | 500 | 25 |

#### 2. Detrended fluctuation analysis (DFA)

^{72}Therefore, there can be DFA-0, DFA-1, and DFA-2 depending on the order of trend filtering.

^{73}A root mean square function $F$ is defined as

Repeating the calculation over segments of different size $n$, if a scaling relationship such as $F\u221d n \alpha $ is obtained , $\alpha $ 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\u2061F(n)=\alpha ln\u2061(n)+C$. Depending on the value of alpha, the following holds:

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

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

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

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

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.

### APPENDIX B: FRACTAL DIMENSION

^{74}According to Refs. 75 and 76, the fractal or Hausdorff dimension quantifies the roughness or smoothness of time series as the observation scale becomes infinitely fine. In other words, the fractal dimension measures a degree of complexity at each scale. A shape with a higher fractal dimension is more intricate or “rough” than one with a lower dimension and fills more space. The fractal dimension is related to the Hurst exponent through the following equation developed by

^{77}

^{78}which is a version of the box-counting estimator that operates on smaller observed scales and avoids the need for general rules in its implementation [see Eq. (B2)]. We also used the semi-periodogram [see Eq. (B3)], wavelets [see Eq. (B4)], variogram, and madogram methods [for both, see Eq. (B5)],

^{79–81}For more details about the methods, refer to Ref. 81. Results for the above fractal dimension determinations are shown in Tables III and IV.

Methods . | D
. |
---|---|

Madograma | 1.341 |

Variograma | 1.340 |

Hall–Wood | 1.349 |

Box–count | 1.302 |

Periodogram | 1.306 |

Wavelet | 1.404 |

Methods . | D
. |
---|---|

Madograma | 1.341 |

Variograma | 1.340 |

Hall–Wood | 1.349 |

Box–count | 1.302 |

Periodogram | 1.306 |

Wavelet | 1.404 |

## REFERENCES

*The Late Cenozoic Glacial Ages*(Yale University Press, 1971), pp. 37–56.

*Cambio Climático: Glaciaciones y Calentamiento Global*, Ciencias Ambientales Series (Fundación Universidad de Bogotá Jorge Tadeo Lozano, 2007).

*Mechanisms of Global Climate Change at Millennial Time Scales*, edited by R. B. Clark, P. U. Webb, and L. D. Keiwin (AGU Geophysical Monograph, 1999), pp. 385–394.

*Nonlinear Time Series Analysis*

*Climatic Determinism*, Meteorological Monographs Vol. 8 (Springer, 1968), pp. 1–3.

*IPCC Special Report on the Ocean and Cryosphere in a Changing Climate*, edited by H.-O. Pörtner, D. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, A. Alegría, M. Nicolai, A. Okem, J. Petzold, B. Rama, and N. M. Weyer (Cambridge University Press, 2022), pp. 203–320.

*Nonlinear Climate Dynamics*

*Scaling Phenomena in Disordered Systems*(Springer, 1986), pp. 1–11.