﻿ Entropic Approach for Emission Rate Estimation of Area Pollutant Sources

American Journal of Environmental Engineering

p-ISSN: 2166-4633    e-ISSN: 2166-465X

2013;  3(1): 56-62

doi:10.5923/j.ajee.20130301.08

### Entropic Approach for Emission Rate Estimation of Area Pollutant Sources

Débora R. Roberti 1, Domenico Anfossi 2, Haroldo F. de Campos Velho 3, Gervásio A. Degrazia 1

1Department of Physics, Federal University of Santa Maria (UFSM), Santa Maria, 97105-900, Brazil

2Institute of Atmospheric Sciences and Climate (ISAC), Italian National Research Council (CNR), Turim, 10133, Italy

3Laboratory for Computing and Applied Mathematics (LAC), National Institute for Space Research (INPE), 12227-010, São José dos Campos, Brazil

Correspondence to: Débora R. Roberti , Department of Physics, Federal University of Santa Maria (UFSM), Santa Maria, 97105-900, Brazil.
 Email:

Abstract

The estimation of the area source pollutant strength is a relevant issue for atmospheric environment. This characterizes an inverse problem in the atmospheric pollution dispersion. In the inverse analysis, an area source domain is considered, where the strength of such area source term is assumed unknown. The inverse problem is formulated as a non-linear optimization approach, whose objective function is given by the square difference between the measured pollutant concentration and the mathematical models, associated with a regularization operator. The forward problem is addressed by a source-receptor scheme, where a regressive Lagrangian model is applied to compute the transition matrix. A quasi-Newton method is employed for minimizing the objective function. The second order maximum entropy regularization is used, and the regularization parameter is calculated by the L-curve technique. This inverse problem methodology is tested with synthetic observational data, producing good inverse solutions.

Keywords: Atmospheric Pollutant, Source Estimation, Regularized Inverse Solution, Entropic Regularization

Cite this paper: Débora R. Roberti , Domenico Anfossi , Haroldo F. de Campos Velho , Gervásio A. Degrazia , Entropic Approach for Emission Rate Estimation of Area Pollutant Sources, American Journal of Environmental Engineering, Vol. 3 No. 1, 2013, pp. 56-62. doi: 10.5923/j.ajee.20130301.08.

### 1. Introduction

A general mathematical theory for dealing with inverse problems is due to the Russian mathematician Andrei Nikolaevich Tikhonov at 1963, introducing the regularization method as a general procedure to solve such mathematical problems. The regularization method looks for the smoothest (regular) inverse solution, where the data model would have the best fitting related to the observation data, subject to constrains. The searching for the smoothest solution is an additional information (or “a priori information), which reads the ill-posed inverse problem a well-posed one[1]. In the mathematical implementation of the method, the inverse problem is formulated as an optimization problem with constrains (a priori information). These constrains can be added to the objective function with the help of a regularization parameter:
 (1)
