## Abstract

The fluctuation statistics of the observed sea-ice extent during the satellite era are compared with model output from CMIP5 models using a multifractal time series method. The two robust features of the observations are that on annual to biannual time scales the ice extent exhibits white noise structure, and there is a decadal scale trend associated with the decay of the ice cover. It is shown that (i) there is a large inter-model variability in the time scales extracted from the models, (ii) none of the models exhibits the decadal time scales found in the satellite observations, (iii) five of the 21 models examined exhibit the observed white noise structure, and (iv) the multi-model ensemble mean exhibits neither the observed white noise structure nor the observed decadal trend. It is proposed that the observed fluctuation statistics produced by this method serve as an appropriate test bed for modelling studies.

This article is part of the theme issue ‘Modelling of sea-ice phenomena’.

## 1. Introduction

Polar amplification and the ice–albedo feedback focus scientific study on the fluctuations in the areal coverage of high latitude ice. By area, on average the Southern Hemisphere ice cover is dominated by the Antarctic ice sheet, whereas in the Northern Hemisphere, the Arctic sea-ice cover dominates. Although we can observe the daily ice cover from space, the substantial changes in Arctic ice mass during recent decades are associated with an ostensibly unmeasurable (approx. 1 W m^{−2}) contribution to the surface energy balance [1]. Hence, given this sensitivity, developing a quantitative understanding of the fluctuations in sea-ice cover is important. According to the Inter-governmental Panel on Climate Change (IPCC) Fifth Assessment Report [2], the most reliably measured characteristic of sea ice is the hemispheric sea-ice extent, to which models are tuned and parametrizations refined. To this end, we compare the statistical structure of satellite observations of ice extent to that from model output.

The global climate models (GCMs) from the IPCC Assessment Report 5 (AR5) project the Arctic to be ice-free by the middle of this century, whereas the previous assessment projected this to occur at the end of the century. The AR5 models project the Arctic to be ice-free as early as 2030 to as late as 2100 ([3], fig. 1 in [4]). The high inter-model variability and the parametrization schemes and tuning used (e.g. the sea-ice albedo, clouds, convective processes) [5–9] constitute key aspects of their veracity at projecting the state of the ice cover.

To quantify the fluctuations in the ice cover from days to decades we have analysed satellite passive microwave data using a multifractal methodology [10]. We studied the Arctic equivalent sea-ice extent (EIE), where the EIE is defined as the total surface area, including land, north of the zonal-mean ice edge latitude and thus is proportional to the sine of the ice edge latitude [11]. (By studying the EIE one minimizes coastal effects.) We find that the EIE is a multifractal in time [10] and thus cannot be explained as an auto-regressive process with a single time scale (a so-called AR-1 process), as is commonly used to characterize Arctic sea ice in GCMs [12–14]. An AR-1 process is inappropriate because (i) the existence of multiple time scales in the data cannot be treated in a quantitatively consistent manner with a single decay autocorrelation. (ii) The strength of the seasonal cycle is such that, if not appropriately removed, retrievals will always produce a single characteristic time of about a year; a time scale at which all moments converge. Hence, the upper bound on the persistence time in any study that assumes an AR-1 process will inevitably be approximately 1 year, as is indeed found for ice area [13,14]. Separate studies with GCMs [15–19] also show that an upper bound exists on the predictability for the Arctic sea-ice cover and that predictability is *skilful* only on the seasonal time scales. Therefore, our approach highlights the dangers of not carefully detrending the seasonality.

## 2. Methods and data

In an effort to improve understanding of Arctic sea-ice behaviour in the AR5 GCMs, we compare the daily satellite observations from 1978–2005 to the model output data with daily frequency for the same period. For parity in this comparison, we analyse the 21 AR5 models (shown in table 1) with daily data in their historical runs, and use the first ensemble member from 11 of these and two from the remaining 10 models. To be consistent in our comparison of all the models and the satellite observations, a common latitude–longitude grid of 0.5° × 0.5° is used and we interpolate the sea-ice concentration (or sea-ice area fraction) from the original model/observation grid onto this common grid. As is done for the observations, the sea-ice extent obtained from the models is then converted to EIE.

