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

TL;DR

This paper presents a straightforward decomposition method for European temperature variability, revealing that long-range correlations are primarily due to short-term relaxation, seasonal cycles, and an interannual oscillatory mode around 7-8 years.

Contribution

It introduces a simple approach to explain temperature variability scaling behavior, emphasizing the role of specific short-term processes and interannual oscillations.

Findings

Interannual variability often exhibits a 7-8 year oscillatory mode.

Scaling behavior can be explained by short-term relaxation and seasonal cycles.

Spatial patterns of parameters support the identified oscillatory mode.

Abstract

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.

Figures

Equations

$$F_x^2(s)={\sigma }_x^2\left(L_q(0,s)+2\sum _{t=1}^{s-1}{C}_{xx}(t)L_q(t,s)\right).$$

$$F_{x+y}^2(s)=F_x^2(s)+F_y^2(s).$$

$$T_t\approx X_t+Y_t+Z_t,$$

$$X_{t_X}=c{X}_{t_X-1}+{ϵ}_{t_X}\text{mbox}andr=-1/\log (c),$$

$$Y_{t_Y}=\frac{{\sigma }_Y}{2}\sin \left(\frac{2\pi t_Y}{1\mathrm{year}}\right),$$

$$Z_{t_Z}=a{Z}_{t_Z-1}+b{Z}_{t_Z-2}+{\eta }_{t_Z},$$

$$P(\omega )=\mathcal{F}[⟨C(t){⟩}_{\Omega }],$$

$$var\left[\log \left(\frac{F_{data}(s)}{F_{model}(s)}\right)\right]$$

$$\sigma =\exp \left(\sum _{s\in S}\log \left(\frac{F_{data}(s)}{F_{model}(s)}\right)\right).$$

$$y_t=\sum _{j=1}^t{x}_j.$$

$$f^2(\nu ,s)=\frac{1}{s}\sum _{t=1+(\nu -1)s}^{\nu s}(y_t-p_t{)}^2.$$

$$F^2(s)=\frac{1}{K}\sum _{\nu =1}^{2K}{f}^2(\nu ,s).$$

$$F^2(s)={\sigma }^2\left(L_q(0,s)+2\sum _{t=1}^{s-1}C(t)L_q(t,s)\right).$$

$$\begin{array}[]{ll}L_{1}(t,s)&=\frac{1}{30(s^{4}-s^{2})}[3t^{5}-5(4s^{2}-1)t^{3}+30(s^{3}-s)t^{2}\\ &-(15s^{4}-35s^{2}+8)t+2(s^{5}-5s^{3}+4s)].\end{array}$$

$$C(t)=c^t,$$

$$F_c^2(s)={\sigma }^2\frac{c^s{J}_c(s)+K_c(s)}{15(c-1{)}^6(s^2-s^4)},$$

$$\begin{array}[]{ll}J_{c}(s)&=60[s^{2}(c^{2}-c)^{2}-3s(c^{3}-c)+2(c^{4}+c^{3}+c^{2})],\\ K_{c}(s)&=s^{5}(c-1)^{5}(c+1)+15s^{4}c(c-1)^{4}\\ &-5s^{3}(c-1)^{3}(1-7c-7c^{2}+c^{3})\\ &-15s^{2}c(c-1)^{2}(1-10c+c^{2})\\ &+2s(c-1)^{3}(2-17c-17c^{2}+2c^{3})\\ &-120c^{2}(1+c+c^{2}).\end{array}$$

$${\sigma }^2={\sigma }_{\eta }^2/(1-c^2),$$

$$C(t)=h_1{g}_1^t+h_2{g}_2^t,$$

$$(1-g_{1}B)(1-g_{2}B)x_t={\eta }_t,$$

$$\tau =\frac{2\pi }{\arctan (\mathrm{Im}(g_1)/\mathrm{Re}(g_1))}$$

$${\sigma }^2=\frac{(1-b){\sigma }_{\eta }^2}{(1+b)(1-a-b)(1+a-b)}$$

$$F_{a,b}^2(s)={\sigma }^2\frac{[c_1{F}_{g_1}^2(1-g_1^2)+c_2{F}_{g_2}^2(1-g_2^2)](1-b)}{(1+b)(1-a-b)(1+a-b)}.$$

Full text

I 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), conditional mean approaches (Jajcay et al., 2016), 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). In (Lovejoy, 2015) the author 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š, 2014; 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; Jajcay et al., 2016), however, to our knowledge there has not been an algorithmic investigation on a similar number of station as we perform it.

We present the foundations of our method in section 3, summarizing previous work in (Meyer and Kantz, 2019). In section 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 section 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.

II Data

We analyze station data of temperature and pressure in Europe (Tank et al., 2002). It can be downloaded from 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 A.

III The method

Our method is based on detrended fluctuation analysis (DFA) (Peng et al., 1994) (see appendix B 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}$

[TABLE]

Here $\sigma_{x}^{2}$ is the variance of the time series and $L_{q}$ is a kernel described in appendix C. 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).

In (Meyer and Kantz, 2019) the authors 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 C.

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 figure 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)

[TABLE]

This will be important for our purpose in the next section.

IV 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

