Finite element simulation of a local scale air quality model over complex terrain

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 e ff ct in the wind field, and the simulation of transport and reaction of pollutants. We apply our methodology to simulate a fictitious pollution episode in La Palma island (Canary Island, Spain).


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 (Hanjalić and Kenjereš, 2005).
Three remarkable uses of unstructured meshes in atmospheric pollution problems are the two-dimensional regionalglobal examples presented in Lagzi et al. (2004); Ahmad et al. (2006), the three-dimensional regional examples, including local refinement with element sizes of 2 km, presented in Tomlin et al. (2000), 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. (2004); Montero et al. (2004).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 areas.In order to simulate it, we have used a mass-consistent model.Several two-dimensional (Winter et al., 1995) and three-dimensional (Montero et al., 1998(Montero et al., , 2005;;Ferragut et al., 2010) adaptive finite element solutions have been developed by the authors.
The convection, diffusion and reaction problem is usually solved using splitting schemes (Martín et al., 2003;Chock et al., 2005) and specific numerical solvers for time integration of photochemical reaction terms (Saylor and Ford, 1995;Sandu et al., 1997;Ahmad et al., 2006).In this paper a stabilized finite element formulation with a Crank-Nicolson temporal integration is proposed to solve the problem (Donea and Huerta, 2003;Rodríguez-Ferran and Sandoval, 2007).The chemistry is simulated with the RIVAD chemical model (Scire et al., 2000).The transport and chemical terms are treated separately with Strang splitting operators (Ropp et al., 2004), and the non-linear chemical part is solved with a second order Rosenbrock method (Verwer et al., 1999).
The paper is organised as follows.In Sect. 2 we describe the main steps of the proposed methodology.Results are shown in Sect.3, and finally the conclusions and future work are presented in Sect. 4.

Algorithm description
In this section a brief description of the different steps of Algorithm 1 is presented.
Published by Copernicus Publications.

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 (Ferragut et al., 1994) 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 (Escobar and Montenegro, 1996) is applied.Finally, the untangling and smoothing procedure described in (Escobar et al., 2003) is used to get a valid mesh and to improve its quality.A detailed description of the mesh generation procedure can be seen in Montenegro et al. (2002).

Wind field simulation
A mass-consistent model (Montero et al., 1998(Montero et al., , 2005;;Ferragut et al., 2010) is used to compute a wind field u in the three-dimensional domain Ω, with a boundary Γ = Γ a ∪ Γ b , that satisfies the continuity equation and the impermeability condition on the terrain Γ a , where n is 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 u 0 = (u 0 , v 0 , w 0 ).The adjusting functional for a field v = ( u, v, w) is defined as where P is a 3 × 3 diagonal matrix with P 1,1 = P 2,2 = 2α 2 1 and P 3,3 = 2α 2 2 .The Lagrange multiplier technique is used to minimise the functional (Eq.2), with the restrictions (Eq.1).Considering the Lagrange multiplier λ, the Lagrangian is de-fined as and the solution u is obtained by finding the saddle point (u, φ) of the Lagrangian Eq. ( 3).This resulting wind field satisfies the Euler-Lagrange equation, where φ is the Lagrange multiplier.As α 1 and α 2 are constant in Ω, the variational approach results in an elliptic problem in φ, by substituting Eq. ( 4) in Eq. ( 1), that is solved by using the finite element method.
The interpolated wind field u 0 can 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 u 0 in the whole computational domain.

Horizontal interpolation
The most common technique of interpolation at a given point, placed at a height z m over the terrain, is formulated as a function of the inverse of the squared distance between that point and the measurement stations, and the inverse of their height differences (Montero et al., 1998) where the value of u n is the velocity measured at station n, N is the number of stations considered in the interpolation, d n is the horizontal distance from station n to the point of the domain where we are computing the wind velocity, |∆h n | is the height difference between station n and 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 (Lalas and Ratto, 1996) in the surface layer, which takes into account the horizontal interpolation (Montero and Sanín, 2001) 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, where u * is the friction velocity, k is von Karman's constant, z 0 is the roughness length (McRae et al., 1982) and z sl is the height of the surface layer.The values of Φ m depend on the Pasquill stability class (Zannetti, 1990), and the friction velocity is obtained from Eq. ( 9) at any point (x, y) by using the horizontal interpolated velocity u 0 (z m ).
The linear interpolation is given by, where u g is the geostrophic wind, z pbl is the height of the planetary boundary layer, and ρ(z) is defined as Finally, this model assumes

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 (Briggs, 1969a,b;Moore, 1974).These works 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 (Olcese and Toselli, 2005) approximate the effective height of a plume z H and the horizontal distance d f from the stack to the point where the plume height reaches z H , 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 d f 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, satisfying certain known conditions.
The vertical motion along the mean trajectory of the plume is defined by an acceleration a 0 (t), a velocity w 0 (t) and z(t), from the initial time t = 0 to the final time t = t f when the plume reaches the effective height, satisfying the following conditions Since there are four conditions on the vertical motion, we propose a cubic approximation of z(t), and therefore a quadratic approximation of w 0 (t), and a linear approximation of a 0 (t).
The horizontal motion is defined by a uniformly accelerated motion, with a constant positive acceleration vector a d = (a dx , a dy ), a velocity u d (t) = (u d (t), v d (t)), and an horizontal relative position vector d(t) = (x(t) − x c , y(t) − x c ) with respect to the centre of the stack, satisfying the following conditions 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 d f .In addition, the final height coincides with the effective height z H . 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 d f is 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.
In this case, we propose a vertical motion along the trajectory of the plume with a constant negative acceleration a 0 , a linear velocity w 0 (t) and a quadratic trajectory z(t).Imposing the conditions ( 12) and ( 13), this vertical motion is completely defined.
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 (González-Yuste et al., 2004) along the Gaussian plume (Green et al., 1980) until all the tetrahedra inside that region fulfill a size criterion.

A. Oliver et al.: FEM simulation of local scale air quality model
Finally, a new ambient wind field u is obtained on the refined mesh with the mass-consistent model described in Sect.2.2.The effect of the gas emission is introduced in this field by modifying its vertical component along the plume.

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: for the spatial coordinates x and time t, (x, t) ∈ Ω × (0, t end ], with initial condition c(x, 0) = c ini (x) on x ∈ Ω, and the following boundary conditions: where ∇ is the gradient with respect to x, and c, u, e and s(c) are respectively the concentration, the perturbed wind velocity, the emission and the chemical vectors with a dimension n c (the number of pollutant species).K is the diffusion matrix of dimension 3 × n c , V d is the deposition diagonal matrix with dimension n c , and n is the outward-pointing normal unit vector, c emi is the concentration of the emission in the top of the stack, and c out the outside concentration at the inlet wind boundaries.Scalar product "•" is applied n c times: the first argument is multiplied by each one of the n c components of the second argument.
The complete description of photochemical reaction of atmospheric species is highly complex (Finlayson-Pitts and Pitts, 1997; Kley, 1997;Andreae and Crutzen, 1997;Ravishankara, 1997).For instance, detailed Volatile Organic Components decomposition involves hundreds of thousand reactions (Atkinson and Arey, 2003;Szopa et al., 2005) 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 (Jimenez et al., 2003;Kirchner, 2005).The most simplified models just involve about ten reactive species (Zlatev, 1995).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 (Scire et al., 2000).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 c involves four species, , and the components of the reaction vector s(c) are: where α 1 (c) = γ 1 /(c 1 + δ 1 c 3 ) and α 3 (c) = γ 3 /(c 3 + δ 3 c 1 ).Note that for values close to zero of the concentration of the primary species c 1 and c 3 , 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 Eq. ( 19).The non-linear case focuses on the development of an splitting method that combines the solution of Eq. ( 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.

Linear chemical problem
In this case the chemical term is linear, that is, s(c) = Ac where A is constant matrix.The resulting Eq. ( 19) is solved with a Crank-Nicolson time integration scheme, and an spatial discretization with a stabilized finite element method, Least-Squares.Concentrations c n and c n+1 at times t n and t n+1 = t n + ∆t are related using a Crank-Nicolson scheme as c n+1 = c n + ∆t 2 ∂c n+1 ∂t + ∂c n ∂t .We define the differential operator L as and a function F in Ω from the known c n Applying the Crank-Nicolson scheme, we can rewrite Eq. (( 19)) as where I is the identity operator.
Using Least-Squares we obtain a symmetric problem such that the weak form of the Eq. ( 25) is where ν is the test function, and (µ, ν) = Ω µν dΩ is the inner product We define the next bilinear forms 28) where (µ, ν) e = Ω e µν dΩ, and e represents the sum over all the mesh elements.Applying the operators L and F in the weak form (Eq. 26), we obtain where e n+ 1 2 = e n+1 +e n 2 .Equation ( 30) can be written as an equation system where c c c n+1 is the concentration vector approximation at t n+1 in the degrees of freedom of the finite element discretization, + e e e n ], and B and F are square matrices with dimension (n c × n dof ), being n dof 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 B is a symmetric positive definite matrix, we have considered a solver based on a conjugate gradient method preconditioned with an incomplete Cholesky factorisation density type (Lin and Moré, 1999).The left-preconditioning is used to improve the convergence of the conjugate gradient method.The original linear system is transformed into T −1 Bc c c n+1 = T −1 f , where T is 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 B can be amortized over many time-steps.More details about the implementation of this system equation solver can be found in Donea and Huerta (2003) and Rodríguez-Ferran and Sandoval (2007).

Non-linear chemical problem
To deal with the non-linearity of the reactive term in the convection-diffusion-reaction Eq. ( 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 Ropp et al. (2004): Once we have split the Eq. ( 19), we solve three equations in different time steps; the reaction Eq. ( 32), the convectiondiffusion Eq. ( 33), and the reaction Eq. ( 34), being finally c n+1 (x) = c * * * (x, ∆t).
The convection-diffusion Eq. ( 33) is solved using the same method proposed in the previous Section, with A = 0.The non-linear chemical Eqs. ( 32) and ( 34) are solved node by node with a second order Rosenbrock method (ROS2) (Verwer et al., 1999).To use the ROS2 method, the Jacobian square matrix of s(c) of dimension n c has to be computed.

Results
The approach has been used to simulate the transport and reaction of pollutants from a fictitious stack in La Palma island.
The studied domain taken under consideration is a rectangular area with dimensions 15 600 × 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 x and y.We add the stack geometry to the topographical data, a stack with a height of 150 m over the terrain and the diameter at its top of 15 m.
The wind field is obtained from four meteorological stations placed in the studied region.Figure 1 represents the interpolated and resulting wind field at a height of 20 m above the terrain.While the interpolated wind field is almost uniform and crosses the terrain surface, the resulting wind field has higher velocities in the peaks of the mountains and follows the terrain and valleys, verifying the wind incompressibility and terrain impermeability conditions.Figure 2 shows the plume rise region where a local mesh refinement has been performed.Figure 3 represents the streamlines of the modified wind field from the top of the stack.Note that streamlines follow the trajectory of the plume rise, so the vertical component of the wind field has been successfully modified.
The plume rise has been calculated using the following values.The stack exit velocity is 5 m s −1 and the gas temperature is 573 K.The velocity of the wind field at the top of the stack is V o = 7.13 m s −1 .With these values, it results that the effective height of the plume is z H = 2347.52m and its horizontal distance is d f = 3700.04m.Finally the transport and reaction of the pollutants has been simulated.The concentration of the primary pollutant at the top of the stack has been fixed to 6 g m −3 .We have considered a horizontal diffusion of 8×10 −6 m 2 s −1 , and a vertical diffusion of 4×10 −6 m 2 s −1 .The time-dependent problem has been simulated with a time step of 10 s.
Figure 4 shows the immission concentration distributions at 30 min for SO 2 and SO 4 , respectively.These Figures are interesting since they give information about the pollutant concentration at the ground-level.Note that primary pollutant tends to have the highest concentrations near the emission source, while the highest concentrations of the secondary pollutant are located further.This is consistent with the chemical reaction effect.
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. 455 953×4 = 1 823 812.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 1 823 812 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 128 GB 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.

Conclusions
The presented methodology promises to be useful to simulate air quality over complex terrains.The modified wind field has more reasonable trajectories and magnitudes than a simple interpolation of the wind data.The local mesh refinement along the Gaussian plume, allows to perturb the ambient wind field to introduce the effect of the pollutant emis-sions.The convection-diffusion-reaction equation obtains the values of concentration for all the pollutants in the whole three-dimensional domain.To really prove the usefulness of the proposed methodology, it will be validated against measured data in the near future.

Figure 1 :
Figure 1: Interpolated (a) and resulting (b) wind field (m s −1 ) at 20 m over the terrain.

Figure 2 :
Figure 2: Refined region along the plume rise.

Figure 3 :
Figure 3: Modified wind field streamlines from the top of stack.

Figure 4 :
Figure 4: Immission concentration distribution (g m −3 ) on the terrain of the primary (a) and secondary (b) pollutants at 30 min.