Contrasting scaling properties of interglacial and glacial climates

Understanding natural climate variability is essential for assessments of climate change. This is reflected in the scaling properties of climate records. The scaling exponents of the interglacial and the glacial climates are fundamentally different. The Holocene record is monofractal, with a scaling exponent H∼0.7. On the contrary, the glacial record is multifractal, with a significantly higher scaling exponent H∼1.2, indicating a longer persistence time and stronger nonlinearities in the glacial climate. The glacial climate is dominated by the strong multi-millennial Dansgaard–Oeschger (DO) events influencing the long-time correlation. However, by separately analysing the last glacial maximum lacking DO events, here we find the same scaling for that period as for the full glacial period. The unbroken scaling thus indicates that the DO events are part of the natural variability and not externally triggered. At glacial time scales, there is a scale break to a trivial scaling, contrasting the DO events from the similarly saw-tooth-shaped glacial cycles.

T he climate system is characterized by complex interactions between atmosphere, oceans, ice, landmasses and the biosphere over a large range of temporal and spatial scales. For understanding natural climate variability and the character of climate change, assessing correlation and persistence times is important. These are reflected in the scaling properties of the climatic records. Scaling was first noted in the seminal work by Hurst 1 on reservoir capacity and runoff in the Nile. Hereafter, it was realized that time series of complex systems such as the climate system are characterized by fractal 2 or even multifractal 3 scaling properties. It has been a long-standing discussion to which extent the fractal nature of the climate dynamics is universal or if it is more specific to the processes and range of scales observed [3][4][5][6] .
To identify the underlying dynamics reflected in scale breaks or robust scaling relations, records covering a large range of temporal scales are necessary. On the atmospheric variability, a range of instrumental records of temperatures and other meteorological parameters have been investigated 3,4,7 . On short time scales the turbulent nature of the atmospheric fields is measured in airplane and drop-sonde campaigns. From these measurements, multifractal scaling of the fields has been reported 8 .
On longer time scales, instrumental records exist for B150 years, typically with daily or twice daily resolution. By filtering out the annual cycle, the scaling properties of temperature variations, covering four to five decades, have been investigated. These indicate universal persistence laws for atmospheric variability 4 , where station data from around the globe shows monofractal (Hurst)-scaling exponents around 0.7. This is significantly different from the value 0.5 characteristic for a trivial white noise process. On even longer time scales we rely on proxy records, where the ice core records are especially suited, as they can be understood as high-resolution sedimentation records from the atmosphere. It was in ice core records that it was first realized that the glacial climate was dominated by millennial scale instabilities, the Dansgaard-Oeschger (DO) events 9 . These events occur all the way from the last inception through to the termination of the last glacial period. One striking feature although is that in the period around the last glacial maximum (LGM, 27-15 kyr) the record only contains a single small event (DO2).
Here we find, by analysing climate periods separately, that the Holocene and the glacial climates have distinctly different scaling properties. The Holocene is monofractal with a scaling exponent HB0.7, whereas the glacial climate is multifractal with H 2 B1.2. The longer persistence time in the glacial period is expected, owing to the presence of the pronounced millennium scale DO events. However, the same scaling is found for the LGM period, indicating that the DO events are not the cause and should be seen as an intrinsic part of the glacial climate. The DO events have a characteristic saw-tooth shape in the records, with rapid warming and slow cooling, similar to the shape of the externally forced glacial cycles with fast terminations and slow inceptions. By analysing the much longer 5 Myr ocean sediment climate record 10 and the 800-kyr Antarctic EPICA (European Project for Ice Coring in Antarctica) ice core record 11 , we find a scale break around 20 kyr, such that on glacial time scales (420 kyr) we have a trivial scaling H ¼ 0.5. We might speculate that this reflects that DO events are internally generated, whereas glacial cycles are externally forced. It is noteworthy that for a trivial red noise process, the socalled Ornstein-Uhlenbeck process, the scaling spectrum shows a cross-over from H ¼ 3/2 for time scales shorter than the correlation time to H ¼ 1/2 for time scales longer than the correlation time. The cross-over time scale indicates the internal time scale of built-up of the large ice sheets.