[TABLE]

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 $\epsilon_{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.

Temperature correlations for short times decay exponentially (Massah and Kantz, 2016). The behaviour is apparent in figure 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)

[TABLE]

where the noise term $\epsilon_{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 $\epsilon_{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 $\sigma^{2}_{\epsilon}$ is a free parameter. For simplicity we will later look at the standard deviation $\sigma_{X}$ of the process instead of $\sigma_{\epsilon}$. They are related by equation (18).

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 equation (16) to $F_{T}$ for small $s$ (see figure 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 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 equation (6)) with $b=-1$, which is equivalent to a sinusoidal signal

[TABLE]

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? In (Plaut et al., 1995) the authors 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 (Jajcay et al., 2016). We want to model interseasonal fluctuations with an AR(2) model

[TABLE]

that - dependent on the parameters $a$, $b$ - can exhibit oscillatory behaviour with period $\tau$ (see equation (21)) or pure relaxations. Again we treat the standard deviation $\sigma_{Z_{t}}$ of the process as a free parameter, which determines $\sigma_{\eta}$ of the noise $\eta$ according to equation (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 figure 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 figure 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).

V 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$, $\sigma_{X}$, $\sigma_{Y}$, $a$, $b$, and $\sigma_{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 $\tau$ 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 figure 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 $\epsilon$ 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 figure 4. The amplitude $\sigma_{Y}$ of seasonality is strongly correlated with $\sigma_{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 $\tau$ between 10 and 22 half years - or not. The set $\Omega_{-}$ of stations without this mode consists of three subsets:

• •

if the slope is $\propto s^{0.5}$ or ever more flat the algorithm indicates an oscillation with a very short period

• •

if the slope is very steep ($\propto 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 A.

For 51% of the analyzed stations we observe an oscillation, see figure 5. We denote the set of these stations by $\Omega_{+}$. The categorization into the two subgroups uncovers clear 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 $\tau$ 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

[TABLE]

where $\mathcal{F}$ is the Fourier transform and the brackets $\langle\rangle$ denote the average over all stations in the set $\Omega_{+}$ or the average over all stations in the set $\Omega_{-}$. 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 figure 5 right).

VI 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.

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.

Appendix A 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 figure 2 if there are at least 4 years recorded. For the analysis of the anomalies in figure 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 $\Omega_{+}$, 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\cdot 75$ and the recursion $s_{i}=sup_{s\in\mathbb{N}}s\leq 0.9s_{i+1}$, which yields a logarithmic scale. For datasets with less than $4\cdot s_{max}$ recorded days the values of $s$ which are larger than $1/4$th of the length of the time series are ignored.

The short time behaviour of the data (4) is fitted for $3\leq s\leq 25$ days. For the analysis in figure 3, the pressure fluctuation function is also fitted for $s\leq 25$ only, unlike in figure 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

[TABLE]

The standard deviation as well as the amplitude of the seasonal cycle is not fitted but rather calculated as

[TABLE]

In figure 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 $\tau=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 B 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

[TABLE]

Then we divide the time axis into $K$ non-overlapping segments of length $s$ and calculate the so-called DFA variance $f^{2}(\nu,s)$ in every segment $\nu$ 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

[TABLE]

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

[TABLE]

Traditionally people look at the asymptotic behaviour ${F}^{2}(s)\sim s^{2\alpha}$, which contains information about the correlation structure. For $\alpha=1/2$ the process is short range correlated and for $\alpha>1/2$ the process is long range correlated.

Appendix C Fluctuation functions of AR(1) and AR(2)

Theoretical fluctuation functions can be calculated from equation (1)

[TABLE]

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

[TABLE]

The correlation function of AR(1) is known to be

[TABLE]

with AR(1) parameter $c$. The fluctuation function for AR(1) can therefore be calculated as

[TABLE]

with $J_{c}(s)$, $K_{c}(s)$ polynomials in $s$

[TABLE]

The variance of the process is given by

[TABLE]

where $\sigma_{\eta}^{2}$ is the variance of the noise.

For AR(2) the correlation function reads

[TABLE]

where $g_{1}$ and $g_{2}$ are real or complex roots obtained from rewriting the definition (6) as

[TABLE]

with the backshift operator $BX_{t}=X_{t-1}$. The constants $h_{1}$ and $h_{2}$ are calculated from $C(1)={a}/{(1-b)}$ and $C(0)=1$. The model might exhibit oscillatory dynamics with a period

[TABLE]

if $g_{1}$ and $g_{2}$ are complex. Its variance is

[TABLE]

The fluctuation function for AR(2) can be expressed in terms of the fluctuation function of AR(1), because of the correlation function (19) being the superposition of two AR(1) fluctuation functions. However, different variances of the two AR(1) processes have to be taken into account

[TABLE]

Here ${F}_{g_{i}}$ are the fluctuation functions of the AR(1) processes with parameters $g_{1}$ and $g_{2}$.

Acknowledgements.

We acknowledge the data providers in the ECA&D project. Klein Tank, A.M.G. and Coauthors, 2002. Daily dataset of 20th-century surface air temperature and precipitation series for the European Climate Assessment. Int. J. of Climatol., 22, 1441-1453. Data and metadata available at http://www.ecad.eu. Map background made with Natural Earth. Free vector and raster map data @ naturalearthdata.com.