Quantification of Non-Linear Dynamics and Chaos of Ambient Particulate Matter Concentrations in Muzaffarabad City

The present study was carried out for quantification of non-linear dynamics and chaos of ambient particulate matter (PM) concentrations in Muzaffarabad city. PM1.0 and PM2.5 concentrations were monitored at six different locations for a continuous 6 h period. The linear behavior of the acquired time series data was analyzed using descriptive statistics. Specifically, the chaotic temporal behavior of the PM concentration was analyzed using phase space reconstruction, the Hurst exponent, and the largest Lyapunov exponent. The average mutual function was used to calculate proper time delay, while the false nearest neighbor method was used to calculate the proper embedding dimension for the phase space reconstruction. No health-protective quantitative standard exists for PM1.0 concentrations. However, the mean concentration of PM2.5 was considerably higher than the standards, developed by WHO, US-EPA and European Union directives, at all six locations. For PM1.0 minimum, 293 ± 149 μg m, and maximum, 544 ± 490 μg m, values were recorded at CMH Chowk and Chehla Bridge locations respectively. For PM2.5 minimum, 394 ± 262 μg m and maximum, 633 ± 426 μg m, values were recorded at CMH Chowk and at old secretariat respectively. The results indicates PM1.0 and PM2.5 concentration time series can be modeled using phase space reconstruction by properly selecting the embedding parameters. The positive largest Lyapunov exponent reveals a strong chaotic signature in the system dynamics for both particulate sizes. Furthermore, Hurst exponent for both particulates was close to 1, showing highly persistent time series pattern.


INTRODUCTION
Particulate matter (PM) is one of the key global pollutants, affecting human and environmental health, due to its ability to travel across the countries and geographic boundaries.The concentration of PM in outdoor air depends on various factors including the amount of local motor traffic, the types of industries present, the dominant approaches to power generation such as domestic coal or biomass fuel combustion, the topography and the weather.Epidemiologic research on long-term exposure to fine particulate matter (FPM) has revealed adversarial health effects, including increased and premature mortality in patients suffering from heart and respiratory disease, lung cancer, aggravated asthma, and adverse reproductive outcomes.Particulate matter (PM) is a heterogeneous and complex mixture of enormously small liquid droplets and solid particles with a cut point smaller than hundred microns (Sioutas et al., 2005;Traversi and Marco, 2008).Evidence suggests that toxicity of PM originates from its physical size rather than from chemical composition (Ferin et al., 1992;Li et al., 1996).Based on the size, particulates are categorized as coarse PM (CPM) with size between 2.5 and 10 µm, or fine PM (FPM) has size 2.5 µm or less.Acute PM 2.5 and PM 10 exposures are specifically connected with increased probability of ventricular ectopy for smokers, suggesting that smokers are more vulnerable to the arrhythmogenic impacts of PM (Liao et al., 2009).Further studies are important to segregate between the acute and chronic impacts.
The FPM are more hazardous than CPM because of their ability to more deeply penetrate human lungs and enter blood, which may increase respiratory and cardiovascular diseases (Schwartz et al., 1996;Hughes et al., 1998;Brook et al., 2004).Fine particulate matter of sizes 0.1-1 µm may reduce visibility, increase surface water acidity, reduce soil nutrient levels, damage forests and crops, diminish ecosystem diversity, alter nutrient balance in river basins and coastal water, damage stone and corrode metals and other materials (Amos et al., 2010).
Motorized transportation vehicles make a very significant increase, 5-80%, in airborne PM concentration not only because of engine exhaust but also from abrasion of tires and brake components and road surface wear (Pallavi and Roy, 2013).FPM dispersion near urban roadways has, in fact, become an issue of serious concern due to the especially hazardous health impacts due to the inhalation of small particulates.Along urban roadways, the determinants of FMP concentration include wind direction and the number of buses and cars with buses with exhausts below the bus having the greatest impact upon pedestrians (Balogh and Mannering, 1992).
In the present study, ambient PM (PM 1.0 and PM 2.5 ) data were collected at six different sites for a continuous period of six hours in the suburb of Muzaffarabad, the Capital of Pakistani held Kashmir.Nonlinear time series analysis techniques were used to characterize chaotic behavior of the time series data.To capture nonlinear dynamics, phase space was reconstructed using an appropriate time delay and embedding dimension.The Lyapunov exponent was computed to determine the evidence of deterministic chaos in the ambient PM time series data.The Hurst exponent is used to explore whether or not the FPM time series data show persistent behavior.

