Research Journal of Chemical Sciences ______________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.Sci. International Science Congress Association 1 Dispersion Modeling of Jaipur Fire, India Biswal D.K., Kumar V. and Barik K. Department of Chemical Engineering, Indian Institute of Technology Roorkee, Roorkee, INDIA Department of Chemical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, INDIAAvailable online at: www.isca.in (Received 13th July 2011, revised 16th January 2012, accepted 19th January 2012)Abstract The continental model equations take the form of partial differential equations which are very difficult to solve analytically. In the present work emphasis is given on the modeling of the horizontal transport of pollutant by taking basic assumptions for the wind velocity, mass burning rate of pollutant from the source etc. Model equations have been developed by taking suitable control volume and are solved by using finite volume method. MATLAB Simulation results give the idea about concentration profile of a pollutant along X and Y co-ordinate of horizontal direction in different time period. Then the model has been used to predict the concentration of SPM at various meteorological stations and is compared with the observed concentration in those stations for the Jaipur Indian Oil Corp. Fire Episode, 2009. ArcGIS is used to determine the geographical locations and distance between meteorological stations for this purpose. Keywords: Continental air pollution,horizontal transport, SPM, ArcGIS, finite volume, Jaipur. Introduction Sometimes traces of air pollutants are found thousands of kilometers away from their sources of generation. The long range transport of air pollutants in the atmosphere sometimes called as continental air-pollution. Long range transport of pollutants is significant in case of forest fire, post-harvest field burning, industrial disaster etc. A typical case of continental air pollution found during the Indian Oil depot fire, Jaipur episode on 29th Oct. 2008 when Earth Observatory satellite captured the picture of smoke flowing pattern shown in figure 1. Figure-1 Direction of Pollutant Transport, IOC depot. Fire, 30thNov. 2009, Sitapura, Jaipur So it is therefore important to develop a mathematical model in order to find out the concentration of the particular pollutants in the atmosphere with respect to time and distance on continental scale. By using the model the approximate concentration of the pollutant may be predicted which may be beneficial if the predicted concentration is higher than the allowable limits of pollutants in the atmosphere set by Pollution Control Board. Predicting the concentration of pollutant before pollutants reach a place may save life and property. In other words mathematical modeling of continental air pollution is necessary. Considerable amount of work has been done in this area of Long range transport of pollutants. The very first model called as Danish Eulerian Model (DEM) has been initiated early of 80’s. The 2-Dimensional dispersion model was developed during the study of the sulphur pollution over Europe. 3-Dimensional advection-diffusion modeling was done for regional scale. Also Modified version of DEM was developed. The three dimensional DEM was split into three sub-models containing vertical exchange operator, a sub-model containing the horizontal advection and the horizontal diffusion operators and a sub-model containing the chemistry, the emissions and the deposition operators. The sub-models were to be treated by appropriate numerical methods. The finite element method is used to discretize the spatial derivatives in the vertical exchange sub-model. The semi-discretized model can be handled by the -method6,7. The chemistry-emission-deposition sub-model can be treated by using partitioning algorithm. The most suitable methods for the solution of Danish Eulerian Model (DEM) are pseudo-spectral algorithm and finite element approach. Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 2 Material and Methods Mathematical Modeling: Five physical and chemical processes have to be described by mathematical terms in the beginning of the development of a continental air pollution model. These processes are: horizontal transport (advection), horizontal diffusion, chemical transformations in the atmosphere combined with emissions from different sources, deposition of pollutants to the surface, and vertical exchange (containing both vertical transport and vertical diffusion).Here the modeling is intended to be done for the horizontal transport of pollutants followed by horizontal diffusion. Chemical transformations and chemical deposition of the pollutants in the atmosphere are ignored while doing the modeling process. There will be some transport of mass and diffusion of pollutants vertically which are ignored in this modeling process. However vertical thickness i.e. thickness of inversion layer (fixed value) is considered while doing the component mass balance along and across the control volume to get the model equations. Basic Assumptions and Model Formulation: The following assumptions apply to the model are described here: i. Transportation of pollutants takes place only in horizontal direction, i.e. in X and Y direction only. ii. Wind velocity has only X component having fixed value in the entire region under consideration, i.e., X axis is taken in the direction of wind. iii. Diffusive mass flux in X-direction is negligible as compared to convection. iv. There is only diffusion and no convection in Y-direction. v. There is an even mass burning rate at the source for all time‘t’. vi. The inversion layer is assumed to be constant having a constant thickness all over the pollutant transport. vii. Model equations are formed doing component mass balance assuming a uniform concentration in a volume element, equal to the concentration at the central point. Component mass balance can be done for the volume element at (i, j) drawn below. \n \r\r \n\r\r  \r\r\r\r   \r\r\r\r ! \r\r# ! 1 Taking the limiting case of equation 1; # " $  % 2 Equation 2 represents the Navier-Stokes equation for horizontal transport of pollutants. Boundary Condition in X-Axis: Doing component mass balance for the above volume element (i, j) where j=1 at x-axis; \n \r\r \n\r\r '  \r\r\r\r  \r\r\r\r(! \r\r# ! 3 Boundary Condition in Y-Axis: Component mass balance at point (i, j) of the above volume element can be written as;  \r\r  \r\r\r\r   \r\r\r\r !"  \r\r# ! 4 Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 3 Computational Technique: The procedure for solving the 2-D partial horizontal continental modeling differential equations has been presented here. Finite volume method is an effective method for finding the solution of Partial differential equations. Formulation of Modeling Equations for Horizontal Transport: Figure-2 shows how the horizontal grids are formed along X-Y plane for discretization process to be performed to solve the equations. Each small solid box represents the top view of each volume element. Each circle at the center of each volume element represents the node point where the value of Conc. can be found out for each time increment. Here ‘i’ and ‘j’ represent the grid points in X and Y direction respectively. In X-direction ‘i’ varies from 1 to m (number of node points in X-direction) and in Y-direction ‘j’ varies from 1 to n (number of node points in Y-direction). Figure-2 Finite volume grid patterns This particular grid represents the variation of conc. in each node points for a particular time step. A new grid of nodes can be obtained for each increment of time step. From equation 1 \n \r\r \n\r\r  \r\r\r\r   \r\r\r\r !" \r\r# ! 5 Dividing by (!) and discretizing and then rearranging;  \n\r\r ") ) ! *  \n\r\r! \n \r\r \n\r \r \n\r \r 6 Formulation of Equation for the Boundary Condition in X-direction: From equation 3 the component mass balance equation along X- axis is given as; \n \r \r \n\r \r '  \r\r\r\r  \r\r\r\r(!=  \r\r# ! 7 Discretizing and rearranging; \n\r \r ") # $ ,  \n\r \r$ \n \r \r \n\r \r 8 In this equation C (i, j-1, k) is same as C (i, j+1, k) due to symmetry along Y-axis. So equation 8 can be written as \n\r \r ") # $ ,  \n\r \r$ \n \r \r, \n\r*\r 9 Formulation of Equation for the Boundary Condition in Y-direction: From equation 4 the component mass balance equation along Y- axis is obtained as;  \r\r  \r\r\r\r   \r\r\r\r !" \r\r# ! 10 Discretizing and rearranging;  \r\r ") # $ ,  \r\r  \r \r  \r \r 11 Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 4 Results and Discussion Finite volume method is used for making and solving the governing model equations. If the initial source concentration of pollutant is known then it is desired to compute the concentration of a particular pollutant on the X-Y plane on the direction of wind at different time intervals. So by taking a fixed value of wind velocity from a definite direction and a fixed value of diffusivity for the specific pollutant during the transport; we can calculate its concentration at different time period on the X-Y plane. It is also important to mention that there are two parameters i.e. wind velocity and diffusivity which play a great role in determining the nature of the concentration profile. Wind velocity can vary depending upon stability conditions of the atmosphere. In general it can take the values from 2 m/s to &#x-3.3;女 6 m/s. Similarly value of diffusivity takes a range of values depending upon the nature of pollutant and atmospheric temperature conditions. Value of molecular diffusivity of pollutants can take values in the order of 10-5 /sin different temperature conditions in the atmosphere whereas the value of Eddy diffusivity can take values in a wide range. It is therefore important to find the wind velocity and the diffusivity of pollutant which can be substituted in the model equations to get the exact concentration distribution. Concentration Distribution of Suspended Particulate Matter (SPM): Figure-3 represents a 3-D plot for a particular pollutant distribution along the Grid [20000, 50]. In x-axis the total distance (20000 m) is divided by 100 grid points and the Y-distance (50 m) is divided by 10 grid points. So each grid point has the value of 200 m in X-direction and 5 m in Y-direction. Time at which the above concentration distribution is obtained is 1 hour. Figure-3Concentration Distribution of Pollutant in X-Y Grid (xmax= 20000 m, ymax=50 m, d=.00013)With a wind velocity of 4 m/s and with time 1 hour, the pollutant would have reached approximately a distance of 14400 m in the direction of wind. So in the figure 3 with an xmax of 20000 m, it is showing a sharp decrease in concentration after 14400 m. Also from figure it is clear that diffusivity has very less effect on the concentration of pollutant as it has got very less value compared to the wind velocity. Model Validation: The Indian Oil Corporation terminal at Sitapura, Jaipur (Rajasthan) was involved in the activity related to receipt, storage and dispatch of petroleum products such as HSD, MS, SKO, and Lubes. A major accident, involving MS vapor cloud explosion and fire, occurred at this facility, on 29th Oct. 2009 between 1830 hrs to 1935 hrs. The fire continued till 0600 hrs on 11th Nov. 2009. After the outbreak of fire in the IOC terminal at Jaipur, the Rajasthan State Pollution Control Board initiated efforts for environmental monitoring and assessment. A quick assessment was carried out by the RSPCB officials that took into consideration the wind direction and the probable impact areas including the population distribution. Accordingly, it was decided to immediately set up 5 air quality monitoring stations at varying distances and in different directions around the IOC terminal. These stations (which were operational from 30th Nov. onwards) were located at Pratap Nagar P.S., STP Delawas, Mathurawala village, Shivdaspura P.S., and Chokhi Dhani. The Central Pollution Control Board also sent a team of officials for the assessment and monitoring of environmental parameters. They have compiled a report entitled ‘The Indian Oil Corporation Ltd. depot (Jaipur terminal of marketing division) fire episode, (29th Oct. 2009 – 11th Nov. 2009), A Report’ (Appendix I) 10. Meteorological Data and Satellite Imagery: An analysis of the satellite observations taken during the fire incident was carried by the State remote Sensing application center, Department of Science and Technology, Govt. of Rajasthan, Jodhpur. Figure 4 is the satellite pictures taken on 30th Oct. showing the direction of pollutant transport. A set of meteorological data and pictures taken by Satellites were collected during and after the incident happened in IOC depot, Sitapur. For analysis purpose, the peak fire incident was taken to be from 30th Oct. to 3rd Nov. so as to be consistent with the air quality data analysis. For Analysis purpose the wind rose of the day 30th Oct. is given by the figure5. From the figure-5, it is clear that the resultant of wind direction on 30th Nov. was towards SE and SSE. Wind speed data analysis shows that on 30th Oct., the average wind speed was around 2.5 m/s. Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 5 Figure-4 Direction of Pollutant transport Figure-5 Wind rose of 30th Nov. 2009, Sitapura Mapping using ArcGIS: Figure-6 represents an ArcGIS figure which shows the geographical location of Sitapura which is around 12 km from Jaipur City. Meteorological stations Chokhi Dhani, Mathurawala Village and Shivdaspura are situated at a distance of around 2 km, 5 km, and 10.5 km, respectively from Sitapura. Concentration Calculation by using Molecular Diffusion: In the present study it is desired to calculate the average concentration of SPM using the continental model in three meteorological stations namely Chokhi Dhani, Mathurawala Village and Shivdaspura stations 24 hour right after the fire broke out at IOC, Sitapur. The reason of choosing these three meteorological stations over other stations for the calculation is due to the way these stations are situated in the map through which maximum transportation of smoke occurred in the peak time of the fire incident. Concentrations of SPM are calculated in these stations assuming the average SPM concentration at Sitapur as 300 µg/m during the peak of the fire incident. Figure-7 is meant to calculate the concentration of SPM at Chokhi Dhani meteorological station. While predicting the concentration, molecular diffusion was taken into account (zero in this case for SPM). By taking a 2000 m grid in wind flow direction and assuming the SPM concentration as 300µg/m, same concentration at Chokhi Dhani was observed as that of Sitapura. Figure-8 represents the concentration profile to calculate concentration of SPM at Shivdaspura meteorological station. By taking a 10500 m grid in wind flow direction and assuming 300 µg/m SPM concentrations at Sitapura, same concentration was observed as shown in figure. The reason why a straight concentration profile was observed both for Chokhi Dhani and Shivdaspura was that while calculating the concentration in these places diffusivity was neglected as molecular diffusivity of SPM is zero. So it should be a straight line if there would be no diffusion. Also same type of concentration profile was observed while calculating it in Mathurawala Village. Concentration Calculation by using Eddy Diffusion: Again by assuming a range of eddy diffusivity values concentrations of SPM were calculated at all the meteorological stations. Eddy diffusivity takes a higher value in case of local air pollution as that of continental air pollution. So a range of values in the order of 10-3 to 10-5were taken and calculations have been done. SPM concentration at Sitapur was assumed 300 µg/m for the calculation. Figure-9 is meant to calculate the concentration of SPM at Chokhi Dhani meteorological station. By taking a 2000 m grid in the wind flow direction and assuming eddy diffusivity to be 0.0013 m/s the concentration at Chokhi Dhani was found to be 280 µg/mFigure-10 is meant to calculate the concentration of SPM at Mathurawala Village meteorological station. By taking a 5000 m grid in wind flow direction and assuming diffusivity to be 0.0013 m/s the SPM concentration was found to be 250 µg/m. From figure-11 concentration of SPM at Shivdaspura was found to be 200 µg/m. Observed Concentrations at Meteorological Stations:Figure-12 shows the observed average concentrations which were calculated at five different meteorological stations using appropriate measuring device during and after fire. From figure the average observed SPM concentration at Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 6 Chokhi Dhani, Mathurawala Village and Shivdaspura were found to be around 280, 260 and 260 µg/m, respectively. Comparision: From figure-13 it is clear that in the absence of any diffusivity the concentration profile along the wind flow direction shows a constant value whereas, in the presence of eddy diffusivity it falls along the length of the pollutant transport because of the diffusion along Y direction. Predicted concentration at Chokhi Dhani and Mathurawala Village matches with the field value of concentration. And at Shivdaspura concentration shows a lesser value as compared to predicted value. It may be due to the sedimentation of pollutants on earth surface during the transport which was not been accounted for during the development of the model. Figure-6 Mapping using ArcGIS Figure-7 Concentration Distribution of SPM at Chokhi Dhani (xmax= 2000 m, ymax=50 m, d=0, u=2.5 m/s) Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 7 Figure-8 Figure-9 Concentration Distribution of SPM at Shivdaspura Concentration Distribution of SPM at Chokhi Dhani (xmax= 10500 m, ymax=50 m, d=0, u=2.5 m/s) (xmax= 2000 m, ymax=50 m, d=0.0013m/s,u=2.5m/s) Figure-10 Figure-11 Concentration Distribution of SPM at Mathurawala Village Concentration Distribution of SPM at Shivdaspura (xmax= 5000 m, ymax=50 m, d=0.0013 m/s, u=2.5 m/s) (xmax= 10500 m, ymax=50 m, d=0.0013 m/s, u=2.5 m/s) Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 8 Figure-12 Average values of Concentrations of SPM, RSPM, NO and SO at each of the five stations during/after the fire episode Figure-13 Comparision between observed and predicted value of Concentration at Chokhi dhani, Mathurawala Village and Shivdaspura Research Journal of Chemical Sciences __________________________________________________________ ISSN 2231-606X Vol. 2(2), 1-9, Feb. (2012) Res.J.Chem.SciInternational Science Congress Association 9 Conclusion Following conclusions may be drawn after the present study: With the help of wind rose and satellite imagery the nature of pollutant transport in the atmosphere can be predicted. Dispersion of pollutants in lateral direction is less than that of predicted by assuming molecular diffusion. The movement of pollutants can be modeled by assuming various parameters under consideration. A much higher value of coefficient of diffusivity as compared to molecular diffusivity indicates presence of micro eddied in apparently laminar flow of wind. Finite volume method is very suitable for describing each physical process inside a control volume unlike other methods where it is very difficult to understand the physical process. Discretization using finite volume method gives significant solutions along the length and breadth of the essential grid under consideration in various time spans. The equations developed by finite difference method are easily programmable and the Simulation of equations gives considerable results for concentration calculation for different pollutants. MATLAB program is written in such a way that it can solve the model equations for horizontal grid of any length with minimum possible time. Geographical locations of the stations can be well viewed by using ArcGIS. Present study is desired to calculate the concentration in horizontal direction of flow of pollutant. In future any other method may be effectively used for the calculation of concentration in vertical direction as well as chemical deposition and transformation of pollutants. ArcGIS may also be used to describe the dispersion pattern in the smoke affected areas. Acknowledgement The authors would like to thank Ministry of Human Resource Development, for funding this research in form of, MHRD postgraduate scholarship. The authors would also like to thank department of Chemical engineering and Institute Computer Center of Indian Institute of Technology Roorkee, India for providing facilities required for carrying out this research work. References 1.Zlatev Z., Mathematical model for studying the sulphur pollution over Europe, Computational and Applied Mathematics, 12, 651-666 (1985)2.Zlatev Z., Berkowicz R. and Prahm L.P., Three dimensional advection-diffusion modeling forregional scale, Atmospheric Environment, 17, 491-499 (1983)3.Zlatev Z., Dimov I. and Georgiev K., Three- dimensional version of the Danish Eulerian Model,Zeitschrift für Angewandte Mathematik und Mechanik.76, 473-476 (1996)4.Zlatev Z., Computer Treatment of Large Air Pollution Models, Kluwer Academic Publishers, Dordrecht, Boston, London (1995)5.Pepper D.W., Kern C.D. and Long P.E., Modeling the dispersion of atmospheric pollution using cubic splines and chapeau functions, Atmospheric Environment, 13, 223-237 (1979)6.Havasi A., Zlatev Z., Trends in Hungarian air pollution levels on a long time-scale, Atmospheric Environment, 36, 4145-4156 (2002)7.Samarskii A.A., Theory of difference schemes, Moscow: Nauka (1977)8.Alexandrov V., Sameh A., Siddique Y. and Zlatev Z.Numerical integration of chemical ODE problems arising in air pollution models Environmental Modeling and Assessment (1997)9.Zlatev Z., Christensen J. and Eliassen A., Studying high ozone concentrations by using the Danish Eulerian Model, Atmospheric Environment., 27, 845-865 (1993)10.A report on Environmental impacts of the fire in Indian Oil Corporation Depot, Sitapura, Jaipur, Submitted to Department of Environment, Govt. of Rajasthan, February (2010)