Three-dimensional finite element modelling of stack pollutant emissions

Abstract In this paper we propose a finite element method approach for modelling the air quality in a local scale over complex terrain. The area of interest is up to tens of kilometres and it includes pollutant sources. The proposed methodology involves the generation of an adaptive tetrahedral mesh, the computation of an ambient wind field, the inclusion of the plume rise effect in the wind field, and the simulation of transport and reaction of pollutants. The methodology is used to simulate a fictitious pollution episode in La Palma island (Canary Island, Spain). Keywords: air quality modelling, finite element method, adaptive three-dimensional mesh, local scale, Eulerian description, mass-consistent model, wind field simulation. 1 Introduction In this paper, we present a new methodology for local scale air quality simulations, summarized in Algorithm 1, by using a non-steady finite element method with unstructured and adaptive tetrahedral meshes. The aim of this proposal is to introduce an alternative to the standard implementation of current models, improving the computational cost of methods that use structured meshes [15]. Three remarkable uses of unstructured meshes in atmospheric pollution problems are the two-dimensional regional–global examples presented in [19, 1], the threedimensional regional examples, including local refinement with element sizes of 2km, presented in [39], and the three-dimensional tetrahedral meshes for local wind field analysis with element sizes ranging from two meters up to two kilometres, see [25, 27]. The ideas of this last approach are considered to include the plume rise effect in the wind field that is used in the air quality problem. The wind field is crucial for the pollutant transport, specially in complex terrain 1 Paper 35 Three-Dimensional Finite Element Modelling of Stack Pollutant Emissions R. Montenegro1, A. Oliver2, E. Rodríguez1, J.M. Escobar1 G. Montero1 and A. Pérez-Foguet2 1 University Institute for Intelligent Systems and Numerical Applications in Engineering (SIANI), University of Las Palmas de Gran Canaria, Spain 2 Laboratori de Càlcul Numèric (LaCàN) Universitat Politècnica de Catalunya--BarcelonaTech, Spain ©Civil-Comp Press, 2012 Proceedings of the Eighth International Conference on Engineering Computational Technology, B.H.V. Topping, (Editor), Civil-Comp Press, Stirlingshire, Scotland areas. In order to simulate it, we have used a mass-consistent model. Several twodimensional [41] and three-dimensional [26, 28, 10] adaptive finite element solutions have been developed by the authors. The convection, diffusion and reaction problem is usually solved using splitting schemes [22, 6] and specific numerical solvers for time integration of photochemical reaction terms [36, 35, 1]. In this paper a stabilized finite element formulation with a Crank-Nicolson temporal integration is proposed to solve the problem [7, 33]. The chemistry is simulated with the RIVAD chemical model [37]. The transport and chemical terms are treated separately with Strang splitting operators [34], and the non-linear chemical part is solved with a second order Rosenbrock method [40]. Algorithm 1 Wind and air pollution modeling algorithm 1: Construct an adaptive tetrahedral mesh of the 3D domain 2: Wind field simulation from experimental data 3: Wind field modification including the plume rise effect 3.1: Compute plume rise trajectory 3.2: Mesh refinement along plume rise trajectory 3.3: Apply step 2 on the refined mesh 3.4: Modify the vertical wind field along the plume rise 4: Air pollution simulation from stack emission data The paper is organised as follows. In Section 2 we describe the main steps of the proposed methodology. Results are shown in Section 3, and finally the conclusions and future work are presented in Section 4. 2 Algorithm description In this section a brief description of the different steps of Algorithm 1 is presented. 2.1 Adaptive tetrahedral mesh The studied domain is limited at the bottom by the terrain and at the top by a horizontal plane. The lateral walls are formed by four vertical planes. A uniform distribution of nodes is defined on the upper boundary. A refinement/derefinement algorithm [11] is applied on this uniform mesh to construct a node distribution adapted to the terrain surface and stacks. Once the node distribution is defined both on the terrain and the upper boundary, we distribute the nodes located between both layers by using a vertical spacing function. Next, a three-dimensional mesh generator based on Delaunay triangulation [8] is applied. Finally, the untangling and smoothing procedure described in [9] is used to get a valid mesh and to improve its quality. A detailed description of the mesh generation procedure can be seen in [24]. 2 2.2 Wind field simulation A mass-consistent model [26, 28, 10] is used to compute a wind field uin the threedimensional domain Ω, with a boundary Γ = Γa∪Γb, that satisfies the continuity equation and the impermeability condition on the terrain Γa, ∇ · u= 0 in Ω n·u= 0 on Γa (1) where nis the outward-pointing normal unit vector. The model formulates a Least-Squares problem in the domain Ωto find a wind field u= (u, v, w), such that it is adjusted as much as possible to an interpolated wind field u0= (u0, v0, w0). The adjusting functional for a field v= (eu, ev, ew)is defined as E(v) = 1 2ZΩ (v−u0)tP(v−u0)dΩ(2) where Pis a 3×3diagonal matrix with P1,1=P2,2= 2α2 1and P3,3= 2α2 2. The Lagrange multiplier technique is used to minimise the functional (2), with the restrictions (1). Considering the Lagrange multiplier λ, the Lagrangian is defined as L(v, λ) = E(v) + ZΩ λ∇ · vdΩ(3) and the solution uis obtained by finding the saddle point (u, φ)of the Lagrangian (3). This resulting wind field satisfies the Euler-Lagrange equation, u=u0+P−1∇φ(4) where φis the Lagrange multiplier. As α1and α2are constant in Ω, the variational approach results in an elliptic problem in φ, by substituting (4) in (1), that is solved by using the finite element method. −∇ · P−1∇φ=∇ · u0in Ω(5) −n·P−1∇φ=n·u0on Γa(6) φ= 0 on Γb(7) The interpolated wind field u0can be constructed from experimental data or meteorological forecasting models. In this paper we consider the first case. Therefore, we consider an horizontal interpolation and a vertical extrapolation of the available measurements to construct u0in the whole computational domain. Horizontal interpolation The most common technique of interpolation at a given point, placed at a height zmover the terrain, is formulated as a function of the inverse of the squared distance 3 between that point and the measurement stations, and the inverse of their height differences [26] u0(zm) = ξ N P n=1 un d2 n N P n=1 1 d2 n + (1 −ξ) N P n=1 un |∆hn| N P n=1 1 |∆hn| (8) where the value of unis the velocity measured at station n,Nis the number of stations considered in the interpolation, dnis the horizontal distance from station nto the point of the domain where we are computing the wind velocity, |∆hn|is the height difference between station nand the studied point, and ξis a weighting parameter (0≤ξ≤1), that allows to give more importance to one of these interpolation criteria. Vertical extrapolation In this work, a log-linear wind profile is considered [20] in the surface layer, which takes into account the horizontal interpolation [29] and the effect of roughness on the wind intensity and the direction. These values also depend on the air stability (neutral, stable or unstable atmosphere) according to the Pasquill stability class. Above the surface layer, a linear interpolation is carried out using the geostrophic wind. The logarithmic profile is given by, u0(z) = u∗ k(log z z0 −Φm)z0< z ≤zsl (9) where u∗is the friction velocity, kis von Karman’s constant, z0is the roughness length [23] and zsl is the height of the surface layer. The values of Φmdepend on the Pasquill stability class [42], and the friction velocity is obtained from (9) at any point (x, y)by using the horizontal interpolated velocity u0(zm). The linear interpolation is given by, u0(z) = ρ(z)u0(zsl) + [1 −ρ(z)]ugzsl < z ≤zpbl (10) where ugis the geostrophic wind, zpbl is the height of the planetary boundary layer, and ρ(z)is defined as ρ(z) = 1 −z−zsl zpbl −zsl 23−2z−zsl zpbl −zsl (11) Finally, this model assumes u0(z) = ugif z > zpbl and u0(z) = 0 if z≤z0. 2.3 Plume rise The plume rise phenomenon is mainly due to the difference of temperature between the released substance and the environment air, and the initial momentum. The trajectory of the plume rise has been widely studied in the past [4, 5, 30]. These works 4 differentiate between two kinds of cases: predominant buoyancy rise and predominant momentum rise. The characterization of these types essentially depends on the ratio between the intensities of the pollutant emission velocity and the wind velocity at the top of the stack. Gaussian plume models [31] approximate the effective height of a plume zHand the horizontal distance dffrom the stack to the point where the plume height reaches zH, depending on the emission characteristics, the ambient wind and the atmospheric stability. The gas elevation mainly depends on the density difference between the emitted gas and the atmospheric air (buoyancy rise) and the emission velocity (momentum rise). Predominant buoyancy rise In all cases with dfdifferent from zero, the driving force is buoyancy, except for stable conditions and calm wind. In order to know the plume rise trajectory, we propose to combine an horizontal and a vertical motion, satisfying certain known conditions. The vertical motion along the mean trajectory of the plume is defined by an acceleration a0(t), a velocity w0(t)and z(t), from the initial time t= 0 to the final time t=tfwhen the plume reaches the effective height, satisfying the following conditions z(0) = z′ cz(tf) = zH(12) w0(0) = wcw0(tf) = 0 (13) Since there are four conditions on the vertical motion, we propose a cubic approximation of z(t), and therefore a quadratic approximation of w0(t), and a linear approximation of a0(t). The horizontal motion is defined by a uniformly accelerated motion, with a constant positive acceleration vector ad= (adx, ady), a velocity ud(t) = (ud(t), vd(t)), and an horizontal relative position vector d(t) = (x(t)−xc, y(t)−xc)with respect to the centre of the stack, satisfying the following conditions |d(tf)|=df(14) ud(0) = u(xc, yc, zc)(15) In order to define the mean trajectory of bent curved plumes considering the influence of complex terrains, we approximate it by a three-dimensional polygonal line taking into account the ambient wind directions, such that the longitude of its projection on the horizontal plane approximates the longitude df. In addition, the final height coincides with the effective height zH. Therefore, this approximation tries to satisfy the main values of the end of the plume considering Briggs’ equations. Predominant momentum rise In all cases where dfis equal to zero, that is when the driving force is momentum or when the driving force is buoyancy with calm wind, the horizontal motion of the plume until reaching the effective height can be considered negligible. Thus the trajectory of the gases is nearly vertical. 5 In this case, we propose a vertical motion along the trajectory of the plume with a constant negative acceleration a0, a linear velocity w0(t)and a quadratic trajectory z(t). Imposing the conditions (12) and (13), this vertical motion is completely defined. tf=2 (zH−z′ c) wc (16) a0=−wc tf (17) w0(z) = wcs1−2 (z−z′ c) wctf (18) In order to modify the ambient vertical wind velocity (w) along the region of the plume rise, we need to have a sufficient mesh resolution in this area. For this reason, we propose to refine locally the mesh [13] along the Gaussian plume [14] until all the tetrahedra inside that region fulfill a size criterion. Finally, a new ambient wind field uis obtained on the refined mesh with the massconsistent model described in Section 2.2. The effect of the gas emission is introduced in this field by modifying its vertical component along the plume. 2.4 Air pollution simulation The air pollution simulation consists of solving the unsteady convection–diffusion– reaction formulation with an stabilized finite element method, specifically Least-Squares method, with a Crank-Nicolson temporal discretization. The equation system governing the problem can be expressed with the following vectorial equation: ∂c ∂t +u· ∇c=∇ · (K∇c) + e+s(c)(19) for the spatial coordinates xand time t,(x, t)∈Ω×(0, tend], with initial condition c(x,0) = cini(x)on x∈Ω, and the following boundary conditions:          c(x, t) = cemi(x)in ΓS: Top of the stack n·K∇c=−Vdcin ΓR: Terrain n· ∇c= 0 in ΓWout : Outlet wind boundaries c(x, t) = cout(x)in ΓWin : Inlet wind boundaries (20) where ∇is the gradient with respect to x, and c,u,eand s(c)are respectively the concentration, the perturbed wind velocity, the emission and the chemical vectors with a dimension nc(the number of pollutant species). Kis the diffusion matrix of dimension 3×nc,Vdis the deposition diagonal matrix with dimension nc, and n is the outward-pointing normal unit vector, cemi is the concentration of the emission in the top of the stack, and cout the outside concentration at the inlet wind boundaries. 6 Scalar product “·” is applied nctimes: the first argument is multiplied by each one of the nccomponents of the second argument. The complete description of photochemical reaction of atmospheric species is highly complex [12, 18, 2, 32]. For instance, detailed Volatile Organic Components decomposition involves hundreds of thousand reactions [3, 38] that needs special methodologies to reduce the number of the modelled reactions and species. Reference models for gaseous phase reactions involve some tens of compounds [16, 17]. The most simplified models just involve about ten reactive species [43]. On the other side, depending on the application, it can be necessary to take into account aqueous phase reactions, that involve several other reactions and species. The RIVAD model is one of the most simplified models that permit to simulate both processes, aqueous and gaseous, involving transport and reaction of four species [37]. In this paper, we have considered the RIVAD model for the chemical term s(c). This model is a pseudo first-order chemical scheme for acid rain simulation, specially calibrated for being used in non-urban areas. The concentration cinvolves four species, c1= [SO2],c2= [SO4],c3= [NO2]and c4= [NO3], and the components of the reaction vector s(c)are: s1(c) = −s2(c) = −α1(c)c1(21) s3(c) = −s4(c) = −α3(c)c3(22) where α1(c) = γ1/(c1+δ1c3)and α3(c) = γ3/(c3+δ3c1). Note that for values close to zero of the concentration of the primary species c1and c3, both α1(c)and α3(c) requires a proper numerical treatment in order to avoid excessively high reaction rates. Next we will treat the linear and the non-linear problem separately. The development of the linear chemical problem focuses on the temporal and spatial discretization of the corresponding linear convection–diffusion–reaction equation (19). The nonlinear case focuses on the development of an splitting method that combines the solution of (19) in a linear case, considering a null chemical term, and the solution of an ordinary differential equation system that approximates the evolution of the chemical reaction separately. 2.4.1 Linear chemical problem In this case the chemical term is linear, that is, s(c) = Ac where Ais constant matrix. The resulting equation (19) is solved with a Crank-Nicolson time integration scheme, and an spatial discretization with a stabilized finite element method, Least-Squares. Concentrations cnand cn+1 at times tnand tn+1 =tn+ ∆tare related using a CrankNicolson scheme as cn+1 =cn+∆t 2h∂cn+1 ∂t +∂cn ∂t i. We define the differential operator Las L=u· ∇ − ∇ · (K∇)−A(23) and a function Fin Ωfrom the known cn F=cn−∆t 2Lcn+∆t 2en+1 +en(24) 7 Applying the Crank-Nicolson scheme, we can rewrite equation (19) as I+∆t 2Lcn+1 =F(25) where Iis the identity operator. Using Least-Squares we obtain a symmetric problem such that the weak form of the equation (25) is hI+∆t 2Liν, hI+∆t 2Licn+1=hI+∆t 2Liν, F(26) where νis the test function, and (µ, ν) = RΩµν dΩis the inner product We define the next bilinear forms a(ν, c) =(ν, u· ∇c−Ac) + (∇ν, K∇c)+ ZΓR Vdcν dΓR (27) b(ν, c) = X e (u· ∇ν−Aν, c)e(28) g(ν, c) = X e (u· ∇ν−Aν, u· ∇c−Ac)e(29) where (µ, ν)e=RΩeµν dΩ, and Perepresents the sum over all the mesh elements. Applying the operators Land Fin the weak form (26), we obtain (ν, cn+1) + ∆t 2a(ν, cn+1) + b(ν, cn+1) + ∆t 2g(ν, cn+1)= (ν, cn)−∆t 2a(ν, cn)−b(ν, cn) + ∆t 2g(ν, cn)+ ∆t(ν, en+1 2) + ∆t 2b(ν, en+1 2)(30) where en+1 2=en+1+en 2. Equation (30) can be written as an equation system Bc c cn+1 =f(31) where c c cn+1 is the concentration vector approximation at tn+1 in the degrees of freedom of the finite element discretization, f= [2F−B]c c cn+∆t 2F[e e en+1 +e e en], and Band F are square matrices with dimension (nc×ndof ), being ndof the number of degrees of freedom. In order to solve this linear system it is necessary to find an efficient solver, using a sparse matrix storage. Since Bis a symmetric positive definite matrix, we have considered a solver based on a conjugate gradient method preconditioned with an 8 incomplete Cholesky factorisation density type [21]. The left–preconditioning is used to improve the convergence of the conjugate gradient method. The original linear system is transformed into T−1Bc c cn+1 =T−1f, where Tis the symmetric positive definite preconditioner obtained with the incomplete Cholesky factorisation. The large fill-in of the complete (i.e. standard) Cholesky factorisation is completely or partially avoided by discarding coefficients along the factorisation process. We have considered an incomplete Cholesky factorisation with no fill-in, such that the incomplete factor L has the same sparsity pattern as the lower triangle of matrix B. The main advantage of the Cholesky method is that the incomplete factorisation of matrix Bcan be amortized over many time-steps. More details about the implementation of this system equation solver can be found in [7] and [33]. 2.4.2 Non-linear chemical problem To deal with the non-linearity of the reactive term in the convection–diffusion–reaction equation (19), we have considered a splitting method that separates this equation into a convection–diffusion equation and a reaction equation. We will make use of the second order splitting operator (Strang splitting) proposed by [34]: dc⋆ dt =s(c⋆) for t∈0,∆t 2and c⋆(x,0) = cn(x) (32) ∂c⋆⋆ ∂t +u· ∇c⋆⋆ =∇ · (K∇c⋆⋆) + e for t∈[0,∆t]and c⋆⋆(x,0) = c⋆(x,∆t 2) (33) dc⋆⋆⋆ dt =s(c⋆⋆⋆) for t∈∆t 2,∆tand c⋆⋆⋆(x,∆t 2) = c⋆⋆(x,∆t) (34) Once we have split the equation (19), we solve three equations in different time steps; the reaction equation (32), the convection–diffusion equation (33), and the reaction equation (34), being finally cn+1(x) = c⋆⋆⋆(x,∆t). The convection–diffusion equation (33) is solved using the same method proposed in the previous Section, with A= 0. The non-linear chemical equations (32) and (34) are solved node by node with a second order Rosenbrock method (ROS2) [40]. To use the ROS2 method, the Jacobian square matrix of s(c)of dimension nchas to be computed. 9 [13] J.M. Gonz´ alez-Yuste, R. Montenegro, J.M. Escobar, G. Montero, E. Rodr´ ıguez, “Local refinement of 3-D triangulations using object-oriented methods”, Advances in Engineering Software, 35(10-11): 693 – 702, 2004, ISSN 0965-9978. [14] A.E.S. Green, R.P. Singhal, R. Venkateswar, “Analytic extensions of the Gaussian Plume Model”, Journal of the Air Pollution Control Association, 30(7): 773–776, 1980. [15] K. Hanjali, S. Kenjere, “Dynamic simulation of pollutant dispersion over complex urban terrains: A tool for sustainable development, control and management”, Energy, 30(8): 1481 – 1497, 2005, ISSN 0360-5442. [16] P. Jimenez, J.M. Baldasano, D. Dabdub, “Comparison of photochemical mechanisms for air quality modeling”, Atmospheric Environment, 37(30): 4179–4194, 2003, ISSN 1352-2310. [17] F. Kirchner, “The chemical mechanism generation programme CHEMATA – Part 1: The programme and first applications”, Atmospheric Environment, 39 (6): 1143–1159, 2005, ISSN 1352-2310. [18] D. Kley, “Tropospheric Chemistry and Transport”, Science, 276(5315): 1043– 1044, 1997. [19] I. Lagzi, D. K´ arm´ an, T. Tur´ anyi, A.S. Tomlin, L. Haszpra, “Simulation of the dispersion of nuclear contamination using an adaptive Eulerian grid model”, Journal of Environmental Radioactivity, 75(1): 59–82, 2004, ISSN 0265-931X. [20] D. Lalas, C. Ratto, Modelling of Atmospheric Flow Fields, World Scientific Publishing, Singapore, 1996. [21] C.J. Lin, J.J. Mor´ e, “Incomplete Cholesky factorizations with limited memory”, SIAM Journal on Scientific Computing, 21(1): 24–45, 1999, CODEN SJOCE3, ISSN 10648275. [22] M.J. Mart´ ın, D.E. Singh, J.C. Mouri˜ no, F.F. Rivera, R. Doallo, J.D. Bruguera, “High performance air pollution modeling for a power plant environment”, Parallel Computing, 29(11-12): 1763–1790, 2003, ISSN 0167-8191. [23] G.J. McRae, W.R. Goodin, J.H. Seinfeld, “Development of a second-generation mathematical model for Urban air pollution–I. Model formulation”, Atmospheric Environment (1967), 16(4): 679–696, 1982, ISSN 0004-6981. [24] R. Montenegro, G. Montero, J.M. Escobar, E. Rodr´ ıguez, J.M. Gonz´ alez-Yuste, “Tetrahedral Mesh Generation for Environmental Problems over Complex Terrains”, in P. Sloot, A. Hoekstra, C. Tan, J. Dongarra (Editors), Computational Science ICCS 2002, Volume 2329 of Lecture Notes in Computer Science, pages 335–344. Springer Berlin / Heidelberg, 2002, ISBN 978-3-540-43591-4. [25] R. Montenegro, G. Montero, J.M. Escobar, E. Rodr´ ıguez, J.M. Gonz´ alez-Yuste, “3-D Adaptive Wind Field Simulation Including Effects of Chimney Emissions”, in Proceedings of WCCM VI/APCOM’04, Beijing, China. Tsinghua University Press and Springer-Verlag, 2004. [26] G. Montero, R. Montenegro, J.M. Escobar, “A 3-D diagnostic model for wind field adjustment”, Journal of Wind Engineering and Industrial Aerodynamics, 74-76(0): 249–261, 1998, ISSN 0167-6105. [27] G. Montero, R. Montenegro, J.M. Escobar, E. Rodr´ ıguez, J.M. Gonz´ alez-Yuste, 16 “Velocity Field Modelling for Pollutant Plume Using 3-D Adaptive Finite Element Method”, in M. Bubak, G. van Albada, P. Sloot, J. Dongarra (Editors), Computational Science - ICCS 2004, Volume 3037 of Lecture Notes in Computer Science, pages 642–645. Springer Berlin / Heidelberg, 2004, ISBN 978-3-54022115-9. [28] G. Montero, E. Rodr´ ıguez, R. Montenegro, J.M. Escobar, J.M. Gonz´ alez-Yuste, “Genetic algorithms for an improved parameter estimation with local refinement of tetrahedral meshes in a wind model”, Advances in Engineering Software, 36 (1): 3–10, 2005, ISSN 0965-9978. [29] G. Montero, N. San´ ın, “3-D modelling of wind field adjustment using finite differences in a terrain conformal coordinate system”, Journal of Wind Engineering and Industrial Aerodynamics, 89(5): 471–488, 2001, ISSN 0167-6105. [30] D.J. Moore, “A comparison of the trajectories of rising buoyant plumes with theoretical/empirical models”, Atmospheric Environment (1967), 8(5): 441 – 457, 1974, ISSN 0004-6981. [31] L.E. Olcese, B.M. Toselli, “Development of a model for reactive emissions from industrial stacks”, Environmental Modelling & Software, 20(10): 1239–1250, 2005, ISSN 1364-8152. [32] A.R. Ravishankara, “Heterogeneous and Multiphase Chemistry in the Troposphere”, Science, 276(5315): 1058–1065, 1997. [33] A. Rodr´ ıguez-Ferran, M.L. Sandoval, “Numerical performance of incomplete factorizations for 3D transient convection–diffusion problems”, Advances in Engineering Software, 38(6): 439–450, 2007, ISSN 0965-9978. [34] D.L. Ropp, J.N. Shadid, C.C. Ober, “Studies of the accuracy of time integration methods for reaction–diffusion equations”, Journal of Computational Physics, 194(2): 544–574, 2004, ISSN 0021-9991. [35] A. Sandu, J. Verwer, M.V. Loon, G. Carmichael, F. Potra, D. Dabdub, J. Seinfeld, “Benchmarking stiff ODE solvers for atmospheric chemistry problems-I. Implicit vs explicit”, Atmospheric Environment, 31(19): 3151–3166, 1997, ISSN 1352-2310. [36] R.D. Saylor, G.D. Ford, “On the comparison of numerical methods for the integration of kinetic equations in atmospheric chemistry and transport models”, Atmospheric Environment, 29(19): 2585–2593, 1995, ISSN 1352-2310. [37] J.S. Scire, D.G. Strimaitis, R.J. Yamartino, A User’s Guide for the Calpuff Dispersion Model (version 5), Earth Tech., Inc, Concord, MA., 2000. [38] S. Szopa, B. Aumont, S. Madronich, “Assessment of the reduction methods used to develop chemical schemes: Building of a new chemical scheme for VOC oxidation suited to three-dimensional multiscale HOx-NOx-VOC chemistry simulations”, Atmospheric Chemistry and Physics, 5(9): 2519–2538, 2005. [39] A. Tomlin, S. Ghorai, G. Hart, M. Berzins, “3-D Multi-scale air pollution modelling using adaptive unstructured meshes”, Environmental Modelling & Software, 15(6-7): 681–692, 2000, ISSN 1364-8152. [40] J.G. Verwer, E.J. Spee, J.G. Blom, W. Hundsdorfer, “A Second-Order Rosenbrock Method Applied to Photochemical Dispersion Problems”, SIAM Journal 17 on Scientific Computing, 20(4): 1456–1480, 1999, CODEN SJOCE3, ISSN 10648275. [41] G. Winter, G. Montero, L. Ferragut, R. Montenegro, “Adaptive strategies using standard and mixed finite elements for wind field adjustment”, Solar Energy, 54 (1): 49–56, 1995, ISSN 0038-092X. [42] P. Zannetti, Air Pollution Modeling, Computational Mechanics Publications, Boston, 1990. [43] Z. Zlatev, Computer Treatment of Large Air Pollution Models, Kluver Academic Publishers, Dordrecht., 1995. 18