A Study of the Dynamic Behaviour of Fine Particulate Matter in Santiago , Chile

We present here a study about the limitations found when trying to develop an accurate atmospheric particulate matter forecasting model based on real data, and evidence that the time series of fine particulate matter concentration exhibit deterministic chaotic behavior. We have calculated the Lyapunov exponents of PM2.5 time series obtained from measurements from four monitoring stations located in the city of Santiago, Chile, in recent years. Values obtained for the largest Lyapunov exponents turned out to be positive and ranging between 0.3 and 0.5 which, according to the theory of chaos, is a condition for the presence of deterministic chaos and random behavior in time series. Given the shape of decay of autocorrelation functions and values of correlation dimension and Hurst exponents, random behavior can be discarded: we therefore conclude that the series are chaotic and very sensitive to initial conditions. The study presented here can be replicated in other mid-sized cities that present similar situations to the city of Santiago, where complexity of topography, meteorology and seasonal trends favor the generation of high concentration episodes of atmospheric particulate matter and where a reliable air quality forecasting model may be important for environmental management.


INTRODUCTION
In urban centers around the world it is observed that on some days the population is exposed to unhealthy concentrations of atmospheric particulate matter.For this reason, it seems convenient that environmental authorities operate a forecasting system which can be used to prevent the population or to enforce emissions restrictions.
Atmospheric particulate matter concentration forecasting models may be classified as statistical or deterministic.Statistical models are based on historical information from the particulate matter (PM) time series and associated meteorology.During a training stage, using samples of an appropriate input-output relation, a number of adjustable parameters are calculated.During operational mode, values of pollutant concentration and meteorological variables measured within the last few hours are introduced in the algorithm that generates an estimation of the future value of the variable of interest.
Particulate matter forecasts using statistical models, in which historical values are used for parameter adjustment, may be of reasonable accuracy when estimation of concentrations a few hours ahead in time is the goal.However with these models in general it is not possible to generate spatial distribution of pollutant concentrations.Besides, errors become considerably large when forecasting more than one day into the future (Pérez et al., 2000;Ordieres et al., 2005;Grivas and Chaloulakou, 2006).Since these models are trained with data from the past, they usually tend to be less accurate for extreme cases (very low or very high pollutant concentration) that are observed with a relative low frequency.From an operational point of view, poor long-range forecasting will be a disadvantage, because high particulate matter episodes are of crucial interest for air pollution management and they are normally observed on few occasions during the year.Among statistical models used for air pollution forecasting we find linear regressions (Thomas and Jacko, 2007), artificial neural networks (Pérez and Reyes, 2006), discriminant analysis (Silva et al., 2001) and Kalman filtering (Van der Wal and Janssen, 2000).
On the other hand, deterministic models are chemistrytransport models (Grell et al., 1994;Byun and Ching, 1999;Tie et al., 2007;Jorba et al., 2008;Stern et al., 2008) for which, given the rate of emission of a variety of pollutants, meteorological fields and chemical boundary conditions, the time evolution of several chemical compounds and aerosols is calculated.In general, these models are composed of a meteorological driver (MD) and a chemical and transport module (CTM).Limitations on the forecasting accuracy of this type of model may be attributed to the difficulty to consider in detail all the relevant sources of air pollution and to take into account all the physical-chemical processes that have an effect on the evolution of particulate matter concentrations (Pool, 1989).Accuracy of forecasting is strongly dependent on the precise estimation of initial conditions (Manders et al., 2009).Konovalov et al (2009) presented a study that combines the properties of a statistical model and a deterministic model, which seems to optimally combine the accuracy of the former with the capability of longer time forecasting of the later.
Average PM concentration forecasting errors with one day anticipation are not significantly lower than 20% with any model, whether statistical or deterministic.It is possible that intrinsic properties of the PM time series do not allow for forecasting accuracies much better than those obtained with existing models, due in part to the evidence that these series are chaotic in nature (Chelani and Devotta, 2006a;Kumar et al., 2008;Lee and Lin, 2008).
According to Lorenz (1963) nonlinearity is inherent to atmospheric systems.Sivakumar et al (2007) points out the existence of nonlinear deterministic behavior in the air quality index in time series such as ozone.Raga and Le Moine (1996) have applied nonlinear dynamic tool on air quality data (NO, CO, SO 2 and O 3 ).
In the decade of the 80's, the analysis of time series through dynamical systems theory (Takens, 1981) has established that relevant information may be obtained by the calculation of a set of parameters (or indexes) as time lag, mutual information (Fraser and Swinney, 1986), embedding dimension, capacity dimension, correlation dimension, entropy of information, Hurst exponent and Lyapunov exponents, among others (Grassberger and Procaccia, 1983).
Among these parameters, the largest Lyapunov exponent is one of the most widely used indicator to detect chaotic behavior in a dynamical system consisting of a time series (Eckmann et al., 1986).
In this paper the PM 2.5 data recorded in four monitoring stations located in the urban area of Santiago, Chile during 2000-2006 is analyzed.The significantly positive Lyapunov exponents obtained are an indication that the evolution of concentrations is very sensitive to initial conditions and the estimation of capacity dimension, correlation dimension and Hurst exponent indicate that they are not random but chaotic.

