## Abstract

We examine the long-term correlations and multi-fractal properties of daily satellite retrievals of Arctic sea ice albedo and extent, for periods of approximately 23 years and 32 years, respectively. The approach harnesses a recent development called multi-fractal temporally weighted detrended fluctuation analysis, which exploits the intuition that points closer in time are more likely to be related than distant points. In both datasets, we extract multiple crossover times, as characterized by generalized Hurst exponents, ranging from synoptic to decadal. The method goes beyond treatments that assume a single decay scale process, such as a first-order autoregression, which cannot be justifiably fitted to these observations. Importantly, the strength of the seasonal cycle ‘masks’ long-term correlations on time scales beyond seasonal. When removing the seasonal cycle from the original record, the ice extent data exhibit white noise behaviour from seasonal to bi-seasonal time scales, whereas the clear fingerprints of the short (weather) and long (approx. 7 and 9 year) time scales remain, the latter associated with the recent decay in the ice cover. Therefore, long-term persistence is re-entrant beyond the seasonal scale and it is not possible to distinguish whether a given ice extent minimum/maximum will be followed by a minimum/maximum that is larger or smaller in magnitude.

## 1. Introduction

Earth's polar oceans are viewed as a bellwether of climate change because their surfaces are covered by a thin (several metres) mosaic of high albedo sea ice floes that modulate the atmosphere/ocean heat flux. Hence, as opposed to the massive meteoric ice sheets that are several kilometres thick, sea ice is considered to be a more sensitive component of the cryosphere to perturbations and feedbacks, particularly the ice-albedo feedback that has driven large-scale climate events over Earth history (Saltzman 2002). Indeed, the retreat of Arctic sea ice coverage during recent decades (figure 1) has captured substantial interest (Kwok & Untersteiner 2011; Serreze 2011). An essential question concerns the nature of the decay in ice coverage; is it a trend associated with the influence of greenhouse forcing, or is it a fluctuation in a quantitative record that is short (approx. 30 years) relative to the dynamics of the cryosphere on climatic epochs ( years)? Our goal is to extract the longest intrinsic time scale from the data to provide a framework for geophysical modelling.

Both past climate data and basic physical arguments indicate that a sufficiently large increase in greenhouse gas concentration will drive the decay of the ice cover, whereas the state-of-the-art global climate models under-project the observed recent decay (Kwok & Untersteiner 2011). The stabilizing response to deleterious perturbations of the ice cover was examined by Tietsche *et al.* (2011), who numerically prescribed ice-free summer states at various times during the projection of twenty-first century climates and found that ice extent typically recovered within several years. Such rapid response times can be captured within the framework of relatively simple theory (Moon & Wettlaufer 2011), but, independent of model complexity, both internal and external forcings and their intrinsic time scales manifest themselves in large-scale observations of the geophysical state of the system. Owing to the fact that we cannot *a priori* exclude the observed decline in the ice cover as being an intrinsic decadal oscillation or non-stationary influence in the climate system, we use the finest temporal resolution in the observed record (days) to examine the action of multiple scales, from that of the weather and beyond. The fingerprints of the noisy dynamics of the system on time scales longer than the seasonal record may reside in that record itself, and our goal here is to extract these to provide a framework for geophysical modelling of the underlying processes that created them. It is hoped that the approach will provide theoretical constraints as well as guidelines for comparative data analysis studies.

Most observational studies of the satellite records of ice coverage extrapolate in time the annual or monthly means (Serreze 2011, fig. 1). While the observed declines over this troika of decades (particularly the last decade) are striking, our goal here is to begin a systematic effort to examine the many processes that influence the ice cover, namely the intrinsic time scales reflected in geophysical data. The general methodological framework we employ allows one to examine time-series data in a manner that can distinguish between long-term correlations and trends. We examine whether there exists a multiplicity of persistent scales in the data that can provide a basis for examining cause and effect in the geophysical scale observables of the system. Therefore, we view several types of large-scale sea ice data as a multi-fractal system in which a spectrum of scaling—the singularity spectrum—characterizes the behaviour of nearby points (Feder 1988; Stanley & Meakin 1988). The basic approach of relevance is the multi-fractal generalization of detrended fluctuation analysis (DFA) advanced and widely applied by Kantelhardt and co-workers (Kantelhardt *et al.* 2002, 2006), aptly called multi-fractal detrended fluctuation analysis (MF-DFA). In the last decade, this approach has been developed in many directions, from studying extreme events with nonlinear long-term memory (Bogachev & Bunde 2011) to examining the influence of additive noise on long-term correlations (Ludescher *et al.* 2011). Here, we use a new extension of this methodology called multi-fractal temporally weighted detrended fluctuation analysis (MF-TWDFA), which exploits the intuition that, in any time series, points closer in time are more likely to be related than distant points and can provide a rather more clear signature of long time scales in the fluctuation function and its moments (Zhou & Leung 2010). In §2, we will summarize MF-TWDFA and describe its connection to the conventional approaches of analysing autocorrelation functions and power spectra. The geophysical data are described in §3, after which we present the principal results from MF-TWDFA in §4 before concluding.

