Power-law behavior in millennium climate simulations

Introduction Conclusions References


Introduction
Power laws for spectra and frequency distributions are commonly found in complex systems from turbulent flows and intensities of earthquakes to some economic variables (Clauset et al., 2009) and are usually intuitively characterised as the system lacking a characteristic length or time scale in the range, where the power law is valid.In the context of climate data, power laws have been found in a wide range of timescales (Bloomfield , 1992;Pelletier, 1998;Wunsch, 2003;Fraedrich and Blender, 2003;Gil-Alana, 2005;Rhines and Huybers, 2011).Huybers and Curry (2006) analysed temperature variability at all timescales from monthly variations to millenial and longer timescales and found different exponents for best-fit power laws in the 1.1-100 yr (smaller βs) and 100-15 000 yr (larger βs) ranges as well as a latitude-dependence Figures of the exponent when analysing the NCEP reanalysis data for periods between two months and 30 yr.Blender et al. (2006) used a long simulation with a general circulation model and ice core data to reach power law fits near Greenland and elsewhere.Blender and Fraedrich (2003) analysed reanalysis data and global warming simulations and in a fit S(f ) ∼ f −β quantified β = 1 over the oceans, β = 0 over the continents and β = 0.3 in coastal areas.More theoretical insight has also been given by, e.g. a simple ocean model (Fraedrich et al., 2004), considerations involving the Hurst exponent (Blender et al, 2006;Vyushin and Kushner, 2009;Vyushin et al., 2009) and by a random-walk framework (Stephenson, 2000).While good fits to power laws have been achieved for many datasets, some scatter around the fits remain.In this study, we utilise long climate simulations belonging to the Millenium experiment (Jungclaus et al., 2010; http://www.mpimet.mpg.de/en/science/internal-projects/millennium.html) to achieve statistical spectral estimates using long time series.

The Earth System Model and the spectral method used
We use the COSMOS network Earth System Model, consisting of the atmospheric model ECHAM5 (Roeckner et al., 2006), the ocean model MPI-OM (Marsland et al., 2003), the land use and vegetation model JSBACH (Raddatz et al., 2007) and the ocean biogeochemistry model HAMOCC (Wetzel et al., 2006).The PRISM/OASIS3 coupler (Valcke et al., 2003) interactively couples the atmospheric and ocean models.
The coupled atmosphere-ocean system in the model features multidecadal variability (Jungclaus et al., 2005) as well as El Nino type variability in the tropical Pacific (Jungclaus et al., 2006).El Nino frequencies tend to be lower than in observations (Jungclaus et al., 2006).The importance of the variability at these timescales is seen in the next section.The Millenium experiment (Jungclaus et al., 2010)  and a 1206-yr simulation with external forcing by volcanoes, solar radiation changes and anthropogenic greenhouse gases and aerosols included.The advantage of the unforced simulation in this context is that we can assume the internal variability to be homogeneous throughout the 1201 simulated years and thus obtain a large sample.Another simulation including the external forcings and thus being more realistic compares better with measurements, as will be shown later.We also briefly analyse simulations including only solar and only volcanic forcing.
The spectral method we used is a discrete Fourier transform (DFT) with varying starting point and length of time window.Previously, we used the same method to analyse the multidecadal (50-80 yr timescales) oscillation observed in measured global mean temperature and in the simulations we analyse here (Henriksson et al., 2012).In that paper we looked at absolute value of Fourier coefficient in order to evaluate the contribution of an oscillation to global mean temperature, but in this study we look at the power spectrum for consistence with earlier studies.By varying the starting point of the time window, we obtain several spectra for a given length of time window and can average over them.Obvious advantages for calculating averages over several spectra are seen in the next section.In Sect.6 in the context of Fig. 11, we also compare our method to an earlier result obtained by multitaper analysis, a method often encountered in climate literature (Yiou et al., 1996;Mann and Park, 1999).

Analysis of global mean temperature
In this section, we fit power laws to power spectra of global mean temperature: The global mean temperature timeseries consists of annual means (except in Fig. 4, where we consider the monthly mean timeseries).Figure 1 shows the power spectrum of the unforced simulation averaged over all possible power spectra from 300-yr DFTs (see Henriksson et al., 2012 for more plots and details).Spectral peaks can be seen at the multidecadal (∼50-80 yr) and El Nino (∼3-6 yr) timescales and smaller peaks at centennial (not shown) and ∼20 yr timescales.Figures

Back Close
Full In the following, we present plots of the logarithm of frequency vs. logarithm of the power spectrum, with the power spectrum calculated using different lengths of time window.We fit power laws for two ranges: the multidecadal to the El Nino timescale (reaching from maximum amplitude at multidecadal scales to minimum before El Nino timescales) and frequencies higher than El Nino (starting at the maximum amplitude at El Nino frequencies).The latter range is quite narrow, corresponding approximately to frequencies between 1/(3.2 yr) and 1/(2 yr).Next, in Fig. 2 we show spectra, where we have not utilised the possibility to average over different spectra and then, in Fig. 3, we show averaged spectra.
The logarithm of the power spectrum obtained by a DFT of the first 64 yr of the simulation is shown in Fig. 2a, together with power-law fits.The best fits are β = 0.24 and β = 7.9.Figure 2b shows the spectrum calculated using the full 1201-yr timeseries.
The best fits are now β = 0.13 and β = 8.6, respectively.The advantage of averaging over several spectra is seen in Fig. 3, showing the power spectrum obtained by averaging over all possible spectra with a time window of 300 yr and 64 yr.Random noise is averaged away and the underlying statistical spectral form is seen more clearly.For the spectrum in Fig. 3a using the 300-yr time window the best fits are β = 0.35 and β = 7.6 and for the spectrum in Fig. 3b using the 64-yr time window the best fits are β = 0.35 and β = 7.5.
When monthly mean values of global mean temperature are used instead of annual mean values, we obtain the result shown in Fig. 4. One observation that can be made is the large amplitude for the annual cycle and its harmonics.The power-law fits are not modified much, with β being 0.30 and 7.9, respectively, for the multidecadal to El Nino and El Nino to interannual frequency ranges and 0.87 for frequencies higher than 1/yr (annual cycle and two first harmonics ignored in fit).detrending is not done) in the multidecadal to El Nino range and 8.3 for higher-than-El Nino frequencies.Figure 5b shows the spectrum using a 64-yr time window (without detrending as the length of time window is so short) and here β is 0.71 and 7.4.β in the multidecadal to El Nino range is thus significantly higher than in the unforced simulation.Most of the difference comes from the low frequencies corresponding to multidecadal and somewhat shorter timescales as can be seen by comparing Figs. 3  and 5.Although we use arbitrary units, the units are the same in these two figures (for the respective subfigures a and b).The minimum between the multidecadal and El Nino frequencies is at the same level, but the multidecadal amplitudes are higher for the forced simulation.
When fitting power laws to a simulation with only volcanic forcing included, best-fit β in the multidecadal to El Nino range is 0.71 (64-yr time window; plot not shown).In a simulation with only solar forcing, best-fit β is 0.52 (64-yr time window; plot not shown), though in this case the fit is not quite as good.Thus it seems that volcanic forcing has a larger effect on β than solar forcing when applied as sole external forcing, but both increase the amplitudes of low frequencies.
These results indicate that the multidecadal and El Nino timescales are the most important when considering statistical properties of spectra on interannual to centennial timescales.As the largest differences between the spectra of the unforced and forced simulations are in the low-frequency range, it seems that the external forcings feed the low frequencies.One question that can be raised is why the spectrum in the multidecadal to the El Nino frequency range has a power-law shape in both cases.It is interesting that internal climate dynamics alone seem to produce a power law, and a much less steep one than that for Milankovitch cycle frequency ranges (Huybers and Curry, 2006).
We can hypothesise that the spectral form has to do with the probability distribution of temperature fluctuations, but we have not managed to develop a satisfactory theoretical explanation along this line as yet.In any case, the spectrum is self-similar in this range.We also looked at spectral dynamics (results not shown here correlations to look for evidence that amplitudes are fed from the multidecadal and El Nino frequencies to higher frequencies, but we did not find statistically significant support for this idea (although weak indication for it existed).Although some earlier studies have been able to explain variability of β in different cases or regions to some extent, we are not aware of a theoretical explanation of the power-law shape of the spectrum in the literature.The question would thus be interesting to study further.

Error and sensitivity analysis
As a sensitivity analysis, we calculate the best fit βs, together with 95 % confidence intervals for different lengths of time window and for different window functions.We have chosen a variety of lengths for the time window that are not integer multiples of each other.We do this for the unforced simulation.The results are shown in Table 1.
As can be seen, all but one of the fits for the multidecadal to El Nino range are nicely within each others error bars, with estimates of beta ranging from 0.34 to 0.39.The single exception is that for the 1201-yr rectangular time window, beta is lower.Also in the frequency range above El Nino frequencies the fits are close to each other, mostly lying in the narrow range of 7.5-8, but indeed the confidence intervals do not overlap.
For the forced simulation results, shown in Fig. 5a Fig. 3a, which would be interesting to study further not the least from a physical point of view, these deviations largely disappear when looking at the spectral range with coarser resolution in Fig. 3b. Figure 6a also shows an important part of the reason behind the scatter within β = 7.5-8.The spectrum has some curvature at the end and also the lowest frequencies included are above the power-law fit as we have defined the fit to start from the El Nino frequency with maximum amplitude.However, we decide to keep the fit for the whole range from El Nino frequency with maximum amplitude to the Nyquist frequency as the distortion is small and the fitted range already narrow as it is, acknowledging that in this range the power law more is merely a numerical fit rather than a fundamental property of the climate system.The curving tail could maybe also be modeled, but it might simply turn out to be mixing with the spectral form of the next higher frequencies shown for instance in Fig. 4.

Local and regional temperatures
In this section, we investigate first the spectra of two well-known modes of internal variability: El Nino and the Atlantic Multidecadal Oscillation (AMO), represented by the mean sea surface temperatures in the Nino 3 region and the North Atlantic (between 0-75 • W and 0-60 • N) in the unforced simulation.We then proceed to study the spatial distribution of β when the fitting is done gridpointwise.The latter is done for the simulation including all external forcings to allow for a better comparison with (Blender and Fraedrich, 2003) using forced simulations and reanalysis data.
Figure 7 shows the spectra of the Nino 3 region and the AMO in the unforced simulation, using a 300-yr time window.The El Nino spectrum in Fig. 7a shows a maximum at frequency ∼1/(3.2yr), decaying thereafter with β = 8.9 ± 0.02.Thus β is larger than in the global average, simply because at ENSO frequencies, temperatures in the Nino 3 vary much more than the global-mean temperature.The AMO spectrum in Fig. 7b shows high amplitudes in the centennial and multidecadal timescales, and local maxima at ∼30, ∼8 and ∼3 yr.Thus a power law does not fit well to any Introduction

Conclusions References
Tables Figures

Back Close
Full frequency range for the AMO, except maybe at the highest frequencies, for which we fitted β = 3.41 ± 0.06. Figure 8 shows the distribution of best-fit β when done gridpointwise for the forced simulation data.As the spectral shape varies regionally, we cannot split up the spectra into frequency ranges and have therefore decided to make a Fourier transform on the full 1206-yr annual mean temperature record and to fit a power law to the full spectrum, which is then varied in Fig. 9.We use a Hann window based on the considerations in Sect. 4. The patterns seen in Fig. 8 are consistent with Blender and Fraedrich (2003), giving βs close to 1 over oceans and close to 0 over land areas, except for the tropical Pacific Ocean with surroundings, where the strong El Nino causes β even to be positive (which can be understood by looking at the previous figure).In Fig. 9, where we restrict the fit to frequencies up until 1/(6 yr), we can see that the βs over ocean generally decrease and the βs over land increase as compared to Fig. 8. Thus the frequency range between 1/(6 yr) and 1/(2 yr) is responsible for a large part of the difference between results over oceans and land in our results.At frequencies higher than 1/(6 yr), over land there is more variability in relation to lower frequencies than in areas over the ocean.These results show the importance of analysing the spectrum piecewise as was done in Sect.3, especially since in a least-squares fit procedure each point receives equal weight and here the interval between 1/(6 yr) and 1/(2 yr) contains two thirds of the resolved frequencies.When considering frequencies below 1/(6 yr), the ocean and land areas seem to exhibit more similar behavior in terms of best fits to power laws, with some regional variation nonetheless.
As a detail from the maps over oceans, we see a local maximum in estimated β south of Greenland as was reported and shown in the CSIRO model results in Fig. 2a and b of Blender et al. (2006)

Comparison to measurements
In this section, we fit power laws to available measurement time series.Although the fits are not as good as for the long simulations, they serve to compare simulation and measurement results and in the case of the Central England temperature (CET), to compare our results with previous results obtained by multitaper analysis.
Figure 10 shows the spectrum of the CRU global mean temperature for the years 1850 to 2011 (Gil-Alana, 2005), detrended quadratically and spectrum averaged over transforms with 132-yr time windows.The El Nino amplitudes are not as strong as in the simulation and at least with the time series of available length, the second power-law range is not distinguishable.However, for multidecadal to El Nino timescales, the fitted β = 0.86 ± 0.09 corresponds quite well to that fitted for the forced simulation (Fig. 5).
Figure 11 shows the spectrum of the CET measurement record (Parker et al., 1992) computed both from measurements and model results (simulation with all external forcings included; as simulation data for the CET we use the gridpoint at 50.1 11a shows the measurement spectrum together with a fit stretching from the multidecadal time scale to the Nyquist frequency with β = 0.35 ± 0.05, very close to what Huybers and Curry presented.This raises confidence in our method as compared to the multitaper method used by Huybers and Curry and many others.Figure 11b shows the same spectrum calculated as an average from the forced simulation.The result β = 0.34 ± 0.02 is consistent with the result of Huybers and Curry.This gives more confidence in the model's ability to simulate regional temperature variations as Central England's β deviates a lot from β for global mean temperature.

Conclusions
In the present study, we applied the method of discete Fourier transform with varying starting point and length of time window to analyze spectra of global-and regional- between multidecadal and El Nino timescales, the spectra for global-mean temperature follow closely a power law: S(f ) ∼ f −β , where β ∼ 0.35-0.4 in the unforced simulation, and β ∼ 0.7-0.8 in the simulation including all external forcings.Volcanic forcing had a larger effect than solar forcing, when applied as the only external forcing in a simulation.
It seems that the external forcings feed the low frequencies more than the high frequencies, which is perhaps not so surprising, but why the power-law form is present in both simulations remains unclear to us.It seems that internal dynamics alone are able to produce a power law, which is less steep than when external forcing is included, which in turn is less steep than power laws obtained in previous studies for Milankovitch cycle time scales.
Comparisons with measurements were reassuring: the power-law fit for measured global mean temperature gave β = 0.86, relatively close to the value in the simulation including external forcing.Similarly, the power-law fit to the simulated Central England temperature spectrum was very close with both measurements and an earlier result obtained by multitaper analysis, also raising confidence in the method we used.
In locally fitted power spectra we made an observation that if a power law is fitted to the full spectrum, there is a clear difference in β over land areas as compared to ocean areas, but that this difference almost disappears if fitting the power law only up until frequency 1/(6 yr).Thus the power-law behavior was similar for these frequencies and only in the highest resolved frequencies was there more variability over land than over the ocean in relation to the lower frequencies.
In all, we managed to achieve good power-law fits to spectra in millenium climate simulations due to the long time series they provide and the fact that they statistically seem to obey these laws.The results were also validated to important parts against measurements and results obtained with other spectral methods.We hope our study will inspire further research in understanding the origin of the spectral form as well as the underlying climate variability.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full Discussion Paper | Discussion Paper | Discussion Paper | consisting of millenniumscale simulations was made to allow for a wide range of studies regarding the climate of the last millenium.We analyse a 1201-yr simulation without external climate forcing Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 5
Figure 5 shows the spectrum from the simulation including all external forcings.Figure 5a shows the spectrum obtained by doing the transform after quadratic detrending for each 300-yr sequence separately to avoid mixing of contributions from longerterm externally-forced trends and oscillatory components.Best-fit β is 0.79 (0.80 if ) and time-lagged Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | , the error bars for β are 0.79 ± 0.06 (0.80 ± 0.06 without detrending) and 8.3 ± 0.01 for the two ranges, and in Fig.5b, they are 0.71 ± 0.11 and 7.4 ± 0.04.Additional checks have been done with Kaiser-Bessel windows (not shown), and results when using them are consistent with results from using the other window functions.The error term in the regressions for the 64-yr time window spectrum from Fig.3bis plotted in Fig.6a.The error term is typical of regressions using any length of time window or window function.It shows what seems like uncorrelated and random scatter in the multidecadal to El Nino range, which is supported by Fig. 6b showing the histogram of the error term from Fig. 3a resembling a Gaussian distribution.While some clustering of amplitudes above the fit line can be seen at periods 20-30 yr and 10-15 yr in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and supported by ice core measurements.In our results, there is a similar maximum also in the northernmost part of the North Pacific gyre.Discussion Paper | Discussion Paper | Discussion Paper | mean temperatures in 1200-yr climate simulations.It is found that for frequencies Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Table 1 .
Power-law estimates for the two frequency ranges using different window functions and different lengths of time window.Introduction