Chaotic Behavior in Air Pollution Time Series
Linear methods for analysis and prediction of time series may be a starting point in the study of behavior of air pollution and meteorological data.However, the most relevant results have been obtained with the introduction of nonlinear tools.The reason for this is that in this type of data, small causes may have large effects in future values, which is the essence of nonlinearity.
Let us consider the following time series: where ∆ is the time difference between measurements of variable x.We can write this series in simpler form by calling x(t + kΔ) = x k .A distance d 0 between two arbitrary samples of the series x i and x j may be defined as: This corresponds to absolute value of difference in one dimension and to Euclidean distance in higher dimensions.
The distance between the respective following values is called d 1 : If we assume that the distance d N = ||x i+N -x j+N || increases exponentially with N, the exponent λ defined by Eq. ( 4) is called the Lyapunov exponent (Hilborn, 2000): Taking natural logarithm: In practice, the Lyapunov exponent is obtained from Eq. ( 5) in the limit of large N, for which saturation of log(d N ) becomes evident (see Fig. 3).What is found is that there are as many different Lyapunov exponents as the embedding dimension D e of the time series, where D e is the integer dimension of phase-space required to reproduce the steady state of the system.
Nonlinear systems are described by multidimensional vectors labeled by a time coordinate.The space in which these vectors lie is known as state space or phase space.What is shown in Eqs.(2)-( 5) corresponds to the case with D e = 1.These equations describe the way to calculate the single Lyapunov exponent for this simple case.In the general case, points must be represented by a D e -dimensional vector: and distances d N must be understood as Euclidian distances.
For a time series of scalar observations, which may be the case of the series of one hour average PM 2.5 concentrations measured at a monitoring station, Pérez and Reyes (2001) have shown a methodology for calculating the embedding dimension and the optimal time lag ∆.In this case, the embedding dimension corresponds to the length of a sequence of observations used as input to a multilayer neural network which are necessary to predict the next step value, where the time separation between inputs is ∆.
Lyapunov exponents may be positive, zero or negative.If the maximal Lyapunov exponent is negative, we are in the presence of a time series representing a dissipative system with a stable fixed point.If the motion settles down onto a limit cycle, the maximal Lyapunov exponent is zero.If a predominantly deterministic system is affected by random noise, the maximal Lyapunov exponent is infinite (Kantz and Schreiber, 2005).A positive finite maximal Lyapunov exponent would be an indication of chaos.For a given time series, the sum of all positive Lyapunov exponents defines its Kolmogorov entropy K (Kolmogorov, 1959), and its inverse quantifies the average time of predictability, T p = 1/K.In general, the largest positive Lyapunov exponent of a system λ L is obtained in a straightforward manner, and then we can establish that providing in this way a strong restriction for modeling (Liu et al., 2004).Lanfredi and Macchiato (1997) have found that random behavior dominates over the chaotic behavior in time series of pollutants like NO x , O 3 and CO.Raga and Le Moyne (1996) have applied nonlinear analysis to O 3 , NO 2 , SO 2 and CO time series obtained from 13 monitoring stations in Mexico City, finding low-dimensional chaotic behavior in those time.Related to total ozone concentration, Chattopadhyay and Chattopadhyay (2008) have found the existence of chaotic behavior in Arosa, using the correlation dimension method.On the other hand, Haase et al. (2002) have indicated that if there is any nonlinearity in the ozone time series it may be very weak.But different to that result, Chelani (2010) has observed chaotic signatures in ozone concentration time series in Delhi, India, for 1 h and 4 h sampling rates.Koçak et al. (2000) applied a local method to predict observed O 3 time series data in Istanbul.The attractor in a multidimensional space was reconstructed, concluding that O 3 concentration is governed by a deterministic chaotic system.
Chelani and Devotta (2006b) implemented a hybrid methodology in order to predict one step ahead air pollutant concentration time series using nonlinear dynamic modeling and applied to NO 2 concentrations observed at a site in Delhi that was more effective than an autoregressive model.Also, using neural networks and nonlinear dynamical systems theory, Chelani et al. (2005) have characterized and predicted ambient nitrogen dioxide time series concentration at three sites in Kolkata, India.At two sites in the Gdansk region in Poland, Khokhlov et al. (2008) have investigated chaotic behavior in the nitrogen dioxide and sulphurous anhydride concentration time series finding low-dimensional chaos (low-dimensional chaos is characterized by only one positive Lyapunov exponent).Using artificial neural networks and a new scheme to predict chaotic time series of air pollutant concentrations (applied to Lorenz map and ozone observations), Gautam et al. (2008) got better predictions as compared to standard backpropagation algorithm.Yu et al. (2011) have found clear evidence of chaotic behavior of PM 10 time series over a time span of 10 years in Lanzhou, China.
An indication of the presence of chaos in a time series is the value of the Lyapunov exponent, for which we describe a calculation method below.