being u the unknowns, A(.) the mathematical model, the observation with δ noise level, Ω(.) the regularization operator, and α the regularization parameter.
This paper is focused on the application of inverse prolem methodology for solving a problem with strong interest in meteorology and environmental science: estimation of the emission rate from the area pollutant sources. The identification of the pollution source is calculated here obtaining a regularized inverse solution.
The problem for identifying the minority gas emission rate for the system ground-atmosphere is an important issue for the bio-geochemical cycle, and it has being intensively investigated. This inverse problem has been solved using regularized solutions[2], Bayes estimation[3, 4], and variational methods[5] – the latter approach coming from the data assimilation studies.
In principle, we have no information about how the area pollution source emission rate is spatially distributed. Seibert[6] and other researchers have suggested a domain decomposition, where each sub-domain has different pollutant emission strength. The direct problem can be expressed as forward-time evolution problem, or even as in a backward-time pollutant diffusion problem. The choice between these approaches (forward-time or backward-time) is given by the proposed scheme that has lower computational effort. The decision is based on the number of sources and number of the measurement points. If the number of receptors is less than the number of sources, the backward-time approach should be adopted; otherwise (if the number of receptors is greater than number of sources) the forward-time approach is suggested.
In this study, the direct mathematical model is described by a Lagrangian particle stochastic code. However, here the direct model is represented by the source-receptor scheme. In such approach, the measured concentration in a given receptor is calculated from the contribution of all sources by means of a transition matrix. The entries of the transition matrix are determined employing the Lagrangian code. Therefore, the direct model is represented by a matrix-vector product, instead of use the very expensive Monte Carlo Lagragian scheme for each iteration of the optimization process. This procedure reduces the computational effort (complexity) of the inverse analysis.
The inverse problem is formulated as a constrained nonlinear optimization problem and solved by a quasi-Newtonian optimization algorithm. The objective function is the norm L2 (see equation (1)) of the differences between the measured concentration and the concentration computed from the mathematical model (Lagrangian code), associated to the regularization operator. The regularization operator for the objective function (1) is the entropy measure of the vector of second-differences of the unknown parameters. This constrains the class of solution candidates into a restricted set of values composed by global smooth sub-domains, similar to the second-order polynomial surface. This regularization approach has been successfully adopted to retrieve the vertical atmospheric profile from satellite data[7], and to identify the turbulent eddy diffusivity[8].
The use of maximum entropy principle to reconstructing sources of inert atmospheric tracers from measurements have been independently introduced by[9] – see also[10] and[11] – an application of the Bocquet’s scheme for reconstruction of the Chernobyl accident source term is described by Davoine and Bocquet[12], for both papers the forward model is given by an Eulerian approach.
In a previous study on higher entropy regularization, it was shown that Morosov’s discrepancy principle was really efficient for the problem of determining the regularization parameter[13]. Here, another method to select the regularization parameter is explored.
The next section is a brief presentation of the formulation of the direct problem. This is followed by a description of the proposed inversion method, and a discussion of the numerical examples. The method is tested over a two-dimensional domain, using synthetic data corrupted with Gaussian white noise.

### 2. LAMBDA: the Lagrangian Stochastic Model

The Lagrangian particle model LAMBDA was developed to study the transport process and pollutant diffusion, starting from the Brownian random walk modeling[14, 15]. In the LAMBDA code, full-uncoupled particle movements are assumed. Therefore, each particle trajectory can be described by the generalized three dimensional form of the Langevin equation for velocity[16]:
 (2a)
and
 (2b)
where , and x is the displacement vector, U is the mean wind velocity vector, u is the Lagrangian velocity vector, ai(x,u,t) is a deterministic term and bij(x,u,t)dWj(t) is a stochastic term, and the quantity dWj(t) is the incremental Wiener process.
The deterministic (drift) coefficient ai(x,u,t) is computed using a particular solution of the Fokker- Planck equation associated to the Langevin equation. The diffusion coefficient bij(x,u,t) is obtained from the Lagrangian structure function in the inertial subrange, (τK << Δt << τL), where τK is the Kolmogorov time scale and τL is the Lagrangian de-correlation time scale. These parameters can be obtained employing the Taylor statistical theory on turbulence[17].
Backward integration can also be applied. This is just to identify which particle arriving in a sensor-j is coming from a source-i.
The drift coefficient ai(x,u,t), for forward and backward integration is given by
 (3a)
with,
 (3b)
and
 (3c)
where cv=1 for forward integration and cv=–1 for backward integration, PE=P(x,u,t) is the non-conditional PDF of the Eulerian velocity fluctuations, and Bij = (1/2)bikbjk. Of course, for backward integration, the time considered is t’=–t, and velocity U’=–U being U the mean wind speed. The horizontal PDFs are considered Gaussians, and for the vertical coordinate the truncated Gram-Charlier type-C of third order is employed[15].
The diffusion coefficients, bij(x,u,t), for both forward and backward integration are given by
 (4)