Study Area
Muzaffarabad is the capital city of Pakistani held Kashmir commonly known as Azad Jammu and Kashmir (AJ&K).It is located at longitude 73.47°E, latitude 34.37°N and 737 m elevation above the sea level at the conjunction of Neelum and Jhelum rivers and is surrounded by lofty mountains.The Muzaffarabad district, with a total area of 6,117 km 2 , is bounded by the Khyber Pakhtunkhawa Province of Pakistan in the west, by the Kupwara and Baramula districts on the Indian side of the Line of Control in the east and Neelum district of Azad Kashmir in the North.The population of Muzaffarabad district is 15% of the overall population of Azad Kashmir, reaching 545,817 in 2011 (Khan, 2012).

Data Acquisition
In this study, the ambient air concentrations of FPM, both PM 1.0 and PM 2.5 , were monitored at six different sites along roadways using an environmental particulate air monitor (Haz-Dust EPAM-5000 Environmental Particulate Air Monitor, Environmental Devices Corporation, 4 Wilder Drive, Building 15 Plaistow, New Hampshire 03865) within the city boundaries.Monitoring locations of FPM are shown in Fig. 1.The measurement device used was a high sensitivity real-time particulate monitoring device that uses an infrared light source positioned at a 90° angle from a photo detector.The airborne particles scatter the infrared beam light so that the amount of light received by the photo detector is directly proportional to the aerosol concentration.It was calibrated using standard Arizona Road Dust (ARD) with NIOSH method 0600 (NIOSH, 1998) with a precision of ± 0.03 mg m -3 and an accuracy of ± 10% for nonrespirable dust.
At each location the device was installed for six consecutive hours to monitor the concentration of fine particulates at a sampling rate of 1 Hz.After acquiring the data, device was returned to laboratory and data transferred to a personal computer with the help of manufacturer provided software.

Non Linear Techniques Basic Review of the Mathematics of Chaotic and Timevarying Systems
Some systems, called chaotic systems, change in very complex ways as a function of time but are deterministic in the sense that if their initial conditions are precisely well known then their state could be exactly predicted.Unfortunately, it is impossible to know their initial or current conditions well enough to perfectly predict their future state because small differences in the initial conditions may cause radically different behaviors in such chaotic systems.The systems, however, are predictable for certain characteristic time periods.The Lyapunov time is a mathematical characteristic of a specific chaotic system under consideration, and predictability is possible for times that are less than a few times the Lyapunov time.If the chaotic systems are very sensitive to the initial conditions then the current or initial conditions must be very well known to make accurate predictions.This sensitivity is measured using the Lyapunov exponent.The Hurst exponent is another property of temporally variable systems that relate to the tendency of the apparent randomness or correlation of the behavior of a temporally variable system over long time periods.These basic have been previously applied to characterized radon concentration, an environmental pollutant sensitive to prevailing local environmental conditions (Shafique et al., 2014), and were applied in the work to the problem of PM in the same portion of the planet.

Phase Space Reconstruction
The method of phase space reconstruction is employed to analyze the physiognomies, or characteristics, of the strange attracters that are functional topological relationships that characterize a chaotic system.The most commonly used method of phase reconstruction is the delay-embedding theorem (Takens, 1981).According to this method, the phase space from the time-series{x(t i ), I = 1,2, 3,…., N} can be reconstructed as follows: where, τ is the time delay and m is the embedding, it is vitally important to suitably select pair m and τ for proper phase space reconstruction.Numerous methods, such as auto correlation function and average mutual information are available in the literature for the estimation of time delay τ.In this study, the average mutual information was used for estimation of τ, which accounts for nonlinear correlations in the time series data.The average mutual information for a time series with consecutive values x i and x i+ τ is given as: where, p(x i ) is the probability of occurrence of i th interval, p(x i+τ )is the probability of occurrence of (i + τ) th interval and p(x i , x i+τ ) is the associated probability that an observation falls into i th and (i + τ) th interval.Finally, I(τ) is plotted against varying τ for estimating suitable value of τ, which is the location of first local minimum when I(τ) decreases non-monotonically with increasing τ (Fraser, 1986).For monotonous decrease of I(t), either the decrease of average mutual information to I(t)/I(0) = 0.2 (Abarbanel, 1993) or (Horak, 2003) can used as the time delay.For a chaotic time series, the embedding dimension can be obtained by several methods (Grassberger and Procaccia, 1983;Abarbanel, 1993;Kantz and Schreiber, 1997).For this work, a false nearest neighbor approach was used for estimating optimal embedding dimension (Kennel et al., 1992).To identify false nearest neighbor, Kennel et al. (1992) suggested following paired criteria be jointly applied: where, δ n is the relative increase in Euclidean distance when dimension of reconstructed state space is increased from m to m + 1.This is computed as: The parameters R tol and A tol are constant thresholds and R A is the standard deviation of a time series.Kennel et al. (1992) suggested a value greater than ten for R tol but provided no suggestion for A tol .