METHODOLOGY
In what follows, the largest Lyapunov exponents for PM 2.5 series are calculated using the algorithm proposed by Wolf et al. (1985).In order to use this algorithm, an estimation of the embedding dimension of the system is necessary.One way to calculate embedding dimension is by analyzing the correlation dimension (Grassberger and Procaccia, 1983) as a function of embedding dimension D e .Correlation dimension is calculated as follows: for each data point generate a hyper dimensional sphere of embedding dimension D e and radius r with the point at its center.The fraction of subsequent points within the sphere is calculated for different values of r.Then the log of this number versus the log of the radius is plotted.The correlation dimension is defined as the average slope of the cumulative curve.As the embedding dimension increases, the correlation dimension also increases, but for most systems it saturates at some value.The smallest integer value of D e within the saturated region corresponds to the correct value of embedding dimension D e .
Once we have estimated the embedding dimension, it is possible to construct the state space from the experimental data record (Small, 2005;Takens, 1981).
Let us consider where y(t) is the reconstructed D e dimensional state vector, x(t) is the observed variable, τ is a time lag, calculated from the first minimum of the Average Mutual Information function (Fraser and Swinney, 1986;Abarbanel, 1996).This quantity is a set theoretic concept and as a function of time delay estimates the mutual information content of two sets of measurements.
It is worth mention that nearly random time series has also large Lyapunov exponents.However, in that case a signature is the linear autocorrelation function, defined as (Abarbanel, 1996): This correlogram is an image of correlation statistics and may be used to investigate if a sequence of data is random or if adjacent observations are related.
For a random series, C(t) goes abruptly to zero, while for chaotic deterministic data it shows slower decay.The calculation of additional statistical parameters as capacity dimension, correlation dimension, Hurst exponent and entropy may be related to the determination of whether the analyzed time series are chaotic.
The capacity dimension is the fractional dimension of any attractor of the time series obtained when data are plotted in a reconstructed phase space with a given embedding dimension.It is calculated by iteratively dividing the phase space with embedding dimension D e into equal hypercubes and plotting the log of the fraction of hypercubes containing data points against the log of the size of the edge of the hypercube.The average slope of the central part of the curve is taken as the capacity dimension.A dimension greater than five is an indication of essentially random data.
The correlation dimension is calculated by counting data points inside hyper-spheres of different radii centered on each data point in a reconstructed phase space with a given embedding dimension.It corresponds to the average slope of the cumulative curve obtained from the log of the record against the log of the radius.Also, a correlation dimension greater than about five implies that the time series represents data from a random process.
The Hurst exponent is a measure of the degree to which the data can be represented by a random walk.Here, the root-mean-square displacement is calculated as a function of time, where each point in the time series is taken as an initial condition.The slope of the curve obtained is the Hurst exponent.For ordinary Brownian motion the exponent is 0.5.Values between 0.5 and 1.0 indicate deterministic behavior and persistence.Hurst exponent significantly less than 0.5 are typical for deterministic behavior and antipersistence (Kantz and Schreiber, 2005).

