Distribution of Air Temperature Multifractal Characteristics Over Greece

In this study, Multifractal Detrended Fluctuation Analysis (MF-DFA) is applied to daily temperature time series (mean, maximum and minimum values) from 22 Greek meteorological stations with the purpose of examining firstly their scaling behavior and then checking if there are any differences in their multifractal characteristics. The results showed that the behavior is the same at almost all stations, i.e., time series are positive long-term correlated and their multifractal structure is insensitive to local fluctuations with large magnitude. Moreover, this study deals with the spatial distribution of the main characteristics of multifractal (singularity) spectrum: the dominant Hurst exponent, the width of the spectrum, the asymmetry and the truncation type of the spectrum. The spatial distributions are discussed in terms of possible effects from various climatic features. In general, local atmospheric circulation and weather conditions are found to affect the shape of the spectrum and the corresponding spatial distributions. Furthermore, the intercorrelation of the main multifractal spectrum parameters resulted in a well-defined group of stations sharing similar multifractal characteristics. The results indicate the usefulness of the non-linear analysis in climate research due to the complex interactions among the natural processes.


Introduction
Complex systems consist of many components that interact with each other in very complicated ways, which are determined by nonlinear laws.A characteristic example of a complex system is the atmosphere and the natural processes that take place in it.The result of such complex interaction among atmospheric processes is the fluctuation of the values of meteorological parameters in an almost random way at many scales in time and space.These nonlinear processes are described by nonlinear partial differential equations [1].Nonlinearity results in linear data analysis being insufficient for a complete analysis of meteorological time series.Besides atmospheric sciences, nonlinearity also occurs in many other fields of science such as economics, biology, physiology, electronics etc. [2].Nowadays, the close cooperation of scientists from various disciplines has helped decisively in the development of new effective methods for the analysis of nonlinear phenomena.Researchers can now analyze processes that obey nonlinearity and they are able to reveal some properties, which could not be detected by linear methods.One such nonlinear method is Detrended Fluctuation Analysis (DFA), which was developed by [3].It should be mentioned that the fluctuations of time series from complex systems obey scaling laws for a broad range of time (and space) scales.By applying DFA, we can reveal the scaling properties and determine long-range correlations, even in nonstationary time series [4].This could not be done by using ordinary methods, such as spectrum analysis.Therefore, researchers from several fields of science applied DFA to their data, taking advantage of its benefits.Thus, DFA has been used in economics [5], stock markets [6], biomedical signals [7], DNA analysis [8], social and natural phenomena [9][10][11], etc.In the atmospheric sciences DFA has also been applied widely to time series of air temperature [12][13][14][15][16], relative humidity [17,18], precipitation [19], drought and flood index [20], ozone [21][22][23][24] and many others.In cases of highly non-linear time series with periodic trends, [25] proposed a variation of DFA (i.e., detrended cross-correlation fluctuation analysis).
Very often, the scaling behavior of a time series cannot be described by only one scaling exponent at all scales in time and space, but a variety of scaling exponents is needed to completely describe the time series.In other words, the time series does not have the same fractal structure at all scales, so we say that this time series is a multifractal signal.In that case, DFA is inadequate to analyze the signal and a generalization of DFA is used, the Multifractal Detrended Fluctuation Analysis (MF-DFA), which was introduced by [26].This method is used in this study and it is described in the next section.MF-DFA has also been applied to various scientific fields.In the field of atmospheric sciences, there are numerous studies which used MF-DFA to analyze time series of meteorological and air quality parameters.MF-DFA has been applied to time series of temperature [27], precipitation [28,29], wind speed [30,31], multiple meteorological parameters [32,33], particulate matter concentrations [34] and many others.Moreover, MF-DFA was used to study the sunspot number fluctuations [35], to investigate the multifractality in the residuals of the connectivity time series of a wind speed monitoring network and to determine how the multifractality degree is related to the correlation threshold [36] or how the spatio-temporal aggregation of the data affects its multifractal properties [37,38].
Greece covers a relatively small geographical area, but with wide climatic variety.For instance, the climate at Western Macedonia in northern Greece is quite different from the climate at Crete in the southernmost area of Greece.Within this context it would be interesting to examine the distribution of multifractality properties in the Greek region.A convenient way to examine multifractality is by using the multifractal spectrum and, more specifically, its basic characteristics such as the spectrum width and the value of the Hölder exponent for which the spectrum takes its maximum value and the asymmetry along with the spectra truncation type.The meaning of these quantities is explained in detail in the methodology section.The application of MF-DFA to temperature time series along with the spatial distribution of the basic characteristics of the resulting spectrums is examined quantitatively and qualitatively in terms of their association with climatic features.In the second section of this paper, the MF-DFA method along with the available air temperature data are presented, while the third section deals with the results and the relevant discussion of the main findings.