where δij is the Kronecker delta, σi2 and τLi are velocity variance at each component and the Lagrangian time scale[17], respectively. With the coordinates and the mass of each particle, the concentration is computed – see equations (5) and (6).
The inverse problem here is to identify the source term S(t). As mentioned, a source-receptor approach is employed for reducing the computer time, instead of running the direct model (equation (2)) for each iteration. This approach displays an explicit relation between the pollutant concentration of the i-th receptor related the j-th sourcers:
 (5)
where the matrix Mij is the transition matrix, and each matrix entry given by
 (6)
where VR,i and VS,j are the volume for the i-th receptor and j-th source, respectively; NS,j and NR,i are the number of particle realised by the j-th source and i-th sensor, respectively; NR,i,j and Ns,i,j are the number of particle released by the j-th source and detected by the i-th receptor.

### 3. Inverse Model

In order to set up the inverse analysis, it is assumed that the concentration obtained with a mathematic model is given by CMod(r,Q), where Q=[Q1(t), Q2(t), … , Qn(t)]T is the vector emission rate, and Qn(t) represents the emission rate of the n-th source and CExp(r) are data from measurement concentrations. The solution of an inverse problem is a function Q hat minimizes the following objective function:
 (7)
with Ω begin a regularization function and α a positive parameter, i.e. a Lagrange multiplier. Each source emission lies in the interval: Qi ∈[li , ui]. The lower bound li and the upper bound ui are chosen to allow the inversion to lie within some physical limits. The least square difference between experimental data and the computed values is represented by
 (8)
The regularization operator can be expressed by an entropic scheme[7, 8, 13]
 (9a)
 (9b)
with
 (9c)
and Q(m) represents the m-th derivative or difference among the unknown function or parameter vector. The application of derivative or difference operators is associated with whether we are dealing with continuous function or a discrete one. The function S(m) attain its global maximum when all k(r) are the same, i.e., a uniform distribution, in contrast, the lowest entropy value kmin=0 is reached when the probability distribution from the k(r) is a Dirac delta. This scheme is based on Jaynes’ criterion of inference[18], called maximum entropy principle (MaxEn-0).
Many schemes have been proposed to find the value of the regularization parameter. This parameter provides a fine balance between the square difference and the regularization terms. Some of these techniques are: Morozov’s discrepancy principle[19, 20], the L-curve[21], and the generalized cross validation[20].

#### 3.1. Optimization Algorithm

The optimization problem is iteratively solved by the deterministic method: quasi-Newtonian optimizer routine E04UCF[22], from the NAG Fortran Library. This routine is designed to minimize an arbitrary smooth function subject to constraints (simple bounds, linear and nonlinear constraints), using a sequential programming method. For the n-th iteration, the calculation proceeds as follows:
1. Solve the direct problem for Sn and compute the objective function J(Sn).
2. Compute by finite differences the gradient ∇J(Sn).
3. Compute a positive-definite quasi-Newton approximation to the Hessian Hn:
 (10)
where
4. Compute the search direction dn as a solution of the following quadratic programming subproblem:
Minimize , subjected to:
5. Set, Sn+1 = Sn + βndn where the step length βn minimizes J(Sn + βndn).

### 4. Numerical Experiment

For testing the procedure to estimate the emission rate procedure, the area pollutant sources are considered as placed in a box volume, where the horizontal domain and vertical height are given by: (1500 m 1000 m) 1000 m. There are two embedded regions R1 and R2 into a computational domain, with the following horizontal domain (600 m 600 m) for each region, and 1 m of height, and they are realising contaminants with two different emission rates. Figure 1 shows the computational scenario in a two-dimensional projection (x,y): the six sensors are placed at 10 m height and spread horizontally with the coordinates presented in Table 1. The domain is divided into sub-domains with 200 m 200 m 1 m. The emission rates for each sub-domain (SAi) are as follows,
R1 = SA2 = SA3 = SA4 = SA7 = SA9 = 10 gm-3s-1 ,
R2 = SA12 = SA13 = SA14 = SA17 = SA18 = SA19 = 20 gm-3s-1 .
For other sub-domians: S=0. The dispersion problem is modelled by the LAMBDA backward-time code, where data from the Copenhagen experiment are used, for the period (12:33 h – 12:53 h, 19/October/1978)[23, 24] It is assumed a logarithmic vertical profile for the wind field
 (11)