Chaos in Santiago PM 2.5 Time Series
PM 2.5 series in several monitoring stations in Santiago, Chile show annual averages that are high relative to national standards.In addition, a seasonal effect is observed, according to which concentrations are significantly higher during fall and winter as compared to the rest of the year.This seasonal increase in concentrations, results from a combination of unfavorable meteorology dominated by severe temperature inversion, the presence of urban emission sources, and a confined geographical area (Ruttland and Garreaud, 1995).This pattern of atmospheric pollution is not unique to Santiago, but is also observed in other cities around the world, like Lohan, Utah, USA (Silva et al., 2007), Tucson, Arizona, USA (Bailey et al., 2011), Graz, Austria (Stadlober et al., 2008), Gothenburg, Sweden (Olofson et al;2000).When episodes of high concentrations occur in the cold period, harmful effects for the population may be prevented in part by the use of a forecasting model.However, when chaos is present in PM 2.5 time series, important limitations to the accuracy of the forecast are unavoidable.
The data from four PM 2.5 monitoring stations in Santiago is analyzed in this study.Their location within the urban area is shown in Fig. 1.Station L is 604 m over sea level and is located near a commercial area.An important contribution to PM 2.5 concentrations measured in this station may be attributed to traffic congestion produced at rush hours.Station M is 785 m over mean sea level (MSL) and it is located in a residential area.Fine particles measured at this station are expected to be produced by vehicle traffic and wood burning (in winter).This station registers PM 2.5 concentrations that in general are lower than those measured by the other three analyzed stations.A reason for this is that its altitude is often higher than altitude of thermal inversions during episodes (between 600 m and 700 m MSL).Station N is 550 m over sea level and is near down town and close to freeway.Station O is 481 m MSL.It is located in a residential area, and being located towards the extreme west and the lowest part of the city, it is expected to measure not only local emissions but also particles transported from other areas.During high concentration episodes this station registers the highest values.The accumulation of air pollution in this area is enhanced by the presence of a night breeze that blows from east to west, which is blocked by a coastal range with elevations over 1000 m that stands between the city and the Pacific Ocean.Given the high PM 2.5 concentrations observed in Santiago, environmental authorities have implemented control strategies, like restrictions to wood burning, vehicle circulation and closing of industries during episodes.We must have in mind that in urban areas an important fraction of the measured PM 2.5 is produced through chemical reactions that occur in the atmosphere (Wang et al., 2008;Seguel et al., 2009).
About 2% of our time series data are missing.For small gaps linear interpolation was applied and for bigger gaps, a cubic spline with single imputation technique was used (Junninen et al., 2004).
One hour average PM 2.5 concentrations during one year are shown in Fig. 2. Seasonal variation is evident from data in the four stations.Being Santiago located in the southern hemisphere, autumn starts in March and winter ends in September.Average concentrations are significantly higher between April and August during which more episodes of high PM 2.5 are observed.During this period, a greater number of episodes are observed in station O as compared with the rest of the year (notice the difference of scale in Fig. 2).Statistical properties of the PM 2.5 series for years from 2000 to 2006 are displayed in Tables 1 and 2. Table 1 shows statistics for the complete year and Table 2 for the period between April and August, which is the season with the worst air quality.From the period averages for all stations and all years it becomes evident the difference between the cold season and the rest of the year.By comparing maximum values in Tables 1 and 2 we can verify that most of them occur between April and August and that the highest values are observed in station O. Large standard deviations suggest that forecasting is not easy.
Skewness indicates the symmetry of the probability density function (PDF) of a time series.A time series with an equal number of large and small values has a skewness of zero (Joanes and Gill, 1998).The positive values of skewness in Tables 1 and 2 is an indication that there are relatively fewer large PM 2.5 concentrations in all series.This will be inconvenient for the construction of precise statistical forecasting models because parameters in them are calculated from historical cases and a low representation of high concentration cases will imply in general a tendency to under-predict extreme values.From an operational point of view, for air quality management in a city, detection of high concentrations situations are particularly important.Average skewness for series including one year data is 2.06, while the average for data from the cold season decreases to 1.54, which is consistent with the observation that higher concentrations occur in the fall and winter periods.
Kurtosis is a measure of the peakedness of the PDF of a time series.A kurtosis value close to three indicates a Gaussian like peakedness.PDFs with relatively sharp peaks have kurtosis greater than three.PDFs with flat peaks have kurtosis less than three.From the values of this parameter obtained from Santiago data, we observe significant differences from one monitoring station to another.It is not unexpected that largest kurtosis corresponds to station O, where higher concentrations within the city are verified, and where under unstable atmospheric conditions, the lowest concentrations are observed.
Since it is of interest to investigate the presence of chaos in PM 2.5 series, we proceeded, as mentioned in the previous section, to calculate the Lyapunov exponent using the Wolff algorithm.As a previous step it was necessary to estimate the time lag τ and the embedding dimension D e .The values obtained for our series of interest remain rather constant, and for all of them we estimated values τ = 11 and D e = 5.
Table 3 shows the largest Lyapunov exponents for PM 2.5 time series (one hour spacing between data) at four monitoring stations in Santiago, between years 2000 and 2006.We can verify that these values are relatively large, and then we may conclude, based on Eq. ( 7), that predictability .69 ( a ) The minimum value for all time series was 1 µg/m 3 .based on information contained in the series is not greater than 3 hours.The average values for all the stations are similar except for station O, which are somewhat smaller.The reason for this may be attributed to the fact that this station is located in an area with the lowest altitude in the city, where concentrations tend to have extreme values, very low under conditions of unstable atmosphere and very high under stable atmosphere, so their values may be easier to predict under specific conditions.One must take into account that this time limit may be increased if we incorporate information from relevant exogenous variables (meteorological for example, which however may also show chaotic behavior).Calculated Lyapunov exponents are between 0.274 and 0.504.These values are of the same order of magnitude as those obtained for PM 10 series in Mexico (Vasquez et al., 2012) and Taiwan (Lee and Lin, 2008).The Lyapunov exponents in Table 3 were calculated using an embedding dimension D e = 5.For the same embedding dimension, the Lyapunov exponent for PM 10 series in Taipei was 0.277 and they found a correlation dimension of 4.32.Both of these values are indications of chaos.Lyapunov exponents for PM 10 series in Mexico city were between 0.4 and 0.6, correlation dimension between 3.6 and 4.4 and a Hurst exponent between 0.2 and 0.25.The values of these three parameters are consistent with attribution of chaotic behavior to the PM 10 time series.Lyapunov exponents for PM 10 time series in Lanzhou are positive but much smaller (Yu et al., 2011), which we believe is due to an erroneous estimation of the time scale.
In our case we can test the exponential divergence by constructing five dimensional vectors, where components are five successive measurements in PM 2.5 series.
For an arbitrary measurement at time t we would have: In order to study the evolution of state vectors we calculated the Euclidian distance between points separated by time n∆.We selected at random points for which the distance to X t+∆ was smaller or of the order of 10 µg/m 3 and then proceeded to calculate the Euclidian distance to successive five dimensional vectors.Fig. 3 shows the logarithm of the distance as a function of time interval for different starting point.The fact that most of the curves show a linear increase may be seen as a verification of exponential divergence that in turn tells us the presence of chaos.
Correlograms provide a useful method of visualizing the spatial dependence between data points in relation to distance (Fortin and Dale, 2005).Furthermore, the correlogram is a tool for checking randomness, rising or declining trend, oscillation, etc., of a time series.If a time series is random, its autocorrelation function should be near zero for all times greater than zero, and if non-random, the autocorrelation function will decay slowly to zero.The correlogram for PM 2.5 series during 2005 is shown in Fig. 4, where significant correlation between measurements close in time is observed.Curves for the other years are not significantly different.This has to do with the presence of a daily cycle, where in general, higher concentrations are observed at times where anthropogenic activity is high and meteorological conditions for pollutant dispersion are unfavorable (around 8 a.m. in the morning and 8 p.m. in the evening).Table 4 shows the result of the calculation of some additional indexes for the case of station N, now for years 2000 through 2008.We observe that for all years capacity dimension and correlation dimension are less than five, which is an indication of non random behavior.Hurst exponent is for all years, except 2001, significantly smaller than 0.5, which is an indication of deterministic behavior and antipersistence (predominance of oscillation between high and low values).
The values of the parameters calculated for station N and shown in Table 4 confirm the conclusion that PM 2.5 series are chaotic.Values for stations L, M, and O are the same order of magnitude and are not shown.
Given the values of the largest Lyapunov exponents displayed in Table 3, the shape of the autocorrelation function shown in Fig. 4 and the values of the additional time parameter characterization presented in Table 4, we can conclude that the analyzed PM 2.5 time series are chaotic.