Experimental Data
In this study, we used data from 22 meteorological stations of the Hellenic National Meteorological Service (HNMS) network.These stations are listed in Table 1, while their spatial distribution is illustrated in Figure 1.The main criteria for the selection of these stations are: a) the length of the time series, b) the time series completeness (i.e., having very few missing values) and c) the spatial coverage (effort was made to cover all climatic and topographic areas of Greece).The experimental data used in this study are the mean, maximum and minimum daily temperature time series (Tmean, Tmax and Tmin, respectively) at 2m above the ground level for the stations and the relevant periods listed in Table 1.The time series are extracted from the Global Summary Of the Day (GSOD) NOAA database [39].In this database, quality assurance checks are    The experimental data used in this study are the mean, maximum and minimum daily temperature time series (T mean , T max and T min , respectively) at 2m above the ground level for the stations and the relevant periods listed in Table 1.The time series are extracted from the Global Summary Of the Day (GSOD) NOAA database [39].In this database, quality assurance checks are performed regularly.In brief, a series of quality tests are applied to the dataset which can detect a variety of data problems, such as inconsistent climatological values with the location of the station, existence of unreal values, comparison with neighboring stations' values and many others (more details can be found in [40]).Furthermore, a manual quality assurance check is performed to increase the reliability of the dataset.
It should be noted that data completeness of all stations is above 95%.The level of data completeness is listed in Table 1.The few missing values are isolated and they are replaced using interpolation of their non-missing neighboring values.This very high data completeness is well over the minimum level of 50% at which the scaling exponent and therefore the scaling behavior of the time series are not seriously affected [41,42].Therefore, MF-DFA can be reliably applied to the time series.
Additionally, daily data of wind speed, relative humidity, pressure and precipitation amount (for the same stations and periods) that are used later in this study have also been extracted from the same database.

Methodology
This study consists of the following steps:

•
The multifractal characteristics of T mean , T max and T min are studied for each station in terms of MF-DFA analysis.

•
Assessment of the spatial distribution of the main multifractal spectrum characteristics is performed.

•
Intercorrelations between multifractal spectrum parameters are examined.
MF-DFA analysis is applied to the temperature time series in order to examine their scaling properties.This method is described in detail by [26].Herein, we give a brief description of the method.

1.
First, the time series has to be deseasonalized.It is obvious that daily temperature time series exhibit periodical trends, which are attributed to the annual seasonal cycle.In general, periodical trends have an influence in the nonlinear properties of the time series [38] and therefore the time series should be deseasonalized before applying the MF-DFA method.An efficient method to deseasonalize a time series is the Seasonal and Trend decomposition using Loess (STL) method, which was introduced by [43].In the STL method, the time series is decomposed into seasonal, trend and remainder components.From the decomposed time series, seasonality is removed and MF-DFA analysis is performed on the deseasonalized time series.A successful utilization of the STL method is presented by [44] with the aim to study the streamflow in the Yellow River basin, using MF-DFA.

2.
Then, we find the 'profile' Y(i) of the deseasonalized time series x k of length N: where <x> is the mean of the time series and i = 1, . . ., N.

3.
Y(i) is then divided into N s ≡ int(N/s) boxes of equal length s, where s is the 'time scale'.It should be noted that N/s must be an integer, otherwise there will be a remaining number of profile points.However, very often N/s is not an integer and this problem is overcome by repeating the same procedure starting from the end.Thus, we get 2N s boxes.4.
In each box of length s, a least squares line is fitted to the data, which represents the trend in that box, i.e., the local trend.By subtracting the local trends, we detrend Y(i) and thus the variance F 2 (ν,s) of each segment (box) ν (ν = 1, . . ., 2N s ) is calculated.

5.
Taking the average over all segments, we find the qth order fluctuation function: 1 q (2) 6.This quantity is calculated repeatedly for all time scales to determine the relationship between F q (s) and s.Typically F q (s) is an increasing function of s.

7.
Making the log-log plots F q (s) versus s for each value of q, we can examine the scaling behavior of F q (s).If the series x k are long-range power-law correlated, F q (s) increases following a power-law: The exponent h(q) usually depends on q.If q = 2, then, for stationary time series, h(2) is the Hurst exponent H.If the time series is monofractal, then h(q) is independent of q.In the case of different scaling between small and large fluctuations, significant dependencies are observed.
From (3), it is clear that h(q) takes a variety of values (for each value of q) for multifractal signals and it is called generalized Hurst exponent.Neglecting for a moment the dependency of Hurst exponent on q and assuming that h(q) = H it is important to note that:

•
If 0 < H < 0.5 then the time series is long-range anticorrelated, that is, an increase in the value is more likely to be followed by a decrease (anti-persistent behavior) and vice versa.

•
If H = 0.5 the time series is uncorrelated (white noise).In this case, the probability that an increase will be followed by an increase or decrease is equal.