being U(z) the main stream, u* the fiction velocity, κ the von Kármán constant, z the height above the ground, and z0=0.06 m the roughness. The wind speed was measured at 10 m, 120 m, and 200 m. The numerical value for u* was obtained from the best fitting with the measured wind speed – see Table 2 – and the equation (11).
 Table 1. Position of sensors in physical domain
For this experiment the wind direction had an angle θ=180°, and the boundary layer height is h=1120 m[24]. The turbulence parameterization follows Degrazia et al.[17], for computing the wind variance (σi2) and the Lagrangian decorrelation time scale ([TL]i). These two turbulent parameters are considered constant for the whole boundary layer, and their numerical values were calculated at a z=10 m level. This characterises a stationary and vertically homogeneous atmosphere.
The sensor’s dimension, where the fictitious particles have arrived, was de 0.1 m 0.1 m 0.1 m, centred in the computational cells (presented in Table 1). Next, the reverse trajectories are calculated for 1000 particles per sensor.
The parameters for numerical simulations were 1800 time-steps with Δt = 1 s, meaning 30 min for the whole simulation. After 10 min, the concentration was computed for each 2 min, for the remaining 20 min of simulation. The mean concentration is found using the following expression:
 (12)
where (xj,yj) is the position of each sensor (the number of sensors is 6), NPES,j=1000 is the number of particles arriving at each sensor. NPVS,I,j,n is the quantity of fictitious particles in the i-th volume source, at the n-th instant.
 Table 2. Measured wind speed for the Copenhagen experiment[24]
 z (m) U (ms-1) 10 2.6 120 5.7 200 5.7
The inverse approach is tested using synthetic measured concentration data They are emulated as:
 (13)
where μ is a random number associated to the Gaussian distribution with zero mean and unitary standard deviation, and η is the level of the noise (Gaussian white noise). In our tests, η=0.05 (Experiment-1) and η=0.1 (experiment-2) were used.
The optimization method is used to find the minimum for the functional:
 (14)
where the vector S=[SA1, SA2,…, SA25 ]T represents the emission rate from 25 sub-domains, presented in Figure 1. However, some sub-domains are considered with null emission rate. Therefore, the number of the unknown parameters to be estimated by the inverse method decreases to 12: S=[SA2, SA4, SA7, SA8, SA9, SA12, SA13, SA14, SA17 , SA18, SA19]T. Another a priori information is smoothness of the inverse solution (regularization operator). As mentioned before in Section 1, the second order maximum entropy is applied (equation (9c), m=2).
Two numerical experiments are performed, using two different noise levels in the synthetic observational data: 5% and 10% of noise.
 Figure 1. (a) Sub-domains representation, in a 2D projection (x,y), with different emission rates, and black points (●) represent the sensor position at 10m height. (b) The gray-scale represents the emission rate
There are several methods to compute the regularization parameter (Section 3). Here, a numerical experimentation was employed. The value found was α=10-6 for Experiment-1, and α=10-5 for Experiment-2[9].
Results for Experiment-2 are presented in Table 4 and Figure 3. The error between the exact values and estimated by the inversions calculated by following form:
 (15)
 Table 3. The values of emission rate estimation, using the Q-N method, with the data from Experiment-1
