ON THE USAGE OF ELLAM TO SOLVE ADVECTION-DIFFUSION EQUATION DESCRIBING THE POLLUTANT TRANSPORT IN PLANETARY BOUNDARY LAYER

The paper deals with the numerical solution of the advection diffusion equation describing the pollutant transport in the atmosphere. The idea is to use the not so common ELLAM framework and to compare the results with the state of the art Walcek method, which is used to solve the advection problems. The real wind model was chosen for the tests in order to get the good performance idea of the two methods. From the performed experiments and the calculation times, one can conclude that ELLAM is suitable to solve the presented problem.


INTRODUCTION
The air pollution modelling is an important and very current topic.It allows predicting the progress of the pollutant dispersion from the specific time to the near future and thus to help to deal with the low-quality air or to prevent from the possible contamination.
The model of the pollutant transport is described by the specific partial differential equation 1.It consists of the several parts that describe the whole process.The first and the most important part from the pollutant behaviour point of view is the advection term.It determines the wind field, commonly changing during time and space, which is the most influencing term in the equation.The diffusion/dispersion is the second important part and also it changes during time and space.The rest of the equation can be summarized into the reaction term which includes the behaviour of the reactant in the atmosphere.All three mentioned parts are in balance with the other side of the equation.It consists of source term the source(s) of the pollution.
Variable c is the concentration, u is the velocity field, D is the diffusion-dispersion model, R is the reaction model and t is time.
The model is usually solved by operator splitting technique, where the equation is separated into the advection, diffusion and reaction models.All the models are then evaluated separately and the calculated concentration is counted together.Although the special simple methods can be used to solve the separate parts, the error due to splitting, which is often neglected, arises [3].
One of such a method was proposed by [8] and it was successfully tested and compared against up to its date best methods that solve the pure advection part of the ADR equation.Consequently, another method, more precisely the framework, for solution of advection dominated pollutant dispersion model was being developed.Although the ELLAM is much more complicated than the Walceks method, it has the advantage of incorporating of the other part, diffusion or reaction and the very important boundary conditions, mostly these near the ground.ELLAM has many variants that were mostly adapted to water fluid environments were the behaviour of the flows is different from the atmospheric turbulent flows.
In this paper the form of ELLAM was tested against the Walceks method.Beside the artificial tests the real wind models with the truly measured data from the performed experiments were chosen for the schemes evaluation.
The paper is organized as follows.In the following section the general and the used ELLAM scheme are outlined.The scheme evaluation including the test inputs and methodology is then described in the ELLAM Evaluation section.Some conclusion and the future work is stated at the end of the paper.

ELLAM
ELLAM method is based on a philosophy of algebraic theory by [2].In this theory, the test functions are used to define the weak form of the governing equation.In the following subsections the general framework as well as the modified version used in this work is presented.

General ELLAM Framework
The procedure of the general ELLAM framework is based on the idea of conversion of equation 1 to the corresponding weak formulation.This is done by multiplying of the governing equation by the test function.Next the whole equation is integrated using Greens formula.
The obtained equation contains several integral terms that have to be evaluated.If we choose the proper form of the test function we can get rid of the one of the integral.This leads from the adjoint equation of equation 1 [6].
The integrals are evaluated analytically for the simple cases, however for practical problems one has to use the numerical approximation.The time discretization is usually done by backward Euler or Runge-Kutta methods.It is an advantage of the scheme that ELLAM can be used generally together within finite difference, finite volume or finite element approaches.The other advantage of the EL-LAM framework is its ability to naturally incorporate the boundary conditions.They simply appear as other integral terms in the equation.
The important part of the ELLAM scheme is the accurate characteristic tracking of the points.The problem of characteristics tracking is described by the ordinary differential equations, thus the solution can be obtained by various numerical methods.The more details on general EL-LAM method can be found for example in [6].

The Presented ELLAM Implementation
Our current implementation of ELLAM framework is in two dimensions (Ω), it is designed for the advectiondiffusion equation and for a rectangular grid.We came out from the work of [4] where the space discretization is based on finite element method.The governing equation is: where f is the function of the source of the pollution and all other variables have the same meaning as in equation 1.The resulting weak formulation for the specified time t n after multiplication by test function w( x,t n ) and applying of Greens formula is: where n is the normal outward unit vector from the element dydt, J n is time domain, Γ n is boundary domain, ).The second integral on the left hand side of the equation 3 is a diffusion term, the third integral is a boundary term and the second integral on the right hand side is a source term.
To evaluate the equation 3, the following procedure is done.The test function was chosen as piecewise linear as is usual for the ELLAM scheme [6].The terms with integration by time are approximated by backward Euler method.The remaining integrals can be evaluated by numerical integration using for example Gaussian quadrature with appropriate integration points.It remains to evaluate the equation w x,t ).This problem leads to the solution of the ordinary differential equation back in time.The common integration methods such as Euler can be used.In our case we decided to use the 4th order Runge-Kutta method.It is a trade-off between speed and accuracy and it behaved very well in cases of the performed experiments.
The last think to explain is the space discretization.We used the rectangular grid of points and the standard FEM process.The equation 3 has to be solved on the whole domain, therefore the elements, on which the approximation of the unknown function c is defined, have to be assembled together.This leads to the system of algebraic equations that has to be solved at each time step.

ELLAM EVALUATION
The most of the tests done for the numerical schemes use the artificial wind velocity fields and then they are evaluated using them.The wind velocities do not have to represent the performance of the scheme in concrete practical applications.Therefore, we have decided to do the tests of the ELLAM scheme against the scheme presented by [8] on the wind velocity fields based on real wind models and real measured data.