•
If H > 0.5 the time series is long range positively correlated, that is, an increase is more likely to be followed by an increase (persistent behavior) and vice versa.
If we use the relationship τ(q) = qh(q) − 1, then (via a Legendre transform) τ'(q) = α and The quantity α is called singularity strength (or Hölder exponent) and f (α) expresses the dimension of the subset of the time series that is characterized by α.The plot of α versus f (α) is called singularity spectrum, or multifractal spectrum.Usually, this plot has the shape of an upside-down parabola as in Figure 2. From the multifractal spectrum we can draw significant conclusions about the multifractality of a time series.The value of α where f (α) becomes maximum is when df (α)/dα = 0 and this happens for q = 0.Then, from (4) we find that f(α) max = 1.This value of α corresponds to the most dominant scaling behavior [45] and it is the dominant Hurst exponent.Hereafter, this quantity will be denoted by α 0 .
The second important property of the spectrum is its width (α max − α min ), where α max and α min are the maximum and the minimum values of α, respectively, for which f (α) = 0 as it is depicted in Figure 2 and it is a measure of the multifractality of the time series.A spectrum with a broad width is indicative of a strong multifractality (i.e., it has a 'fine' structure).If the width becomes smaller, then the time series tends to be a monofractal one.The spectrum width of a pure monofractal time series is equal to zero.More details about the multifractal spectrum can be found in [46].A measure of the width of the multifractal spectrum can be obtained by fitting a second-order equation to the curve of the spectrum around α 0 as recommended by [47]: Among the three constants A, B, C, the most important is B, which is an asymmetry parameter.When B = 0, the shape of the spectrum is symmetrical, whereas for positive B values, the spectrum is right-skewed and for negative left-skewed [48].A right-skewed spectrum is related to relatively strongly weighted high fractal exponents, while a left-skewed spectrum is indicative of low fractal exponents (a more 'regular' time series).According to [46], a time series with a high value of α 0 , broad width and right-skewed spectrum is considered to be more 'complex' than a time series with the opposite characteristics.Furthermore, the shape of the multifractal spectra can be also characterized by its truncation type.Considering the range of the q values, the shape of the spectra can be symmetrical, left or right truncated.In our case, the truncation type is grouped into 4 categories, LL, L, S and R where:  LL: there is a high degree of truncation on the left side.(i.e., when the left leg of a spectrum is truncated more than its half-way point).
 R: the spectrum is right-truncated.
 S: the spectrum is symmetrical (there is no significant truncation).
To make this notation clearer, a representative plot for each case using experimental data is depicted in Figure 3.The left end of the curve of the spectrum corresponds to the greatest value of q (in our case q = 6) and the right end corresponds to the smallest q value (here q = −6), while the maximum value of f(α) corresponds to q = 0. Therefore, the spectrum is left-truncated then the change of f(α) is smaller for positive values of q than for negative ones, which means that the multifractal structure of the time series shows an insensitivity to local fluctuations with large magnitudes.Working in a similar way for the right-truncated spectrum, we find that there is insensitivity to local fluctuations with small magnitudes.If the spectrum is symmetrical, then the multifractal structure of the time series exhibits the same sensitivity regardless the magnitude of local fluctuations [45].
Schematic illustration of a multifractal spectrum and its main characteristics.
Furthermore, the shape of the multifractal spectra can be also characterized by its truncation type.Considering the range of the q values, the shape of the spectra can be symmetrical, left or right truncated.In our case, the truncation type is grouped into 4 categories, LL, L, S and R where: • L: the spectrum is left-truncated.

•
LL: there is a high degree of truncation on the left side.(i.e., when the left leg of a spectrum is truncated more than its half-way point).• R: the spectrum is right-truncated.• S: the spectrum is symmetrical (there is no significant truncation).
To make this notation clearer, a representative plot for each case using experimental data is depicted in Figure 3.The left end of the curve of the spectrum corresponds to the greatest value of q (in our case q = 6) and the right end corresponds to the smallest q value (here q = −6), while the maximum value of f (α) corresponds to q = 0. Therefore, the spectrum is left-truncated then the change of f (α) is smaller for positive values of q than for negative ones, which means that the multifractal structure of the time series shows an insensitivity to local fluctuations with large magnitudes.Working in a similar way for the right-truncated spectrum, we find that there is insensitivity to local fluctuations with small magnitudes.If the spectrum is symmetrical, then the multifractal structure of the time series exhibits the same sensitivity regardless the magnitude of local fluctuations [45].

189
 S: the spectrum is symmetrical (there is no significant truncation).

190
To make this notation clearer, a representative plot for each case using experimental data is 191 depicted in Figure 3.The left end of the curve of the spectrum corresponds to the greatest value of q 192 (in our case q = 6) and the right end corresponds to the smallest q value (here q = −6), while the 193 maximum value of f(α) corresponds to q = 0. Therefore, the spectrum is left-truncated then the   According to [27], multifractality is caused due to either a broad probability density function for the time series values or to different long-range correlations for fluctuations of large and small magnitudes.In order to examine which of the two cases causes multifractality in a time series, its values are ordered randomly (shuffled) and MF-DFA is applied to the shuffled time series.In the case where multifractal characteristics are preserved, then multifractality is caused by a broad probability density function.In the other case (i.e., the multifractality of the shuffled time series appears to be significantly weaker than that of the original time series), the multifractality is mainly due to different long-range correlations for different magnitudes of the time series fluctuations.
However, in nature, multifractality can be driven by one of the sources or simultaneously by both of them [49,50].
Finally, the effects of seasonality and trends on the analyzed time series are subject to uncertainties regarding the DFA and therefore, MF-DFA results [41,51].