In the Experiment-1, a global error around 2% is found from to inverse solution – see Table 3. For the region-1, the error in the estimation is 10%, and for region-2 the error is 3%. The biggest error estimation are verified for sub-domains A9 and A14, with error of 55% and 42%, respectively.
Table 4 shows the estimated emission rate for each sub-domain. The global estimation error enhances to 6.7%. The error for region-1 from the inverse analysis decreases to 1%, but the error for region-2 increases to 11%. The worse estimation is related to region-2, where the sub-domains A12 and A14 present approximately 6% of error.
 Table 4. The values of emission rate estimation, using the Q-N method, with the data from Experiment-2
The error for the inverse solution for experiment-2 is greater than experiment-1, as expected, since the level of noise in the measured data is higher in experiment-2. The results suggest that the inversion approach is effective, and the entropy regularization of second order is also appropriated.
One feature to point out is the use of source-receptor strategy. Indeed, this scheme allows a reduction of the computational complexity to compute the forward problem. As the inverse analysis is formulated as a non-linear optimization problem, this implies an iterative process requiring to execute several (thousands) times the forward problem.
 Figure 2. Representation of the of emission rate estimation in the domain with data from Experiment-1
 Figure 3. Representation of the of emission rate estimation in the domain with data from Experiment-2

### 5. Conclusions

The methodology proposed for estimating the emission rate term and localization of the source/sink was effective to produce good estimations for pollutant emission rate with second order maximum entropy. Good results were obtained even for a high level of noise, showing that the inverse model is robust related to the noise in experimental data. The determination of the regularization parameter by L-curve allowed finding an appropriated value for regularization operator, following the procedure presented by Hansen[21]. The source-receptor scheme is employed in the inverse analysis, in order to reduce the computational effort.
Future works will include the application of the methodology to real gas fluxes with different types of soil covering. Further we will employ other entropic regularization operators[25], and different optimization schemes for the minimizing functional (1), using stochastic schemes and/or approaches from artificial intelligence (neural networks, or fuzzy systems).

### ACKNOWLEDGEMENTS

The authors would like to thank to the Brazilians (CNPq, FAPESP) and Italian agencies for research financial support. The first author acknowledges CNPq-Brasil (Proc. 20.0046/04-7) for the sandwich doctorate fellowship in the ISAC-CNR, Turin – Italy.

