A simple decomposition of European temperature variability capturing the variance from days to a decade

We analyze European temperature variability from station data with the method of detrended fluctuation analysis. This method is known to give a scaling exponent indicating long range correlations in time for temperature anomalies. However, by a more careful look at the fluctuation function we are able to explain the emergent scaling behaviour by short time relaxation, the yearly cycle and one additional process. It turns out that for many stations this interannual variability is an oscillatory mode with a period length of approximately 7–8 years, which is consistent with results of other methods. We discuss the spatial patterns in all parameters and validate the finding of the 7–8 year period by comparing stations with and without this mode.


Introduction
Variability of temperature time series has been described in two ways. Some people concentrate on scaling behavior of correlations and its seeming power law behaviour (long range correlations). Popular methods are detrended fluctuation analysis (DFA) (Fraedrich and Blender 2003), and multifractal versions (Moon et al. 2018), power spectral densities (Fredriksen and Rypdal 2016) and wavelet analysis (Lovejoy 2018). Such results consider the complete power of climate fluctuations. They are useful for estimations of prediction errors (Massah and Kantz 2016;Ludescher et al. 2016) and the validation of more complicated general circulation models (Govindan et al. 2002).
Other projects have tried to identify characteristic timescales in the data (Ghil et al. 2002). Here people use for example singular spectrum analysis (SSA) (Plaut et al. 1995), versions of Monte Carlo SSA (Paluš and Novotná 2004), and again wavelet methods (Baliunas et al. 1997). Results of these methods also help to validate general circulation models and contribute to the physical understanding of the atmosphere as well as to make predictions if the component is significant like the El Nino-Southern Oscillation (Chekroun et al. 2011). They can explain observed slowdown (or fastening) of global warming (Steinman et al. 2015;Hu and Fedorov 2017).
The disadvantage of the first approach is that it neglects properties that are already known and reduces the information to one property which is difficult to interpret (Maraun et al. 2004). The disadvantage of the second approach is that it does not take the full power of the fluctuations into account (Lovejoy 2015).
While parts of the literature concentrating on characteristic timescales completely ignore the findings of power law scaling, others mention the so called continuum variability (Huybers and Curry 2006) and treat it as one of the features contained in temperature fluctuations. Recently interest in this continuum variability is growing again (Lovejoy 2018;Fredriksen and Rypdal 2017). Lovejoy (2015) suggest to filter the known frequency modes, just like it is normally done with the seasonal cycle, in order to understand the nature of the continuum better. Our approach in fact follows this spirit. We will apply a novel method based on DFA that in contrast to traditional usage of DFA concentrates on the characteristic timescales. Therefore we show connections between the two approaches mentioned at the beginning. However, we conclude that all the power of the fluctuations can be explained by a superposition of few simple processes for the time range under consideration.
We want to concentrate on climate data from Europe where a large number of stations with several decades of daily recordings exists. For this region there is already a long history of investigations on both individual stations or grid data. Most importantly, people have repeatedly reported a 7-8 year cycle in various stations and other records (Plaut et al. 1995;Paluš and Novotná 1998;Sen and Ogrin 2016), the origin of which is still unclear. Spacial patterns were described in Grieser et al. (2002), Pišoft et al. (2009) andJajcay et al. (2016), however, to our knowledge there has not been an algorithmic investigation on a similar number of stations as we perform it.
We present the foundations of our method in Sect. 3, summarizing previous work in Meyer and Kantz (2019). In Sect. 4 we postulate our model for weather and macroweather variability in Europe and extend our method to processes with more than one characteristic timescale. In Sect. 5 we show results of a large scale algorithmic application of our method to station data of European temperatures. We are able to reproduce the 7-8 year cycle in part of this data. We show that our method is able to separate data with and without this frequency mode.