Air temperature multifractal characteristics
Prior to MF-DFA analysis, the STL decomposition method is applied to the temperature time series.The resulting components along with the original time series are illustrated in Figure 4.According to [27], multifractality is caused due to either a broad probability density function for the time series values or to different long-range correlations for fluctuations of large and small magnitudes.In order to examine which of the two cases causes multifractality in a time series, its values are ordered randomly (shuffled) and MF-DFA is applied to the shuffled time series.In the case where multifractal characteristics are preserved, then multifractality is caused by a broad probability density function.In the other case (i.e., the multifractality of the shuffled time series appears to be significantly weaker than that of the original time series), the multifractality is mainly due to different long-range correlations for different magnitudes of the time series fluctuations.However, in nature, multifractality can be driven by one of the sources or simultaneously by both of them [49,50].
Finally, the effects of seasonality and trends on the analyzed time series are subject to uncertainties regarding the DFA and therefore, MF-DFA results [41,51].

Air Temperature Multifractal Characteristics
Prior to MF-DFA analysis, the STL decomposition method is applied to the temperature time series.The resulting components along with the original time series are illustrated in Figure 4.
Upon removing the seasonal component, MF-DFA analysis is applied to the deseasonalized time series.The plots of: i) Fluctuation function F q (s) versus the time scale s (log-log plots), ii) Generalized Hurst exponent h(q) versus q and iii) Multifractal spectrum f (α) versus α are obtained for each time series for q ranging from −6 up to +6.Due to space limitations, only selected characteristic plots are presented.In Figure 5, the plots of F q (s) vs s, h(q) vs q and the multifractal spectrum f (α) vs α are illustrated for Kerkira station (the same shape of the relevant plots was found for almost all stations under study).The values of s used for the plotting of F q (s) are approximately from 30 days (s ≈ 10 1.5 ) to 8.5 years (s ≈ 10 3.5 ).It should be noted that the upper limit of the time scale is N/5, where N is the time series length and for stations with different time series lengths (i.e., different observation periods), the upper limit changes accordingly.