## 2. Multi-fractal temporally weighted detrended fluctuation analysis

Here, we harness a new variant of MF-DFA to study the long-range correlations and fractal scaling properties of time-series data of basin-wide Arctic sea ice extent and satellite albedo retrievals. The new variant, called MF-TWDFA, was recently developed by Zhou & Leung (2010) and uses weighted regression to express the intuition that points closer in time are more likely to be similar than those that are more distant. To ensure that this paper is reasonably self-contained, we motivate the approach by first discussing the general rationale for multi-fractals versus the treatment of time series using autocorrelations and spectral analysis. We then describe MF-DFA followed by a development of the aspects that distinguish MF-TWDFA from MF-DFA.

### (a) Rationale for a multi-fractal approach

Consider a time series of length *i*=1,…,*N* describing a quantity *X*_{i}. The linear two-point autocorrelation function *C*(*s*) for realizations of *X*_{i} separated by an increment *s* is given by
2.1in which is the variance and is the mean. *Long-term persistence* in a time series occurs when on average the linear correlations become arbitrarily long (Feder 1988). Therefore, *C*(*s*)∝*s*^{−γ} and the mean correlation time diverges as for 0<*γ*<1, whereas, when *γ*≥1, we have finite and hence we say that the *X*_{i} are *short-term correlated*. Finally, when *C*(*s*)=0 for *s*>0 then the *X*_{i} are *uncorrelated*. Therefore, there may be a time *s*^{★} that delimits correlated from uncorrelated *X*_{i}, namely *C*(*s*<*s*^{★})>0 and *C*(*s*>*s*^{★})=0.

Despite the simplicity of characterizing a system using a single scaling exponent *γ*, it has long been understood that a wide range of natural and laboratory systems have correlations that cannot be captured by such an approach (Feder 1988; Stanley & Meakin 1988). It is particularly important to extract the scaling behaviour for the long-term correlations, wherein the variability on long time scales may be poorly represented by a single-scaling exponent and long-term trends or periodicity may obscure the calculation of the autocorrelation function. Moreover, the distinction between trends and long-term correlations in stationary time series can become blurred. These constitute some of the reasons for the introduction of a multi-fractal description of the correlations which replaces a single value with a continuous spectrum of scaling exponents.

### (b) Multi-fractal detrended fluctuation analysis