DISCUSSION
The evidence that the time series of atmospheric concentrations of particulate matter four monitoring stations in Santiago, Chile, exhibit deterministic chaotic behavior imposes severe restrictions about the predictability of air quality in the city.Although a positive finite maximal Lyapunov exponent is not by itself an unambiguous indication of chaos.In our case, the results of calculation of additional statistical parameters confirm this behavior, which means that the series are deterministic and very sensitive to initial conditions.Statistical forecasting models, which constitute a pragmatic approach, adjusting parameters on the basis of historical data, would appear as a reliable tool for operational forecasting of particulate matter concentrations in urban regions only for a few hours into the future, provided relevant meteorological information is included as input (Klingner and Sähn, 2008;Pérez and Salini, 2008;Stadobler et al., 2008;Stern et al., 2008).Deterministic photo-chemical forecasting models (Konovalov et al., 2009) would require detailed consideration of topography, emission sources and atmospheric processes.Any spatial or temporal averaging oriented to simplify the analysis may be critical if they leave out conditions that are relevant for estimating future patterns.A combination of deterministic and statistical approaches may be an option for more accurate results and longer time forecasting (Konovalov et al., 2009).For high concentration episodes, which in general are of great interest for air pollution management, special effort must be put towards overcoming shortcomings of prognostic weather prediction models for the simulation of stagnant weather situations.The results of this study show that particulate matter forecasting in a large city like Santiago will exhibit intrinsic limitations on  accuracy.However, this conclusion may not discourage us from the necessary task of modeling in order to provide the environmental authorities with a tool that nevertheless will help to protect the health of citizens.The situation described here may apply in other places such as mid-sized cities, where the complexity of topography, multiple emission sources, seasonal trends and meteorology impose limitations to the estimation of future pollutant concentrations.

Fig. 1 .
Fig. 1.Location of PM 2.5 monitoring stations in Santiago urban area.

Fig. 2 .
Fig. 2. One hour average of PM 2.5 concentrations measured in Santiago monitoring stations during 2006.(a) Station M, (b) station L, (c) station N, (d) station O.

Table 1 .
Statistics of PM 2.5 series for years 2000-2006 as measured in four monitoring stations.Data from January through December.

Table 2 .
Statistics of PM 2.5 series for years 2000-2006 as measured in four monitoring stations.Data from April through August.

Table 3 .
Lyapunov exponents for PM 2.5 annual and winter time series in four monitoring station in Santiago, Chile.
Fig. 3. Logarithm of Euclidian distance between successive five dimensional vectors constructed from 2007 PM 2.5 time series in station O.

Table 4 .
Statistics of PM 2.5 series for years 2000-2006 as measured in four monitoring stations.Data from April through August are considered.