224
Generalized Hurst exponent h(q) versus q and iii) Multifractal spectrum f(α) versus α are obtained for 225 each time series for q ranging from −6 up to +6.Due to space limitations, only selected characteristic 226 plots are presented.In Figure 5, the plots of Fq(s) vs s, h(q) vs q and the multifractal spectrum f(α) vs α (F q (s)) q=-6 q=-3 q=0 q=3 q=6 q = 6: y = 0.67x − 0.22 (R 2 = 0.995) q = 3: y = 0.69x − 0.33 (R 2 = 0.996) q = 0: y = 0.71x − 0.45 (R 2 = 0.996) q = −3: y = 0.73x − 0.56 (R 2 = 0.995) q = −6: y = 0.76x − 0.68 (R 2 = 0.993)  Generalized Hurst exponent h(q) versus q and iii) Multifractal spectrum f(α) versus α are obtained for each time series for q ranging from −6 up to +6.Due to space limitations, only selected characteristic plots are presented.In Figure 5, the plots of Fq(s) vs s, h(q) vs q and the multifractal spectrum f(α) vs α are illustrated for Kerkira station (the same shape of the relevant plots was found for almost all stations under study).The values of s used for the plotting of Fq(s) are approximately from 30 days (s ≈ 10 1.5 ) to 8.5 years (s ≈ 10 3.5 ).It should be noted that the upper limit of the time scale is N/5, where N is the time series length and for stations with different time series lengths (i.e., different observation periods), the upper limit changes accordingly.
(a) From the Fq(s) vs s log-log plots for all temperature time series of all stations, it is observed that they have the form of an almost straight line (i.e., constant slope).It is noteworthy that slope values decrease when q has greater values, i.e., in Figure 5 the slope is 0.76 for q = −6, while for q = 6 the slope is 0.67.Moreover, the coefficient of determination (Pearson correlation coefficient) for all q values is higher than 0.98, which means that these scaling exponents are considered to be representative of the underlying scaling phenomenon [52].
It should be mentioned that, in the small segments (s small), it is easier to distinguish local periods with small fluctuations (q < 0) from periods with large fluctuations (q > 0).On the other hand, the large segments include periods of small and large fluctuations and thus the differences in magnitude are cancelled.Therefore, for large time scales, the behavior of Fq(s) is more similar to that of monofractal time series [46].This can be noted in Figure 5a, where for large time scales the lines of Fq(s) are coming closer together.It can also be observed that the log-log fit of Fq(s) is not very good for small values of s and especially for q = −6.This could be attributed to the fact that the root mean square value is very sensitive to local fluctuations with small magnitudes (q < 0, the smaller the q, the higher the sensitivity), especially for small segments (which correspond to small s values).This is reflected in the existence of a long right tail at the multifractal spectrum, as it has previously been explained in methodology section.
Calculating the slope of Fq(s) for a multitude of q values, the values of generalized Hurst exponent h(q) are obtained.Here, h(q) is calculated for q values from −6 up to +6 in steps of 0.1 and therefore h(q) is found for 120 values of q.The multifractality of the time series can be concluded From the F q (s) vs s log-log plots for all temperature time series of all stations, it is observed that they have the form of an almost straight line (i.e., constant slope).It is noteworthy that slope values decrease when q has greater values, i.e., in Figure 5 the slope is 0.76 for q = −6, while for q = 6 the slope is 0.67.Moreover, the coefficient of determination (Pearson correlation coefficient) for all q values is higher than 0.98, which means that these scaling exponents are considered to be representative of the underlying scaling phenomenon [52].
It should be mentioned that, in the small segments (s small), it is easier to distinguish local periods with small fluctuations (q < 0) from periods with large fluctuations (q > 0).On the other hand, the large segments include periods of small and large fluctuations and thus the differences in magnitude are cancelled.Therefore, for large time scales, the behavior of F q (s) is more similar to that of monofractal time series [46].This can be noted in Figure 5a, where for large time scales the lines of F q (s) are coming closer together.It can also be observed that the log-log fit of F q (s) is not very good for small values of s and especially for q = −6.This could be attributed to the fact that the root mean square value is very sensitive to local fluctuations with small magnitudes (q < 0, the smaller the q, the higher the sensitivity), especially for small segments (which correspond to small s values).This is reflected in the existence of a long right tail at the multifractal spectrum, as it has previously been explained in methodology section.
Calculating the slope of F q (s) for a multitude of q values, the values of generalized Hurst exponent h(q) are obtained.Here, h(q) is calculated for q values from −6 up to +6 in steps of 0.1 and therefore h(q) is found for 120 values of q.The multifractality of the time series can be concluded from the h(q) versus q plots.The dependency of the generalized Hurst exponent h(q) on q is evident from Figure 5b and this is a characteristic of multifractal time series, whereas in monofractal time series, the Hurst exponent is constant.Moreover, the fact that h(q) > 0.5 for all temperature time series leads to the conclusion that these time series exhibit long-range positive correlations.This means that an increase in temperature values is more likely to be followed by another increase.Furthermore, the slope of h(q) is greater for negative q (i.e., for small fluctuations) than for positive values of q (large fluctuations), which means that there is a higher degree of multifractality for the smaller fluctuations, in agreement with the aforementioned relevant conclusion from F q (s) plots.As it was previously mentioned, the multifractal spectrum reveals significant properties about the multifractality of the time series.Therefore, we focused on its main characteristics, the value of α 0 , the width α max − α min , the asymmetry parameter B and the truncation type of the spectrum.The results are presented in Table 2. Regarding the α 0 results, it is observed that the values are within a narrow zone ranging from 0.662 up to 0.775 for all time series.Thus, the dominant Hurst exponent of the examined temperature time series is of the order of 0.7, which means that the temperature time series are long range positively correlated.It should also be noted that for T max the majority of α 0 values is between 0.65 and 0.70, while for T mean and T min the more α 0 values are in the range 0.70-0.75.Regarding the multifractal (singularity) spectrum width values, (Table 2) it is observed that there is a significant variation, ranging from 0.342 to 0.759 for T max , from 0.380 to 0.727 for T mean and from 0.286 to 0.716 for T min and their distributions are illustrated in the histograms of Figure 6.It is observed that for all temperature time series (T mean , T max and T min ), the majority of the spectrum width values are within 0.4 and 0.7.The asymmetry parameter B values are in the range from −0.484 up to 0.685.It is easily noticeable that the majority of the B values are positive and that negative values are observed only for T min .Thus, most of the spectra are right-skewed, which indicates that, in most cases, there are relatively strongly weighted high fractal exponents.Focusing on the truncation type, it is clear that the majority of the spectra are left-truncated, leading to the conclusion that, in general, the multifractal structure of temperature time series is robust to local fluctuations with large magnitudes.
Concerning the type of multifractality, MF-DFA application to the shuffled time series reveals that in all cases the multifractality is very weak and it is mainly due to different long-range correlations for small and large fluctuations.A characteristic case is presented in Figure 7, where the multifractal spectrum for Kerkira daily maximum temperature is depicted for both the original and the shuffled time series.It was found that in all temperature time series, the spectra for the original and shuffled time series are similar to those of Figure 7.In shuffled time series we observe that the width of the spectrum for the shuffled time series is much smaller than the original time series and the value of α 0 is very close to 0.5 and, therefore, diminished multifractality (compared to the unshuffled series) shows that long-range correlations play the main role in the multifractality of the data [26].
It should be noted that in order to avoid statistical errors, data shuffling of each time series was repeated 100 times, where the initial data from the time series were randomly shuffled at each repetition, and the mean values of multifractal parameters were calculated [49].The corresponding mean value of the spectrum width for all stations and temperature time series is of the order of 0.2, i.e., in our case a time series with spectrum width less than 0.2 can be considered as monofractal [53].Considering the values of spectrum width (Table 2), it is concluded that all temperature time series exhibit a multifractal behavior.
the majority of the spectra are left-truncated, leading to the conclusion that, in general, the multifractal structure of temperature time series is robust to local fluctuations with large magnitudes.Concerning the type of multifractality, MF-DFA application to the shuffled time series reveals that in all cases the multifractality is very weak and it is mainly due to different long-range correlations for small and large fluctuations.A characteristic case is presented in Figure 7, where the multifractal spectrum for Kerkira daily maximum temperature is depicted for both the original and the shuffled time series.It was found that in all temperature time series, the spectra for the original and shuffled time series are similar to those of Figure 7.In shuffled time series we observe that the width of the spectrum for the shuffled time series is much smaller than the original time series and the value of α0 is very close to 0.5 and, therefore, diminished multifractality (compared to the unshuffled series) shows that long-range correlations play the main role in the multifractality of the data [26].
It should be noted that in order to avoid statistical errors, data shuffling of each time series was repeated 100 times, where the initial data from the time series were randomly shuffled at each repetition, and the mean values of multifractal parameters were calculated [49].The corresponding mean value of the spectrum width for all stations and temperature time series is of the order of 0.2, i.e., in our case a time series with spectrum width less than 0.2 can be considered as monofractal [53].
Considering the values of spectrum width (Table 2), it is concluded that all temperature time series exhibit a multifractal behavior.