Results
The climate record. To determine to which extent the scaling behaviour extends to scales beyond the length of the instrumental temperature records, we thus rely on paleoclimatic proxies. The stable isotope (d 18 O and dD) records obtained from the Greenland and Antarctic ice sheets constitute such temperature proxies, with a sufficiently linear relationship to the average atmospheric temperature mixed with an independent noise, which the scaling properties of the atmospheric temperature can be assumed to be reflected in the record 12 .
The issue of nonlinearity and multifractality on multi-glacial cycle time scales has been addressed before in the analysis of the Antarctic Vostok record 13 and comparison with Greenland Icecore Project (GRIP), Greenland Ice Sheet Project (GISP) (Greenland) and Taylor dome (Antarctica) ice core records. The analysis show slightly different scaling exponents between the different ice core records, indicating that dating accuracy and period analysed are important for the results (GRIP and GISP cores should, due to close proximity, give the same result). The influence of chronology on the scaling properties is confirmed by an analysis of the North Greenland Icecore Project (NGRIP) record 14,15 (Fig. 1c), where the scaling properties of the recent 2,000-year d 18 O record is different from the properties of the record older than 2,000 years 16 . This should be expected, as the climate is influenced by different processes operating at different time scales. The NGRIP ice core 14 has been dated with unprecedented accuracy over the past 122 kyr (ref. 15) This enables us to accurately calculate the scaling properties for the different climate states separately.
On multi-millennial time scales, it is known that glacial cycles are linked, in a nonlinear manner 13 , to periodic and quasiperiodic changes in the insolation (incoming solar radiation) from variations in Earth's orbit around the Sun. The climate response is B100 kyr glacial cycles since the middle Pleistocene transition around 1 Myr ago, where the climate has shifted regularly between the glacial and interglacial climate states.
The temperature and proxy records ranging from the instrumental record to five million years, obtained from the stacked deep ocean sediment record 10 , over a huge range of scales are shown in Fig. 1. As there is a connection between spatial and temporal scales, the records shown will be more local in the top panels and more global in the bottom panels. The Greenland ice core record represent a composite of a local and a Northern Atlantic climate signal. To eliminate the influence of long-term trends on the scaling properties, we employ the multifractal detrended fluctuation analysis (MF-DFA) method 17 for the analysis of the records. We have employed the DFA analysis of both first and second order (DFA1 and DFA2), and found that our results are robust in the sense that there is very little difference between the two for the analysed data (see Supplementary Figs 1 and 2). Here we report for DFA1 (see Methods section for description of the MF-DFA). The MF-DFA is complementary to a power spectrum analysis, where the discrete components such as the diurnal and seasonal cycle, as well as the orbital periods on Milankovitch timescales, are mixed with the continuous part of the spectrum 18,19 . If the continuous spectrum scales with frequency P(f)Bf À b , there is a direct linear relation between the scaling exponents for the spectrum and for the fluctuations 20 ; H ¼ (b þ 1)/ 2. Thus, the trivial white noise power spectrum b ¼ 0 corresponds to H ¼ 0.5, whereas the trivial red noise spectrum b ¼ 2 corresponds to H ¼ 1.5. This simply follows from the power spectrum being the Fourier transform of the autocorrelation function. The multifractality is not captured in the powerspectrum.
Multifractal detrended fluctuation analysis. We first analyse the Holocene climate (0-11.7 kyr B2K). The scaling of the Holocene represented by the ice core compares well with the instrumental record (Fig. 2); this shows a remarkable range of scaling over more than five decades, from a day to a few thousand years. In the insert, the corresponding power spectra are shown. Here the pronounced spectral peaks, the year and the Milankovitch periods, are mixed with the continuous part of the spectra. The Holocene climate shows monofractal scaling, with a scaling exponent HB0.7, significantly different from the trivial value H ¼ 0.5. Figure 3a shows the spectra F q (s) for q ¼ ±2, for both Holocene and the glacial records. For Holocene, the black line corresponds to H ¼ 0.7, which fits for both values of q. This is in contrast to the glacial record that shows multifractal scaling with H À 2 ¼ 1.4 and H 2 ¼ 1.2. The Holocene record is tested against a Monte Carlo reshuffling, which preserved the probability density of the data, see Fig. 3b and blue histogram in Fig. 3c. We have also tested the record against an autoregressive process with identical lack-one autocorrelation, with similar results (not shown). To further investigate the reliability, we have generated time series of same length as the record from a fractional Brownian motion, with Hurst exponent 0.7 (Supplementary

Figs 3 and 4)
. The probability density of measured exponents is shown in Fig. 3c (orange histogram). It is noteworthy that the consistency with the record (red bar) is a much weaker result than rejection of the null hypothesis above. It indicates the uncertainty in the estimated scaling exponent due to the limited length of the record. As we do not have a full theory of the underlying climate processes generating the fractal structure, we cannot rule out that the Holocene is weakly multifractal but the time series is to short to detect the multifractality. To obtain some indications on this possibility, we have simulated a known weakly multifractal process of similar length to the record ( Supplementary Figs 5-7), where the multifractallity can indeed be detected. This we interpret as further supporting the observation that the Holocene record is monofractal. The findings for the Holocene is in agreement with findings from climate model millennial simulations 21 . In Supplementary  Fig. 8 we present an analysis of two of the CESM1-CAM5 Last Millennium Ensemble runs.
For longer time scales we enter the glacial climate state and on even longer time scales we observe a scale break at the Milankovitch time scales (420 kyr). Analysis of the 800-kyr Antarctic Epica Dome C (EDC) core 11 with 3 kyr resolution and the 5.3-Myr stacked ocean benthic foraminiferal isotope record 10 (green and orange curves in Fig. 1d,e) does indeed show that these records have a trivial scaling with a Hurst exponent close to 0.5 (green and orange dots in Fig. 2). As a further consolidation of the robustness we also analysed GRIP ice core on the GICC05 time scale, in agreement with the results for NGRIP (Supplementary figs 9 and 10). Furthermore, we have analysed the 420-kyr Vostok core 22 , which agreed with the results for EDC (not shown).
The ice core record (Fig. 1c) shows that the warmer climate of the Holocene period (0-11.7 kyr B2K) is more stable than that of the last glaciation (12-120 kyr B2K) 23 . The difference in climate states is reflected in the scaling properties of the proxy temperature signal; thus, we split the NGRIP d 18 O signal into two parts covering the Holocene warm period and the last glacial period. The glacial period was characterized by millennium scale sudden climate shifts, the DO events 9 , the cause of which is still unknown. The predominant assumption for the cause of the DO s (years) events is abrupt changes in the Atlantic meridional overturning circulation 24,25 perhaps triggered by (negative) changes in freshwater forcing. Many mechanisms have been proposed as a trigger, from oscillations in the ice sheets 26 , ice shelf breakup 27 or sea ice switching 28 , or changes in solar output 29 . The occurrence of DO events influences the correlation time and thus the scaling properties of the record. One could speculate that the larger glacial Hurst exponent is a consequence of the presence of the DO events alone. This is not the case, as the LGM period, 15-27 kyr B2K, with only one short DO event (DO2) does also show a scaling, which is significantly different from the Holocene climate (black crosses in Fig. 3d). To assess the robustness, we furthermore analyse non-overlapping 12 kyr glacial periods (same duration as the Holocene). For all these periods and for the full glacial period, we get multifractal scaling with H ¼ H 2 B1.2 (thin black lines in Fig. 3d). As the DO events are a part of the glacial record scaling, this indicates that they are part of the internal variability and not externally caused, in contrast to the glacial cycles, which are forced by the Milankovitch cycles and show trivial scaling. Confirming the robustness of the results is rather delicate, as there is only a limited set of truly independent paleoclimatic records, each influenced by independent noise processes, which might mask the scaling properties inferred for the climate. One possible test of the observed difference between the interglacial and the glacial climate is to investigate the previous interglacial and glacial periods separately. This we have done by splitting the EDC record into interglacial and glacial periods, and analysing them separately. The result is shown in Fig. 4 and are consistent with the results found for the Greenland ice cores. Furthermore, we have confirmed the results by repeating the same analysis for the EDML and Vostok cores ( Supplementary Figs 11 and 12).
The Greenland stable isotope ice core records are proxies for some large-scale North Atlantic mean temperature. Nearby North Atlantic sea surface temperatures has been reconstructed for the Holocene from ocean sediment cores 30,31 Despite proximity to the Greenland ice core sites, these show much more variability in the Holocene climate. The scaling exponents are HB1.1 for LO09-14 from the North Atlantic current at the Reykjanes Ridge south of Iceland and HB1.4 for MD95-2011 in the Norwegian Atlantic current (Supplementary Fig. 13). The two records are slightly anticorrelated (corr ¼ À 0.25), indicating that they locally monitor fractions of the northward Atlantic heat transport, which can be much more variable than the mean climate of the North Atlantic.
To investigate further their multifractality and difference, Fig. 3e shows the multifractal scaling exponent t q (see Methods section) for Holocene and the glacial period, where a change in slope from negative to positive values of q is observed. This means that there is an asymmetry in the scaling properties of large and small fluctuations. This has been quoted a multifractal phase transition 32 . In multifractal phase transition, two critical q orders defined as q s þ and q s À exist 16 . When q4q s þ or qoq s À , t q will be linear as a function of q, as the largest fluctuations dominate the empirical moments 6 . We observe that q s þ E3 and q s À E À 2 for the glacial record; thus, we can limit the plotting interval for q from À 2 to 3 in Fig. 3e. Figure 3f shows the multifractal spectrum f(a) versus a¼dt=dq. For monofractals the (a, f(a)) spectrum collapses to the point (1, 1). The small red wedge indicates that the Holocene record is monofractal, whereas the black curve shows strong multifractality in the glacial climate.
The tail of the probability distribution. The difference of multifractality in the records of the two periods reflects that the climate of the interglacial is quite different from the glacial climate: The DO events introduce long-range correlations related to the waiting times of several thousand years for jumping between the stadial and the interstadial states. This jumping between states is absent in the Holocene climate. As the DO events do not lead to a scale break in the scaling of the glacial climate signal, we speculate that this is an indication that they do not have a trigger, which is disconnected (such as changing solar radiation) from the climate dynamics giving rise to the scaling properties. The occurrence of the DO events seems random in nature 33 , which agrees well with internal fluctuations in the Atlantic meridional overturning circulation as the cause of these events. The reason why DO variability is absent in the Holocene climate could be attributed to the absence of the large ice sheets. However, this does not necessarily imply ice sheet instability to be the trigger. It could well be that the larger short-term variability in the glacial climate strongly enhance the triggering. The short-term fluctuations can be represented by , where x(t) is the given evenly spaced time series. A generic character of multifractal processes is that they have fat-tailed probability distributions, that is, P (DX4Dx)B(Dx) À g , for large Dx (and similar for the negative tail: P (DXo À Dx)), where g is the corresponding probability exponent 6 . Figure 5 shows the comparison of probability distributions P (|DX|4D) versus Dx for the the Holocene (red) and last glacial (black) (positive and negative tails have identical distributions). The sizes of the fluctuations are larger in the glacial period but interestingly the probability distributions for both the Holocene and the glacial show fat-tail scaling (straight black line) for large Dx with gE À 8.3 (Holocene) and gE À 6.3 (glacial). One could speculate that the less extreme climate of the Holocene prevents triggering of DO events.