Data
We analyze station data of temperature and pressure in Europe (Tank et al. 2002). It can be downloaded from http:// www.ecad.com. The resolution of the time series is one day. The datasets strongly vary in length. We restrict ourselves to long time series as described in Appendix 1.

The method
Our method is based on detrended fluctuation analysis (DFA) (Peng et al. 1994) (see Appendix 2 for implementation) and its theoretical understanding (Hoell and Kantz 2015b). DFA produces the so called fluctuation function F of the time series, which is in fact a nonlinear transformation of the autocorrelation function C Here 2 x is the variance of the time series and L q is a kernel described in Appendix 3. It turns out, that looking at F(s) has two advantages over looking at C(t): on the one hand it is numerically more stable and on the other hand DFA is able to remove polynomial trends from the signal, so they do not contribute to the fluctuation function (Höll et al. 2019). We only use DFA1, because the trend due to climate change is sufficiently small compared to the fluctuations (Massah and Kantz 2016). Meyer and Kantz (2019) have shown that this method has the ability to uncover characteristic timescales in time series. The idea is to fit the theoretical fluctuation function of autoregressive models of different orders: AR(1) models for relaxation times, and AR(2) for noisy oscillations. The method works well for data that is dominated by one characteristic timescale, like approximations to yearly global mean temperatures (Meyer et al. 2018) or an 11-month averaged El Nino signal. We show the theoretical fluctuation functions in (16) and (23) in Appendix 3.
As a first example we want to show the application of the method to the time series of air pressure measured at Potsdam, Germany. We deseasonalize both the pressure and its variance in order to get a homogeneous time series. Therefore we calculate the long time climatological average pressure for each calender day and subtract it from each day. We also divide each value by the average variance for the calender day.
In the long time limit F(s) scales like 1/2, implying short range correlations. Figure 1 shows that the data can be described by an AR(1) model in a fairly good approximation. Obviously, for a complex system like the atmosphere, we would not claim that pressure is purely autoregressive and no other effects happen at characteristic timescales. However, we can see that those other timescales are not significant for a smoothed measure such as the DFA-fluctuation function. The method shows us a data model that can be used for stochastic predictions and qualitative understanding of the dynamics.
In Fig. 1 we also show temperatures and temperature anomalies (T with climatological average temperatures for each calender day subtracted) from Potsdam. In DFA the temperature anomalies exhibit anomalous scaling proportional to s 0.65 , which could be described by a model with long range correlations (Massah and Kantz 2016). In this article we want to show a different approach. Instead of looking at the anomalies, we look at F T itself. We want to identify characteristic timescales and therefore extend our understanding beyond the description of emergent scaling. If the time series is indeed a superposition of several independent processes we can make use of the superposition principle of DFA (Hu et al. 2001) This will be important for our purpose in the next section.