Spatial Distributions
The spatial distributions of spectral width, a0 and asymetry parameter for the Tmean, Tmax and Tmin over Greece are illustrated in Figures 8, 9 and 10, respectively.For a better interpretation of those figures, isolines are also shown.

Spatial Distributions
The spatial distributions of spectral width, a 0 and asymetry parameter for the T mean , T max and T min over Greece are illustrated in Figures 8-10, respectively.For a better interpretation of those figures, isolines are also shown.

Spatial Distributions
The spatial distributions of spectral width, a0 and asymetry parameter for the Tmean, Tmax and Tmin over Greece are illustrated in Figures 8, 9 and 10, respectively.For a better interpretation of those figures, isolines are also shown.Regarding the spatial distribution of spectral width, two distinct regions are observed.A region of relatively higher multifractal spectrum width values, which includes a great part of continental Greece and Aegean Sea and a part of Ionian islands, and a region with lower spectrum values, which includes southeastern areas, a part of south mainland and parts of northwestern Greece.The first region has a richer multifractal spectrum, which means that weather conditions-and consequently the temperature-in this region are affected in a more 'complex' way than the rest of Greece.
Possible causes of local day-to-day temperature variations are the depression tracks over Greece from western directions [54].The first region includes the coastal areas of eastern mainland and some parts of the western Greece.Taking also into account the topography of Greece, temperature at these areas could be affected by other factors, such as cold outbreaks in the winter mainly in eastern Greece [55], the onset of Etesians winds in the summer mainly over Aegean Sea [56], katabatic winds causing Föhn effect, thermal instability in summer at continental areas, or probably other local circulations.On the other hand, the more 'stable' climatic conditions at the marine and coastal parts of the second region and the less irregular topography contribute towards a relatively simpler temperature variability.
The spatial distribution of α0 values (Figure 9) exhibits a rather smooth field of the order of 0.7 Regarding the spatial distribution of spectral width, two distinct regions are observed.A region of relatively higher multifractal spectrum width values, which includes a great part of continental Greece and Aegean Sea and a part of Ionian islands, and a region with lower spectrum values, which includes southeastern areas, a part of south mainland and parts of northwestern Greece.The first region has a richer multifractal spectrum, which means that weather conditions-and consequently the temperature-in this region are affected in a more 'complex' way than the rest of Greece.Possible causes of local day-to-day temperature variations are the depression tracks over Greece from western directions [54].The first region includes the coastal areas of eastern mainland and some parts of the western Greece.Taking also into account the topography of Greece, temperature at these areas could be affected by other factors, such as cold outbreaks in the winter mainly in eastern Greece [55], the onset of Etesians winds in the summer mainly over Aegean Sea [56], katabatic winds causing Föhn effect, thermal instability in summer at continental areas, or probably other local circulations.On the other hand, the more 'stable' climatic conditions at the marine and coastal parts of the second region and the less irregular topography contribute towards a relatively simpler temperature variability.
The spatial distribution of α 0 values (Figure 9) exhibits a rather smooth field of the order of 0.7 (α 0 values are in the range 0.66-0.77 in the Greek area).Therefore, the prevailing behavior of all temperature time series is persistent, which means that an increase in temperature is more likely to be followed by another increase, in agreement with the values of h(q), as in Figure 5b.In more detail, relatively larger values are noticed at two geographical areas: i) At southeastern Greece and ii) At some parts of coastal western Greece.Considering that α 0 is the dominant Hurst exponent, it is deduced that temperature time series for these areas appear to have relatively more persistent behavior than the other areas of Greece.A possible explanation could be that these areas have more 'stable' climatic conditions throughout the year, which agrees with [14] results for Turkey.Additionally, the absence of rainfall in the summer period is more pronounced at coastal and marine areas, whereas in continental Greece, the precipitation amount exhibits a more uniform interannual distribution, due to atmospheric instability especially in the warm period [57,58].Furthermore, it is also interesting that continental areas in Greece appear to have greater temperature variability than the marine and coastal areas [59] which is in partial agreement with the α 0 spatial distribution findings.It is noteworthy that higher α 0 values for the T min time series-compared to the corresponding values for T max and T mean -are observed for a significant part of continental Greece.This could be attributed to the fact that summer thermal instability is observed at noon and in the afternoon whereas the minimum temperature is observed in the early morning.
Regarding the spectral asymmetry spatial distribution, it is observed that the lower values of parameter B are found, in general, mainly in southeastern Greece.Spectra with higher B values contain a greater number of high fractal components and therefore, according to [47,48] are related to more 'complex' time series and it could be assumed that the pre-mentioned greater climate stability at the southeastern area leads to less 'complex' time series.
Finally, it is noteworthy to mention that in order to achieve a coverage of the entire Greek area for the spatial distribution of the main multifractal characteristics, we used the relevant temperature data from 13 neighboring stations.These stations together with their main multifractal spectrum characteristics (which were found using the same analysis as for the Greek stations) are listed in Table 3.A brief remark about the stations of Table 3 is that the values of the singularity spectrum width are in the same ranges as the Greek stations.