### References

 [1] A.N. Tikhonov and V.I. Arsenin, Solutions of Ill-posed Problems, John Wiley & Sons, 1977. [2] P. Kasibhatla et al. (Editors), Inverse Methods in Global Biogeochemical Cycles, American Geophysical Union, USA, 2000. [3] I.G. Enting, Inverse Problems in Atmospheric Constituent Transport, Cambridge: University Press, UK. 2002. [4] N. R. Gimson, M. Uliasz, “The Determination of Agricultural Methane Emissions in New Zealand Using Inverse Modeling Techniques”, Atmospheric Environment, vol. 37, no. 28, pp. 3903-3912, 2003. [5] H. Elbern, A. Strunk, H. Schmidt, O. Talagrand, “Emission rate and chemical state estimation by 4-dimensional variational inversion”, Atmospheric Chemistry and Physics, vol. 7, no. 2, pp. 3749–3769, 2007. [6] P. Seibert, “Inverse Modelling of Sulfur Emissions in Europe Based on Trajectories”, in Inverse Methods in Global Biogeochemical Cycles, Geophysical Monograph Series, vol. 114, edited by P. Kasibhatla et al., pp. 147–154, AGU, Washington, D. C., 2000. [7] F. M. Ramos and A. Giovannini, “Résolution d'un problème inverse multidimensionnel de diffusion de la chaleur par la méthode des éléments analytiques et par le principe de l'entropie maximale. International Journal of Heat and Mass Transfer, vol. 38, no.1, pp. 101-111, 1995. [8] H.F. Campos Velho, M.R. de Moraes, F.M. Ramos, G.A. Degrazia, D. Anfossi, “An Automatic Methodology for Estimating Eddy Diffusivity from Experimental Data”, Nuovo Cimento, vol. 23-C, no. 1, pp. 65-84, 2000. [9] D.R. Roberti, Inverse Problems in Atmospheric Physics, Doctor of Science thesis, Physics Department, Federal University of Santa Maria, 2005. [10] Roberti, D.R., Anfossi, D., H.F. de Campos Velho, Degrazia, G.A., “Estimation of Location and Strength of the Pollutant Sources”, Ciencia e Natura, pp. 131-134, 2005. [11] Marc Bocquet, "Reconstruction of an atmospheric tracer source using the principle ofnmaximum entropy. I: Theory", Quaterly Journal of the Royal Meteorological Society, vol. 131, no. 610, pp. 2191–2208, 2005. [12] X. Davoine, M. Bocquet, "Inverse modelling-based reconstruction of the Chernobyl source term available for long-range transport", Atmospheric Chemistry and Physics, vol. 7, no. 6, pp. 1549-1564, 2007. [13] W. B. Muniz, F.M. Ramos, H.F. Campos Velho, “Entropy- and Tikhonov-based Regularization Techniques Applied to the Backwards Heat Equation”, Computers & Mathematics with Applications, 40(8/9), pp. 1071-1084, 2000. [14] E. Ferrero, D. Anfossi, G. Brusasca, G. Tinarelli, "Lagrangian Particle Model LAMBDA: Evaluation against Tracer Data". International Journal of Environment and Pollution, vol. 5, no. 4-5, pp. 360-374, 1955. [15] D. Anfossi, E. Ferrero, “Comparison among Empirical Probability Density Functions of the Vertical Velocity in the Surface Layer Based on Higher Order Correlations”, Boundary-Layer Meteorology, vol. 82, no. 2, pp. 193-218, 1997. [16] D. J. Thomson, "Criteria for the Selection of Stochastic Models of Particle Trajectories in Turbulent Flows", Journal of Fluid Mechanics, vol. 180, no. 1, pp. 529-556, 1987. [17] G. A. Degrazia, D. Anfossi, J. C. Carvalho, T. Tirabassi, H. F. Campos Velho, "Turbulence Parameterization for PBL Dispersion Models in All Stability Conditions", Atmospheric Environment, vol. 34, no. 21, pp. 3575–3583, 2000. [18] E. T. Jaynes, "Information Theory and Statistical Mechanics", Physical Review, vol. 106, no. 4, pp. 620-630, 1957. [19] V. A. Morosov, "On the Solution of Functional Equations by the Method of Regularization", Sovietic Mathematic Doklady, vol. 7, pp. 414-417, 1966. [20] M. Bertero, P. Boccacci, Introduction to Inverse Problems in Imaging, Taylor & Francis, USA, 1998. [21] Christian Hansen, "Analysis of Discrete Ill-posed Problems by means of the L-curve", SIAM Review, vol. 34, no. 4, pp. 561-580, 1992. [22] E04UCF routine: NAG Fortran Library. Oxford: Mark 17. 1995. [23] S.E. Gryning, and E. Lyck, “Atmospheric dispersion from elevated source in an urban area: comparison between tracer experiments and model calculations”, Journal of Climate Applied Meteorology, vol. 23, no. 4, pp.651-660, 1984. [24] S.E. Gryning, and E. Lyck, The Copenhagen tracer experiments: reporting of measurements. RIS National Laboratory, ISBN 87-550-2395-9. 54 p (version available in the internet:http://www.risoe.dk/rispubl/VEA/ris-r-1054.htm), 1998. [25] H.F. de Campos Velho, F.M. Ramos, E.H. Shiguemori, J.C. Carvalho (2006): A Unified Regularization Theory: The Maximum Non-extensive Entropy Principle, Computational and Applied Mathematics, 25(2-3), 1-24 – also available in the internet:www.scielo.br/scielo.php?script=sci_arttext&pid=S1807-03022006000200011