Next, we remove the mean seasonal cycle (the mean EIE on each calendar day of the year) from the model output. Figures 1 and 2 show the mean seasonal cycle for the Arctic sea-ice extent and the EIE for both the AR5 models and the observations. The inter-model variability is remarkable. It is noted in the IPCC report *Evaluation of climate models* [2] that models are tuned to match the observed climate system, which in turn helps to give skilful predictions. However, as is evident in figure 2, where although the models may be tuned to match the Arctic sea-ice extent (figure 1), if one calculates the EIE from the model output the variation from the observed seasonal cycle is substantial, not only in the magnitude of the mean seasonal cycle, but also in the qualitative shape of the curves. Namely, some models have a *cycloidal structure*, with a plateau shaped region in winter and a deep well structure in summer. We emphasize again that the EIE not only mitigates the effects of land on the ice extent, but it also appropriately characterizes regional sea-ice characteristics. Thus, even if a model is able to faithfully reproduce the total ice extent, it is possible that regional differences between the models and the observations are large and summing these acts to minimize these differences [20], thereby leading to specious conclusions.

When viewing a time series as a temporal multifractal, there is no constraint on the number of allowable time scales present within the data. Thus, depending on the temporal resolution of the time series, one can extract all time scales corresponding to the physical processes governing the system. Kantelhardt *et al.* [21] developed a multifractal generalization of detrended fluctuation analysis (a modification of the rescaled range method used by Hurst [22] to study the dynamics of river discharge) called multifractal detrended fluctuation analysis (MF-DFA). We use a new variant of MF-DFA called multifractal temporally weighted detrended fluctuation analysis (MF-TW-DFA), which exploits the intuition that points closer in time are more likely to be related than points distant in time, providing a clearer signature of long time scales present in a time series [23]. The detailed algorithm for MF-TW-DFA can be found in ([10,24,25] and references therein) and this is the approach we use to analyse and compare the EIE (without a seasonal cycle) from AR5 models to that from the satellite observations. Because of the fact that there are no *a priori* assumptions made regarding the physical processes in the system, this method has also been used to explain the dynamics of Arctic sea-ice velocity fields [24], and to detect exoplanets [25] and their atmospheres [26]. We note again that when the seasonal cycle is not removed from the sea-ice data multifractality is inhibited and hence all the time scales longer than annual are masked, which leads to erroneous interpretation of data or model output as obeying an AR-1 process.

The fluctuations in a time series are quantified with respect to a smooth *profile*, defined as the cumulative sum produced by time-weighting the data. These fluctuations are then combined to construct the fluctuation function, *F*_{q}(*s*), for each time scale, *s*, under examination, where *q* denotes the *q*th moment of this function. For a given moment, the key behaviour examined is the *s*-dependence of *F*_{q}(*s*), which is characterized by a generalized Hurst exponent *h*(*q*), namely
2.1For example, in a monofractal time series *h*(*q*) is independent of *q*, and thus equivalent to the classical Hurst exponent *H*. The exponent *h*(2) is related to the decay of the power spectrum, *S*(*f*). Hence, if *S*(*f*)∝*f*^{−β}, with frequency *f*, then *h*(2) = (1 + *β*)/2 (e.g. [27]). Thus, for white noise *β* = 0, which gives , whereas for red noise *β* = 2, giving . Therefore, the dominant time scales found using MF-TW-DFA also capture the corresponding temporal dynamics, such as white noise, red noise and correlation structure, which allows us to construct stochastic models for the observed processes, such as the statistical structure and dynamics of sea-ice velocity fields [24]. The ‘crossover’ or ‘dominant’ time scale is defined as that at which the fluctuation function log_{10}*F*_{2}(*s*) changes slope with respect to log_{10}*s*. These occur when slope of the curve exceeds a threshold of *C*_{th} = 0.01.

## 3. Discussion

Here we compare the fluctuation functions from the AR5 models to the observations. The nature and number of figures are such that the majority of the results are presented in appendix A, and here we summarize the key features, which are (i) large inter-model variability in the crossover time scales extracted from the models, (ii) none of the models exhibit the decadal time scales found in the satellite observations, and (iii) only five (CNRM-CM5, MRI-CGCM3 , MRI-ESM1, FGOALS-g2, MIROC5 (second ensemble member)) of the 21 models exhibit the white noise structure (shown in figure 3) observed in the satellite record on annual to biannual time scales. Additionally, only four models exhibit time scales approaching decadal (7.2 years in BCC1-1, 4.8 years in FGOALS-g2, 7.9 years in NorESM1, 4.5 years in the second ensemble member of NorESM1 and 6 years in GFDL-ESM2G), whereas in all of the other models the longest time scale is of the order a year.