Multifractal Spectrum Parameters Intercorrelations
Regarding the multifractal spectrum parameter intercorrelations, we examined initially the relationship of a 0 and width with B and subsequently their relationship with the truncation type.Although no clear relationship was found between α 0 and spectral width, the analysis revealed a well-defined group of stations with LL truncation type, for α 0 > 0.68 and spectral width >0.6 for all temperature time series (T mean , T max and T min ) as illustrated in Figure 11.These cases, considering that a left truncated spectrum has mainly high fractal exponents (i.e., 'fine' structure), can be considered more complex than the remaining time series.It should also be noted that no significant geographical correlations were revealed; a similar finding was found in another Mediterranean region [48].

Conclusions
The analysis reveals that air temperature over the study area exhibits long-range temporal correlations that cannot be fully described by a single scaling exponent.According to the Fq(s) vs. s plots, the multifractal character of the examined temperature time series is verified for all sites.The generalized Hurst exponent h(q) is dependent on q in all cases, thus confirming the multifractal nature of the examined time series.It is also found that h(q) is greater than 0.5, which reveals that the time series are long-term positive correlated, meaning that an increase in temperature is more likely to be followed by another increase.Emphasizing the multifractal spectrum (i.e., f(α) vs. α plots), significant multifractal features were deduced from examining the dominant Hurst exponent (α0), the spectrum width, the asymmetry parameter B and the truncation type.The spatial distribution of the above parameters indicates that larger multifractal spectrum widths could be associated with regions where weather conditions and consequently the temperature are affected in a more 'complex' way than in other areas.Furthermore, relatively higher α0 values as well as higher values of the parameter B may be associated with regions of more stable climatic conditions.The dominant truncation type is left truncation, and thus for the vast majority of the time series, their multifractal Finally, we examined the correlation between the multifractal parameters of temperature time series and the main descriptive statistical parameters (i.e., mean, standard deviation, skewness and kurtosis) of the main climatic parameters (i.e., relative humidity, wind speed and precipitation amount).The results did not exhibit a clear association, verifying the complexity of the factors that affect temperature fluctuations.