Decomposition of one time series
In the examples above fitting fluctuation functions was exclusively applied to systems with one dominant characteristic timescale. What are we supposed to do if the fluctuation function implies a more complex dynamics? In this section, by analysing historic temperature data from Potsdam station, we show that our method still works for systems with more than one characteristic timescale if these timescales are sufficiently well separated. We assume that the time series can be written as a superposition of three processes namely, a short range part X t , the seasonal cycle Y t and one unspecified component, which describes interseasonal variability. We will see that these components are clearly separated in time, such that Y t and Z t might as well be interpreted as additional components of the noise input t of X t . We will describe X t and Z t by simple autoregressive models, while the seasonal cycle Y t is a regular oscillation.
Our model contains some simplifications. This is not only because it only consists of three processes described by few parameters, but also our assumption of independence of the three processes. This is not always true as shown in Paluš (2014) and Jajcay et al. (2016). However our focus in this text are not the cross-scale interactions, but only the observation of the different timescales. We also do not discuss spatial correlations.
Temperature correlations for short times decay exponentially (Massah and Kantz 2016). The behaviour is apparent in Fig. 2 in the fluctuation function F T (s) for small s. Other subseasonal timescales are not pronounced in the fluctuation function. This lack of characteristic timescales is known and the reason why predictions on these timescale are so difficult (Ghil et al. 2019). Hence X t is best described by an autoregressive model of order one AR (1) where the noise term t X is the output of some process on a much smaller timescale. The relaxation time r is determined by the AR-parameter c. The shape of the distribution of t X is not important for our model (Meyer and Kantz 2019), however, it seems to be well approximated by a Gaussian distribution. The variance 2 is a free parameter. For simplicity we will later look at the standard deviation X of the process instead of . They are related by Eq. (18).
(4) X t X = cX t X −1 + t X and r = −1∕ log(c), Due to the temporal separation of the three processes X t , Y t and Z t , and because autoregressive models are low pass filters, the impact of the slower processes can be neglected for short times. For Potsdam the fit of Eq. (16) to F T for small s (see Fig. 2) yields a relaxation time of 3.8 days. We subtract the fitted AR(1) fluctuation function F 2 X from F 2 T . Since we now covered the fluctuations for short times we change the timescale to 10 days by dividing F T−X and s by 10 and remove small values of s. The by far most pronounced mode in the fluctuation function is the seasonal cycle (Deng et al. 2018). Although it is not perfectly periodic (Paluš et al. 2005) it has a phase-locked frequency of 1/365.25 days, which is enforced by the periodic driving of the sun. We do not remove the cycle by subtracting the climatological average since this popular procedure does not account for the whole dynamics on the yearly timescale. It does not take into account that the phase of the oscillation fluctuates. Instead we approximate the fluctuation function F Y by an AR(2) process (see Eq. (6)) with b = −1 , which is equivalent to a sinusoidal signal which does not describe the phase fluctuations, but accounts for their full power. Figure 2 shows that this is a good approximation to the remaining fluctuation function. Again we subtract F 2 Y from F 2 T−X , change the timescale to half a year (182.625 days) by dividing both s and F by 18.2625 and remove points with small s.
So what is left if the first two components are removed from the fluctuation function? Plaut et al. (1995) claim to have found an 7-8 year cycle in the record of central England temperatures. Since then this mode was found in several other time series of European temperatures (Paluš and Novotná 1998;Pišoft et al. 2009;Sen and Ogrin 2016). We want to model interseasonal fluctuations with an AR(2) model that-dependent on the parameters a, b-can exhibit oscillatory behaviour with period (see Eq. (21)) or pure relaxations. Again we treat the standard deviation Z t of the process as a free parameter, which determines of the noise according to Eq. (22) For Potsdam and other stations we looked at, our fit seems to show good agreement with F 2 T−X−Y , although the fluctuation function is not as smooth as for shorter times. We obtain a period of 7.5 years. The complete fluctuation function F 2 X+Y+Z is shown in Fig. 1(right panel). It shows excellent agreement with the measured F 2 T . It has been suggested that climate can be described by a continuum variability with anomalous scaling (Lovejoy 2015) and quasi-periodic perturbations. For the time-range under consideration, which is admittedly shorter than in most studies that concentrate on this issue, we do not see any indicators that such a model would be necessary. The scale invariance of the anomalies visible in Fig. 1(right) is explained in our model by superposition of white noise driving at each identified timescale. This is because AR(2) is a low pass filter which 'hides' the fluctuations of its driver for short times. Several AR(2) processes with different periods therefore lead to a successive increase of the background noise for longer time. This emergent scaling is an alternative interpretation to dynamical long range correlations, the origin of which is still not fully understood (Fraedrich et al. 2004;Fredriksen and Rypdal 2017;Meyer and Kantz 2017).