Discussion
In summary, the interglacial climate shows scaling over a remarkable range of scales from daily to millennial. The versus Dx between the Holocene (red) and last glacial (black) records.
For both parts of the record, the extreme tail is a power law, P(DX 0 4D)BD g , represented by the straight lines with slopes gE À 8.3 (Holocene) and gE À 6.3 (glacial).

s (years)
s (years)  The glacial record also shows a clear multifractal scaling, with an asymmetry between small and large fluctuations. The DO events are a part of the scaling process, indicating that they are part of the internal variability, and not externally caused, in contrast to the glacial cycles, which are forced by the Milankovitch cycles and show trivial scaling.

Methods
Multifractal detrended fluctuation analysis. The MF-DFA 17 is a robust and easily implemented analysis of the scaling properties in strongly fluctuating or non-stationary time series. It is performed in five steps as follows: (i) determine the cumulated data series where hxi is the mean value of the time series x k (k ¼ 1,...,N). (ii) The profile is divided into N s ¼ int(N/s) disjoint segments with same size s. As the congruence between N and s is often not zero, a part will remain after division. To preserve this part, the same dividing procedure is repeated from the opposite end. As a consequence, 2N s segments are generated. (iii) The variance is calculated as where v ¼ 1, ..., N s and where v ¼ N s þ 1,..., 2N s . y v (i) is the least square-fitting line in segment v. Next, the q-th order fluctuation function is ; for q ¼ 0: (iv) The above steps (ii) and (iii) are repeated as the segment size s increases. (v) The scaling exponent is determined by fitting log-log plots of F q (s) versus s as where H q is the generalized Hurst exponent. Next, the multifractal spectrum (f(a) versus a) can be derived using the following relationship 16,17 : where a is the singularity strength, t q is the multifractal scaling exponent and f (a) is the dimension as a function of the a. x ¼ H 1 À 1 is the scaling exponent of the mean fluctuations. For monofractal time series, the (a, f(a)) spectrum collapses to the point (1, 1). In practice, for finite time series f (a) versus a will be a tiny arc solely due to random fluctuation. For multifractal time series, the multifractal spectrum will typically have a single-humped parabolic shape 17 .