Lyapunov Exponents
The Lyapunov exponent (λ) characterizes the nature of time evolution of close trajectories in the phase space and is considered as the key component chaotic behavior.Various approaches have been proposed to compute Lyapunov exponent (Rosenstein et al., 1993;Wolf et al., 1985).In the present study, Wolf's algorithm is used to calculate the largest Lyapunov exponent (Wolf et al., 1985).This algorithm involves embedding the time series in m-dimensional space and monitoring the divergence of trajectories initially close to each other in that phase space.The Lyapunov exponent can be defined as: Negative Lyapunov exponents are characteristic of dissipative or non-conservative systems, whose trajectories approach a common fixed point.A Lyapunov exponent of zero indicates that the system is in some sort of steady state mode and the trajectories maintain their relative positions.The positive Lyapunov exponent is usually taken as an indication that the system is chaotic.

Hurst Exponent
The Hurst exponent provides a quantitative measure of the persistence of a time series data.In this study rescaled range (R/S) analysis is used to calculate the Hurst exponents (Jens, 1998).The R/S statistic is the range of partial sums of deviations of times series from its mean, rescaled by its standard deviation.The R/S analysis statistics for a discrete time series data x i is given as: where The Hurst exponent, H, was evaluated by plotting the log(R/S) against log(N) and measuring the slope of the straight line.The value for H in the range 0.5-1.0carries the signatures of persistence behavior.The Hurst exponent is related to the fractal dimension D of the time series curve by the formula D = 2 -H.

RESULTS AND DISCUSSIONS
The statistical parameters including mean value, standard deviation (STD) and 95% confidence interval (CI) for PM 1.0 and PM 2.5 , at selected sites in the downtown city of Muzaffarabad, are listed in Table 1.Each of values listed for PM 1.0 and PM 2.5 , at specific locations, are averages of 21,600 data points.For PM 1.0 minimum value equal to 293 ± 149 µg m -3 was recorded at CMH Chowk having a 95% CI of 291-295 µg m -3 , while the maximum, 544 ± 490 µg m -3 , corresponded to the Chehla Bridge location with a 95% CI of 537-550 µg m -3 .

Phase Space Reconstruction
The selection of time delay τ and optimum minimum embedding dimension m is vitally important for phase space reconstruction.The average mutual information, AMI, a generalization of autocorrelation function is used for determination of time delay τ.In Fig. 2, the AMI is plotted against varying τ, for optimal determination of its value for PM 1.0 and PM 2.5 .It evident from the figure, AMI decreased monotonically with increasing τ, so I(t)/I(0) = 1/e (Horak, 2003) was used as the criteria for determination of optimal time delay.The optimal time delay, τ, for particulates PM 1.0 and PM 2.5 using this criterion are 17 and 18 respectively.
The false nearest neighbors (FNN) approach is used to find the minimum embedding dimension m.For any given embedding dimension m, the proportion of the identified FNN to all the neighbors was computed for the given time delay τ.The percentages of the FNN are plotted as function of the embedding dimension.A zero FNN percentage indicates the optimum, i.e., minimum, embedding dimension.In Fig. 3   the results of the FNN approach for determining the optimum embedding dimension of N = 21,600 data points of PM 1.0 and PM 2.5 using various values for the threshold parameters R tol and A tol are shown.In this study, the value of the parameter R tol was varied from 20 to 100 with a step size of 20 and A tol = 0.2 × R tol was used.The results summarized in Fig. 3 showed that the optimum embedding dimension obtained by the FNN approach varies from 5 to 8 depending on the values of thresholds R tol and A tol .The higher values of optimal embedding showed that the concentration time series of PM 1.0 and PM 2.5 have dominant degrees of freedom, indicating dominance of some complexity in their concentrations respectively.The complexity in PM concentrations may be due to the various influencing factors such as temporal variations, humidity, temperature, wind directions and speed, smoke produced from automobiles.