Analysis for European station data
The decomposition of the fluctuation function as described in the previous section was incorporated into an algorithm and applied to temperature data of stations around Europe. We obtain six parameters, r, X , Y , a, b, and Z , for each station, which model the fluctuations of timescales up to a bit more than a decade, describing weather and macroweather. The aim is to investigate the spatial dependence of these parameters and especially for which regions the AR(2) model (6) indicates an oscillatory mode on the timescale of several years. The period time of this oscillation can be calculated from the parameters a and b.
The short time dynamics is described by the relaxation time r. We obtain it by fitting F X to F T . The results are displayed in Fig. 3. r shows clear regional patterns, which are much stronger than the uncertainty of our method. The correlation time is short in the west and the north of Europe close to the Atlantic. It is longer in central and eastern Europe and around the Baltic sea. For most stations r is just a bit smaller than 4 days.
The dynamics of pressure data for short times turns out to be similar to that of temperature. Both describe the typical timescale of weather patterns. However, there seems to be no connection between r T and r p for one specific station. On the contrary the regional patterns of pressure strongly differ from those of temperature. Pressure values relax faster in the south compared to the north.
The variance of the noise in the short time dynamics is very low close to the coasts and has its highest values in Russia and Scandinavia. The variances of all three processes X, Y and Z are displayed in Fig. 4. The amplitude Y of seasonality is strongly correlated with X .
When fitting F Z of the AR(2) model to F T−X−Y we only distinguish two cases. Either the model has an oscillatory mode with a period between 10 and 22 half years-or not. The set − of stations without this mode consists of three subsets: • if the slope is ∝ s 0.5 or ever more flat the algorithm indicates an oscillation with a very short period • if the slope is very steep ( ∝ s 1.5 ) the algorithm indicates an oscillation with a very long period • for stations in between the AR(2) model has real roots indicating no oscillation.
All three cases have in common, that there is no sharp crossover in the fitted s-range and therefore statements about periods can only be made with high uncertainty. Examples are shown in Appendix 1. For 51% of the analyzed stations we observe an oscillation, see Fig. 5. We denote the set of these stations by + . The categorization into the two subgroups uncovers regional patterns with few exceptions which might be interpreted as false classifications. The frequency mode is significant in an area including England, Belgium, the south of Scandinavia, central Europe north of the Alps and parts of eastern Europe. The values for that we obtain do not show clear regional patterns, which indicates that we indeed only see the known 7-8 year cycle and the deviations are the uncertainty of our method. On average the period length is 7.6 years with large standard deviation of 1.8 years.
We test the ability of our method to distinguish between stations with and without an observable period by looking at the conditionally averaged power spectrum where F is the Fourier transform and the brackets ⟨⟩ denote the average over all stations in the set + or the average over all stations in the set − . The validations is successful as the stations where the period was found show a clear peak at 8 years while the others do not (see Fig. 5 right).