Used Wind Model
The used wind model is in the stationary form and it depends on the z (height) variable.
The special coefficients are needed in order to calculate the wind speed at given position.These are Monin-Obukhov length (L), roughness length (z 0 ), von Krmn constant (k) and friction velocity (u * ) [7,9].Then, the wind speed is given by: The coefficient z b = min (|L| , 0.1 • h), where h is a height of the unstable boundary layer.
The function Ψ m is of the form: (5)

Experiment Data
The above explained wind model has the specific parameters that differ in cases of different experiments.The experiments done in Copenhagen [1] were chosen to completely describe the wind model.
There were 9 experiments performed in Copenhagen, in which all of the required parameters were measured.The all parameters of the experiments that were used for calculations are shown in Table 1.

The Experiments
Experiments were done for the all cases of the performed Copenhagen experiments.The space was discretized to 100 × 100 points with the 40m spacing.The wind velocity field was pre-calculated using the equation 4 and the data in Table 1.The initial conditions for the test cases were defined by the initial concentration profile at zero time.
Two cases were tested.The first one has the shape of cone and the second one has the initial shape in form of cylinder.The cylinder case simulates the very sharp edges of concentrations which is very critical for the numerical schemes.
Table 1 The Parameters of the Performed Experiments in Copenhagen On the other hand the cone shape has the single maximum value of concentration and the numerical schemes have often tendency to smooth it.The diffusion in these test was set to zero, therefore the concentration profile should remain the same at the end of the experiments.
Fig. 1 The concentration profiles at the beginning of the simulation (red) and at the end of the simulation (blue) in case of no artificial diffusion addition

Oscillations
For the very steep concentration profiles like the cylinder, the numerical oscillations of the solution in case of EL-LAM scheme appear.The situation is shown in Figure 1.
To avoid the oscillations the artificial diffusion was added to the model.As a result, the oscillations disappeared afterwards.The new situation is shown in Figure 2.
Fig. 2 The concentration profiles at the beginning of the simulation (red) and at the end of the simulation (blue) in case of artificial diffusion addition

Error Measurement and Results
Each of the nine experiments was evaluated by several error measurements.The first one is the peak error (Equation 6).It demonstrates how the scheme is able to preserve the level of the concentration with respect to certain profile.The oscillations have the very negative effect to this measure.
peak c resp.peak 0 is the peak and min c resp.min 0 is the minimum of calculated resp.precise concentration profile.
The second measure is the distribution error (Equation 7).It describes the difference between the concentration distribution at the beginning and at the end of the simulation.The position of the profile has no influence to the error size.
Ω i and Ω j refers to domain where calc i, j and exact i, j differs from min 0 .
The last error measurement is referred to the mass conservation principle.In other words the sum of the concentration in the domain at the beginning of the simulation should be the same as at the end of the simulation.For the completeness of the evaluation we measured the calculation times of the schemes as well.
The results of the first set of tests are summarized in Table 2. ELLAM scheme gained the smaller peak error than the Walcek scheme, the similar distribution and final mass error.The calculation times was nearly the same, however one should be noted that the used ELLAM scheme that is described in section 2.2 includes the calculation of the diffusion term which was set to zero for the purposes of the error evaluation.The second test includes the cylindrical profile of the concentration.As was mentioned in section 3.4, the undesired oscillation due to the very sharp edges can occur during the simulations.Therefore, we performed the tests for the two versions of the ELLAM scheme.The first was with zero diffusion and the second one was with the artificial diffusion added to the model.The results are shown in Table 3, where the scheme with diffusion is marked as ELLAM D .It can be seen that in case of the zero diffusion ELLAM scheme gained relatively big peak error but very low distribution error and vice versa.It can be seen that the trade-off between the two cases should be appropriate.For further work on this problem we can be inspired by the [10] where the authors developed the technique to prevent the oscillation for 1D case.The calculation times are nearly the same as in previous test case.

CONCLUSION
We used the specific ELLAM scheme to solve the advection-diffusion equation with the parameters of atmospheric wind velocity fields.We performed the experiments and compared the scheme with the used advection integration scheme proposed in [8].From the performed experiments it is obvious that the ELLAM scheme is more than equal competitor to the state of the art method.It comes from the fact that the errors as well as the calculation time were similar and the ELLAM framework contained the diffusion computation as well.
For the further development we plan firstly to try to avoid oscillation in the scheme together with preserving the low distribution error.The next step would be the testing of the ELLAM scheme with the complete advection-diffusion atmospheric model and its evaluation against the real collecting data measurement.From the effectiveness point of view we are planning to rewrite the ELLAM code to the parallel platform using the very popular OpenCL framework [5].and computer graphics.For more information see please http://www.fit.vutbr.cz/idvorak.
František Zbořil is an Associate Professor of Computer Science at the Faculty of Information Technology, Brno University of Technology, Czech Republic.He received his M.Sc. in 1968 and Ph.D. in 1978 (both in Computer Science) at the same university.He started his research activities on analogue and hybrid computers with simulation of continuous systems, namely of systems described by partial differential equations.His next research was focused on classical artificial intelligence, robotics and neural networks.Now, the main objects of his professional interests are the simulation of combined dynamic systems and soft computing problems.He is the author of more than 100 papers and several lecture notes and he has supervised about 200 bachelors, masters and doctoral theses.He is a member of the board of the Czech and Slovak Simulation Society (CSSS) and a member of several other educational, research and academic boards or societies.

Table 2
The results of the performed experiments in case of cone profile

Table 3
The results of the performed experiments in case of cylinder profile