Lyapunov Exponents
The Lyapunov exponent, λ, characterizes the nature of time evolution of close trajectories in the phase space and is considered as the key component of chaotic behavior.In the present study, Wolf's algorithm was used to calculate the largest Lyapunov exponent (Wolf et al., 1985) using the same embedding parameters as before.In Table 2, the value of the largest Lyapunov exponent λ of PM 1.0 and PM 2.5 at different sites along with their mean and standard deviation PM 1.0 and PM 2.5 values are presented.The positive values of λ indicates an exponential divergence of the trajectories and hence a strong signature of chaos for both particulate sizes.

Hurst Exponent
In this study rescaled range (R/S) analysis is used to calculate the Hurst exponents (Jens, 1998).The R/S statistic is the range of partial sums of deviations of times series from its mean, rescaled by its standard deviation.The Hurst exponent H is evaluated by plotting the log(R/S) against log(N) and measuring the slope of the straight line and then fractal dimension D for the for both particulates was obtained by the formula D = 2 -H.In Fig. 4, the log(R/S) plotted against log(N), for particulates PM 1.0 and PM 2.5 at six sites are shown.
The results for Hurst exponent and fractal dimension are summarized in Table 3. Hurst exponent obtained at different locations, in time series data for particulate matter PM 1.0 , ranged in the interval 0.85 ≤ H ≤ 0.92 and those for PM 2.5 ranged in the interval 0.83 ≤ H ≤ 0.92.It is evident from the findings that values of obtained Hurst exponent for both PM 1.0 and PM 2.5 are very close to '1' showing highly persistent time series pattern.

CONCLUSIONS
In this study, the concentrations of FPM, specifically PM 1.0 and PM 2.5 , in ambient air were monitored at six different sites along roadways in the downtown city of Muzaffarabad.The spatial behavior of PM 1.0 and PM 2.5 was analyzed using statistical and nonlinear times series analysis techniques.The statistical measures showed alarmingly high concentration of both PM 1.0 and PM 2.5 in the vicinity of Muzaffarabad.The results indicate that PM 1.0 and PM 2.5 concentration time series can be modeled using phase space reconstruction by properly selecting the embedding parameters m and τ.The concentration time series of both particulates indicated dominance of some complexity in their concentrations as depicted by higher values of optimal embedding dimension.The positive largest Lyapunov exponent depicts strong chaotic signature in the system dynamics of both particulates.Furthermore, Hurst exponent for both particulates was close to 1, showing highly persistent time series pattern.The study can have implications for estimating spatial or temporal future patterns of particulates in the ambient and indoor environment, which can be helpful for long term forecasting of these particulates.Furthermore, the study can be helpful for the Government and environmental protection agencies to protect the citizen from hazardous effect of atmospheric pollutants and for air pollution management.We have not considered the effects of wind speed and its direction, temperature and chemical characterization of PM and gases on the ambient concentrations of Particulates.It is suggested to carry out an exhaustive study in future that can address above mentioned limitations of current study for better understanding of PM concentrations over long period of time.

Fig. 1 .
Fig. 1.Site map of six monitoring locations taken from Google earth.

Fig. 2 .
Fig.2.The optimal time delay τ, for particulates PM 1.0 and PM 2.5 using AMI method.

Fig. 3 .
Fig. 3. False nearest neighbor curves with various thresholds for determining the embedding dimension of particulates PM 1.0 and PM 2.5 .

Table 1 .
, Statistical parameters (Mean ± STD and 95% CI) of the particulate matter concentrations at selected sites.

Table 2 .
LLE Values, N = 21600 for all the sites.

Table 3 .
Hurst and FD values for all the selected sites.