Conclusions
We introduce a method for the detection of characteristic timescales in time series based on detrended fluctuation analysis that works for data where the dominant timescales are sufficiently well separated. It is an intuitive tool for the derivation of approximate dynamics that accounts for the full power of the fluctuations in the system. In this way we bridge the gap between spectral methods that are applied for the detection of oscillations and scaling methods that detect colored noise. The method does not detect nonlinear effects like cross-scale interactions but only covers the characteristic timescales individually.
While pressure data usually only shows exponential decay of correlations, temperature data shows a much slower decay, often interpreted as long range correlations. We model these time series with a simple model, a superposition of short range decay, a seasonal cycle and a potentially oscillating noisy component for longer times. This model is able to describe the observed fluctuation functions accurately. DFA is a smoothing method that is not designed to find every tiny structure in the time series. For timescales up to a decade it is not necessary to claim a colored background noise, as many authors do, who look at longer timescales. The increase of the background noise can be explained by the input noise of the interseasonal variability. By fitting the parameters we obtain the power of each component for each station. more importantly we are able to reproduce a previously found 7-8 period. Our method is with only few false classifications able to distinguish between stations where this component is significant, and stations, where it is not.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Technical details
We restrict ourselves to datasets which cover at least 60 years of (non-blended) daily temperature measurements and have at most 6.7% of missing days. We only consider one dataset for each station in cases where there are more. In an area around Germany with very good coverage we even require 75 years recordings. We use pressure datasets from the same stations in Fig. 2 if there are at least 4 years recorded. For the analysis of the anomalies in Fig. 1 we remove the 29th of February in leap years. In our main analysis, this is not necessary, since we are dealing with the real values. In total 336 stations were considered, out of which 171 belong to + , i.e. exhibit the interseasonal oscillation. The fluctuation function is calculated for a set S of s-values, constructed from the maximal value s max = 365 × 75 and the recursion s i = sup s∈ℕ s ≤ 0.9s i+1 , which yields a logarithmic scale. For datasets with less than 4 × s max recorded days the values of s which are larger than 1/4th of the length of the time series are ignored.
The short time behaviour of the data (4) is fitted for 3 ≤ s ≤ 25 days. For the analysis in Fig. 3, the pressure fluctuation function is also fitted for s ≤ 25 only, unlike in Fig. 1, where it is fitted for the whole fluctuation function as displayed. Figure 2 shows the full procedure for temperatures. The seasonal cycle (5) is fitted for 8 < s < 38 times 10 days. The AR(2) model (6) is fitted for 8 < s times 6 months. In each case fitting is done by minimizing the variance The standard deviation as well as the amplitude of the seasonal cycle is not fitted but rather calculated as In Fig. 6 we show the decomposition of F 2 T for stations in Arkhangelsk, Bologna, and Oslo Blindern. While Oslo Blindern is another example for a station where an oscillation could be identified with a measured value of = 8.2 years, the others are examples for stations where the frequency mode is not strong enough to be significant. For Arkhangelsk our algorithm yields AR(1)-like interseasonal dynamics with a relaxation time of approximately 2.5 years. For Bologna it yields a slow oscillation due to the steep increase of the fluctuation function. The measured period length is 39 years, which is outside our investigated time range and therefore not a reliable result.

Appendix 2: Detrended fluctuation analysis
Here we recall the basics of DFA (Peng et al. 1994;Hoell and Kantz 2015a). Given a time series {x n } N n=1 we first calculate the integral Then we divide the time axis into K non-overlapping segments of length s and calculate the so-called DFA variance x j . f 2 ( , s) in every segment which is given by the squared error sum of y t and a fitting polynomial p t of order q. This order gives the name of the method DFAq We repeat the procedure for an other K non-overlapping segments of length s, but this time starting from the end of the series. Finally, the square of the DFA fluctuation function is the average of all the DFA variances over all segments Traditionally people look at the asymptotic behaviour F 2 (s) ∼ s 2 , which contains information about the correlation structure. For = 1∕2 the process is short range correlated and for > 1∕2 the process is long range correlated.

Appendix 3: Fluctuation functions of AR(1) and AR(2)
Theoretical fluctuation functions can be calculated from Eq.
(1) as proposed in Hoell and Kantz (2015b) and performed in Hoell and Kantz (2015a) and Meyer and Kantz (2019). The kernel L 1 in DFA-1 is blue 'x'), however, here we show the full decomposition (all three steps) in one panel. The fits are drawn in the same color as the data they a fitted to. In order to increase clarity the fluctuation functions are not rescaled after each step. The dashed line in each panel represents the result F 2 X+Y+Z of the decomposition. Left: Arkhangelsk, where F 2 Z is best described by an AR(1) model; Center: Bologna, where F 2 T−X−Y is very steep, which means that the next oscillatory mode is outside the presented s-range; right: Oslo Blindern, where an oscillatory mode with period length = 8.2 years was detected The correlation function of AR(1) is known to be with AR(1) parameter c. The fluctuation function for AR(1) can therefore be calculated as with J c (s) , K c (s) polynomials in s The variance of the process is given by where 2 is the variance of the noise.