The high-frequency processes governing shorter time scales are associated with the influence of synoptic systems, the intermediate time scales are seasonal and the longer—decadal—time scales are generally ascribed to climate processes [10]. For example, Agarwal & Wettlaufer [24] demonstrated that on annual–biannual time scales, the sea-ice extent is largely controlled by the wind fields over the Arctic. By comparing the sea-ice volumes from the models with the PIOMAS [28] dataset, Shu *et al.* [29] showed that the sea-ice thickness produced in the CMIP5 models is much less than what is observed. Therefore, although analysis of the Arctic sea-ice velocity fields from the CMIP5 models shows that most models reproduce the statistical characteristics of the observed velocity fields, if the sea ice produced by the models is too thin, the spectrum of its response to wind forcing will change. For example, Colony & Thorndike [30] showed that local low-frequency ice motion is linearly correlated with local synoptic wind forcing, but the local inertial motion is not, despite a high coherency between the inertial motion of the ice on the same length (approx. 100 km) scales. This highlights the role of mechanical interactions between ice floes on short time scales, which would not be captured in models producing thinner ice.

Wang & Overland [4] analysed the sea-ice extent in the CMIP5 models and concluded that the Arctic would be ice-free in the next 14–36 years with a median of 28 years, based on the model spread. Massonnet *et al.* [31] tried to reduce the CMIP5 uncertainty for the Arctic to become ice-free, finding a spread of about 40–60 years. However, the absence of the longer time scales in the models that we have shown here provides a challenge for their predictive capability on decadal time scales. The origin of this lack of parity between satellite observations and the models will of course differ from model to model, but should provide an observational and statistical constraint for model physics, which has evidently not been met by the tuning and parametrization schema used thus far. Indeed, Swart *et al.* [32] showed how a deliberately-biased picking of short-term trends to predict the long-term behaviour of sea-ice extent can lead to false conclusions. Similarly, our study demonstrates both quantitatively and qualitatively the reason behind such high inter-model spread on decadal time scales and longer.

## 4. Conclusion

The observed Arctic sea-ice extent can be thought of as the ‘output’ of the highly complex nonlinear interactions governing the air/sea/ice system. It is a central goal of modelling to try and reproduce this output. A model can be tuned to match some of the observations, but if one is interested in making predictions about the future state of the system, one seeks to know the key underlying physical and dynamical processes and how they govern its central statistical response. Many studies compare the multi-model ensemble mean (MMEM) to observations, and argue that a sufficient condition for predictability is that the observations are within 1 s.d. of the MMEM. However, here we find that even the MMEM does not exhibit the robust observational features of white noise characteristics on annual–biannual time scales or the characteristic decadal time scales. Figure 4 compares the fluctuation functions for the EIE from the satellite observations (blue) to the MMEM. The longest crossover time scale present in the MMEM data is 1.8 years and that from the observations is 8.6 years. Therefore, simply summing different models to produce the MMEM does not provide the observed statistical structure of the ice pack.

A central role of observations is to produce robust statistical metrics that can serve as a target for modelling. Whereas, even with incorrect physics, it is possible to match features of the observations, capturing the appropriate statistical properties should provide important litmus tests for models. Our method is agnostic with regard to the number of crossover time scales, and hence associated processes, that may be present in the system. Therefore, we hope that the fact that we capture both the shortest and longest observed processes in the ice pack will provide a useful and robust test bed for modelling studies.

## Data accessibility

This article has no additional data.

## Competing interests

We declare we have no competing interests.

## Funding

The support of NASA grant NNH13ZDA001N-CRYO is acknowledged by both authors. J.S.W. acknowledges Swedish Research Council grant no. 638-2013-9243 and a Royal Society Wolfson Research Merit Award. Data for the CMIP5 models is publicly available and the sea-ice satellite observations are available from NSIDC.

## Appendix A

Below we provide additional figures 5–8 of the fluctuation functions, *F*_{q}(*s*), versus time scale, *s*, for the CMIP5 models labelled therein and described in table 1. Changes in the slope of *F*_{q}(*s*) occur at the particular crossover time scales denoted in the caption. The main text describes the dynamics associated with these time scales.

## Footnotes

One contribution of 15 to a theme issue ‘Modelling of sea-ice phenomena’.

- Accepted June 22, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.