There are four stages in the implementation of MF-DFA (Kantelhardt *et al.* 2002). First, one constructs a non-stationary *profile* *Y* (*i*) of the original time series *X*_{i}, which is the cumulative sum
2.2Second, the profile is divided into *N*_{s}=int(*N*/*s*) segments of equal length *s* that do not overlap. Excepting rare circumstances, the original time series is not an exact multiple of *s* leaving excess segments of *Y* (*i*). These are dealt with by repeating the procedure from the end of the profile and returning to the beginning and hence creating 2*N*_{s} segments. Thirdly, within each of the *ν*=1,…,2*N*_{s} segments the variance Var (*ν*,*s*) of the profile relative to a local least squares polynomial fit *y*_{ν}(*i*) of *n*th order is
2.3Finally, the generalized fluctuation function is formed as
2.4and the principal tool of MF-DFA*n* is to examine how *F*_{q}(*s*) depends on the choice of time segment *s* for a given degree of polynomial fit *n* and the order *q* of the moment taken. The scaling of the generalized fluctuation function is characterized by a generalized Hurst exponent *h*(*q*), namely
2.5When the time series is monofractal then *h*(*q*) is independent of *q* and is thus equivalent to the classical Hurst exponent *H*. For the case of *q*=2, MF-DFA and DFA are equivalent (Kantelhardt *et al.* 2002). Therefore, a time series with long-term persistence has *h*(2)=1−*γ*/2 for 0<*γ*<1. However, short-term correlated data, decaying faster than 1/*s*, have *γ*>1 and finite leading to a change at *s*=*s*^{★}, and asymptotic behaviour defined by *h*(2)=1/2. Moreover, the connection between *h*(2) and the decay of the power spectrum *S*(*f*)∝*f*^{−β}, with frequency *f*, is *h*(2)=(1+*β*)/2 (Rangarajan & Ding 2000). Therefore, one sees that, for classical white noise, *β*=0 and hence *h*(2)=1/2, whereas for Brownian or red noise *β*=2 and *h*(2)=3/2.

In general multi-fractal time series exhibit a scaling dependence on the moment *q*. One can relate the generalized Hurst exponents that we focus on here to other multi-fractal exponents (Hentschel & Procaccia 1983; Halsey *et al.* 1986; Feder 1988) and this is the subject of a separate publication. Such complementary exponents are useful because (i) different exponents are more easily extracted from different datasets and (ii) multi-fractality can originate both in a broad probability density as well as in large and small fluctuations having a different long-term persistence. The MF-DFA procedure can distinguish between these (Kantelhardt *et al.* 2002) and thus different exponents can provide tests of distinct multi-fractal origins. Here, in §4*c*, we focus on the generalized Hurst exponents, which make the most transparent contact with the observations that are the topic of our study.

### (c) Moving windows and MF-TWDFA

The generalization of MF-DFA by Zhou & Leung (2010) applies a variant of the weighted least-squares approach to fit the polynomial *y*_{ν}(*i*) to the profile *Y* (*i*) on each interval *ν*. Here, rather than using *n*th order *y*_{ν}(*i*)'s to estimate *Y* (*i*) *within* a fixed window, without reference to points in the profile outside that window, a moving window which is smaller than *s* but determined by the distance between points is used to construct a point-by-point approximation to the profile. Thus, instead of equation (2.3), we compute the variance up (*ν*=1,…,*N*_{s}) and down (*ν*=*N*_{s}+1,…,2*N*_{s}) the profile as
2.6Therefore, we replace the global linear regression of fitting the polynomial *y*_{ν}(*i*) to the data as appears in equation (2.3), with a weighted local estimate determined by the proximity of points *j* to the point *i* in the time series such that |*i*−*j*|≤*s*. A larger (or smaller) weight *w*_{ij} is given to according to whether |*i*−*j*| is small (large).

In the ordinary least-squares method, one minimizes the sum of the squared differences between the predicted and the original function, whereas, in weighted least squares, a weight factor is applied to the squared differences before performing the minimization. There are many varieties of weighted least squares originating in time-series analysis for evenly spaced observations (Macauley 1931) and ‘robust’ locally weighted regression for unevenly spaced data (Cleveland 1979). The standard weighted linear regression scheme fits observed variables *Y* (*i*) and assuming a linear relationship as
2.7where *ϵ*(*i*) is the error, which may (but need not) be homoscedastic. Here, *a*_{k}(*i*) is the *k*th fitting parameter at time *i*, the vector **a** of which is determined from
2.8where **w** is a diagonal matrix having weight elements *w*_{ij} with magnitudes that depend on the proximity to *i*, thereby estimating an *a*_{k}(*i*) according to a prescribed criterion of proximity. Hence, is determined by a weighted least-squares approach using
2.9Here again, the fitting parameters *a*(*i*)=[*a*_{0}(*i*),*a*_{1}(*i*)]^{T} depend on the place *i* in the series and are determined from equation (2.8) in which the *w*_{ij}, with *j*=1,…,*n*, are the diagonals of an *n*×*n* diagonal matrix, and is an *n*×2 matrix with the first column being unity and the second running from 1 to *n*. We use the same bisquare function for the weights
2.10as did Zhou & Leung (2010), although one can envisage a range of possibilities for dealing with the temporal weighting.

We conclude this section by emphasizing that a principal advantage of MF-TWDFA over conventional MF-DFA is the robustness of extracting the scaling of the fluctuation function, and hence the crossover points between scalings, which indicate the underlying processes reflected in the data. Of particular importance is the fidelity of extracting the long-term scaling behaviour in geophysical datasets, which can often be obscured in MF-DFA and, as we shall see below, compromises MF-TWDFA when a strong seasonal cycle remains in the original time series. Moreover, in MF-DFA, the profile of the time series is fitted using discontinuous polynomials, which can introduce errors in the determination of crossover times for new scalings, and can be particularly questionable at long time scales. Finally, for time series of length *N*, while MF-DFA is typically informative only up to *N*/4, MF-TWDFA can be carried out to *N*/2.

## 3. Sea ice geophysical data

We use MF-TWDFA to examine the multi-scale structure of two satellite-based geophysical datasets for Arctic sea ice; the equivalent ice extent (EIE) and albedo retrievals from the Advanced very high-resolution radiometer (AVHRR) Polar Pathfinder (APP) archive. The EIE data derive from retrievals of satellite passive microwave radiances over the Arctic converted to daily sea ice concentration using the NASA Team sea ice algorithm in the same manner as described by Eisenman (2010), who focused on the origin of the difference between ice extent and EIE. We refer the reader to Eisenman's paper for a detailed description of the determination of EIE. The mean EIE seasonal cycle from 1978 to present is shown in figure 2. Daily satellite retrievals of the directional–hemispheric apparent albedo are determined from the APP archive as described in a separate publication dedicated to a different form of analysis from that presented here (Agarwal *et al.* 2011). The apparent albedo is what would be measured by upward and downward looking radiometers and thus varies with the state of the atmosphere and the solar zenith angle.

The APP dataset has been refined for use in a wide range of polar studies and is described in detail in Agarwal *et al.* (2011) and references therein. In brief, the AVHRR channels range from the visible to the thermal infrared (0.58–12.5 μm) and measure top of the atmosphere reflectances and brightness temperatures. With 5×5 km resolution, we analyse albedo retrievals from 1 January 1982 to 31 December 2004, taken daily at 1400 h. Sea ice is distinguished from land and open water using microwave brightness temperatures and filtering the surface type mask data with the NASA Team sea ice algorithm that distinguishes between first-year (FYI) and multi-year ice (MYI) concentrations. The approach ascribes a MYI concentration flag to a region containing at least 50 per cent of this ice type, with uncertainties depending on the season (e.g. melt pond fraction) and region (near the ice edge), along with the surface type categories. On physical grounds, the albedo data are filtered to remove any values greater than 1 or less than 0.2. Each pixel is assessed every day for the presence of ice and then the albedo is averaged for that pixel (Agarwal *et al.* 2011). Examples of histograms for mid-March and mid-September are shown in figure 3.

## 4. Results and discussion

Here, we focus most of our discussion on the nature of the fluctuation functions *F*_{q}(*s*) that emerge from MF-TWDFA and return later to the generalized Hurst exponents. First, we note that we have eliminated the possibility of spurious multi-fractality (Schumann & Kantelhardt 2011) by checking the usual measures Δ*h*_{20}≡*h*(−20)−*h*(20) and Δ*α*. Second, while we show the fluctuation functions for a range of *q*, for clarity, we also show plots solely for *q*=2 and remind the reader that, for *h*(2)=1/2, the system exhibits a completely uncorrelated white noise dynamics, and, for *h*(2)=3/2, there are red (or Brownian) noise correlations. However, for 0<*h*(2)<1/2, the dynamics are anticorrelated. Moreover, as noted above, for analysis by power spectra *β*=2*h*(2)−1. It is evident from equation (2.4) that small temporal fluctuations are characterized by *h*(*q*)'s for *q*<0, whereas large temporal fluctuations are characterized by *h*(*q*)'s for *q*>0. As we discuss below, the fluctuation functions *F*_{q}(*s*) do indeed distinguish the scaling of small and large fluctuations.

On all segments of the *F*_{q}(*s*) plots, the slopes are computed and the ‘crossovers’ refer to the times *s* where the curve changes slope. As a test of the fidelity of the many crossovers in slope (associated with particular intrinsic time scales) detected with MF-TWDFA, we applied *both* MF-DFA and MF-TWDFA. We found that the crossovers for time scales of 2 years or longer would not have been captured by MF-DFA because of large-amplitude fluctuations in this range. Note that, for the EIE data with a seasonal cycle, this result was confirmed for polynomial fitting in the detrending step of MF-DFA of up to order 9, and for the other data up to order 3. Moreover, and that which we focus on here, the MF-TWDFA analysis used on the time series after removing the seasonal cycle leads to the extraction of crossovers associated with long-term persistence that are ‘masked’ when the seasonal cycle is not removed. Next, we discuss this finding in more detail.

### (a) Masking long-term correlations by strong seasonal cycle

A prominent feature of our analysis is the role played by the seasonal cycle. We lay bare the distinction between profiles with and without the seasonal cycle in figure 4, which is the origin of the striking distinction between the upper and lower panels of figures 5 and 6 for the EIE and albedo, respectively. It is seen that the strength of the seasonal cycle is such that it ‘masks’ the dynamics on time scales longer than seasonal, thereby suppressing the fingerprints of long-term persistence. This is clearly displayed in the fluctuation functions. First, in figures 5*a*,*b* and 6*a*,*b*, which have a seasonal cycle, we see that the curves for all *q*'s converge at the 1 year time scale. This convergence is removed when the seasonal cycle is removed as seen in figures 5*c*,*d* and 6*c*,*d*. Second, when the original time series still possesses the seasonal cycle, the slope abruptly drops below *h*(2)=1/2 at the 1 year time scale to transition to an anticorrelated structure. This behaviour is distinct from the clear transitions and positive slopes at longer time scales that are seen when the seasonal cycle is removed from the original time series. While there is a finite union of underlying processes that influence the ice extent and the ice albedo, that union is clearly not one to one (see discussion in Agarwal *et al.* (2011)). Thus, although we find similar qualitative structure in the analysis of the EIE and albedo data with and without a seasonal cycle, there are important quantitative distinctions. For example, the high-frequency weather time scales of the order of weeks are clearly seen in all of the data. However, figures 5*a*,*b* and 6*a*,*b* clearly show masking beginning at 1 year, the transition to anticorrelated behaviour and then a steeper—white noise—slope of *h*(2)=1/2 returning at 5.25 and 2.2 years in the EIE and albedo, respectively. Therefore, all long-term persistence is destroyed by the strength of the seasonal cycle, which is associated with the strong periodicity remaining in the profiles of figure 4*a*,*b*. We see that the general steep/flat/steep behaviour of the slopes of the fluctuation function with increasing time is seen in all panels of figures 5 and 6. However, the essential distinction is that, when the seasonal cycle is removed, long-term persistence is re-entrant. Indeed, time scales of approximately 7 and 9 years are revealed in *both* the EIE and albedo data.

### (b) White noise on seasonal time scales

The region between 1 and 2 years in the EIE data without a seasonal cycle shows *h*(2)≈1/2. In other words, although the strength of the seasonal cycle is such that it dominates the power spectrum (not shown), when removing the seasonal scale from the original EIE data, the system exhibits a white noise behaviour from the seasonal to the bi-seasonal time scales. However, the clear fingerprints of the short (weather) and long (approx. 7 and 9 year) time scales remain. The implications of this finding are, unfortunately for the goal of forecasting, rather clear; a given ice minimum (maximum) can be followed by a minimum (maximum) with a larger or smaller magnitude. The same logic follows for minima/maxima separated by 2 years and the EIE states on time scales between. Most of the discussion in the literature focuses on the extremes in the observations (maxima and minima), which we find here to exhibit a nearly white time distribution. If we interpret the results in the limit *h*(2)=1/2, then the autocorrelation will have strong peaks at the seasonal maxima, minima and bi-seasonal maxima and minima with rapid decays for times from 1 to 2 years. As discussed below, this is consistent with aspects of other studies in which rapid decorrelations are found.

It is important to be clear regarding what these results *do not* tell us. We are not saying that the minima or maxima (or states between) are uncorrelated in time, but solely that the magnitudes of those states cannot be predicted to be larger or smaller from these data alone. Moreover, the lack of an extended autocorrelation does not *a priori* constrain the probability distribution(s) associated with the process(es) from which the value of the observable originates, but rather solely refers to the temporal distribution of the observable. In other words, while a particular probability distribution may underlie the nature of how a given magnitude of the EIE or albedo is reached, we find that the temporal distribution of the sequence of observed magnitudes is white on seasonal to bi-seasonal time scales. Thus, on such time scales, these data are *short-term correlated*, and are associated with a finite and an autocorrelation that decays faster than ∼*s*^{−1}. Therefore, on these short time scales, the processes controlling the EIE or albedo states in the Arctic are the familiar regular strong seasonal scale radiative and advective forcings. Indeed, this is explained by the robustness of the response times associated with the seasonal control of the ice state through the competition between the ice-albedo feedback in summer and the long-wave forcing in winter (e.g. Moon & Wettlaufer 2011; Tietsche *et al.* 2011). However, on longer time scales, we can associate the retreat of the ice cover with the observation of the upswing in *h*(2) towards a more ‘Brownian’ dynamics; *h*(2)∼3/2. Such trends, when accumulating future data over decadal time scales, could well be compromised by a higher degree of noise in the state of the system, but to speculate more moves beyond the level of firm evidence.

We test the idea that the rapid decline in the ice coverage during the most recent decade is dominating the longest crossovers found in our analysis as follows. The EIE record is approximately a decade longer than that for the albedo and covers the decade of the most rapid decline. Thus, we ask what crossovers are seen when we reanalyse this record between 1982 and 2004, the time range corresponding to the albedo dataset? It is found that the two long-term crossovers decrease from 7 to 6.2 and from 9 to 8.4 years. Moreover, the seasonal to bi-seasonality of the white noise persists. Therefore, the increase in the crossover times as we include the recent data appears to be the most plausible reflection of the recent declines. However, we leave for future study what would probably be a less quantitative endeavour of correlating these longer time scales with climate indices that have time scales ranging from 2–7 years (for the El Niño Oscillation) to a decade or many decades (such as for the Pacific Decadal, North Atlantic or Arctic Oscillations).

### (c) Generalized Hurst exponents

Figure 7 shows the trends in the generalized Hurst exponents, *h*(*q*), associated with the strong seasonal cycle. We calculate these exponents, which describe the slope of the fluctuation curve, for each *q* as the linear approximation of the entire curve. It is principally used here to demonstrate the strong influence of the seasonal cycle in suppressing *h*(*q*) for all *q*. There are several important effects immediately evident. First, the general reduction in the magnitudes of the *h*(*q*) associated with the seasonal cycle is synonymous with the suppression of long-term persistence, restricting the crossovers to short-term dynamics and masking behaviour on time scales longer than seasonal, as was seen in the fluctuation functions themselves in figures 5 and 6. Second, this suppression of long-term persistence by definition suppresses multiple crossovers as well as trends on long time scales. Third, the combination of the reduced magnitude and the general flatness of the curves with the seasonal cycle present indicates a suppression of the ‘dynamic range’ of multi-fractality. The analysis and interpretation of related multi-fractal exponents are the subject of a separate publication. The main point of figure 7 is to display the behaviour and interpretation of figures 5 and 6 in a more compact manner.

### (d) Discussion

In his EIE observational analysis, Eisenman (2010) noted that the overall retreat of the ice cover is governed by a noisy signal and particularly so at the September minimum. Observations using approximately 30 years of monthly ice extent and area data show rapid decorrelations of a few months (Blanchard-Wrigglesworth *et al.* 2011) using a lag-correlation method and assuming that the underlying process is a first-order autoregressive process. A similar approach was applied to fit hindcasting simulation output, which showed a persistence time scale for the September ice area of 1.2 years (Armour *et al.* 2011*a*). In a study of the forecast skill of a linear empirical model, Lindsay *et al.* (2008) found no skill in predicting detrended data for time scales of three or more months. Finally, the rapid decay in ice area anomalies is also found in climate models (see Armour *et al.* 2011*b* and references therein) and this is argued to underlie the lack of hysteresis associated with the loss of summer sea ice in both global models and theoretical treatments such as in Eisenman & Wettlaufer (2009).

The corpus of these studies led to the general conclusion of rapidly decaying correlations and hence compromised predictability. Thus, in this respect they are heuristically consistent with our interpretation here. However, there remain important distinctions. First, even a single moment of the fluctuation function demonstrates the existence of multiple time scales in the data which, as noted in §2*a*, cannot be treated in a quantitatively consistent manner with a single decay autocorrelation. Moreover, we showed in §4*a* that if the seasonal cycle is not removed one will always observe a single crossover time (longer than the synoptic time scale) of 1 year; a time scale at which all moments converge. Thus, the upper bound on the persistence time in any study that assumes an autoregressive process will inevitably be approximately 1 year, as is indeed found for ice area (Armour *et al.* 2011*a*; Blanchard-Wrigglesworth *et al.* 2011). Apart from the single moment evidence of multiple scales in the system, the approach here highlights the dangers of not carefully detrending the seasonality and dealing with stationarity. Second, both datasets as used for this analysis and many of the studies mentioned above do not explicitly distinguish between ice types. Thus, this analysis applies to the aggregate of the ice cover, whereas, subject to a wide range of caveats, many models can diagnose features such as thickness or type. The present work focuses on extending this analysis into various ice types.

## 5. Concluding remarks

We have examined the long-term correlations and multi-fractal properties of daily satellite retrievals of Arctic sea ice albedo and EIE, for periods of approximately 23 years and 32 years, respectively, with and without the seasonal cycle removed. A recent development (MF-TWDFA), which exploits the intuition that in any time series points closer in time are more likely to be related than distant points, was adapted for use with these data. Points in the records nearer each other are weighted more than those farther away in order to determine a polynomial used to fit the time-series *profile*, which is a cumulative sum of the time series. As a methodology, the approach offers several advantages over the more generally applied MF-DFA. The profile in MF-DFA is fitted using discontinuous polynomials, which can introduce errors in the determination of crossover times for new scalings, and can be particularly questionable for long time scales. Whereas MF-DFA is typically informative only up to *N*/4 for time series of length *N*, MF-TWDFA can be carried out to *N*/2. Additionally, the generalized fluctuation functions *F*_{q}(*s*) for all moments *q* as a function of time scale *s* are substantially smoother for all *s* and this is particularly so for large values. This facilitates clear extraction of crossover times from one scale to another.

The generalized Hurst exponents and multiple crossover timescales were found to range from the synoptic or weather time scale to decadal, with several in between. Such multiple time scales were exhibited in both datasets and hence the approach provides a framework to examine ice dynamical and thermodynamical responses to climate forcing that goes beyond treatments that assume a process involving a single autocorrelation decay, such as a first-order autoregressive process. Indeed, the method shows that single decay autocorrelations cannot be meaningfully fitted to these geophysical observations. One of our most important findings is that the strength of the seasonal cycle is such that it dominates the power spectrum and ‘masks’ long-term correlations on time scales beyond seasonal. When detrending the seasonality from the original record, the EIE data exhibit a white noise behaviour from seasonal to bi-seasonal time scales, but the clear fingerprints of the short (weather) and long (approx. 7 and 9 year) time scales remain, demonstrating a *re-entrant long-term persistence*. Therefore, it is not possible to distinguish whether a given EIE minimum (maximum) will be followed by a minimum (maximum) that is larger or smaller. This means that while it is tempting to use an anomalous excursion associated with a low ice year to predict the following year's minimum, or that of 2 years henceforth, the present data do not justify such a prediction. It is unfortunate, but it is so.

We tested the idea that the rapid decline in the ice coverage during the most recent decade is associated with the longest crossovers found in our analysis of the EIE. This was done by reanalysing this record between 1982 and 2004, the time range corresponding to the albedo dataset. In so doing, we found that the exclusion of the most recent decade led to the two long-term crossovers decreasing from 7 to 6.2 and from 9 to 8.4 years. Moreover, the white noise structure on seasonal to bi-seasonal times is the same in both records. Hence, the increase in the crossover times associated with including the last decade of data appears to be the most plausible reflection of the recent declines. However, we have not attempted to correlate these longer time scales with other climate indices.

Finally, other methods of analysing such data find solely a rapid decorrelation, whereas we find multi-year and decadal transitions as well as the origin of the dominance of the seasonal cycle in long-term persistence. Hence, we believe that combining such multi-fractal studies of model output and other observations will substantially improve the acuity with which one can disentangle the strength of the seasonal cycle in this highly forced system from the longer term trends.

## Acknowledgements

W.M. thanks NASA for a graduate fellowship. J.S.W. thanks the Wenner-Gren and John Simon Guggenheim Foundations, the Swedish Research Council and Yale University for support. The authors thank N. Untersteiner, A. J. Wells and the referees for their comments and suggestions.

- Received December 14, 2011.
- Accepted February 20, 2012.

- This journal is © 2012 The Royal Society