VIEWS: 1 PAGES: 18 CATEGORY: Emerging Technologies POSTED ON: 3/20/2013
International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com , Volume 2, Issue 2, February 2013, ISSN 2319 – 4847, ISRA Journal Impact Factor: 2.379
International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 Advection - Diffusion numerical model of air pollutants emitted from an urban area source with removal mechanisms by considering point source on the boundary Lakshminarayanachari K.1, *Sudheer Pai K.L.2 , Siddalinga Prasad M.3 and Pandurangappa, C.4 1 Department of Mathematics, Sai Vidya Institute of Technology, Bangalore - 560 064, INDIA 2 Principal, R N S First Grade College, Bangalore - 560 098, INDIA 3 Department of Mathematics, Siddaganga Institute of Technology, Tumkur - 572 103, INDIA 4 Department of Mathematics, Raja Rajeswari College of Engineering, Bangalore - 560 074, INDIA *corresponding author ABSTRACT A two-dimensional mathematical model of primary and secondary pollutants of an area source with chemical reaction and dry deposition by considering point source on the boundary is presented. One of the important atmospheric phenomena is the conversion of air pollutants from gaseous to particulate form. The primary pollutants which are emitted directly to the atmosphere are converted into secondary pollutants by means of chemical reaction. The governing partial differential equation of primary and secondary pollutants with variable wind velocity and eddy diffusivity is solved by using Crank-Nicolson implicit finite difference technique. Consistency, stability and convergence criteria have been tested for the numerical scheme used in this model. Concentration contours are plotted and the results are analyzed for primary and secondary pollutants in stable and neutral atmospheric conditions for various meteorological parameters and transformation processes. Keywords: Primary and secondary pollutants; chemical reaction; variable wind velocity; eddy diffusivity; Crank Nicolson implicit scheme. 1. INTRODUCTION Rapid industrialization and urbanization have posed a serious threat to human life and environment in recent years. The toxic gases and small particles could accumulate in large quantities over urban areas, under certain meteorological conditions. This is one of the serious health hazards in many of the cities in the world. Over the past three or four decades, there have been important advances in the understanding of the actions, exposure-response characteristics, and mechanisms of action of many common air pollutants. A multidisciplinary approach using epidemiology, animal toxicology, and controlled human exposure studies has contributed to the database. This review will emphasize studies of humans but will also draw on findings from the other disciplines. Air pollutants have been shown to cause responses ranging from reversible changes in respiratory symptoms and lung function, changes in airway reactivity and inflammation, structural remodeling of pulmonary airways, and impairment of pulmonary host defenses, to increased respiratory morbidity and mortality [1]. Tropospheric ozone is known to be an important pollutant and it has been established that ozone has significant impact on human health globally [2]. The problem of air pollution in cities has become so severe that there is a need for timely information about changes in the pollution level. Today forecasting of air quality is one of the major topics of air pollution studies due to the health effects caused by these airborne pollutants in urban areas during pollution episodes [3]. Atmospheric particulate matter (PM) in the micrometer size range follows a bi–modal distribution, with a coarse mode including particles produced by mechanical processes, such as soil dust, cloud droplets and many biological particles [4] and a fine mode dominated by both primary anthropogenic pollution from combustion processes and gas–to–particle conversion [5]-[7]. The primary pollutants which are emitted directly to the atmosphere is converted into secondary pollutants by means of chemical reaction. Study of such conversion process in which sulphur dioxide is converted to particulate sulfate, nitrogen oxide to particulate nitrate and hydrocarbons to particulate organic material would reveal a lot on urban plume characteristics [8]. The study of the secondary pollutants is important as the life period of the secondary pollutants is longer than that of the primary pollutants and it is more hazardous to human life and environment protection. Several studies have reported in which the downwind measurements of large urban complexes were carried out in order to obtain material balances on gaseous and particulate pollutants [9]-[12]. Acid precipitation can occur when the Volume 2, Issue 2, February 2013 Page 251 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 particles and the gases are removed from the atmosphere by (i) rain and snow (wet deposition) (ii) by impaction on water, soil and vegetation surfaces (dry deposition) and (iii) gravitational settling velocity. The transformation of the gases to the particles in the atmosphere and the deposition rates of the gases and the particles are a function of meteorological conditions over a time and the distribution of emission sources. In order to justify controls on the emissions of acid precursors, the relationship between a source and the deposition pattern that it produces needs to be understood [13]-[14]. The prediction of the concentration of pollutants is difficult to accomplish by field monitoring. Mathematical modeling is the only way to estimate the relative contribution of sources to the total deposition at a particular receptor over a long period of time. There is an interesting Lagrangian finite difference model, which has been developed using fractional step method by [15] to compute time dependent advection of air pollutants. In this model, the Eulerian grid used for the diffusion part of the pollutant transport equation remains unchanged. The finite difference scheme used avoids the numerical pseudo diffusion and it is unconditionally stable. The numerical solution obtained is compared with the corresponding analytical solution of steady state case and a reasonable agreement has been found [15]. In this paper he has given a brief discussion of related numerical schemes like puff-in-cell methods by [16] and conservation of momentum model by [17]. However, the work of [15] does not deal with the secondary pollutant or any kind of removal processes. Later on, [18] have presented an intermediate range grid model with point source for atmospheric sulpher dioxide and sulphate concentrations. This model deals with the atmospheric transport, diffusion, and wet and dry depositions over a region of several hundred kilometers. In this long range analytical model, some of the restrictive assumptions like constant velocity and eddy diffusivity profiles are imposed. [19] presented a two-dimensional analytical model for turbulent dispersion of pollutants in stable atmospheric layer with a quadratic exchange coefficient and a linear velocity profile. This model does not take into account any removal mechanisms. A numerical model is developed for the dispersion with chemical reaction and dry deposition from area source, which is steady state in nature [20]. All these area source models deal with only primary inert air pollutants. Subsequently, [21] presented a time-dependent numerical model for both primary and secondary pollutants in order to obtain time dependent contours of pollutant concentration in urban area. This model has been solved by using the fractional step method taking into account the specified functional form of vertical eddy diffusivity and velocity ( ) profiles. [22] presented a time dependent area source mathematical model of chemically reactive air pollutants and their byproduct in a protected zone above the surface layer with rainout/ washout and settling. But the horizontal homogeneity of pollutants and constant eddy diffusivity are assumed in this model. [23] has presented a time dependent two-dimensional air pollution model for both the primary pollutant (time dependent emission) and the secondary pollutant with instantaneous (dry deposition) and delayed (chemical conversion, rainout/washout and settling) removals. However his model being analytical, deals with the uniform velocity and eddy diffusivity profiles. Both [24] and [25] have presented a modeling of the urban heat island in the form of mesoscale wind and of its effect on air pollution dispersal. These numerical models deal with only primary inert air pollutants. A numerical model for the primary and the secondary pollutants with more realistic wind velocity and eddy diffusivity profiles by considering point source on the boundary is proposed in this paper. The effect of removal mechanism i.e, dry deposition velocity on primary and secondary pollutants is analysed. The model has been solved by using Crank- Nicolson implicit finite difference technique. The finite difference scheme used in this model avoids numerical pseudo- diffusion and it is unconditionally stable. Consistency, stability and convergence criteria have been tested for the numerical scheme used in this model. We have computed concentration residue, obtained after every time step against the number of time steps and analyzed. Accuracy depends on the fall in residue. It is seen that the residue settles to around 10-6. For the grid-independence study we have computed concentration for , , and grids and analyzed. The analysis reveals that concentration of the primary and secondary pollutants for and grids differ considerably against those on grids. Further, there is no perceptible change of values occurring on grids from that of grids. It is therefore reasonable to assume that the solution obtained on grids is an independent solution. The above computational cycle is then repeated for each of the next time levels. The steady state solution is obtained when the convergence criteria for the residual defined as is satisfied. Here, is the concentration which stands for both and , n refers to time, i and j refer to the space coordinates. 2. MODEL DEVELOPMENT The dispersion of chemically reactive pollutant concentration in a turbulent atmospheric medium using K-theory approach is usually described by the following equation [26]. Volume 2, Issue 2, February 2013 Page 252 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 where C is the pollutant concentration in air at any location (x, y, z) and time t; and are the coefficients of eddy diffusivity in the x, y and z directions respectively; U, V and W are the wind velocity components in x, y and z directions respectively and S is the source or sink of the air pollutants. The physical problem consists of an area source of an urban city with finite downwind and infinite crosswind dimensions with the point source (industrial stack) on the z-axis with and m. The grid points may miss the source because the source is at an arbitrary point. In this case the grid points have to be taken on the source point. To overcome this difficulty, one can think of the following two methods. One is to use Gaussian distribution for pollutants source at the initial line which is equivalent to the above point source and the other is distributing the point source to its neighboring two grid points. We have used the second procedure in this numerical model for air pollutants to take into account of a source at an arbitrary point on the z-axis. We have equally distributed the point source to its neighboring two grid points on z-axis. We assume that the pollutants are emitted at a constant rate from uniformly distributed area source. The Physical description of the model is shown schematically in Figure 1. Figure 1. Physical layout of the model We intend to compute the concentration distribution both in the source region and source free region till the desired distance = 12000 meters in the downwind. We have taken the primary source strength at ground level from an area source and the mixing height is selected as 624 meters. We have considered the point source strength as in this model. 2.1 Primary pollutant The governing equation of primary pollutant can be written as (1) where is the ambient mean concentration of pollutant species, is the mean wind speed in x-direction, is the turbulent eddy diffusivity in z-direction and k is the first order chemical reaction rate coefficient of primary pollutant . Equation (1) is derived under the following assumptions: The lateral flux of pollutants along the crosswind direction is assumed to be small i.e., where is the velocity in the y - direction and is the eddy-diffusivity coefficient in the y - direction. Horizontal advection is greater than the horizontal diffusion for not too small values of wind velocity, i.e., meteorological conditions are far from stagnation. The horizontal advection by the wind dominates over the horizontal diffusion, i.e., where and are the horizontal wind velocity and horizontal eddy diffusivity along x - direction respectively. Vertical diffusion is greater than the vertical advection since the vertical advection is usually negligible compared to the diffusion owing to the small vertical component of the wind velocity. We assume that the region of interest is free from the pollution at the beginning of the emission. Thus, the initial condition is (2) Where is the length of desired domain of interest in the wind direction and is the mixing height. We assume that there is a point source (industrial stack) on , i.e either at the beginning or at the end of the urban city. The wind direction is assumed towards the city from the point source. Thus Volume 2, Issue 2, February 2013 Page 253 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 (3) where z is taken as height. We assume that the chemically reactive air pollutants are emitted at a constant rate from the ground level and they are removed from the atmosphere by ground absorption. Hence the corresponding boundary condition takes the form (4) where is the emission rate of primary pollutant species, is the source length in the downwind direction and is the dry deposition velocity. The pollutants are confined within the mixing height and there is no leakage across the top boundary of the mixing layer. Thus (5) The term in equation (1) represents conversion of gaseous pollutants to particulate material as long as the process can be represented approximately by first-order chemical reaction. Bimolecular, Three-body, Photolysis and Thermal Decomposition are some the chemical reactions which take place in the atmosphere to form secondary pollutants and secondary aerosols. 1.2 Secondary pollutant The governing equation for the secondary pollutant is (6) where, is the mass ratio of the secondary particulate species to the primary gaseous species which is being converted. In deriving equation (6) we have made similar assumptions as in the case of primary pollutant. The appropriate initial and boundary conditions on C s are: (7) (8) Since there is no direct source for secondary pollutants, we have (9) (10) where is the dry deposition velocity of the secondary pollutant . 3. METEOROLOGICAL PARAMETERS The treatment of equations (1) and (6) mainly depends on the proper estimation of the diffusivity coefficient and the velocity profile of the wind near the ground/or the lower layers of the atmosphere. The meteorological parameters influencing eddy diffusivity and velocity profiles are dependent on the intensity of turbulence, which is influenced by atmospheric stability. Stability near the ground is dependent primarily upon the net heat flux. In terms of the boundary layer notation, the atmospheric stability is characterized by the parameter L [27], which is also a function of net heat flux among several other meteorological parameters. It is defined by (11) where u* is the friction velocity, the net heat flux, the ambient air density, the specific heat at constant pressure, the ambient temperature near the surface, the gravitational acceleration and the Karman’s constant 0.4. and consequently represents stable atmosphere, and represent unstable atmosphere and and represent neutral condition of the atmosphere. The friction velocity is defined in terms of geostrophic drag coefficient and geostrophic wind such that Volume 2, Issue 2, February 2013 Page 254 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 (12) where is a function of the surface Rossby number , where is the Coriolis parameter due to earth’s rotation and is the surface roughness length. [28]Lettau (1959) gave the value of the drag coefficient for a neutral atmosphere in the form (13) The effect of thermal stratification on the drag coefficient [28] can be accounted through the relations: for unstable flow, (14) for slightly stable flow and (15) for stable flow. (16) In order to evaluate the drag coefficient, the surface roughness length z0 may be computed according to the relationship developed by [29] i.e., , where is the effective height of roughness elements, a is the frontal area seen by the wind and is the lot area (i.e., the total area of the region divided by the number of elements). Finally, in order to connect the stability length L to the Pasquill stability categories, it is necessary to quantify the net radiation index. The following values of H f (Table 1) are used for urban area [30]. Table 1. Net heat flux Net radiating index : 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 Net heat flux : 0.24 0.18 0.12 0.06 0.0 -0.03 -0.06 3.1 Eddy diffusivity profiles Following gradient transfer hypothesis and dimensional analysis, the eddy viscosity KM is defined as (17) By using [27] similarity theory, the velocity gradient may be written as (18) Substituting equation (18) in equation (17), we have (19) The function depends on, where is Monin-Obukhov stability length parameter. It is assumed that the surface layer terminates at for neutral stability. For stable conditions, surface layer extends to For neutral stability with (within surface layer) and (20) For stable flow with , (21) (22) For stable flow with , (23) It has been shown that by [31]. In the PBL (planetary boundary layer), where is greater than the limits considered above and , we have, the following expressions for For neutral stability with , Volume 2, Issue 2, February 2013 Page 255 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 (24) For stable flow with , up to , the mixing height, . (25) Equations (19) to (25) give the eddy viscosity for the conditions needed for the model. However, the model deals with the transport of mass rather than the transport of momentum, as implied by the use of viscosity. Since both the mass and the momentum are transported by turbulent eddies, it is physically reasonable to assume that the turbulent viscosity coefficient is numerically equivalent to the eddy diffusivity coefficient . Also, there is some experimental evidence that the ratio remains constant and equal to unity, at least in the surface layer as shown by [31]. The common characteristic of is that it has a linear variation near the ground, a constant value at mid mixing depth and a decreasing trend as the top of the mixing layer is approached. An expression based on theoretical analysis of neutral boundary layer [32], in the form (26) where H is the mixing height. For stable condition, [33] used the following form of eddy-diffusivity, , (27) The above form of was derived from a higher order turbulence closure model which was tested with stable boundary layer data of Kansas and Minnesota experiments. Eddy-diffusivity profiles given by equations (26) and (27) have been used in this model developed for neutral and stable atmospheric conditions. 1.3 Wind velocity profiles In order to incorporate more realistic form of velocity profile in our model which depends on Coriolis force, surface friction, geosrtophic wind, stability characterizing parameter L and vertical height z, we integrate equation (13) from to for neutral and stable conditions. So we obtain the following expressions for wind velocity: In case of neutral stability with , we get (28) In case of stable flow with , we get (29) In case of stable flow with , we get (30) In the planetary boundary layer, above the surface layer, power law scheme has been employed. (31) where, is the geostrophic wind, the wind at , the top of the surface layer, the mixing height and is an exponent which depends upon the atmospheric stability. The values for the exponent , obtained from the measurements made from urban wind profiles [34], as follows: Volume 2, Issue 2, February 2013 Page 256 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 Wind velocity profiles given by equations (28) - (31) due to [30] are used in this model. 4. NUMERICAL METHOD It is to be noted that it is difficult to obtain the analytical solution for the coupled equations (1) and (6) because of the complicated form of wind speed and eddy diffusivity profiles considered in this model. Hence, we have used the numerical method based on the Crank-Nicolson finite difference scheme to obtain the solution. The dependent variable is a function of the independent variables and i.e., . First, the continuum region of interest is overlaid with or subdivided into a set of equal rectangles of sides and , by equally spaced grid lines, parallel to axis, defined by , and equally spaced grid lines parallel to axis, defined by respectively. Time is indexed such that where is the time step. At the intersection of grid lines, i.e. grid points, the finite difference solution of the variable is defined. The dependent variable is denoted by , where and indicate the value at a node point and value at time level respectively. We employ the implicit Crank-Nicolson scheme to discretize the equation (1). The derivatives are replaced by the arithmetic average of their finite difference approximations at the and time steps. Then equation (1) at the grid points and time step can be written as (32) This analog is actually the same as the first order correct analog used for the forward difference equation, but is now second order-correct, since it is used to approximate the derivative at the point . We use the backward differences for advective term in the primary pollutant equation. i.e (33) (34) (35) Also, for the diffusion term, we use the second order central difference scheme Hence, Volume 2, Issue 2, February 2013 Page 257 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 (36) Similarly, (37) Substituting equations (34) to (37) in equation (32) and rearranging the terms we get the finite difference equations for the primary pollutant in the form (38) for each , for each where, and are the values at and respectively and x is the value of at . Equation (38) is true for interior grid points. At the boundary grid points we have to use the boundary conditions (2) to (5). The initial condition (2) is The condition (3) becomes The boundary condition (4) can be written as (39) for and (40) for The boundary condition (5) can be written as (41) The above system of equations (38) to (41) has a tridiagonal structure and is solved by [35]. Volume 2, Issue 2, February 2013 Page 258 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 Similarly the finite difference equations for the secondary pollutant C s can be obtained as (42) for , where, The initial and boundary conditions on secondary pollutant Cs are for for , (43) - (44) Vg is the mass ratio of the secondary particulate species to the primary gaseous species which is being converted. The system of equations (42) to (44) also has tridiagonal structure but is coupled with equations (38) to (41). First, the system of equations (38) to (41) is solved for , which is independent of the system (42) to (44) at every time step n. This result at every time step is used in equations (38) to (41). Then the system of equations (42) to (44) is solved for at the same time step n. Both the systems of equations are solved using Thomas algorithm for tri-diagonal equations (38) to (41) and (42) to (44). Thus, the concentration for primary and secondary pollutants is obtained. An important question concerning computational solutions is what guarantee can be given that the computational solution will be close to the exact solution of partial differential equation and under what circumstances the computational solution will coincide with the exact solution. The second part of this question can be answered by requiring that the approximate or computational solution should converge to the exact solution as the grid spacing x , z and t shrink to zero. However, convergence is very difficult to establish directly and therefore an indirect route as indicated below can be adopted. Figure 2 Solution of partial differential equations The indirect route requires that the system of algebraic equations formed by the discretization process should be consistent with the governing partial differential equation. Consistency implies that the discretization process can be reversed through Taylor series expansion to recover the governing partial differential equations. In addition, the algorithm is said to be convergent, if the approximate solution approaches the exact solution of the partial differential equations, for each value of the independent variable as the grid spacing tends to zero. Convergence can be established by Lax Equivalence theorem “Given a properly posed linear initial value problem and a finite difference approximation to it that satisfies the consistency condition, stability is the necessary and sufficient condition for the convergence” Volume 2, Issue 2, February 2013 Page 259 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 Consistency is necessary if the approximate solution is to converge to the exact solution of the partial differential equations under consideration. Stability is the tendency for any spontaneous perturbation (such as round-off error) in the solution of the system of algebraic equations to decay. The concept of stability is concerned with the growth or decay of errors introduced at any stage of the computation. In this context, errors referred to are not those produced by incorrect logic but by those which occur because the computer cannot give answers to an infinite number of decimal places. In practice each calculation made on the computer is carried out to a finite number of significant figures which introduces a round-off error at every step of the computation. A particular method is stable if the cumulative effect of the entire round off errors produced in the application of the algorithm is negligible. Therefore we will now carry out consistency and stability analysis of the finite difference scheme to solve the partial differential equation. 4.1 Consistency analysis of the governing equation We assume that the pollutants do not undergo removal mechanisms, chemically inert, constant eddy diffusivity and wind velocity. Under these assumptions the governing partial differential equation (1) reduces to (45) We have applied a forward difference scheme for time derivative and backward difference for advective term. For the diffusion term we have used second order central difference scheme for the above equation. Thus we get (46) We have (47) (48) (49) (50) (51) Substituting equations (47) to (51) in equation (46) and as we have We get the original partial differential equation. Hence the implicit finite-difference Crank-Nicolson scheme used for this model to obtain the solution is consistent. 4.2 Stability analysis of the governing equation Now we shall discuss the stability analysis of the governing equation Volume 2, Issue 2, February 2013 Page 260 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 We have used the implicit Crank-Nicolson scheme to discretize the equation (45). The derivatives are replaced by the arithmetic average of their finite difference approximations at the nth and (n+1)th time steps. Then equation (45) at the grid points (i, j) and time step n+1/2 can be written as We have used forward difference for time derivative, backward difference for advection term and second order central difference scheme for the diffusion term. (52) Each Fourier component of the solution is written as where, is the amplitude function at time-level n of the particular component whose wave numbers are (wave length ), (wave length ) and . Define the phase angles, and , . Then Similarly, Substituting in equation (52) we have I i j Canceling the common term e , we get On simplification we obtain (53) We evaluate the amplification factor G from the equation where i.e (54) Volume 2, Issue 2, February 2013 Page 261 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 Note that, i.e., the amplitude factor varies for each Fourier components. If solutions are to remain bounded, we must have for all and . (55) From equations (54) and (55), we have, which is satisfied for all possible and i.e for all possible Fourier mode. Thus the scheme is unconditionally stable. We have analysed the above numerical scheme for consistency and stability. Consistency is investigated for the implicit discretization of the governing diffusion equation. The derivatives are replaced by the arithmetic average of its finite difference approximations at the nth and (n+1)th time steps. The resulting equation coincides with governing diffusion equation as tends to zero. Hence the implicit finite difference scheme used for the solution of this model is consistent. We have used the Von Neumann’s method to study the stability analysis. Using Fourier mode analysis it is shown that the employed numerical scheme is unconditionally stable. Therefore the whole scheme is unconditionally stable for equations (1) and (6). 4.3 Residual The difference between the exact solution and the approximate solution is called the residual. We discuss the residual when the concentration of the pollutant reaches the steady state. When the system has reached the steady state, the time derivative of the physical quantity tends to zero. When the numerical solution obtained is not exactly steady, the time discrete derivatives will not be precisely zero. The non zero value is called residual. The magnitude of the residual indicates the accuracy of the method. When Computational Fluid Dynamics experts are comparing the relative merits of two or more different algorithms for a time marching solution to the steady state, the magnitude of the residuals and their rate of decay are often used as figures of merit. The algorithm which gives the fastest decay of the residuals to the smallest value is usually looked upon most favorably. In this problem we obtain the steady state residual to indicate the accuracy of the Crank - Nicolson method. The concept of residual can be understood from the following procedures. Consider the governing equation, (56) When an upwind version of the Crank - Nicolson method is applied to this equation, we get (57) When the steady state is reached the time derivative of the physical quantity should approach zero if the solution is exactly steady. Since the numerical values of the derivative are not precisely zero, the non zero value of the time derivative is called the residual. This is the left hand side of the equation (56) which is computed from (58) As time progress and as the steady state is approached , the time derivative (58) should approach zero. Since the numerical value of the left hand side of equation (57) are not precisely zero they are called residuals [36]. We have computed the residuals obtained after every time step against the number of time steps and analyzed in figure 3. It is seen that the residuals settle to around 10-6. 1 0 .1 0 .0 1 Residual 1 E -3 1 E -4 1 E -5 1 E -6 1 E -7 0 2 0 0 0 4 0 0 0 6 0 0 0 8 0 0 0 1 0 0 0 0 T im e Figure 3. Plot of residual versus time. 5. RESULTS AND DISCUSSION A numerical model for the computation of the ambient air concentration of the pollutant along the down-wind and the vertical directions emitted from an area source along with the point source on the boundary with the removal Volume 2, Issue 2, February 2013 Page 262 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 mechanism and the transformation process has been presented. The numerical model permits the estimation of the concentration distribution for more realistic meteorological conditions. An area source is an emission source which is spread out over the surface of the city with finite downwind and infinite crosswind dimensions where major source being vehicular exhausts due to traffic flow. We have taken a point source arbitrarily at the beginning of the city. The grid lines are not passing through the point source and it is difficult to adopt point source in numerical method. Therefore we have distributed the concentration of point source equally to its two neighboring grid points. We have considered grid size as along -direction and along -direction. Figures 4 and 5 demonstrates the concentration of the primary pollutant with respect to height for different values of distance with varying values of chemical reaction rate coefficient for the stable and the neutral atmospheric conditions. Since we have considered the point source at and at , the concentration of the pollutant is maximum around 20.5 m height in both the atmospheric conditions. As the distance increases, the concentration of the pollutant decreases due to removal mechanisms or transformation processes. In the stable atmospheric condition, near the ground level the initial concentration of the pollutant is around 30 and then it slowly decreases. It reaches zero around 15 m height and then again increases up to 20.5 m height and then rapidly decreases in the stable case. From Figures 4(a) and 4(b), it is understood that as the value of chemical reaction rate coefficient increases there is no much change in the concentration of primary pollutant near the point source for stable atmospheric condition This behavior is because the increase in the value of the chemical reaction rate coefficient is too smaller when compared to the continuous emission of the point source strength. The similar effect is observed in the neutral atmospheric condition. In the neutral case, the concentration of the primary pollutant is less comparing to the stable case for the same values of distance. The concentration of the pollutant reaches zero in stable case around height and 35 0 P r im a r y p o llu ta n t P r im a r y p o llu t a n t 3 5 0 30 0 3 0 0 k = 0 .0 0 0 0 8 k = 0 .0 1 x= 7 5 x = 7 5 V d= 0 .0 0 5 V d= 0 .0 0 5 25 0 2 5 0 Concentration Concentration x = 1 5 0 20 0 2 0 0 x = 1 5 0 x = 2 2 5 15 0 1 5 0 x = 2 2 5 10 0 1 0 0 5 0 5 0 0 0 (a ) (b ) 0 5 1 0 1 5 2 0 2 5 3 0 3 5 4 0 0 5 10 1 5 2 0 2 5 3 0 35 4 0 H e ig h t H e ig h t Figure 4. Concentration of primary pollutant with respect to height for different distances with varying values of chemical reaction rate coefficient for stable case. 1 6 0 1 40 P r im a r y p o llu t a n t P r i m a r y p o l lu t a n t 1 4 0 k = 0 .0 0 0 0 8 1 20 V d= 0 .0 0 5 k = 0 .0 1 1 2 0 V d= 0 .0 0 5 x = 7 5 1 00 x= 7 5 1 0 0 Concentration 80 Concentration 8 0 x = 1 5 0 60 6 0 x = 15 0 x = 2 2 5 40 4 0 x= 2 2 5 20 2 0 0 0 (a ) (b ) -2 0 0 2 0 4 0 6 0 8 0 1 0 0 1 2 0 0 20 4 0 60 80 100 1 20 H e ig h t H e ig h t Figure 5. Concentration of primary pollutant with respect to height for different distances with varying values of chemical reaction rate coefficient for neutral case. 350 P r im a r y p o l l u t a n t 350 P r im a r y p o llu t a n t 300 V d= 0 .0 1 300 V d= 0 .2 5 k = 0 .0 0 0 0 8 x=75 k = 0 .0 0 0 0 8 x=75 250 250 Concentration Concentration 200 x=150 200 x=1 50 x=225 150 150 x=225 100 100 50 50 0 0 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 H e ig h t H e ig h t Figure 6. Concentration of primary pollutant with respect to height for different distances with varying values of dry deposition velocity for stable case. around height in the neutral case. This behavior is because the neutral atmospheric condition enhances the vertical diffusion of the pollutant to the greater heights. In Figures 6 and 7, the concentration of the primary pollutant with respect to height for different distances with varying values of the dry deposition is analysed for the stable and the neutral atmospheric conditions. As the point source is considered at the beginning of the city region and at 20.5 m height, the concentration is maximum at that height and then the concentration decreases as the distance increases. It reaches zero at 35 m height in the stable case. In the stable Volume 2, Issue 2, February 2013 Page 263 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 case, the ground level concentration of the primary pollutant at is around 35 for . But when , the ground level concentration of the primary pollutant becomes zero at . This clearly indicates that as the value of the dry deposition velocity increases, the ground level concentration of the pollutant decreases. Similar effect is observed in the neutral case. But the concentration of the primary pollutant becomes zero at 100 m height in the neutral atmospheric condition. The concentration of the primary pollutant is less in the neutral case comparing to the stable case for the same values of distance. Also as the value of the dry deposition velocity increases there is no much variation in the concentration of the primary pollutant near the point source for the stable and the neutral atmospheric conditions. This behavior is because the increase in the value of the dry deposition velocity is negligible when compared to the continuous emission of the point source strength and we have considered the dry deposition velocity at the ground level. 160 1 6 0 P r i m a r y p o ll u t a n t P r im a r y p o llu ta n t 140 1 4 0 V = 0 .0 1 d V d = 0 .2 5 k = 0 .0 0 0 0 8 1 2 0 120 k = 0 .0 0 0 0 8 x= 75 x = 7 5 1 0 0 100 Concentration Concentration 8 0 80 x = 1 5 0 x=150 6 0 60 x = 2 2 5 x= 225 4 0 40 2 0 20 0 0 0 2 0 4 0 6 0 8 0 1 0 0 1 2 0 0 20 40 60 80 100 120 H e ig h t H e ig h t Figure 7. Concentration of primary pollutant with respect to height for different distances with varying values of dry deposition velocity for neutral case. 0 .0 0 2 8 S e c o n d a r y p o llu ta n t 0 .0 0 2 8 S e c o n d a ry p o llu ta n t 0 .0 0 2 4 V d = 0 .0 1 0 .0 0 2 4 V d = 0 .2 5 k = 0 .0 0 0 0 8 x= 75 x=75 0 .0 0 2 0 k 0 .0 0 0 0 8 0 .0 0 2 0 x= 150 Concentration Concentration 0 .0 0 1 6 x=150 0 .0 0 1 6 0 .0 0 1 2 0 .0 0 1 2 x=225 x=225 0 .0 0 0 8 0 .0 0 0 8 0 .0 0 0 4 0 .0 0 0 4 0 .0 0 0 0 0 .0 0 0 0 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 H e ig h t H e ig h t Figure 8. Concentration of secondary pollutant with respect to height for different distances with varying values of dry deposition velocity for stable case. 0 .0 0 2 4 S e c o n d a r y p o ll u ta n t S e c o n d a r y p o l lu t a n t V d = 0 .0 1 k = 0 .0 0 0 0 8 0 .0 0 1 6 0 .0 0 2 0 V d= 0 .2 5 x=75 k = 0 .0 0 0 0 8 x=75 0 .0 0 1 6 0 .0 0 1 2 Concentration x=150 Concentration 0 .0 0 1 2 x=150 0 .0 0 0 8 x=225 0 .0 0 0 8 x=225 0 .0 0 0 4 0 .0 0 0 4 0 .0 0 0 0 0 .0 0 0 0 0 20 40 60 80 100 120 0 20 40 60 80 100 120 H e ig h t H e ig h t Figure 9. Concentration of secondary pollutant with respect to height for different distances with varying values of dry deposition velocity for neutral case. Figures 8 and 9 represent the concentration of the secondary pollutants versus height for varied distance with different values of the dry deposition velocity for the stable and the neutral cases. In the stable case, when , initially as the height increases the concentration of the secondary pollutant decreases. Again the concentration increases up to the height 20.5 m and then it decreases and reaches zero at 38 m height. When the value of the dry deposition velocity increases to 0.25, the ground level concentration of the secondary pollutant becomes zero at the beginning. The concentration increases in the downwind distance up to 20.5 m height and then it decreases. It is observed that the concentration of the secondary pollutant decreases as the distance increases in both the stable and the neutral atmospheric conditions. In the neutral case, the ground level concentration of the secondary pollutant is low as compared to that in the stable case when . When , the concentration increases up to 20.5 m height Volume 2, Issue 2, February 2013 Page 264 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 then it decreases and reaches zero around 110 m height. Therefore comparing to the stable case, the concentration of the secondary pollutant reaches zero at the greater heights in the neutral case. Also as the value of dry deposition velocity increases there is no much effect of the secondary pollutant near the point source for the stable and the neutral atmospheric conditions. This behavior is because the dry deposition is introduced at the ground level and increase in the value of the dry deposition velocity is too smaller when compared to the continuous emission of the point source strength. From Figures 10 and 11, the concentration of the primary and the secondary pollutants with respect to distance for different heights for the stable and the neutral atmospheric conditions is studied. In the stable case the concentration of the primary pollutant increases up to and it decreases outside the city because there is no area source beyond the city region. As the height increases, the concentration of the primary pollutant decreases. The concentration of the primary pollutant is maximum near the point source and at the end of the city region. The concentration of the secondary pollutant is zero at the beginning of the city region and it increases with respect to distance. The concentration of the secondary pollutant is more in the outskirt of the urban city at the ground level when compared to inside the city region. In the neutral case the similar effect is observed as in the case of the stable atmosphere for both the primary and the secondary pollutants. But the concentration of the pollutants is less in the neutral case when compared to the stable case for the same values of height. P r im a r y p o llta n t 0 .1 4 S e c o n d a r y p o llu t a n t 30 0 0 .1 2 z= 2 25 0 0 .1 0 z= 2 Concentration 20 0 z= 1 2 Concentration 0 .0 8 z= 1 2 15 0 0 .0 6 z= 2 2 10 0 0 .0 4 z= 3 0 0 .0 2 5 0 z= 22 0 .0 0 z= 3 0 0 - 0 .0 2 0 20 0 0 4 0 00 60 0 0 80 0 0 10 0 00 1 20 00 0 2 0 0 0 4 0 0 0 60 0 0 8 0 0 0 1 0 00 0 1 2 0 0 0 D is ta n c e D is ta n c e Figure 10. Concentration of primary and secondary pollutants with respect to distance for different heights for stable case. 0 .0 2 0 110 P r i m a r y p o l lu t a n t S e c o n d a r y p o llu ta n t 100 0 .0 1 8 z=22 0 .0 1 6 90 z=2 0 .0 1 4 z=2 80 z=12 Concentration 0 .0 1 2 Concentration 70 z=12 z=22 60 0 .0 1 0 0 .0 0 8 z=30 50 z=30 40 0 .0 0 6 30 0 .0 0 4 20 0 .0 0 2 10 0 .0 0 0 0 2000 4000 6000 8000 10000 12000 0 2000 4000 6000 8000 10000 12000 D is t a n c e D is t a n c e Figure 11. Concentration of primary and secondary pollutants with respect to distance for different heights for neutral case. The concentration contours of the primary and the secondary pollutants are drawn in Figures 12 and 13 for both the stable and the neutral atmospheric conditions. It is noted that the concentration of the primary pollutant attains peak value at the downwind end of the source region. Whereas, the concentration of the secondary pollutant attains its peak value at the outside of the city region in the downwind direction. The secondary pollutants are more concentrated away from the ground and are spread out evenly throughout the region. This is due to the chemical reaction taking place in the atmosphere converting the primary into secondary pollutants. The magnitude of the concentration of the primary and the secondary pollutants is higher in the stable case and lower in the neutral case. This behavior is because the neutral atmosphere carries the pollutants to the greater heights and thus the concentration is less. 2 0 0 6 0 P r i m a r y p o l lu t a n t P r im a r y p o ll u t a n t (b ) 1 7 5 6 .5 7 5 5 0 (a ) 1 5 .6 8 1 5 0 7 .8 4 1 4 0 1 2 5 2 3 .5 2 Height 1 0 0 Height 3 0 3 9 .2 0 7 5 1 3 .6 5 2 2 01 9 6 0 1 1 .. 7 4 7 .0 5 7 0 .6 27 . 7 3 5 5 4 .8 9 5 0 1 7 .1 9 6 2 .7 3 3 1 .3 6 2 0 .7 3 1 0 2 5 1 1 6 .2 2 4 .2 6 1 0 .1 1 3 8 . 4 13 1 . 3 4 2 7 .8 0 2 0 0 0 4 0 0 0 6 0 0 0 8 0 0 0 1 0 0 0 0 2 0 0 0 4 0 0 0 6 0 0 0 8 0 0 0 1 0 0 0 0 D is t a n c e D is ta n c e Figure 12. Concentration contours of primary pollutant for (a) stable and (b) neutral atmospheric conditions Volume 2, Issue 2, February 2013 Page 265 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 5 0 S e c o n d a r y p o llu ta n t 3 0 0 S e c o n d a r y p o llu ta n t 2 7 5 (a ) 4 0 1 .4 E -0 2 2 5 0 (b ) 2 .4 E - 0 3 1 .2 E - 0 2 2 2 5 3 2 .1 E - 0 2 .9 E - 0 3 2 0 0 3 0 1 .4 E - 0 3 5 .0 E -0 3 2 .4 E -0 2 2 .6 E -0 2 1 7 5 Height Height 1 5 0 3 .6 E - 0 3 5 .7 E - 0 3 2 .2 E -0 2 1 2 5 2 0 7 .2 E -0 3 1 .9 E -0 2 1 0 0 9 .6 E -0 3 3 6 .4 E -07 .2 E -0 3 1 .7 E -0 2 7 5 1 0 5 0 2 5 7 .9 E - 0 3 4 .3 E -0 3 2 0 0 0 4 0 0 0 6 0 0 0 8 0 0 0 1 0 0 0 0 2 0 0 0 4 0 0 0 6 0 00 80 0 0 1 0 0 0 0 D is t a n c e D is ta n c e Figure 13. Concentration contours of secondary pollutant for (a) stable and (b) neutral atmospheric conditions. 6. CONCLUSION It is well documented that the raised ambient concentrations of the air pollution can have an effect on susceptible people; it is very difficult to predict the likely size of such effects from a single point source and area source. The reality is that the current epidemiological methods cannot conclusively relate effects with the speciﬁc emissions and many studies suffer from the poor estimates of the exposure or small numbers. Any effect at a local level is likely to be very small. Therefore a time dependent two dimensional mathematical model of air pollution due to area source along with a point source on the boundary is presented to simulate the dispersion processes of the gaseous air pollutants in an urban area. The numerical model computes the ambient air concentration emitted from the urban area source and point source, undergoing removal mechanism and transformation process. This model takes into account the realistic form of variable wind velocity and eddy diffusivity profiles, which are the functions of vertical height, frictional velocity, terrain categories, geostrophic wind and several other stability dependent parameters. Reliability of the implicit finite difference scheme is ensured by studying the consistency, stability and convergence through residue fall. The results of this model have been analysed for the dispersion of the air pollutants in the urban area downwind and vertical direction for the stable and the neutral conditions of the atmosphere. The ground level concentration of primary pollutants attains peak value at the downwind end of the source region. Whereas, the concentration of the secondary pollutant attains its peak value at the source free region in the downwind direction. The model also predicts that the ground level concentration of the secondary pollutant at a particular downwind distance is always higher in the stable atmosphere condition than that in the neutral atmospheric condition. There is negligible effect of chemical reaction rate coefficient and dry deposition velocity near the point source. These removal mechanisms plays an important role in reducing the concentration of the pollutants everywhere in the city region except near the point source. Also in the case of the stable atmospheric condition the maximum concentration of the pollutants is observed at the ground level. Same phenomenon is noticed in the neutral condition but the magnitude of the concentration in the neutral condition is comparatively less than that in the stable atmospheric condition. The present study shows that the stable condition is an unfavorable condition for the animals and plants from air pollution point of view. REFERENCES [1] Folinsbee, L.J., 1993. Human health effects of air pollution. Environmental Health Perspectives 100, 45-56 [2] Om P. Tripathi, Stephen G. Jennings, Colin O’Dowd1, Barbara O’Leary, Keith Lambkin, Eoin Moran, Simon J. O'Doherty, T. Gerard Spain, 2012. An assessment of the surface ozone trend in Ireland relevant to air pollution and environmental protection. Atmospheric pollution research 3, 341-351. [3] Anikender Kumar, Pramila Goyal, 2011. Forecasting of air quality in Delhi using principal component regression technique. Atmospheric pollution research 4, 436-444. [4] Seinfeld, J.H., Pandis, S.N., 1997. Atmospheric chemistry and physics from air pollution to climate change, 1st Edition; November, John Wiley & Sons. [5] Friedlander, S.K., 2000. Smoke, dust, and haze: fundamentals of aerosol dynamics, 2nd Edition, Oxford University Press. [6] Almeida, S.M., Pio, C.A., Freitas, M.C., Reis, M.A., Trancoso, M.A., 2005. Source apportionment of fine and coarse particulate matter in a sub urban area at the Western European Coast. Atmospheric Environment 39, 3127 3138. [7] Antonio Febo, Fabio Guglielmi, Maurizio Manigrasso, Valerio Ciambottini, Pasquale Avino, 2010. Local air pollution and long-range mass transport of atmospheric particulate matter: A comparative study of the temporal evolution of the aerosol fractions. Atmospheric pollution research 3, 141-146. [8] Manju Agarwal, Abhinav Tandon, 2009. Interconversion of primary air pollutants to secondary air pollutants, along with their removals through delayed mechanisms. International Journal of Ecological Economics and Statistics 15, 68-92. Volume 2, Issue 2, February 2013 Page 266 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 [9] Haagensen, P.L, Morris, A.L.,1974. Forecasting the behavior of the St. Louis, Missouri pollutant plume. Journal of Applied Meteorology 13, 901-909 [10] Stampfer, J.F. Jr., Anderson, J.A., 1975. Locating the St. Louis urban plume at 80 and 120 km and some of its characteristics. Atmospheric Environment 6, 743-757. [11] Breeding, R.J., Haagensen, P.L., Anderson, J.A., Lodge, J.P. Jr., Stampfer, J.F. Jr.,1975. The urban plume as seen at 80 and 120 km by five different sensors. Journal of Applied Meteorology 14, 204-216. [12] Breeding, R.J., Klonis, H.B., Lodge, J.P. Jr., Pate, J.B., Sheesley, D.C., Englert, T.R., Sears, D.R., 1976. Measurement of atmospheric pollutants in the St. Louis area. Atmospheric Environment 10, 181-194. [13] Lakshminarayanachari,.K., Pandurangappa,.C., Venkatachalappa, M., 2011. Mathematical model of air pollutant emitted from a time defendant area source of primary and secondary pollutants with chemical reaction, International Journal of Computer Applications in Engineering, Technology and Sciences 4, 136-142. [14] Pandurangappa, C., Lakshminarayanachari,.K., Venkatachalappa, M., 2011. Effect of mesoscale wind on the pollutant emitted from a time dependent area source of primary and secondary pollutants with chemical reaction, International Journal of Computer Applications in Engineering, Technology and Sciences 4, 143-150. [15] Hinrichsen, K., 1986. Comparision of four analytical dispersion models for near-surface releases above a grass surface. Atmospheric Environment 20, 29-40. [16] Sheih, C..M., 1977. Application of a statistical trajectory model to the simulation of sulfur pollution over northeastern states. Atmospheric Environment 11, 173-186. [17] Eagan, B.A., Mahoney, J.R., 1972. Numerical modeling of advection and diffusion of urban area source pollutants. Journal of Applied Meteorology 11, 312-322. [18] Ragland, K.W., Dennis, R.L., 1975. Point source atmospheric diffusion model with variable wind and diffusivity profiles. Atmospheric Environment 9, 175-183. [19] Robson, R.E., 1987. Turbulent dispersion in a stable layer with a quadratic exchange coefficient. Boundary Layer Meteorology 39, 207-218. [20] Pal, D., Sinha, D.K., 1987. An area source numerical model on atmospheric dispersion with chemical reaction and deposition. International Journal of Environmental Studies 29, 269-280. [21] Arora, U., Gakhar, S., Gupta, R.S., 1991. Removal model suitable for air polluatants emitted from an elevated source. Applied Mathematical Modeling 15, 386-389. [22] Rudraiah, N., Venkatachalappa, M., Sujit Kumar Khan, 1997. Atmospheric diffusion model of secondary pollutants. International Journal of Environmental Studies 52, 243-267. [23] Khan, S.K., 2000. Time dependent mathematical model of secondary air pollutant with instantaneous and delayed removal. Association for the advancement of modeling and simulation techniques in enterprises 61, 1-14. [24] Manju Agarwal, Abhinav Tandon., 2010. Modeling of the urban heat island in the form of mesoscale wind and its effect on air pollution dispersal. Applied Mathematical modeling 34, 2520-2530. [25] Lakshminarayanachari, K., Sudheer Pai, K.L., Siddalinga Prasad, M., Pandurangappa, C., 2013. A two dimensional numerical model of primary pollutant emitted from an urban area source with mesoscale wind, dry deposition and chemical reaction. Atmospheric Pollution Research, doi:10.5094/APR.2013.011. [26] Pasquill, F., Smith, F.B., 1983. Atmospheric Diffusion. John Wiley & Sons. 89-115. [27] Monin, A.S., Obukhov, A.M., 1954. Basic laws of turbulent mixing in the ground layer of the atmosphere. Doklady Akademii SSSR 151, 163-172. [28] Lettau, H.H., 1959. Wind surface stress and Geostrophic Drag Coefficients in the Atmospheric Surface Layer in Advances in Geophysics Atmospheric Diffusion and Air pollution Academic Press New York 6, 241-254. [29] Lettau, H.H., 1970. Physical and Meteorological Basis for Mathematical Models of Urban diffusion processes Proceedings of Symposium on Multiple Source Diffusion Models USEPA publication AP-86. [30] Ragland, K.W., 1973. Multiple box model for dispersion of air pollutants from area sources. Atmospheric Environment 7, 1071-1089. [31] Webb, E.K., 1970. Profile relationships: the long–linear range, and extension into strong stability Quarterly Journal of Royal Meteorological Society 96, 67-90. [32] Shir, C.C., 1973. A preliminary numerical study of a atmospheric turbulent flows in the idealized planetary boundary layer. Journal of Atmospheric Sciences 30, 1327-1341. [33] Ku, J.Y., Rao, S.T., Rao, K.S., 1987. Numerical simulation of air pollution in urban areas; model development. Atmospheric Environment, 21(1) 201-214. [34] Jones, P.M., Larrinaga, M.A.,,Wilson, C.B., 1971. The urban wind velocity profile. Atmospheric Environment 5, 89-102. Volume 2, Issue 2, February 2013 Page 267 International Journal of Application or Innovation in Engineering & Management (IJAIEM) Web Site: www.ijaiem.org Email: editor@ijaiem.org, editorijaiem@gmail.com Volume 2, Issue 2, February 2013 ISSN 2319 - 4847 [35] Akai, T.J., 1994. Applied Numerical Methods for Engineers. John Wiley and Sons Inc : 86-90. [36] Venkatachalappa, M., Sujit Kumar Khan, Khaleel Ahmed G Kakamari, 2003. Time dependent mathematical model of air pollution due to area source with variable wind velocity and eddy diffusivity and chemical reaction, Proceedings of Indian National Science Academy, 69, A, No.6, 745-758. Volume 2, Issue 2, February 2013 Page 268