Adaptive finite element simulation of stack pollutant emissions over complex terrains

Adaptive finite element simulation of stack pollutant emissions over complex terrains A. Olivera,∗, G. Monterob, R. Montenegrob, E. Rodr´ıguezb, J.M. Escobarb, A. P´erez-Fogueta aLaboratori de C`alcul Num`eric (LaC`aN), www-lacan.upc.edu Departament de Matem`atica Aplicada III Universitat Polit`ecnica de Catalunya - BarcelonaTech Jordi Girona 1-3, 08034 Barcelona, Spain bUniversity Institute for Intelligent Systems and Numerical Applications in Engineering (SIANI), www.dca.iusiani.ulpgc.es/proyecto2012-2014, University of Las Palmas de Gran Canaria, 35017 Las Palmas, Spain Abstract A three-dimensional finite element model for the pollutant dispersion is presented. In these environmental processes over a complex terrain, a mesh generator capable of adapting itself to the topographic characteristics is essential. The first stage of the model consists on the construction of an adaptive tetrahedral mesh of a rectangular region bounded in its lower part by the terrain and in its upper part by a horizontal plane. Once the mesh is constructed, an adaptive local refinement of tetrahedra is used in order to capture the plume rise. Wind measurements are used to compute an interpolated wind field, that is modified by using a mass-consistent model and perturbing its vertical component to introduce the plume rise effect. Then, we use an Eulerian convection-diffusion-reaction model to simulate the pollutant dispersion. In this work, the transport of pollutants is considered and dry deposition is formulated as a boundary condition. The discretization of the stack geometry allows to define the emissions as boundary conditions. The proposed model uses an adaptive finite element space discretization, a Crank-Nicolson time scheme, and an splitting operator. This approach has been applied in La Palma island. Finally, numerical results and conclusions are presented. Keywords: Power plant emissions, Air quality modelling, Eulerian description, Wind field simulation, Finite element method, Adaptive tetrahedral mesh 1. Introduction Numerical simulation of pollutant transport and reaction on atmosphere has been the result of important advances in the last thirty years. However, nowadays it remains as a scientific challenge. Key analysis are related to acid rain, ozone, particle matter and toxic emissions [1]. Air quality modelling systems mainly involve three components: emissions, meteorology and transport–chemistry. The first component characterizes different emission sources of chemical compounds [2], the second one determines the atmospheric phenomena as wind and temperature fields, and the third one simulates the transport of pollutants (convection and diffusion) and their chemical reactions. Environmental Protection Agency (www.epa.gov) classifies Air Quality Modelling Systems as Dispersion, Photochemical and Receptor models. Dispersion models estimate pollutant concentrations at ground level ∗Corresponding author Email addresses: [email protected] (A. Oliver), [email protected] (G. Montero), [email protected] (R. Montenegro), [email protected] (E. Rodr´ıguez), [email protected] (J.M. Escobar), [email protected] (A. P´erez-Foguet) Preprint submitted to Energy October 19, 2012 Accepted for publication in Energy, Oct. 2012 near a punctual source. Photochemical models consider all the sources in a large area. Receptor models identify and characterize the emission sources using receptor measures. Main dispersion models include Gaussian plume models [3], particle tracking models [4], and puff models [5, 6, 7]. The two latters models are based on a Lagrangian approach and usually include first order chemistry reactions or linearised photochemical-reaction models. Versions with complex non-linear photochemical-reaction models have been also developed [8]. In contrast with dispersion models, photochemical ones follow an Eulerian description of the coupled transport and hydrodynamic problem. Most of these Eulerian models approximate the solution to this problem by using finite difference schemes. The modelling domains usually vary from about few thousands to tens of kilometres, and the grid spatial discretization varies from about tens of kilometres to one kilometre using nested sub-grids [9, 10, 11]. The height of the computational domain usually varys from 4 up to 10 km with a non-uniform vertical discretization between 10 and 20 levels. Vertical spatial resolution ranges from few tens meters close to the ground level to one thousand meters above two kilometres over the terrain. The number of grid points in this kind of problems can vary from tens to hundreds thousand. The most advanced photochemical models consider local emissions and are known as Plume in Grid (PIG). For this purpose, CAMx [12] and CMAQ [13] include a puff model, and UAM-V [14] uses a plume model. In environmental management, dispersion and photochemical models have very different applications. The first ones are usually applied to local emission impact assessments, and the second ones to regional planning and monitoring. Both models have a clearly different application scale and photochemical reaction complexity. However, because of the great coupling among key components as ozone, nitrogen compounds, and Volatile Organic Components [15, 16], and due to the awareness about the socio-economic impacts of their immission [17, 18, 19], some references about the need of coupling local emissions using regional planning with Plume in Grid photochemical models can be found [20, 9, 10]. Although Plume in Grid models can couple local scales (up to a resolution of one kilometre) with regional scales, several limitations has been reported [9]. For example, it is clear that this local high resolution can be insufficient for complex terrains. Thus, the search of alternatives is justified. In this paper, we present a new methodology for local scale air quality simulations 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 [21]. Three remarkable uses of unstructured meshes in atmospheric pollution problems are the two-dimensional regional–global examples presented in Lagzi et al. [22], Ahmad et al. [23], the three-dimensional regional examples, including local refinement with element sizes of 2 km, presented in Tomlin et al. [24], and the three-dimensional tetrahedral meshes for local wind field analysis with element sizes ranging from two meters up to two kilometres, see Montenegro et al. [25], Montero et al. [26]. The ideas of this last approach are considered to determine the wind field, that includes the effect of the plume rise, used in the air quality problem. The wind field is crucial for the pollutant transport, specially in complex terrain areas. In order to simulate it, we have used a mass-consistent model. Several two-dimensional [27] and three-dimensional [28, 29, 30] adaptive finite element solutions have been developed by the authors. The convection, diffusion and reaction problem is usually solved using splitting schemes [31, 32] and specific numerical solvers for time integration of photochemical reaction terms [33, 34, 23]. A non-steady and non-linear transport model is presented in this paper. A stabilized finite element formulation, specifically Least-Squares, with a Crank-Nicolson temporal integration is proposed to solve the problem [35, 36]. The chemistry is simulated by using the RIVAD/ARM3. This is a very simplistic empirical chemistry model used in codes such as CALPUFF [37] and is employed here as a proof of concept. More complete chemistry models will be included in future work. The transport and chemical terms are treated separately with Strang splitting operators [38]. The non-linear chemical part is solved node by node with a second order Rosenbrock method [39]. A previous description of the proposed procedure can be found in P´erez-Foguet et al. [40]. In the second Section of the paper we describe the proposed methodology. It mainly involves the following steps: an automatic tetrahedral mesh generator, the wind field simulation, the plume rise approximation 2 Accepted for publication in Energy, Oct. 2012 and the air pollution simulation. In the third Section we apply the proposed air quality model in La Palma island (Canary Islands, Spain). 2. Wind and air pollution modelling In this section we introduce the models to simulate the wind field and the transport and reaction of pollutants. The evaluation of the wind field is based on a mass-consistent model, while the transport of pollutant is calculated by a convection–diffusion–reaction PDE by using a non-linear chemical model. The proposed methodology is summarized in Algorithm 1 Algorithm 1 Wind and air pollution modelling 1: Construct an adaptive tetrahedral mesh of the domain 1.1: Adaptive discretization of the terrain surface 1.2: Vertical spacing function and 3D distribution of points 1.3: Three-dimensional mesh generation 1.4: Mesh optimisation 2: Wind field simulation from experimental or forecasting data 2.1: Construction of the initial interpolated wind field 2.2: Approximation of the wind field with a mass-consistent model 3: Wind field modification including the plume rise effect 3.1: Compute a plume rise trajectory 3.2: Mesh refinement along the plume rise trajectory 3.3: Simulate a wind field in the new mesh (applying step 2) 3.4: Modify vertical components of the 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.1. Automatic tetrahedral mesh generation 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 [41] 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 [42] is applied. Finally, the untangling and smoothing procedure described in [43] is used to get a valid mesh and to improve its quality. A detailed description of the mesh generation procedure can be seen in [44, 45]. 2.2. Wind field simulation Once the tetrahedral mesh is constructed, we consider a mass-consistent model [28, 29, 30] to compute a wind field uin the three-dimensional domain Ω, with a boundary Γ = Γa∪Γb, that verifies the continuity equation (mass conservation) for constant density 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 wind field u verifies the Euler-Lagrange equation, u=u0+P−1∇φ(2) 3 Accepted for publication in Energy, Oct. 2012 where φis the Lagrange multiplier, Pis a 3 ×3 diagonal matrix with P1,1=P2,2= 2α2 1and P3,3= 2α2 2, being α1and α2constant in Ω. The variational approach results in an elliptic problem in φ, by substituting (2) in (1), that is solved by using the finite element method. −∇ · P−1∇φ=∇ · u0in Ω (3) −n·P−1∇φ=n·u0on Γa(4) φ= 0 on Γb(5) The interpolated wind field u0can be constructed from experimental data or meteorological forecasting models. In this paper we consider the first case. Therefore, we use an horizontal interpolation and a vertical extrapolation of the available measurements to construct u0in the whole computational domain. The horizontal interpolation is formulated as a function of the inverse of the squared distance to the measurement stations, and the inverse of the height differences [28]. For the vertical extrapolation, a log-linear wind profile is considered. It takes into account the horizontal interpolation and the effect of roughness on the wind velocity [46, 47, 48, 49, 50]. For proof of concept this work adopts the simpler Pasquill stability class. Future work can employ a more preciselly planetary boundary layer theory. Note that the log law is recommended only for bi-dimensional fluid flow. In three-dimensional fluid flows there is a significant change on the near wall behavior. The mean velocity and the shear stress vector may shift directions with the distance from the wall; furthermore, the shear stress and the mean vector are not, generally, aligned [51]. However, we use the log law because there is no information about the near three-dimensional wall behavior, and it is a simple method to construct a reasonable interpolated wind field that will be modified by the mass-consistent model taking into account the topography. 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 [52, 53, 54]. These works differentiates two kind 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. The previous studies are based on an important assumption such as a constant ambient wind field with null vertical component. This assumption produce that the mean trajectory of the plume is placed on a vertical plane, that contains the stack. Therefore, the resulting plume rise trajectory cannot be directly applied in the case of complex terrains. In this case, the plume rise trajectory follows a bent curve due to the variation of the ambient wind field. The first objective of this Section is to construct a reasonable plume rise mean trajectory that takes into account the Briggs’ equations and the ambient wind field obtained in the previous Section. In the future work we will implement a more precise numerically integrated approach, as the one used in the PRIME algorithm [55]. Once the plume trajectory is defined, the effect of the pollutant emission will be introduced by modifying the ambient vertical wind velocity (w) along the region of the plume rise trajectory, taking into account the vertical velocity of the plume (w0). As stated by Briggs [53] w0is a convenient variable with which to identify the plume. In this Section we introduce a method to give an approximation of w0at any point in the surroundings of the bent plume trajectory. 2.3.1. Plume rise mean trajectory Using the Briggs formulation, we know the predominant force in the plume rise, and the values of the effective height of the plume (zH) and the horizontal distance (df). In this section a formulation to know the trajectory of the plume rise and the vertical component of the perturbed wind field along this trajectory has to be defined. Firstly, we will consider separately the two different types of trajectories, bent-over or vertical, in the vertical plane defined by the top of the stack and the end of the plume. 4 Accepted for publication in Energy, Oct. 2012 Secondly, we will extend the previous results to three-dimensional trajectories in the case of a complex terrain. Bent-over plume mean trajectory in a vertical plane In all cases with df different 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, verifying 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, verifying the following conditions z(0) = z0 cz(tf) = zH(6) w0(0) = wcw0(tf) = 0 (7) 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) given by the following expressions z(t) = z0 c+wct+−2wctf+ 3(zH−z0 c) t2 f t2+wctf−2(zH−z0 c) t3 f t3(8) w0(t) = wc+−4wctf+ 6(zH−z0 c) t2 f t+3wctf−6(zH−z0 c) t3 f t2(9) a0(t) = −4wctf+ 6(zH−z0 c) t2 f +6wctf−12(zH−z0 c) t3 f t(10) To satisfy the conditions (7), and considering that w0(t) is a decreasing function in [0, tf], we impose a0(t)≤0 along the trajectory, obtaining the following condition for the unknown tf 3 2(zH−z0 c)≤wctf≤3 (zH−z0 c) (11) To determine the value of tf, an additional assumption has to be done introducing the horizontal motion. 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, verifying the following conditions |d(tf)|=df(12) ud(0) = u(xc, yc, zc) (13) Note that at the top of the stack the ambient wind u(xc, yc, zc) has only horizontal components because of the imposition of the impermeability condition (1) at the top of the stack when solving the mass-consistent problem. The horizontal trajectory can be expressed as x(t) = xc+u(xc, yc, zc)t+1 2adxt2 y(t) = yc+v(xc, yc, zc)t+1 2adyt2 (14) From the previous conditions, we obtain an additional restriction for the unknown tf tf=1 ad− |u(xc, yc, zc)|+q|u(xc, yc, zc)|2+ 2addf(15) 5 Accepted for publication in Energy, Oct. 2012 where ad=|ad|is still an unknown. Imposing the condition (11) we obtain the following expression for ad ad= (1 + δ)2wc 3(zH−z0 c)(1 + δ)wc 3(zH−z0 c)df− |u(xc, yc, zc)|(16) being 0 ≤δ≤1. For δ= 0, the value of tfis related to the upper bound in (11) and, for δ= 1, to the lower bound. The case δ= 1/2 corresponds to a value of tfwhich produces a constant vertical acceleration a0, a linear vertical component of velocity w0(t) and a quadratic vertical position z(t). Note that if the parameter δis fixed, the horizontal and vertical motion is completely defined verifying strictly the conditions (12), (13), (6) and (7). In this paper we will consider the value δ= 1/2 due to the considerations presented at the end of the following Section. Vertical plume trajectory In all cases where df is equal to zero, the horizontal motion of the plume until reaching the effective height can be considered negligible. Thus the trajectory of the gases is nearly vertical. 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 (6) and (7) this vertical motion is completely defined. tf=2 (zH−z0 c) wc (17) a0=−wc tf (18) w0(z) = wcs1−2 (z−z0 c) wctf (19) We remark that using the value of tffrom (17) in the equation (8), the coefficient of its cubic term is null. As we have commented in the previous Section this cancellation also occurs when δ= 1/2 is chosen in (16). It is also important to note that using the value of tffrom (17) in the equation (10) results (18). Thus, fixing δ= 1/2, the vertical motion in the bent-over plume case, given by equations (8), (9) and (10), and in the vertical plume case, given by equations (17), (18) and (19), are the same. This result makes compatible the vertical motion in both cases and justifies the election of δ= 1/2. Extension of the plume mean trajectory to 3D over complex terrains Several approximations can be considered to define the mean trajectory of bent curved plumes considering influence of complex terrains. Extension of vertical plumes is straightforward. In this paper, we propose a new solution that takes into account the trajectories obtained in the vertical plane case, and the non uniform three-dimensional ambient wind obtained in Section 2.2. A three-dimensional bent trajectory by an iterative process starting from the emission point (xc, yc, zc). We uniformly subdivide the interval of time [0, tf] defining n+1 time instants t0, t1, . . . ,tn, such that ti=i∆t for i= 0,1, . . . , n. For each time instant tia position of the mean trajectory is calculated. Starting in time t0at position x0= (xc, yc, zc), and known the point xi= (xi, yi, zi), xi+1can be computed from (xiwith the following formulae xi+1 =xi+u(xi, yi, zi) |(u(xi, yi, zi), v(xi, yi, zi))|u(xc, yc, zc)∆t+1 2ad∆t2(20) yi+1 =yi+v(xi, yi, zi) |(u(xi, yi, zi), v(xi, yi, zi))|v(xc, yc, zc)∆t+1 2ad∆t2(21) zi+1 =z(ti) (22) where z(ti) is defined by (8), and the horizontal movement is deduced from (14). The mean trajectory obtained by this method is a three-dimensional polygonal line, such that the longitude of its projection on the horizontal plane approximates the longitude df. In addition, the final height zncoincides with the effective height zH. Therefore, this method tries to verify the main values of the end of the plume considering Briggs’ equations. 6 Accepted for publication in Energy, Oct. 2012 2.3.2. Local mesh refinement along a Gaussian plume 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 along the Gaussian plume [56] until all the tetrahedra inside that region verify a size criteria. The region can be defined in function of the diffusion coefficients of the Gaussian plume (σy,σz). To decide if a sample point xof a tetrahedron Tis inside this region, the first step is to find its closest point xp∈x0,x1,...,xn. The second step is to calculate the distance d(tp) = p(x(tp)−xc)2+ (y(tp)−yc)2 where the time tp=p∆t, and x(tp), y(tp) are evaluated with equation (14). Then, we can compute the values σz(d(tp)) and σy(d(tp)) in function of the Pasquill stability class. These values define vertical ellipses (centred in xp) whose major and minor semi-axes are proportional to σyand σzrespectively. Finally, fixing a value K1>1, we consider that the sample point xof the tetrahedron Tis inside the Gaussian plume region if |x−xp|< K1σy. For the refinement of the marked tetrahedra, we use a local refinement algorithm [57] based on the 8-subtetrahedron subdivision developed in [58]. 2.3.3. Wind field perturbation along the plume rise Finally, a new ambient wind field uis obtained on the refined mesh with the mass-consistent model described in Section 2.2. The effect of the gas emission is introduced in this field by modifying its vertical component. In the case of bent curve trajectory, the vertical component of the ambient wind velocity uis modified at any point xof the domain Ω located inside the region that was defined in the previous Section. More precisely, the vertical component wof the velocity uin xis imposed by using equation (9) as w0(tp), being xpthe closest point to xin the discrete trajectory x0,x1,...,xnof the plume. In the case of vertical trajectory, the vertical component of the ambient wind velocity is modified inside a standard cylinder defined in the previous Section. That is, given a point x= (x, y, z) inside the cylinder, we impose the vertical velocity wof uby using equation (19) as w0(z). 2.4. Air pollution simulation The air pollution simulation consists in 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) (23) 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 (24) where ∇is the gradient with respect 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 nis 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. Scalar product “·” is applied nc times: 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 [59, 60, 61, 62]. For instance, detailed Volatile Organic Components decomposition involves hundreds of thousand reactions [63, 64] that needs special methodologies to reduce the number of the modelled reactions and 7 Accepted for publication in Energy, Oct. 2012 species. Reference models for gaseous phase reactions involve some tens of compounds [65, 66]. The most simplified models just involve about ten reactive species [67]. 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/ARM3 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/ARM3 model for the chemical term s(c). 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 (23) is solved with a Crank-Nicolson time integration scheme, and an spatial discretization with an stabilized finite element method, Least-Squares. Concentrations cnand cn+1 at times tnand tn+1 =tn+∆t are related using a Crank-Nicolson scheme as cn+1 =cn+∆t 2h∂cn+1 ∂t +∂cn ∂t i. And using Least-Squares we obtain a symmetric problem that can be written as the following equation system Bc c cn+1 =f(25) where c c cn+1 is the concentration vector approximation at tn+1 in the degrees of freedom of the finite element discretization, fdepends on c c cnand pollutant emissions, and Bis an square matrix 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 an sparse matrix storage. Since Bis a symmetric positive definite matrix, we have considered a solver based on a conjugate gradient method preconditioned with an incomplete Cholesky factorisation density type [68]. More details about the implementation of this system equation solver can be found in [35] and [36]. 2.4.2. Non-linear chemical problem To deal with the non-linearity of the reactive term in the convection–diffusion–reaction equation (23), we have considered an 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 Ropp et al. [38]. Using this method the complete equation is splited into two equations. The convection– diffusion equation that is solved using the same method proposed in the previous Section, being A= 0. And the non-linear chemical equations that are solved node by node with a second order Rosenbrock method (ROS2) [39]. To use the ROS2 method, the Jacobian square matrix of s(c) of dimension nchas to be computed. 3. Results In this section a realistic example is presented by using the methodology described in this paper. We propose to study a region in La Palma island, where several wind measurements data are given, and SO2and NO2emissions from an stack are considered. The topography of the island is real, from a digital elevation model, but the wind field measurements, and the stack location and emissions, are simulated. 3.1. Adaptive tetrahedral mesh The studied domain taken under consideration is a rectangular area with dimensions 15 600 m ×22 803 m. The topography of the terrain is highly complex ranging from the sea level up to a maximum height of 2279 m with several deep valleys. The upper boundary of the domain has been placed at h=9000 m. The digital elevation model of the area is defined over a uniform grid with a spacing step of 200 m in directions xand y. We add the stack geometry to the topographical data. Let us consider a stack with a height of 150 m over the terrain and the diameter at its top of 15 m, the location of the stack is UTM 28N 224760 3178743, and is shown in Figure 4. The mesh must be able to detect the details of the stack. Therefore, we impose an element size of about 2 m ×2 m in the stack. If we start from a uniform mesh τ1of the rectangular area with an element size of about 2000 m ×2000 m, we must make ten refinement steps close to the stack. However, 8 Accepted for publication in Energy, Oct. 2012 we make only five global refinement steps over the whole domain and five additional local refinements in the surroundings of the stack. We have considered the derefinement parameter of ε=40 m. Thus, the adapted mesh nears the terrain surface with an error less than that value. On the other hand, the node distribution of τ1is the one considered on the upper boundary of the domain. Figure 1 (a) shows the resulting tetrahedral mesh, with 198146 nodes and 1123012 tetrahedra. Figure 2 (a) shows a detail of the surroundings of the stack. Figure 3 (a) shows the distribution of the mesh element quality using the mean ratio algebraic shape quality metric. In order to have a more suitable mesh for the air quality problem, all the elements located below the height 3000 m have been refined. Figure 1 (b) shows the refined mesh with 455953 nodes and 2622454 tetrahedra. Figure 2 (b) shows a detail of the surroundings of the stack and Figure 3 (b) shows the distribution of the mesh element quality. (a) (b) Figure 1: La Palma tetrahedral mesh before (a) and after (b) the refinement under 3000 m (a) (b) Figure 2: Mesh detail around the stack before (a) and after (b) the refinement 3.2. Ambient wind field The wind field is constructed from the measurements of 4 stations. The locations of these stations are shown in Figure 4, and the values of the corresponding wind velocities are listed in Table 1. In the 9 Accepted for publication in Energy, Oct. 2012 (a) (b) (c) (d) (e) (f) 0.0e+00 2.5e-06 5.0e-06 7.5e-06 1.0e-05 Figure 10: Evolution of the immission concentration distribution (g m−3) of the primary pollutant SO2after (a) 5 min, (b) 10 min, (c) 15 min, (d) 20 min, (e) 25 min and (f) 30 min 16 Accepted for publication in Energy, Oct. 2012 (a) (b) (c) (d) (e) (f) 0.0e+00 2.5e-06 5.0e-06 7.5e-06 1.0e-05 Figure 11: Evolution of the immission concentration distribution (g m−3) of the secondary pollutant SO4 after (a) 5 min, (b) 10 min, (c) 15 min, (d) 20 min, (e) 25 min and (f) 30 min 17 Accepted for publication in Energy, Oct. 2012 (a) (b) (c) (d) (e) (f) Figure 12: Evolution of the isosurface corresponding to a concentration of 10 µ g m−3for the primary pollutant SO2after (a) 5 min, (b) 10 min, (c) 15 min, (d) 20 min, (e) 25 min and (f) 30 min 18 Accepted for publication in Energy, Oct. 2012 (a) (b) (c) (d) (e) (f) Figure 13: Evolution of the isosurface corresponding to a concentration of 10 µ g m−3for the secondary pollutant SO4after (a) 5 min, (b) 10 min, (c) 15 min, (d) 20 min, (e) 25 min and (f) 30 min 19 Accepted for publication in Energy, Oct. 2012 Figure 14: Wind field at 10 m over the terrain, and location of the stack and the points A-E 0 1 2 3 4 5 6 7 8 9 10 0 5 10 15 20 25 30 35 40 45 concentration [microgram/m3] time [min] A B C D E Figure 15: Immission concentration evolution of the primary pollutant at points A-E 20 Accepted for publication in Energy, Oct. 2012 (C,D,E) that are located at different heights behind the crest of a mountain. The approximated horizontal distance from the stack to this three points are 3000 m,4000 m and 5000 m respectively. Figure 15 shows the evolution of the immission concentration at the points (A-E). It can be observed the influence of the distance to the stack and the orography. The numerical results verifies the qualitative standard behaviour in the first group of points (concentration decreases with the distance). Moreover, reasonable values are reached behind the crest of the mountain. Note that the value of the steady concentration at the deepest point (C) is lower than the values at points D and E. A final comment about the computational complexity of the evolution process should be done. For each time step we have to solve a finite element problem with a number of degrees of freedom about the number of nodes multiplied by the number of species, i.e. 455953 ×4 = 1823812. The number of time steps in the simulation period (about 30 min) is 30 ×60/10 = 180. Therefore, in the whole evolution process about 180 linear equation systems with 1823812 unknowns must be solved. The computational cost corresponding to the mesh generation, wind simulation, and the resolution of the ordinary differential equations in the splitting method are insignificant with respect to the resolution of the unsteady convection–diffusion equation. In a computer with 128GB of RAM memory and 2.34 GHz, the total computing time is about 40 min. In a future work, the present computational complexity will be significantly reduced by using a refinement/derefinement strategy that follows the front of the pollutant plume, minimising the number of degrees of freedom in each time step. 4. Conclusions We have introduced a new methodology for solving air pollution problems over a complex terrain. The adaptive three-dimensional mesh generation discretize domains defined over complex terrains which include several stack geometries, with a minimal user intervention and low computational cost. The mass-consistent model obtains an ambient wind field that takes into account the complex orography. The local mesh refinement along the Gaussian plume, allows to perturb the ambient wind field to introduce the effect of the pollutant emissions. The convection–diffusion–reaction equation obtains the values of concentration for all the pollutants in the whole three-dimensional domain. The proposed methodology uses the necessary mesh resolution to consider the main effects of the physical phenomena that are involved in air pollution process. The numerical experiment shows a reasonable behaviour of the proposed method. The present manuscript introduce a solid framework to develop more accurate models in the future. Acknowledgements This work has been supported by the Spanish Government, “Ministerio de Econom´ıa y Competitividad”, and FEDER, grant contracts: CGL2011-29396-C03-01, CGL2008-06003-03-01, CGL2008-06003-0302, UNLP08-3E-010, CSD2006-00032C; and by CONACYT-SENER (“Fondo Sectorial CONACYT SENER HIDROCARBUROS”, grant contract: 163723). Finally, the authors wish to express their sincere thanks to the three anonymous reviewers, whose comments and suggestions were very helpful. References [1] Seigneur C. Air pollution: Current challenges and future opportunities. AIChE Journal 2005;51(2):356–64. doi: 10.1002/aic.10458. [2] Nazari S, Shahhoseini O, Sohrabi-Kashani A, Davari S, Paydar R, Delavar-Moghadam Z. Experimental determination and analysis of co2, so2 and nox emission factors in irans thermal power plants. Energy 2010;35(7):2992–8. doi: 10.1016/j.energy.2010.03.035. [3] Olcese LE, Toselli BM. Development of a model for reactive emissions from industrial stacks. Environmental Modelling & Software 2005;20(10):1239–50. doi:10.1016/j.envsoft.2004.08.008. [4] Hernandez J, Cremades L, Baldasano J. Dispersion modelling of a tall stack plume in the spanish mediterranean coast by a particle model. Atmospheric Environment 1995;29(11):1331–41. doi:10.1016/1352-2310(94)00346-M. [5] Bourque CA, Arp P. Simulating sulfur dioxide plume dispersion and subsequent deposition downwind from a stationary point source: A model. Environmental Pollution 1996;91(3):363–80. doi:10.1016/0269-7491(95)00041-0. 21 Accepted for publication in Energy, Oct. 2012 [6] Souto M, Souto J, P´erez-Mu˜nuzuri V, Casares J, Berm´udez J. A comparison of operational lagrangian particle and adaptive puff models for plume dispersion forecasting. Atmospheric Environment 2001;35(13):2349–60. doi:10.1016/S13522310(00)00537-9. [7] Walcek CJ. Effects of wind shear on pollution dispersion. Atmospheric Environment 2002;36(3):511–7. doi:10.1016/S13522310(01)00383-1. [8] Hurley PJ, Physick WL, Luhar AK. TAPM: a practical approach to prognostic meteorological and air pollution modelling. Environmental Modelling & Software 2005;20(6):737–52. doi:10.1016/j.envsoft.2004.04.006. [9] Morris R, Yarwood G, Emery C, Wilson G. Recent advances in photochemical air quality modeling using the CAMx Model: Current update and ozone modeling of point source impacts. In: Proceedings of The Air & Waste Management Association’s 95th Annual Conference & Exhibition. A&WMA; 2002,Paper #43180. [10] Taghavi M, Cautenet S, Arteta J. Impact of a highly detailed emission inventory on modeling accuracy. Atmospheric Research 2005;74(1-4):65–88. doi:10.1016/j.atmosres.2004.06.007. [11] San Jos´e R, P´erez JL, Gonz´alez RM. An operational real-time air quality modelling system for industrial plants. Environmental Modelling & Software 2007;22(3):297–307. doi:10.1016/j.envsoft.2005.07.030. [12] ENVIRON International Corporation . CAMx Comprehensive Air Quality Model with extensions. User’s Guide v. 5.40. ENVIRON International Corporation. Novato, California.; 2011. [13] Byun D, Ching J. Science Algorithms of the EPA Models–3 Community Multiscale Air Quality (CMAQ) Modelling System. EPA/600/R–90/030. Atmospheric Modeling Division. U.S. Environmental Protection Agency. Washington DC.; 1999. [14] ICF Consulting . User’s Guide to the variable–grid Urban Airshed Model (UAM-V). Systems Applications International, Inc., ICF Consulting, Inc., San Rafael, California.; 1999. [15] Ryerson TB, Trainer M, Holloway JS, Parrish DD, Huey LG, Sueper DT, et al. Observations of ozone formation in power plant plumes and implications for ozone control strategies. Science 2001;292(5517):719–23. doi:10.1126/science.1058113. [16] Lee SM, Fernando HJ, Princevac M, Zajic D, Sinesi M, McCulley JL, et al. Transport and diffusion of ozone in the nocturnal and morning planetary boundary layer of the phoenix valley. Environmental Fluid Mechanics 2003;3:331–62. [17] Mauzerall DL, Sultan B, Kim N, Bradford DF. NOx emissions from large point sources: variability in ozone production, resulting health damages and economic costs. Atmospheric Environment 2005;39(16):2851–66. doi: 10.1016/j.atmosenv.2004.12.041. [18] He K, Lei Y, Pan X, Zhang Y, Zhang Q, Chen D. Co-benefits from energy policies in china. Energy 2010;35(11):4265 –72. doi:10.1016/j.energy.2008.07.021. [19] Czarnowska L, Frangopoulos CA. Dispersion of pollutants, environmental externalities due to a pulverized coal power plant and their effect on the cost of electricity. Energy 2012;41(1):212–9. doi:10.1016/j.energy.2011.08.004. [20] Hidy GM. Multiscale impact of fuel consumption on air quality. Energy & Fuels 2002;16(2):270–81. doi:10.1021/ef0101659. [21] Hanjali´c K, Kenjereˇs S. Dynamic simulation of pollutant dispersion over complex urban terrains: A tool for sustainable development, control and management. Energy 2005;30(8):1481 –97. doi:10.1016/j.energy.2004.05.001. [22] Lagzi I, K´arm´an D, Tur´anyi T, Tomlin AS, Haszpra L. Simulation of the dispersion of nuclear contamination using an adaptive eulerian grid model. Journal of Environmental Radioactivity 2004;75(1):59–82. doi:10.1016/j.jenvrad.2003.11.003. [23] Ahmad NN, Bacon DP, Hall MS, Sarma A. Application of the multidimensional positive definite advection transport algorithm (mpdata) to environmental modelling on adaptive unstructured grids. International Journal for Numerical Methods in Fluids 2006;50(10):1247–68. doi:10.1002/fld.1113. [24] Tomlin A, Ghorai S, Hart G, Berzins M. 3-D multi-scale air pollution modelling using adaptive unstructured meshes. Environmental Modelling & Software 2000;15(6-7):681–92. doi:10.1016/S1364-8152(00)00038-4. [25] Montenegro R, Montero G, Escobar JM, Rodr´ıguez E, Gonz´alez-Yuste JM. 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] Montero G, Montenegro R, Escobar JM, Rodr´ıguez E, Gonz´alez-Yuste JM. Velocity field modelling for pollutant plume using 3-d adaptive finite element method. In: Bubak M, van Albada G, Sloot P, Dongarra J, editors. Computational Science - ICCS 2004; vol. 3037 of Lecture Notes in Computer Science. Springer Berlin / Heidelberg. ISBN 978-3-54022115-9; 2004, p. 642–5. [27] Winter G, Montero G, Ferragut L, Montenegro R. Adaptive strategies using standard and mixed finite elements for wind field adjustment. Solar Energy 1995;54(1):49–56. doi:10.1016/0038-092X(94)00100-R. [28] Montero G, Montenegro R, Escobar JM. A 3-D diagnostic model for wind field adjustment. Journal of Wind Engineering and Industrial Aerodynamics 1998;74-76(0):249–61. doi:10.1016/S0167-6105(98)00022-1. [29] Montero G, Rodr´ıguez E, Montenegro R, Escobar JM, Gonz´alez-Yuste JM. Genetic algorithms for an improved parameter estimation with local refinement of tetrahedral meshes in a wind model. Advances in Engineering Software 2005;36(1):3–10. doi:10.1016/j.advengsoft.2004.03.011. [30] Ferragut L, Montenegro R, Montero G, Rodr´ıguez E, Asensio M, Escobar JM. Comparison between 2.5-D and 3-D realistic models for wind field adjustment. Journal of Wind Engineering and Industrial Aerodynamics 2010;98(10-11):548–58. doi: 10.1016/j.jweia.2010.04.004. [31] Mart´ın MJ, Singh DE, Mouri˜no JC, Rivera FF, Doallo R, Bruguera JD. High performance air pollution modeling for a power plant environment. Parallel Computing 2003;29(11-12):1763–90. doi:10.1016/j.parco.2003.05.018. [32] Chock DP, Whalen MJ, Winkler SL, Sun P. Implementing the trajectory-grid transport algorithm in an air quality model. Atmospheric Environment 2005;39(22):4015–23. doi:10.1016/j.atmosenv.2005.03.037. [33] Saylor RD, Ford GD. On the comparison of numerical methods for the integration of kinetic equations in atmospheric chemistry and transport models. Atmospheric Environment 1995;29(19):2585–93. doi:10.1016/1352-2310(95)00187-4. 22 Accepted for publication in Energy, Oct. 2012 [34] Sandu A, Verwer J, Loon MV, Carmichael G, Potra F, Dabdub D, et al. Benchmarking stiff ODE solvers for atmospheric chemistry problems-I. Implicit vs explicit. Atmospheric Environment 1997;31(19):3151–66. doi:10.1016/S13522310(97)00059-9. [35] Donea J, Huerta A. Finite Element Methods for Flow Problems. John Wiley and Sons Ltd, West Sussex.; 2003. [36] Rodr´ıguez-Ferran A, Sandoval ML. Numerical performance of incomplete factorizations for 3D transient convection– diffusion problems. Advances in Engineering Software 2007;38(6):439–50. doi:10.1016/j.advengsoft.2006.09.003. [37] Scire JS, Strimaitis DG, Yamartino RJ. A User’s Guide for the Calpuff Dispersion Model (version 5). Earth Tech., Inc, Concord, MA.; 2000. [38] Ropp DL, Shadid JN, Ober CC. Studies of the accuracy of time integration methods for reaction–diffusion equations. Journal of Computational Physics 2004;194(2):544–74. doi:10.1016/j.jcp.2003.08.033. [39] Verwer JG, Spee EJ, Blom JG, Hundsdorfer W. A second-order Rosenbrock method applied to photochemical dispersion problems. SIAM Journal on Scientific Computing 1999;20(4):1456–80. doi:DOI:10.1137/S1064827597326651. [40] P´erez-Foguet A, Oliver A, Escobar JM, Rodr´ıguez E. Finite element simulation of chimney emissions: A proposal for near field impact assessment in highly complex terrains. In: Topping B, Montero G, Montenegro R, editors. Procedings of the Fifth International Conference on Engineering Computational Technology. 2006,Paper #101. [41] Ferragut L, Montenegro R, Plaza A. Efficient refinement/derefinement algorithm of nested meshes to solve evolution problems. Communications in Numerical Methods in Engineering 1994;10(5):403–12. doi:10.1002/cnm.1640100506. [42] Escobar JM, Montenegro R. Several aspects of three-dimensional Delaunay triangulation. Advances in Engineering Software 1996;27(1-2):27–39. doi:10.1016/0965-9978(96)00006-3. [43] Escobar JM, Rodr´ıguez E, Montenegro R, Montero G, Gonz´alez-Yuste JM. Simultaneous untangling and smoothing of tetrahedral meshes. Computer Methods in Applied Mechanics and Engineering 2003;192(25):2775–87. doi:10.1016/S00457825(03)00299-8. [44] Montenegro R, Montero G, Escobar JM, Rodr´ıguez E. Efficient strategies for adaptive 3-D mesh generation over complex orography. Neural, Parallel & Scientific Compututation 2002;10:57–76. [45] Montenegro R, Montero G, Escobar JM, Rodr´ıguez E, Gonz´alez-Yuste JM. Tetrahedral mesh generation for environmental problems over complex terrains. In: Sloot P, Hoekstra A, Tan C, Dongarra J, editors. Computational Science ICCS 2002; vol. 2329 of Lecture Notes in Computer Science. Springer Berlin / Heidelberg. ISBN 978-3-540-43591-4; 2002, p. 335–44. [46] Panofsky H, Dutton J. Atmospheric turbulence. Models and methods for engineering applications. New York: Wiley; 1984. [47] McRae GJ, Goodin WR, Seinfeld JH. Development of a second-generation mathematical model for urban air pollution–I. Model formulation. Atmospheric Environment (1967) 1982;16(4):679–96. doi:DOI: 10.1016/0004-6981(82)90386-9. [48] Zannetti P. Air Pollution Modeling. Boston: Computational Mechanics Publications; 1990. [49] Lalas D, Ratto C. Modelling of Atmospheric Flow Fields. Singapore: World Scientific Publishing; 1996. [50] Montero G, San´ın N. 3-D modelling of wind field adjustment using finite differences in a terrain conformal coordinate system. Journal of Wind Engineering and Industrial Aerodynamics 2001;89(5):471–88. doi:10.1016/S0167-6105(00)000751. [51] Degani AT, Smith FT, Walker JDA. The structure of a three-dimensional turbulent boundary layer. Journal of Fluid Mechanics 1993;250:43–68. doi:10.1017/S0022112093001375. [52] Briggs GA. Optimum Formulas for Buoyant Plume Rise. Philosophical Transactions of the Royal Society of London Series A, Mathematical and Physical Sciences (1934-1990) 1969;265(1161):197–203. [53] Briggs GA. Plume rise. Tech. Rep.; U.S. Atomic Energy Commission, Division of Technical Information; 1969. [54] Moore DJ. A comparison of the trajectories of rising buoyant plumes with theoretical/empirical models. Atmospheric Environment (1967) 1974;8(5):441 –57. doi:10.1016/0004-6981(74)90060-2. [55] Schulman LL, Strimaitis DG, Scire JS. Development and evaluation of the prime plume rise and building downwash model. Journal of the Air & Waste Management Association 2000;50(3):378–90. doi:10.1080/10473289.2000.10464017. [56] Green AES, Singhal RP, Venkateswar R. Analytic extensions of the gaussian plume model. Journal of the Air Pollution Control Association 1980;30(7):773–6. [57] Gonz´alez-Yuste JM, Montenegro R, Escobar JM, Montero G, Rodr´ıguez E. Local refinement of 3-D triangulations using object-oriented methods. Advances in Engineering Software 2004;35(10-11):693 – 702. doi: 10.1016/j.advengsoft.2003.07.003. [58] L¨ohner R, Baum JD. Adaptive h-refinement on 3d unstructured grids for transient problems. International Journal for Numerical Methods in Fluids 1992;14(12):1407–19. doi:10.1002/fld.1650141204. [59] Finlayson-Pitts BJ, Pitts JN. Tropospheric air pollution: Ozone, airborne toxics, polycyclic aromatic hydrocarbons, and particles. Science 1997;276(5315):1045–51. doi:10.1126/science.276.5315.1045. [60] Kley D. Tropospheric chemistry and transport. Science 1997;276(5315):1043–4. doi:10.1126/science.276.5315.1043. [61] Andreae MO, Crutzen PJ. Atmospheric aerosols: Biogeochemical sources and role in atmospheric chemistry. Science 1997;276(5315):1052–8. doi:10.1126/science.276.5315.1052. [62] Ravishankara AR. Heterogeneous and multiphase chemistry in the troposphere. Science 1997;276(5315):1058–65. doi: 10.1126/science.276.5315.1058. [63] Atkinson R, Arey J. Atmospheric degradation of volatile organic compounds. Chemical Reviews 2003;103(12):4605–38. doi:10.1021/cr0206420. [64] Szopa S, Aumont B, Madronich S. 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 2005;5(9):2519–38. doi:10.5194/acp-5-2519-2005. [65] Jimenez P, Baldasano JM, Dabdub D. Comparison of photochemical mechanisms for air quality modeling. Atmospheric 23 Accepted for publication in Energy, Oct. 2012 Environment 2003;37(30):4179–94. doi:10.1016/S1352-2310(03)00567-3. [66] Kirchner F. The chemical mechanism generation programme CHEMATA – Part 1: The programme and first applications. Atmospheric Environment 2005;39(6):1143–59. doi:10.1016/j.atmosenv.2004.09.086. [67] Zlatev Z. Computer Treatment of Large Air Pollution Models. Kluver Academic Publishers, Dordrecht.; 1995. [68] Lin CJ, Mor´e JJ. Incomplete Cholesky factorizations with limited memory. SIAM Journal on Scientific Computing 1999;21(1):24–45. doi:DOI:10.1137/S1064827597327334. 24 Accepted for publication in Energy, Oct. 2012