Conclusions
The analysis reveals that air temperature over the study area exhibits long-range temporal correlations that cannot be fully described by a single scaling exponent.According to the F q (s) vs. s plots, the multifractal character of the examined temperature time series is verified for all sites.The generalized Hurst exponent h(q) is dependent on q in all cases, thus confirming the multifractal nature of the examined time series.It is also found that h(q) is greater than 0.5, which reveals that the time series are long-term positive correlated, meaning that an increase in temperature is more likely to be followed by another increase.Emphasizing the multifractal spectrum (i.e., f (α) vs. α plots), significant multifractal features were deduced from examining the dominant Hurst exponent (α 0 ), the spectrum width, the asymmetry parameter B and the truncation type.The spatial distribution of the above parameters indicates that larger multifractal spectrum widths could be associated with regions where weather conditions and consequently the temperature are affected in a more 'complex' way than in other areas.Furthermore, relatively higher α 0 values as well as higher values of the parameter B may be associated with regions of more stable climatic conditions.The dominant truncation type is left truncation, and thus for the vast majority of the time series, their multifractal structure is insensitive to larger local variations of temperature.In general, atmospheric circulation and weather conditions affect the shape of the spectrum, and the corresponding spatial distributions results highlight the complexity with which the natural processes interact and the importance of finding associations among the multifractal characteristics and different atmospheric processes.Moreover, the intercorrelation results among α 0 and spectrum width with the truncation type revealed a well-defined group of stations that exhibit similar multifractal characteristics.In particular, left truncation spectra are, in general, associated with higher spectrum width values.In addition, the analysis presented in this study could be useful for assessing the performance of climate models in terms of their ability to simulate reasonably well the observed air temperature non-linear dynamics.More specifically, as future work, the proposed methodology will be performed at larger or even global spatial scales for experimental temperature datasets and climate model simulations.

Figure 1 .
Figure 1.Meteorological station spatial distribution.The numbers denote the stations as listed in

Figure 1 .
Figure 1.Meteorological station spatial distribution.The numbers denote the stations as listed in Table1.

Figure 2 .
Figure 2. Schematic illustration of a multifractal spectrum and its main characteristics.

18 181Figure 2 .
Figure 2. Schematic illustration of a multifractal spectrum and its main characteristics.
spectrum is right-truncated.
194 change of f(α) is smaller for positive values of q than for negative ones, which means that the 195 multifractal structure of the time series shows an insensitivity to local fluctuations with large 196 magnitudes.Working in a similar way for the right-truncated spectrum, we find that there is 197 insensitivity to local fluctuations with small magnitudes.If the spectrum is symmetrical, then the 198 multifractal structure of the time series exhibits the same sensitivity regardless the magnitude of 199 local fluctuations [45].

Figure 3 .
Figure 3. Representative multifractal spectra plots (circles) related to the truncation type for: (a) T min for Mitilini (case LL); (b) T max for Alexandroupolis (case L); (c) T mean for Kos (case S); (d) T min for Herakleion (case R).

Atmosphere 2019 , 18 220Figure 4 .
Figure 4. Original time series (at the top of each plot), remainder, seasonal and trend components

227
are illustrated for Kerkira station (the same shape of the relevant plots was found for almost all stations 228 under study).The values of s used for the plotting of Fq(s) are approximately from 30 days (s ≈ 10 1.5 ) to 8.5 229 years (s ≈ 103.5 ).It should be noted that the upper limit of the time scale is N/5, where N is the time 230 series length and for stations with different time series lengths (i.e., different observation periods),231the upper limit changes accordingly.

Figure 4 .
Figure 4. Original time series (at the top of each plot), remainder, seasonal and trend components (from top to the bottom in each plot), for T max at Kerkira station.

Figure 4 .
Figure 4. Original time series (at the top of each plot), remainder, seasonal and trend components (from top to the bottom in each plot), for Tmax at Kerkira station.

Figure 5 .
Figure 5. Plots of (a) Log-log plot of Fq(s) vs s; (b) Generalized Hurst exponent h(q) vs q; (c) Multifractal spectrum f(α) vs α.All the plots are for daily maximum temperature time series at Kerkira station.

Figure 7 .
Figure 7. Multifractal spectra for daily maximum temperature at Kerkira station, for the original (right) and for the shuffled time series (left).

Figure 7 .
Figure 7. Multifractal spectra for daily maximum temperature at Kerkira station, for the original (right) and for the shuffled time series (left).

Figure 7 .
Figure 7. Multifractal spectra for daily maximum temperature at Kerkira station, for the original (right) and for the shuffled time series (left).

Figure 9 . 18 Figure 9 .Figure 10 .
Figure 9. Spatial distribution of α 0 values for (a) T mean ; (b) T max ; (c) T min in grayscale (higher values are illustrated by lighter shades in grayscale).

Figure 10 .
Figure 10.Spatial distribution of spectral asymmetry values for (a) T mean ; (b) T max ; (c) T min in grayscale (higher values are illustrated by lighter shades in grayscale).

Atmosphere 2019 , 18 Figure 11 .
Figure 11.Scatter plots of spectrum width vs dominant Hurst exponent (i.e., α0 values) for (a) Tmean; (b) Tmax; (c) Tmin daily temperature time series, with respect to truncation type (LL: *, L : •, S : ■, R : ▲).Finally, we examined the correlation between the multifractal parameters of temperature time series and the main descriptive statistical parameters (i.e., mean, standard deviation, skewness and kurtosis) of the main climatic parameters (i.e., relative humidity, wind speed and precipitation amount).The results did not exhibit a clear association, verifying the complexity of the factors that affect temperature fluctuations.

Table 1 .
Meteorological stations and observation period.

Table 1 .
Meteorological stations and observation period.

Table 2 .
Values of the main multifractal characteristics.

Table 3 .
Values of the main multifractal characteristics for neighboring stations.