In conclusion, 6th – 8th February 2006 is a period which is representative for one of the highest O3 episodes during the dry season of the year in HCMC. The primary pollutants show highest values in centre of city where the highest density of traffic is found. Therefore, the huge part of population in the centre of HCMC is living with unfavorable conditions due to high concentration of primary pollutants. However, in the case of secondary pollutant we can see that it has most favorable conditions for the population living in the centre of HCMC and unfavorable for the population living in the north and north-west of HCMC. Once the models were shown able to reproduce and understand the principal characteristics of pollution in HCMC, it is very useful to study different strategies to reduce pollution for the city in the future.

#### **4.4 Air pollution abatement strategy**

For over 15 years, a lot of studies to evaluate air pollution abatement strategies have been carried out by using air quality models (Metcalfe et al. 2002; Palacios et al. 2002; Zarate et al. 2007 and to cite a few). The previous section shows that it is urgent to establish emission control scenarios for HCMC. Over some recently years, the results of air quality monitoring have shown that the pollutants exceeded regularly the standard limits in HCMC due to the emissions from the traffic source (HEPA, 2005; HEPA, 2006). The local government has started to design some emission control plans for traffic in HCMC. The plans are designed for the year of 2015 and 2020. The two emission control plans are named: (1) Emission reduction scenario for 2015 and (2) Emission reduction scenario for 2020.

The main ideas are that in 2015 the HCMC government will perform many activities to control air pollution concerning the road traffic source (Trinh, 2007): (1) controlling the emission of all vehicles (Thang, 2004), (2) the first metro line will be finished at the end of 2014 (Bao du lich, 2008) and (3) HCMC government will add 3000 new buses during 2006-2015 (Tuong, 2005).

For the year of 2020, four metro lines will be constructed (metro system will replace 50% of total motorcycle) (Bao du lich, 2008) and the number of buses will be increased to 4500 during 2006-2020 (Bang et al., 2010).

#### **Spatial distribution of O3**

The impacts of two strategies on the levels of troposphere O3 in HCMC are shown in Fig .16(a) and Fig .16(b). If HCMC follows the reduction plan: (i) For the year of 2015, the reductions for the HCMC grid area as a whole are 3.3% for CO, SO2, and CH4. There is an increase in NOx emissions of 8%. Mean values of O3 are reduced from 28.5 ppb to 28.0 ppb and the maximum from 150 ppb to 136 ppb on 7th Feb. (ii) For the year of 2020, the reductions for HCMC grid area as a whole are 8.6% for CO and 12.5 % for CH4, there are increases in NOx, SO2 and NMVOC emissions of 20.1%, 7.6% and 6.2%, respectively. Mean values of O3 are reduced from 28.5 ppb to 27.6 ppb and the maximum from 150 ppb to 120 ppb for 7th Feb.

The highest reduction of O3 concentration is found at the same place of the principal O3 plume for both abatement strategies. A deeper analysis of this reduction will be discussed by plotting the O3 concentration variable with respect to time and its uncertainties. In this research, we select some stations where we find the maximum of O3 reduction, the medium of O3 reduction and the minimum of O3 reduction. The MA, D2 and HB stations are chosen for the maximum, medium and minimum reduction zones, respectively (Fig .16 (a)). Their results are shown in the following section.

Urban Air Pollution 33

Fig. 17. Reduction of O3 (in ppb) in 2015 (strategy in 2015; Fig .17 (a, b, c)\_upper panels) and in 2020 (strategy in 2020; Fig .17 (d, e, f)\_lower panels) from the O3 concentration in 2006 and

its uncertainties. The results are plotted for the first layer near the ground during the

O3, ppb

O3, ppb

O3, ppb

Compare at D2 station

Delta O3

Standard deviation O3

6 12 18 24 6 12 18 24 6 12 18 24 FEB 6 FEB 7 FEB 8

Delta O3

Standard deviation O3

6 12 18 24 6 12 18 24 6 12 18 24 FEB 6 FEB 7 FEB 8

Delta O3

Standard deviation O3

6 12 18 24 6 12 18 24 6 12 18 24 FEB 6 FEB 7 FEB 8

Compare at HB station

Compare at M A station

Compare at M A station

Delta O3

Standard deviation O3

6 12 18 24 6 12 18 24 6 12 18 24 FEB 6 FEB 7 FEB 8

6 12 18 24 6 12 18 24 6 12 18 24 FEB 6 FEB 7 FEB 8

6 12 18 24 6 12 18 24 6 12 18 24 FEB 6 FEB 7 FEB 8

selected episode (on 6th - 8th Feb) at MA, D2 and HB stations.

Compare at D2 station

Delta O3

Standard deviation O3

Compare at HB station

Delta O3

Standard deviation O3

O3, ppb

O3, ppb

O3, ppb

Fig. 16. Effect of two strategies on O3 concentration (ppb) fields for the 7th Feb 2006 at 1300LT ground level. Figure 16(a) represents the reduction of O3 concentration in 2015 from the ozone concentration in 2006. Figure 16(b) represents the reduction of O3 concentration in 2020 from the ozone concentration in 2006. The measurement stations are shown in Fig. 13.

#### **Analysis of O3 at different measuring stations**

(b)

Strategy in 2015: Figure 17(a, b, c)-upper panels shows the reduction of O3 concentrations ( *O*3 in equation 4) in 2015 from the O3 concentration in 2006 and its uncertainties ( *O*3 in equation 5) in different stations (we consider the standard deviation as the uncertainty).

The Delta O3 ( *O*3 ) and uncertainties of O3 ( *O*3 ) presented on Fig .17 are calculated using equations (4) and (5). These values are calculated from 100 MC simulations for the base case and 100 MC simulations for each strategy.

**DO**

MA

**HB**

**TS**

**TD**

**TD**

in

) presented on Fig .17 are calculated using

**D2 ZO**

**D2 ZO**

**DO**

MA

**HB**

Fig. 16. Effect of two strategies on O3 concentration (ppb) fields for the 7th Feb 2006 at 1300LT ground level. Figure 16(a) represents the reduction of O3 concentration in 2015 from the ozone concentration in 2006. Figure 16(b) represents the reduction of O3 concentration in 2020 from the ozone concentration in 2006. The measurement stations are shown in Fig. 13.

Strategy in 2015: Figure 17(a, b, c)-upper panels shows the reduction of O3 concentrations ( *O*3 in equation 4) in 2015 from the O3 concentration in 2006 and its uncertainties ( *O*3

equations (4) and (5). These values are calculated from 100 MC simulations for the base case

equation 5) in different stations (we consider the standard deviation as the uncertainty).

**TS**

**QT**

**TN**

**QT**

**TN**

(a)

(b)

**Analysis of O3 at different measuring stations**

The Delta O3 ( *O*3 ) and uncertainties of O3 ( *O*3

and 100 MC simulations for each strategy.

**BC**

**BC**

Fig. 17. Reduction of O3 (in ppb) in 2015 (strategy in 2015; Fig .17 (a, b, c)\_upper panels) and in 2020 (strategy in 2020; Fig .17 (d, e, f)\_lower panels) from the O3 concentration in 2006 and its uncertainties. The results are plotted for the first layer near the ground during the selected episode (on 6th - 8th Feb) at MA, D2 and HB stations.

Urban Air Pollution 35

main source of pollution in cities. The environmental managers strongly need for air quality management systems to help them control air quality in the cities. The available air quality management system is discussed. Emission inventory is the most important part in air quality management system as it is the mostly used tool to identify the pollution source. However there are a lot of difficulties to be faced when we generate EIs. In this chapter, we also present an example of generation an EI for developing cities where there are not enough available data. The results of EI are used as input parameters for air quality modeling. Then we present many Air pollution abatement strategies over HCMC as a case study of developing cities.

Ermakov, X.M., 1997. Monte Carlo methodology and its relations. Translated by Pham,T.N.

Molina, L., Molina, M., 2002. Air Quality in the Mexico Megacity. An integrated assessment. Kluwer Academic Publishers, Dordrecht, The Netherlands. ISBN 1-4020-045204.

Moussiopoulos, Nicolas., 2003. Air Quality in Cities (book). Springer, Heidelberg, Germany.

Bang Q. HO, Clappier, A., 2011. Road traffic emission inventory for air quality modelling

Bang Q, HO, Clappier, A., Zarate, E., Hubert, V.D.B., Fuhrer, O., 2006. Air quality meso-

Bougeault, P., Lacarrere, P., 1989. Parameterization of orography-induced turbulence in a

Clappier, A., Perrochet, P., Martilli, A., Muller, F., Krueger, B.C., 1996. A new non-hydrostatic

Collella, P., Woodward, P., 1984. The piecewise parabolic method (PPM) for gas dynamical

Gurjar, B.R., J. Lelieveld., 2005. New Directions: Megacities and global change. Atmospheric

Generoso, S and Bey, I., 2007: A satellite and model-based assessment of the 2003 Russian

fires: Impact on the Arctic region. J. Geophys.Res.,vol 112, D15302

'96. Computational Mechanics Publications, Southampton, pp. 527–553. Clappier, A., 1998. A correction method for use in multidimensional time-splitting advection

mesobeta-scale model. Monthly Weather Review 117, 1872–1890.

simulations. Journal of Computational Physics 54, 174–201.

and to evaluate the abatement strategies: a case of Ho Chi Minh City, Vietnam. Atmospheric Environment Vol 45, Issue 21 (2011) pp. 3584-3593. ISSN: 1352-2310. Bang Q. HO, Clappier, A., Golay F., 2011. Air pollution forecast for Ho Chi Minh City,

Vietnam in 2015 and 2020. Air Quality, Atmosphere & Health, Volume 4, Number

scale modeling in Ho Chi Minh City: evaluation of some strategies' efficiency to reduce pollution. Journal Science and Technology Development. Vol 9, N° 5, 2006. Belalcazar, L., Fuhrer, O., Ho. D., Zarate, E., Clappier, A., 2009. Estimation of road traffic

emission factors from a long term tracer study in Ho Chi Minh City (Vietnam),

mesoscale model using a control volume finite element (CVFE) discretisation technique. In: Borrell, P.M., et al. (Ed.), Proceedings of the EUROTRAC Symposium

algorithms: application to two- and three-dimensional transport. Monthly Weather

Rainer, Friedrich., Stefan, Reis., 2004. Emissions of Air Pollutants, Springer. 333p

atmospheric environment, vol. 43 (2009) 5830–5837.

**6. References**  - Print Books

Dinh, X.T., 2003. Air pollution, VNU-HCM edition. 399p.

ISBN 3-540-00842-x. 298 p.

2, p.145-158. ISSN: 1873-9318.

Review 126, 232–242.

Environment 39 (2), 391-393.

Sciences and Technology publishers. 271p.

$$\text{Debta } \mathbf{O}\_3 = \overline{\boldsymbol{\Lambda O}}\_3 = \frac{\sum\_{i=1}^{100} \left( \mathbf{O}\_{3\\_Base\\_case}^i - \mathbf{O}\_{3\\_2015}^i \right)}{100} \tag{4}$$

where: i is the number of simulation

3 \_ *i O Base case* is the O3 concentration of base case (2006) for the ith simulation

3 \_ 2015 *i O* is the O3 concentration of strategy in 2015 for the ith simulation,

$$\text{Standard deviation of } \mathbf{O}\_3 = \sigma\_{\overline{\underline{\alpha O}}} = \sqrt{\frac{\sum\_{i=1}^{100} (\Delta O\_3^i - \overline{\Delta O}\_3)^2}{99}} \tag{5}$$

Where: 3 *i O* is the difference in O3 concentration between the base case and the strategy in 2015 for the ith simulation.

The highest uncertainties of O3 reduction appear at the same time of the highest O3 reduction at 1200LT - 1400LT of each day (Fig .17). The highest reduction of O3 concentration at MA, D2 and HB stations during 6th - 8th Feb are 14ppb, 6.7ppb and 3.9ppb, respectively. However, the highest uncertainties of O3 reduction at MA, D2 and HB stations are 3ppb, 9ppb and 3.5ppb respectively. Therefore, the uncertainties of O3 reduction are in general similar to the O3 reduction. We can not conclude that the change in O3 concentration is due to the impact of the emission control plan, because the change can probably be due to the impact of uncertainties of input parameters.

For the evolution of primary pollutants of the strategy in 2015, the concentrations of NOx in simulations will increase 7% than those were in 2006. However, the concentrations of CO and CH4 decrease around 10% and 8%, respectively than those were in 2006. In conclusion, there is very little impact of the emission control plan in 2015.

Strategy in 2020: Figure 17 (d, e, f)-lower panels shows that the highest uncertainties of O3 reduction ( *O*3 ) also appear at the same time of the highest O3 reduction ( *O*3 ) at 1200LT - 1400LT of each day. In general, the O3 concentration in 2020 is lower than the O3 concentration in 2006. The highest reduction of O3 concentration at MA, D2 and HB stations during 6th – 8th Feb are 23.5ppb, 13.4ppb and 7.8ppb, respectively. While, the highest uncertainties of O3 reduction at MA, D2 and HB stations are 8.2ppb, 11ppb and 5.1ppb, respectively. The O3 reduction is higher than the uncertainty of O3 reduction in all stations. It means that the change in O3 concentration is due to the change in emissions from the emission control plan.

For the evolution of primary pollutants of the strategy in 2020, the concentrations of NOx in simulations will increase 16% than those were in 2006. However, the concentrations of CO decrease around 7% than those were in 2006. In conclusion, this emission control plan in 2020 has strong impact on the primary and secondary pollutants in HCMC.

#### **5. Summary and key recommendations**

In this chapter, the state of the art for studying urban air pollution was reviewed and discussed. Air in the large cities is polluted, especially in developing cities. Road traffic is the

3 \_ 3 \_ 2015

100

1

*i*

(4)

(5)

2

3 3

( )

*O O*

*i*

99

( )

100

*i i Base case*

*O O*

3

*O*

*i O* is the difference in O3 concentration between the base case and the strategy in

100

3

The highest uncertainties of O3 reduction appear at the same time of the highest O3 reduction at 1200LT - 1400LT of each day (Fig .17). The highest reduction of O3 concentration at MA, D2 and HB stations during 6th - 8th Feb are 14ppb, 6.7ppb and 3.9ppb, respectively. However, the highest uncertainties of O3 reduction at MA, D2 and HB stations are 3ppb, 9ppb and 3.5ppb respectively. Therefore, the uncertainties of O3 reduction are in general similar to the O3 reduction. We can not conclude that the change in O3 concentration is due to the impact of the emission control plan, because the change can probably be due to

For the evolution of primary pollutants of the strategy in 2015, the concentrations of NOx in simulations will increase 7% than those were in 2006. However, the concentrations of CO and CH4 decrease around 10% and 8%, respectively than those were in 2006. In conclusion,

Strategy in 2020: Figure 17 (d, e, f)-lower panels shows that the highest uncertainties of O3

1400LT of each day. In general, the O3 concentration in 2020 is lower than the O3 concentration in 2006. The highest reduction of O3 concentration at MA, D2 and HB stations during 6th – 8th Feb are 23.5ppb, 13.4ppb and 7.8ppb, respectively. While, the highest uncertainties of O3 reduction at MA, D2 and HB stations are 8.2ppb, 11ppb and 5.1ppb, respectively. The O3 reduction is higher than the uncertainty of O3 reduction in all stations. It means that the change

For the evolution of primary pollutants of the strategy in 2020, the concentrations of NOx in simulations will increase 16% than those were in 2006. However, the concentrations of CO decrease around 7% than those were in 2006. In conclusion, this emission control plan in

In this chapter, the state of the art for studying urban air pollution was reviewed and discussed. Air in the large cities is polluted, especially in developing cities. Road traffic is the

in O3 concentration is due to the change in emissions from the emission control plan.

2020 has strong impact on the primary and secondary pollutants in HCMC.

) also appear at the same time of the highest O3 reduction ( *O*3 ) at 1200LT -

1

*i*

*O*

3 3

*i O Base case* is the O3 concentration of base case (2006) for the ith simulation

*i O* is the O3 concentration of strategy in 2015 for the ith simulation,

Standard deviation of O

the impact of uncertainties of input parameters.

**5. Summary and key recommendations**

there is very little impact of the emission control plan in 2015.

Delta O

where: i is the number of simulation

3 \_

3 \_ 2015

Where: 3

reduction ( *O*3

2015 for the ith simulation.

main source of pollution in cities. The environmental managers strongly need for air quality management systems to help them control air quality in the cities. The available air quality management system is discussed. Emission inventory is the most important part in air quality management system as it is the mostly used tool to identify the pollution source. However there are a lot of difficulties to be faced when we generate EIs. In this chapter, we also present an example of generation an EI for developing cities where there are not enough available data. The results of EI are used as input parameters for air quality modeling. Then we present many Air pollution abatement strategies over HCMC as a case study of developing cities.

#### **6. References**

Dinh, X.T., 2003. Air pollution, VNU-HCM edition. 399p.

Rainer, Friedrich., Stefan, Reis., 2004. Emissions of Air Pollutants, Springer. 333p

Urban Air Pollution 37

Streets.D.G, Bond.T.C, Lee.T and Jang.C., 2004: On the future of carbonaceous aerosol

Varinou, M., Kallos. G, Tsiligiridis, G and Sistla., 1999. The role of anthropogenic and

Stockwell, W.R., Kirchner, F., Kuhn, M., Seefeld, S., 1997. A new mechanism for regional

Zarate, E., Belalcazar, L.C., Clappier, A., Manzi, V and Hubert V, D, B. 2007. Air quality

ADB (Asian Development Bank) and Clean Air Initiative for Asian Cities Center. Country Synthesis Report on Urban Air quality Management: Vietnam. Dec., 2006. Belalcazar, Luis.C. Alternative Techniques to Assess Road Traffic Emissions. Ph.D Thesis of

Clappier, Alain. Modélisation numérique des polluants atmosphériques. 2001. 98p. Cours de

DOSTE (Department of Science, Technology and Environment of Ho Chi Minh city). Urban

European Environment Agency (EEA). 1999. EMEP/CORINAIR. Emission inventory

HEPA (Ho Chi Minh environmental protection agency). Report 2006 on air quality in Ho

HEPA (Ho Chi Minh environmental protection agency)., 2005. Last report of 2005 on

IER., 2006. Report annual. Environmental monitoring in South of Vietnam, Zone III. Air

Junier, Martin., 2004. Gas phase chemistry mechanisms for air quality modeling: generation and application to case studies. (2004). 112 p. Thèse Doctorat EPF Lausanne, no 2936 Krpo, A., 2009. Development and Application of a Numerical Simulation System to Evaluate

Kuentz Burchi, C., 1996. Les polluants atmosphériques. Approche toxicologique de

Kreinovich, J., Beck, C., Ferregut, A., Sanchez, G. R., Keller, M., Averill., Starks, S.A. Monte-

Mics-Asia project. International Institute for Applied Systems Analysis, Laxenburg, Austria.

quality monitoring program in south of Vietnam, Institute of Environment and

the Impact of Anthropogenic Heat Fluxes on Urban Boundary Layer Climate. Ph.D

l'évaluation des risques. Thèse de Doctorat. Université Louis Pasteur de Strasbourg

Carlo-Type Techniques for Processing Interval Uncertainty, and Their Potential Engineering Applications College of Engineering and NASA Pan-American Center for Earth and Environmental Studies (PACES), University of Texas, El Paso, TX

inventory of emissions sources for HCMC. December 2005.

DOSTE. Project on Energy – air pollution in HCMC. 2000. DOSTE - HCMC and ADEME. Duong. T,M,H. 2004. Proposed technical solutions for diminishing air pollution level by

transport energy demand and emission analysis – Case study of HCM city. N° 1

traffic and for analyzing the solution to overcome barriers in application. MSc/IER

emission inventories. . Atmospheric Environment, 41, 6302–6318.

biogenic emissions on tropospheric ozone formation over Greece. Phys. Chem.

atmosphere rich chemistry modeling. Journal of Geophysical Research 102, 25847–

modeling over Bogota, Colombia: Combined techniques to estimate and evaluate

emission. J. Geophys.Res.,vol 109,D24212

Earth (C), Vol. 24, No. 5, pp. 507-513, (1999).

25879.

EPFL. N° 4504 (2009).

(phase II). 2001.

guidebook.

thesis. EPFL.

79968, USA .

2000.

troisième année EPF Lausanne.

Chi Minh City. December 2006.

Resources (IER), (2006).

Gong, W., Cho, H., 1993. A numerical scheme for the integration of the gas phase chemical

Hanna, S.R., Chang, J.C., Fernau, M.E., 1998. Monte Carlo estimates of uncertainties in

Hanna, S.R., Lu, Z., Frey, C.H., Wheeler, N., Vukovich, J., Arunachalam, S., Fernau,

Hwang, D., Karami, H.A., Byun, D.W., 1998. Uncertainty analysis of environmental models within GIS environments. Computers & Geosciences 24 (2), 119–130. Ho M.D., Dinh X.T. Estimation of emission factors of air pollutants from the road traffic in

HOUTRANS, 2004. The Study on Urban Transport Master Plan and Feasibility Study in Ho

Junier, M., Kirchner, F., Clappier A. and Hubert, V.D.B., 2004. The chemical mechanism

in a three-dimensional mesoscale simulation, Atmos. Environ. 39, 1161-1171 Kuhlwein, J and Friedrich, R., 2000. Uncertainties of modeling emissions from road

Louis, J.F., 1979. A parametric model of vertical eddies fluxes in the atmosphere. Boundary-

Martilli, A., Clappier, A., Rotach, M.W., 2002. An urban surface exchange parameterization

Martilli, A., Roulet, Y.-A., Junier, M., Kirchner, F., Rotach, M.W and Clappier, A., 2003. On

Metcalfe, S.E., Whyatt, J.D., Derwent, R.G., O'Donoghue, M., 2002. The regionaldistribution

Nenes, A., Pandis, S., Pilinis, C., 1998. ISORROPIA: A new thermodynamic equilibrium

Palacios, M., Kirchner, F., Martilli, A., Clappier, A., Martin, F., Rodriguez, M.E., 2002. Summer

Abatement Strategies by Modeling. Atmospheric Environment 36, 5323–5333. Rappenglück, B., Oyola, P., Olaeta, I., Fabian, P., 2000. The evolution of Pho-tochemical

Ranjeet S. Sokhi and Nutthida Kitwiroon, 2008. World Atlas of atmospheric pollution.

Streets.D.G, Bond.T.C, Carmichael. G.R, Fernandes.S.D, Fu.Q, He.D, Klimont.D, Nelson.S.M,

the impact of urban exchange parameterization on air quality simulations: the

of ozone across the British Isles and its response to control strategies. Atmospheric

model for multiphase multi-component inorganic aerosols. Aquatic Geochemistry

ozone episodes in the Greater Madrid area. Analyzing the Ozone Response to

Smog in the Metropolitan Area of Santiago de Chile. Journal of Applied

Tsai.N.Y, Wang.M.Q, Woo.J.H and Yarber.K.F., 2003: An inventory of gaseous and primary aerosol emissions in Asia in the year 2000. J. Geophys. Res. 108 (D21)

for mesoscale models. Boundary-Layer Meteorology 104, 261–304.

Environment. 27A, 2147–2160.

ALMEC Corporation

Layer Meteorology 17, 187–202

Environment 36, 4045–4055.

Meteorology 39, 275-290.

Anthem Press, first edition. 144p.

(2003), p. 8809 doi:10.1029/2002JD003093.

4, 123-152.

variables. Atmos Environ 1998; 32:3619-3628.

transport, Atmos. Environ., (2000) 34, 4603-4610.

Athens case, Atmos.Environ.37, 4217-4231.

OTAG domain. Atmospheric Environment 35 (5), 891–903.

HCMC. VNU Journal of Science, Earth Sciences 24 (2008) 184-192

rate equations in three-dimensional atmospheric models. Atmospheric

predictions by a photochemical grid model (UAM-IV) due to uncertainties in input

M.,Hansen, A., 2000. Uncertainties in predicted ozone concentration due to input uncertainties for the UAM-V photochemical grid model applied to the July 1995

Chi Minh Metropolitan Area. In: Vol. 6, No. 1 Traffic and Transport Surveys.

generation program CHEMATA, part II: Comparison of four chemical mechanisms

**0**

**2**

1,2*Brazil* 3*Italy*

**Analytical Model for Air Pollution in the**

Daniela Buske1, Marco Tullio Vilhena2, Bardo Bodmann2

The worldwide concern due to the increasing frequency of ecological disasters on our planet urged the scientific community to take action. One of the possible measures are analytical descriptions of pollution related phenomena, or simulations by effective and operational models with quantitative predictive power. The atmosphere is considered the principal vehicle by which pollutant materials are dispersed in the environment, that are either released from the productive and private sector or eventually in accidental events and thus may result in contamination of plants, animals and humans. Therefore, the evaluation of airborne material transport in the Atmospheric Boundary Layer (ABL) is one of the requirements for maintenance, protection and recovery of the ecological system. In order to analyse the consequences of pollutant discharge, atmospheric dispersion models are of need, which have to be tuned using specific meteorological parameters and conditions for the considered region. Moreover, they shall be subject to the local orography and shall supply with realistic information concerning environmental consequences and further help reduce impact from potential accidents such as fire events and others. Moreover, case studies by model simulations may be used to establish limits for the escape of pollutant material from specific

The theoretical approach to the problem essentially assumes two basic forms. In the Eulerian approach, diffusion is considered, at a fixed point in space, proportional to the local gradient of the concentration of the diffused material. Consequently, it considers the motion of fluid within a spatially fixed system of reference. These kind of approaches are based on the resolution, on a fixed space-time grid, of the equation of mass conservation of the pollutant chemical species. Lagrangian models are the second approach and they differ from Eulerian ones in adopting a system of reference that follows atmospheric motions. Initially, the term Lagrangian was used only to refer to moving box models that followed the mean wind trajectory. Currently, this class includes all models that decompose the pollutant cloud into discrete "elements", such as segments, puffs or pseudo-particles (computer particles). In Lagrangian particle models pollutant dispersion is simulated through the motion of computer

**1. Introduction**

sites into the atmosphere.

**Atmospheric Boundary Layer**

2*Federal University of Rio Grande do Sul, Porto Alegre, RS* 3*Institute ISAC, National Research Council, Bologna*

1*Federal University of Pelotas, Pelotas, RS*

and Tiziano Tirabassi3

http://www.ineris.fr/aida/?q=consult\_doc/consultation/<IP_ADDRESS>.8.4339

	-

### **Analytical Model for Air Pollution in the Atmospheric Boundary Layer**

```
Daniela Buske1, Marco Tullio Vilhena2, Bardo Bodmann2
                                      and Tiziano Tirabassi3
                       1Federal University of Pelotas, Pelotas, RS
       2Federal University of Rio Grande do Sul, Porto Alegre, RS
             3Institute ISAC, National Research Council, Bologna
                                                       1,2Brazil
                                                          3Italy
```
#### **1. Introduction**

38 Air Pollution – Monitoring, Modelling and Health

Nguyen, T.T.M., 2000. Proposed methods to reduce air pollution by traffic in HCMC:

Nguyen, T, 2002: Asian Regional Research Program in Energy, Environment and Climate

Sathya, V., 2003 Uncertainty analysis in air quality modeling - the impact of meteorological

Statistical yearbook Ho Chi Minh City., 2003. Ho Chi Minh Statistical Office [réf. du Juin

Thang, Q. D., 2004. A Vision for Cleaner Emissions from Motorcycles in Viet Nam. Paper

Tremback, C.J., Kessler, R., 1985. A surface temperature and moisture parameterization for

World Health Organization (WHO), 2005. WHO Air Quality guidelines for particulate

Zarate, E., 2007. Understanding the Origins and Fate of Air Pollution in Bogotá, Colombia.

Le, T.G.; Dan, G and Nao, I., 2008. Clean Air Initiative, "Air Pollution Blamed as Study Finds

Nguyen, D.T., Pham, T.T. Air pollution in HoChiMinh City, Vietnam. Conference on: "Better Air quality in Asian and Pacific Rim Cities (BAQ 2002), Dec.2002, Hong Kong.

EU directive N° 96/62/CE of 27/09/96 concernant l'évaluation et la gestion de la qualité de

Madronich, S., 1998. TUV troposphere ultraviolet and visible radiation model, from the

Tuong, L., 2005. Go together by bus. Conference in: "The solutions for reducing the transport congestion for HCMC City in 2020", HCMC. 2005. Available at http://vietbao.vn/Phong-su/Di-xe-buyt-thoi-xang-tang-gia/30066850/263/.

 http://www.ineris.fr/aida/?q=consult\_doc/consultation/<IP_ADDRESS>.8.4339 Gale group., 2007. Ho Chi Minh City becomes one of 100 rapid economic growth cities.

March, 2007. Ipr strategic business information database – articles.

Website: http://acd.ucar.edu/models/open/tuv/tuv.html/.

Respiratory Illness Hitting HCMC's Children," March 26, 2008.

Numerical Weather Prediction, Montreal, Quebec, Canada, June 17–20. World Health Organization (WHO), 2000. Guidelines for Air Quality. World Health

input uncertainties. Thesis N°2318 (2003). EPFL

Vietnam", HCMC. August 2007.

Organization, Geneva.

l'air ambiant.

Health Organization, Geneva.

Doctoral thesis N° 3768, EPFL.

http://www.encyclopedia.com/doc/1G1-160479731.html

Statistical yearbook Ho Chi Minh City., 2006. Ho Chi Minh Statistical Office.

VNU.

2004].

developing public transportation and using lead free gasoline. (2000), MSc/IER-

(ARRPEEC). Project phase III, supported by SIDA organization, from 2002 to 2004.

presented at the Cleaner Vehicles and Fuels in Viet Nam Workshop, 13–14 May 2004, Hanoi, Viet Nam. Vietnamese Ministry of Transportation and US-EPA. Trinh, N.G., 2007. Motorcycles do not meet emissions standards should be upgraded or

replaced. Conference in: "Control emission from motorcycles in major cities of

use in mesoscale numerical models, Proceedings of Seventh conference on

matter, ozone, nitrogen dioxide and sulphur dioxide. Global update 2005. World

The worldwide concern due to the increasing frequency of ecological disasters on our planet urged the scientific community to take action. One of the possible measures are analytical descriptions of pollution related phenomena, or simulations by effective and operational models with quantitative predictive power. The atmosphere is considered the principal vehicle by which pollutant materials are dispersed in the environment, that are either released from the productive and private sector or eventually in accidental events and thus may result in contamination of plants, animals and humans. Therefore, the evaluation of airborne material transport in the Atmospheric Boundary Layer (ABL) is one of the requirements for maintenance, protection and recovery of the ecological system. In order to analyse the consequences of pollutant discharge, atmospheric dispersion models are of need, which have to be tuned using specific meteorological parameters and conditions for the considered region. Moreover, they shall be subject to the local orography and shall supply with realistic information concerning environmental consequences and further help reduce impact from potential accidents such as fire events and others. Moreover, case studies by model simulations may be used to establish limits for the escape of pollutant material from specific sites into the atmosphere.

The theoretical approach to the problem essentially assumes two basic forms. In the Eulerian approach, diffusion is considered, at a fixed point in space, proportional to the local gradient of the concentration of the diffused material. Consequently, it considers the motion of fluid within a spatially fixed system of reference. These kind of approaches are based on the resolution, on a fixed space-time grid, of the equation of mass conservation of the pollutant chemical species. Lagrangian models are the second approach and they differ from Eulerian ones in adopting a system of reference that follows atmospheric motions. Initially, the term Lagrangian was used only to refer to moving box models that followed the mean wind trajectory. Currently, this class includes all models that decompose the pollutant cloud into discrete "elements", such as segments, puffs or pseudo-particles (computer particles). In Lagrangian particle models pollutant dispersion is simulated through the motion of computer

equation in each layer is then solved by the Laplace Transform technique. For more details

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 41

A further mile-stone in the direction of the present approach is given in references Costa et al. (2006; 2011) by the Generalized Integral Advection Diffusion Multilayer Technique (GIADMT), that presents a general solution for the time-dependent three-dimensional ADE again assuming the stepwise approximation for the eddy diffusivity coefficient and wind profile and proceeding further in a similar way as the multi-layer approach (Moreira et al. (2006a)), however in a semi-analytical fashion. In this chapter we report on a generalisation of the afore mentioned models resulting in a general space-time (3 ⊕ 1) solution for this problem,

To this end, we start from the three-dimensional and time dependent ADE for the ABL as the underlying physical model. In the subsequent sections we show how the model is solved analytically using integral transform and a spectral theory based methods which then may provide short, intermediate and long term (normalized) concentrations and permit to assess the probability of occurrence of high contamination level case studies of accidental scenarios. The solution of the 3 ⊕ 1 space-time solution is compared to other dispersion process approaches and validated against experiments. Comparison between different approaches may help to pin down computational errors and finally allows for model validation. In this line we show with the present discussion, that our analytical approach does not only yield an acceptable solution for the time dependent three dimensional ADE, but further predicts tracer concentrations closer to observed values compared to other approaches from the literature. This quality is also manifest in better statistical coefficients including model validation.

Differently than in most of the approaches known from the literature, where solutions are typically valid for scenarios with strong restrictions with respect to their specific wind and vertical eddy diffusivity profiles, the discussed method is general in the sense that once the model parameter functions are determined or known from other sources, the hierarchical procedure to obtain a closed form solution is unique. More specifically, the general analytical solution for the advection-diffusion problem works for eddy diffusivity and wind profiles, that are arbitrary functions with a continuous dependence on the vertical and longitudinal

A few technical details are given in the following, that characterise the principal steps of the procedure. In a formal step without physical equivalence we expand the contaminant concentration in a series in terms of a set of orthogonal eigenfunctions. These eigenfunctions are the solution of a simpler but similar problem to the existing one and are detailed in the next section. Replacing this expansion in the time-dependent, three-dimensional ADE in Cartesian geometry, we project out orthogonal components of the equation, thus inflating the tree-dimensional advection-diffusion equation into an equation system. Note, that the projection integrates out one spatial dimension, so that the resulting problem reduces the problem effectively to 2 ⊕ 1 space-time dimensions. Such a problem was already solved by the Laplace transform technique as shown in Moreira et al. (2009a), Moreira et al. (2009b), Buske et al. (2010) and references therein. In the next sections we present the prescription for constructing the solution for the general problem, but for convenience and without restricting generality we resort to simplified models. Further, we present numerical simulations and

about this methodology see reference Moreira et al. (2006a).

given in closed form.

spatial variables, as well as time.

indicate future perspectives of this kind of modelling.

particles whose trajectories allow the calculation of the concentration field of the emitted substance. The underlying hypothesis is that the combination of the trajectories of such particles simulates the paths of the air particles situated, at the initial instant, in the same position. The motion of the particles can be reproduced both in a deterministic way and in a stochastic way.

In the present discussion we adopt an Eulerian approach, more specifically to the K model, where the flow of a given field is assumed to be proportional to the gradient of a mean variable. K-theory has its own limits, but its simplicity has led to a widespread use as the mathematical basis for simulating pollution dispersion. Most of of Eulerian models are based on the numerical resolution of the equation of mass conservation of the pollutant chemical species. Such models are most suited to confronting complex problems, for example, the dispersion of pollutants over complex terrain or the diffusion of non-inert pollutants.

In the recent literature progressive and continuous effort on analytical solution of the advection-diffusion equation (ADE) are reported. In fact analytical solutions of equations are of fundamental importance in understanding and describing the physical phenomena. Analytical solutions explicitly take into account all the parameters of a problem, so that their influence can be reliably investigated and it is easy to obtain the asymptotic behaviour of the solution, which is usually more tedious and time consuming when generated by numerical calculations. Moreover, like the Gaussian solution, which was the first solution of ADE with the wind and the K diffusion coefficients assumed constant in space and time, they open pathways for constructing more sophisticated operative analytical models. Gaussian models, so named because they are based on the Gaussian solution, are forced to represent real situations by means of empirical parameters, referred to as "sigmas". In general Gaussian models are fast, simple, do not require complex meteorological input, and describe the diffusive transport in an Eulerian framework exploring the use of the Eulerian nature of measurements. For these reasons they are widely employed for regulatory applications by environmental agencies all over the world although their well known intrinsic limits. For the further discussion we restrict our considerations to scenarios in micro-meteorology were the Eulerian approach has been validated against experimental data and has been found useful for pollution dispersion problems in the ABL.

A significant number of works regarding ADE analytical solution (mostly two-dimensional solutions) is available in the literature. By analytical we mean that starting from the advection-diffusion equation, which we consider adequate for the problem, the solution is derived without introducing approximations that might simplify the solution procedure. A selection of references considered relevant by the authors are Rounds (1955), Smith (1957), Scriven & Fisher (1975), Demuth (1978), van Ulden (1978), Nieuwstadt & Haan (1981), Tagliazucca et al. (1985), Tirabassi (1989), Tirabassi & Rizza (1994), Sharan et al. (1996), Lin & Hildemann (1997), Tirabassi (2003). It is noteworthy that solutions from citations above are valid for specialized problems. Only ground level sources or infinite height of the ABL or specific wind and vertical eddy diffusivity profiles are admitted. A more general approach that may be found in the literature is based on the Advection Diffusion Multilayer Method (ADMM), which solves the two-dimensional ADE with variable wind profile and eddy diffusivity coefficient (Moreira et al. (2006a)). The main idea here relies on the discretisation of the Atmospheric Boundary Layer (ABL) in a multilayer domain, assuming in each layer that the eddy diffusivity and wind profile take averaged values. The resulting advection-diffusion 2 Will-be-set-by-IN-TECH

particles whose trajectories allow the calculation of the concentration field of the emitted substance. The underlying hypothesis is that the combination of the trajectories of such particles simulates the paths of the air particles situated, at the initial instant, in the same position. The motion of the particles can be reproduced both in a deterministic way and in a

In the present discussion we adopt an Eulerian approach, more specifically to the K model, where the flow of a given field is assumed to be proportional to the gradient of a mean variable. K-theory has its own limits, but its simplicity has led to a widespread use as the mathematical basis for simulating pollution dispersion. Most of of Eulerian models are based on the numerical resolution of the equation of mass conservation of the pollutant chemical species. Such models are most suited to confronting complex problems, for example, the

In the recent literature progressive and continuous effort on analytical solution of the advection-diffusion equation (ADE) are reported. In fact analytical solutions of equations are of fundamental importance in understanding and describing the physical phenomena. Analytical solutions explicitly take into account all the parameters of a problem, so that their influence can be reliably investigated and it is easy to obtain the asymptotic behaviour of the solution, which is usually more tedious and time consuming when generated by numerical calculations. Moreover, like the Gaussian solution, which was the first solution of ADE with the wind and the K diffusion coefficients assumed constant in space and time, they open pathways for constructing more sophisticated operative analytical models. Gaussian models, so named because they are based on the Gaussian solution, are forced to represent real situations by means of empirical parameters, referred to as "sigmas". In general Gaussian models are fast, simple, do not require complex meteorological input, and describe the diffusive transport in an Eulerian framework exploring the use of the Eulerian nature of measurements. For these reasons they are widely employed for regulatory applications by environmental agencies all over the world although their well known intrinsic limits. For the further discussion we restrict our considerations to scenarios in micro-meteorology were the Eulerian approach has been validated against experimental data and has been found useful

A significant number of works regarding ADE analytical solution (mostly two-dimensional solutions) is available in the literature. By analytical we mean that starting from the advection-diffusion equation, which we consider adequate for the problem, the solution is derived without introducing approximations that might simplify the solution procedure. A selection of references considered relevant by the authors are Rounds (1955), Smith (1957), Scriven & Fisher (1975), Demuth (1978), van Ulden (1978), Nieuwstadt & Haan (1981), Tagliazucca et al. (1985), Tirabassi (1989), Tirabassi & Rizza (1994), Sharan et al. (1996), Lin & Hildemann (1997), Tirabassi (2003). It is noteworthy that solutions from citations above are valid for specialized problems. Only ground level sources or infinite height of the ABL or specific wind and vertical eddy diffusivity profiles are admitted. A more general approach that may be found in the literature is based on the Advection Diffusion Multilayer Method (ADMM), which solves the two-dimensional ADE with variable wind profile and eddy diffusivity coefficient (Moreira et al. (2006a)). The main idea here relies on the discretisation of the Atmospheric Boundary Layer (ABL) in a multilayer domain, assuming in each layer that the eddy diffusivity and wind profile take averaged values. The resulting advection-diffusion

dispersion of pollutants over complex terrain or the diffusion of non-inert pollutants.

for pollution dispersion problems in the ABL.

stochastic way.

equation in each layer is then solved by the Laplace Transform technique. For more details about this methodology see reference Moreira et al. (2006a).

A further mile-stone in the direction of the present approach is given in references Costa et al. (2006; 2011) by the Generalized Integral Advection Diffusion Multilayer Technique (GIADMT), that presents a general solution for the time-dependent three-dimensional ADE again assuming the stepwise approximation for the eddy diffusivity coefficient and wind profile and proceeding further in a similar way as the multi-layer approach (Moreira et al. (2006a)), however in a semi-analytical fashion. In this chapter we report on a generalisation of the afore mentioned models resulting in a general space-time (3 ⊕ 1) solution for this problem, given in closed form.

To this end, we start from the three-dimensional and time dependent ADE for the ABL as the underlying physical model. In the subsequent sections we show how the model is solved analytically using integral transform and a spectral theory based methods which then may provide short, intermediate and long term (normalized) concentrations and permit to assess the probability of occurrence of high contamination level case studies of accidental scenarios. The solution of the 3 ⊕ 1 space-time solution is compared to other dispersion process approaches and validated against experiments. Comparison between different approaches may help to pin down computational errors and finally allows for model validation. In this line we show with the present discussion, that our analytical approach does not only yield an acceptable solution for the time dependent three dimensional ADE, but further predicts tracer concentrations closer to observed values compared to other approaches from the literature. This quality is also manifest in better statistical coefficients including model validation.

Differently than in most of the approaches known from the literature, where solutions are typically valid for scenarios with strong restrictions with respect to their specific wind and vertical eddy diffusivity profiles, the discussed method is general in the sense that once the model parameter functions are determined or known from other sources, the hierarchical procedure to obtain a closed form solution is unique. More specifically, the general analytical solution for the advection-diffusion problem works for eddy diffusivity and wind profiles, that are arbitrary functions with a continuous dependence on the vertical and longitudinal spatial variables, as well as time.

A few technical details are given in the following, that characterise the principal steps of the procedure. In a formal step without physical equivalence we expand the contaminant concentration in a series in terms of a set of orthogonal eigenfunctions. These eigenfunctions are the solution of a simpler but similar problem to the existing one and are detailed in the next section. Replacing this expansion in the time-dependent, three-dimensional ADE in Cartesian geometry, we project out orthogonal components of the equation, thus inflating the tree-dimensional advection-diffusion equation into an equation system. Note, that the projection integrates out one spatial dimension, so that the resulting problem reduces the problem effectively to 2 ⊕ 1 space-time dimensions. Such a problem was already solved by the Laplace transform technique as shown in Moreira et al. (2009a), Moreira et al. (2009b), Buske et al. (2010) and references therein. In the next sections we present the prescription for constructing the solution for the general problem, but for convenience and without restricting generality we resort to simplified models. Further, we present numerical simulations and indicate future perspectives of this kind of modelling.

(0, *y*0, *HS*) is located at the boundary of the domain **r***S* ∈ *δ*Γ. Note, that in cases where the source is located in the domain, one still may divide the whole domain in sub-domains, where the source lies on the boundary of the sub-domains which can be solved for each sub-domain separately. Moreover, a set of different sources may be implemented as a superposition of independent problems. Since the source term location is on the boundary, in the domain this term is zero everywhere (*S*(**r**) ≡ 0 for **r** ∈ Γ\*δ*Γ), so that the source influence may be cast in form of a condition, where we assume that our coordinate system is oriented such that the *x*-axis is aligned with the mean wind direction. Since the flow crosses the plane perpendicular

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 43

height of the source (in units of *m*), *Lx* and *Ly* are the horizontal domain limits (in units of *m*)

In this section we first introduce the general formalism to solve a general problem and subsequently reduce the problem to a more specific one, that is solved and compared to

In order to solve the problem (3) we reduce the dimensionality by one and thus cast the problem into a form already solved in reference Moreira et al. (2009a). To this end we apply the integral transform technique in the *y* variable, and expand the pollutant concentration as

where **R** = (*R*1, *R*2,...)*T* and **Y** = (*Y*1,*Y*2,...)*T* is a vector in the space of orthogonal

∇*T***K** + (**K**∇)*T*

ˆ *∂T*

*y* **K** + (**K**ˆ

(∇2**R***T***Y**) +

For convenience we introduce some shorthand notations, ∇2 = (*∂x*, 0, *∂y*)*T* and ˆ

*dy***Y**[**F**] = *Ly*

0

here **F** is an arbitrary function and ∧ signifies the dyadic product operator, and making use of orthogonality renders equation (8) a matrix equation. The appearing integral terms are

eigenfunctions, given by *Ym*(*y*) = cos(*λmy*) with eigenvalues *λm* = *m π*

*∂y***Y** =

2 **K** + (**K**∇2)*T*

 *Ly* 0

*u*¯*c*¯(0, *y*, *z*, *t*) = *Qδ*(*y* − *y*0)*δ*(*z* − *HS*) , (6)

*c*¯(*x*, *y*, *z*, *t*) = **R***T*(*x*, *z*, *t*)**Y**(*y*), (7)

∇2**R***T***Y** + **R***T* ˆ

*∂y*)*T* (**R***T* ˆ

**F***T* ∧ **Y** *dy* (9)

*Ly* for *m* = 0, 1, 2, . . ..

*∂y***Y**  *∂y* =

*∂y***Y**) . (8)

*s* ), *h* the height of the ABL (in units of *m*), *HS* the

to the propagation (here the *y*-*z*-plane) the source condition reads

where *Q* is the emission rate (in units of *g*

(0, *∂y*, 0)*T*, so that equation (3) reads now,

Upon application of the integral operator

= ∇*T*

∇2**R***T***Y** + **R***T* ˆ

(*∂t***R***T*)**Y** + **U**¯

**3. A closed form solution**

experimental findings.

**3.1 The general procedure**

and *δ*(*x*) represents the Cartesian Dirac delta functional.

#### **2. The advection-diffusion equation**

The advection-diffusion equation of air pollution in the atmosphere is essentially a statement of conservation of the suspended material and it can be written as (Blackadar (1997))

$$\frac{D\overline{c}}{Dt} = \partial\_t \overline{c} + \mathbf{U}^T \nabla \overline{c} = -\nabla^T \overline{\mathbf{U}^\prime c^\prime} + \mathcal{S} \tag{1}$$

where *c*¯ denotes the mean concentration of a passive contaminant (in units of *g m*3 ), *D Dt* is the substantial derivative, *∂t* is a shorthand for the partial time derivative and **U**¯ = (*u*¯, *v*¯, *w*¯)*T* is the mean wind (in units of *m s* ) with Cartesian components in the directions *x*, *y* and *z*, respectively, ∇ = (*∂x*, *∂y*, *∂z*)*T* is the usual Nabla symbol and *S* is the source term. The terms **U**�*c*� represent the turbulent flux of contaminants (in units of *g sm*2 ), with its longitudinal, crosswind and vertical components.

Observe that equation (1) has four unknown variables (the concentration and turbulent fluxes) which lead us to the well known turbulence closure problem. One of the most widely used closures for equation (1), is based on the gradient transport hypothesis (also called K-theory) which, in analogy with Fick's law of molecular diffusion, assumes that turbulence causes a net movement of material following the negative gradient of material concentration at a rate which is proportional to the magnitude of the gradient (Seinfeld & Pandis (1998)).

$$\overline{\mathbf{U}^{\prime}c^{\prime}} = -\mathbf{K}\nabla\overline{c} \tag{2}$$

Here, the eddy diffusivity matrix **K** = diag(*Kx*, *Ky*, *Kz*) (in units of *m*2 *s* ) is a diagonal matrix with the Cartesian components in the *x*, *y* and *z* directions, respectively. In the first order closure all the information on the turbulence complexity is contained in the eddy diffusivity.

Equation (2), combined with the continuity equation of mass, leads to the advection-diffusion equation (see reference Blackadar (1997)).

$$
\partial\_t \boldsymbol{\tilde{c}} + \bar{\mathbf{U}}^T \nabla \boldsymbol{\tilde{c}} = \nabla^T \left( \mathbf{K} \nabla \boldsymbol{\tilde{c}} \right) + \mathbf{S} \tag{3}
$$

The simplicity of the K-theory of turbulent diffusion has led to the widespread use of this theory as a mathematical basis for simulating pollutant dispersion (open country, urban, photochemical pollution, etc.), but K-closure has its known limits. In contrast to molecular diffusion, turbulent diffusion is scale-dependent. This means that the rate of diffusion of a cloud of material generally depends on the cloud dimension and the intensity of turbulence. As the cloud grows, larger eddies are incorporated in the expansion process, so that a progressively larger fraction of turbulent kinetic energy is available for the cloud expansion.

Equation (3) is considered valid in the domain (*x*, *y*, *z*) ∈ Γ bounded by 0 *< x < Lx*, 0 *< y < Ly* and 0 *< z < h* and subject to the following boundary and initial conditions,

$$\mathbf{K}\nabla\overline{\boldsymbol{\varepsilon}}|\_{\left(0,0,0\right)} = \mathbf{K}\nabla\overline{\boldsymbol{\varepsilon}}|\_{\left(L\_{\times}L\_{g},h\right)} = \mathbf{0} \tag{4}$$

$$
\vec{\varepsilon}(x, y, z, 0) = 0 \,. \tag{5}
$$

Instead of specifying the source term as an inhomogeneity of the partial differential equation, we consider a point source located at an edge of the domain, so that the source position **r***S* = (0, *y*0, *HS*) is located at the boundary of the domain **r***S* ∈ *δ*Γ. Note, that in cases where the source is located in the domain, one still may divide the whole domain in sub-domains, where the source lies on the boundary of the sub-domains which can be solved for each sub-domain separately. Moreover, a set of different sources may be implemented as a superposition of independent problems. Since the source term location is on the boundary, in the domain this term is zero everywhere (*S*(**r**) ≡ 0 for **r** ∈ Γ\*δ*Γ), so that the source influence may be cast in form of a condition, where we assume that our coordinate system is oriented such that the *x*-axis is aligned with the mean wind direction. Since the flow crosses the plane perpendicular to the propagation (here the *y*-*z*-plane) the source condition reads

$$
\overline{u}\overline{c}(0, y, z, t) = Q\delta(y - y\_0)\delta(z - H\_\mathcal{S})\ , \tag{6}
$$

where *Q* is the emission rate (in units of *g s* ), *h* the height of the ABL (in units of *m*), *HS* the height of the source (in units of *m*), *Lx* and *Ly* are the horizontal domain limits (in units of *m*) and *δ*(*x*) represents the Cartesian Dirac delta functional.

#### **3. A closed form solution**

4 Will-be-set-by-IN-TECH

The advection-diffusion equation of air pollution in the atmosphere is essentially a statement

substantial derivative, *∂t* is a shorthand for the partial time derivative and **U**¯ = (*u*¯, *v*¯, *w*¯)*T*

respectively, ∇ = (*∂x*, *∂y*, *∂z*)*T* is the usual Nabla symbol and *S* is the source term. The

Observe that equation (1) has four unknown variables (the concentration and turbulent fluxes) which lead us to the well known turbulence closure problem. One of the most widely used closures for equation (1), is based on the gradient transport hypothesis (also called K-theory) which, in analogy with Fick's law of molecular diffusion, assumes that turbulence causes a net movement of material following the negative gradient of material concentration at a rate

with the Cartesian components in the *x*, *y* and *z* directions, respectively. In the first order closure all the information on the turbulence complexity is contained in the eddy diffusivity. Equation (2), combined with the continuity equation of mass, leads to the advection-diffusion

The simplicity of the K-theory of turbulent diffusion has led to the widespread use of this theory as a mathematical basis for simulating pollutant dispersion (open country, urban, photochemical pollution, etc.), but K-closure has its known limits. In contrast to molecular diffusion, turbulent diffusion is scale-dependent. This means that the rate of diffusion of a cloud of material generally depends on the cloud dimension and the intensity of turbulence. As the cloud grows, larger eddies are incorporated in the expansion process, so that a progressively larger fraction of turbulent kinetic energy is available for the cloud expansion.

Equation (3) is considered valid in the domain (*x*, *y*, *z*) ∈ Γ bounded by 0 *< x < Lx*, 0 *< y <*

Instead of specifying the source term as an inhomogeneity of the partial differential equation, we consider a point source located at an edge of the domain, so that the source position **r***S* =

*Ly* and 0 *< z < h* and subject to the following boundary and initial conditions,

*Dt* = *∂tc*¯ + **U**¯ *T*∇*c*¯ = −∇*T***U**�*c*� + *S* (1)

*s* ) with Cartesian components in the directions *x*, *y* and *z*,

**U**�*c*� = −**K**∇*c*¯ (2)

*∂tc*¯ + **U**¯ *T*∇*c*¯ = ∇*T* (**K**∇*c*¯) + *S* (3)

**K**∇*c*¯|(0,0,0) = **K**∇*c*¯|(*Lx*,*Ly*,*h*) = **0** (4)

*c*¯(*x*, *y*, *z*, 0) = 0 . (5)

*m*3 ), *D*

*sm*2 ), with its longitudinal,

*s* ) is a diagonal matrix

*Dt* is the

of conservation of the suspended material and it can be written as (Blackadar (1997))

where *c*¯ denotes the mean concentration of a passive contaminant (in units of *g*

which is proportional to the magnitude of the gradient (Seinfeld & Pandis (1998)).

Here, the eddy diffusivity matrix **K** = diag(*Kx*, *Ky*, *Kz*) (in units of *m*2

*Dc*¯

terms **U**�*c*� represent the turbulent flux of contaminants (in units of *g*

**2. The advection-diffusion equation**

is the mean wind (in units of *m*

crosswind and vertical components.

equation (see reference Blackadar (1997)).

In this section we first introduce the general formalism to solve a general problem and subsequently reduce the problem to a more specific one, that is solved and compared to experimental findings.

#### **3.1 The general procedure**

In order to solve the problem (3) we reduce the dimensionality by one and thus cast the problem into a form already solved in reference Moreira et al. (2009a). To this end we apply the integral transform technique in the *y* variable, and expand the pollutant concentration as

$$\vec{\boldsymbol{c}}(\mathbf{x}, \boldsymbol{y}, \boldsymbol{z}, \mathbf{t}) = \mathbf{R}^{\mathsf{T}}(\mathbf{x}, \boldsymbol{z}, \mathbf{t}) \mathbf{Y}(\boldsymbol{y}), \tag{7}$$

where **R** = (*R*1, *R*2,...)*T* and **Y** = (*Y*1,*Y*2,...)*T* is a vector in the space of orthogonal eigenfunctions, given by *Ym*(*y*) = cos(*λmy*) with eigenvalues *λm* = *m π Ly* for *m* = 0, 1, 2, . . .. For convenience we introduce some shorthand notations, ∇2 = (*∂x*, 0, *∂y*)*T* and ˆ *∂y* = (0, *∂y*, 0)*T*, so that equation (3) reads now,

$$\begin{aligned} \left(\partial\_l \mathbf{R}^T\right) \mathbf{Y} + \bar{\mathbf{U}} \left(\nabla\_2 \mathbf{R}^T \mathbf{Y} + \mathbf{R}^T \partial\_{\bar{y}} \mathbf{Y}\right) &= \left(\nabla^T \mathbf{K} + (\mathbf{K} \nabla)^T\right) \left(\nabla\_2 \mathbf{R}^T \mathbf{Y} + \mathbf{R}^T \partial\_{\bar{y}} \mathbf{Y}\right) \\ &= \left(\nabla\_2^T \mathbf{K} + (\mathbf{K} \nabla\_2)^T\right) \left(\nabla\_2 \mathbf{R}^T \mathbf{Y}\right) + \left(\partial\_{\bar{y}}^T \mathbf{K} + (\mathbf{K} \partial\_{\bar{y}})^T\right) \left(\mathbf{R}^T \partial\_{\bar{y}} \mathbf{Y}\right) . \end{aligned} \tag{8}$$

Upon application of the integral operator

$$\int\_{0}^{L\_{y}} dy \, \mathbf{Y}[\mathbf{F}] = \int\_{0}^{L\_{y}} \mathbf{F}^{T} \wedge \mathbf{Y} \, dy \tag{9}$$

here **F** is an arbitrary function and ∧ signifies the dyadic product operator, and making use of orthogonality renders equation (8) a matrix equation. The appearing integral terms are

where **Λ** = diag(*λ*2

1, *λ*2

which is equivalent to the problem

by virtue of **B** being a diagonal matrix.

dependence, which is novel in the literature.

time averaged solution for *Kz*.

recursion step on.

2,...). The simplified equation system to be solved is then,

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 45

*z***R***T***B** − *Ky***R***T***ΛB** (14)

*z***R** − *Ky***ΛR** (15)

*∂t***R***T***B** + *u*¯*∂x***R***T***B** = (*∂zKz*)*∂z***R***T***B** + *Kz∂*2

*∂t***R** + *u*¯*∂x***R** = (*∂zKz*)*∂z***R** + *Kz∂*2

The specific form of the eddy diffusivity determines now whether the problem is a linear or non-linear one. In the linear case the **K** is assumed to be independent of *c*¯, whereas in more realistic cases, even if stationary, **K** may depend on the contaminant concentration and thus renders the problem non-linear. However, until now no specific law is known that links the eddy diffusivity to the concentration so that we hide this dependence using a phenomenologically motivated expression for **K** which leaves us with a partial differential equation system in linear form, although the original phenomenon is non-linear. In the example below we demonstrate the closed form procedure for a problem with explicit time

The solution is generated making use of the decomposition method (Adomian (1984; 1988; 1994)) which was originally proposed to solve non-linear partial differential equations, followed by the Laplace transform that renders the problem a pseudo-stationary one. Further we rewrite the vertical diffusivity as a time average term *K*¯ *z*(*z*) plus a term representing the variations *κz*(*z*, *t*) around the average for the time interval of the measurement *Kz*(*x*, *z*, *t*) = *K*¯ *z*(*z*) + *κz*(*z*, *t*) and use the asymptotic form of *Ky*, which is then explored to set-up the

The function **R** = ∑*j* **R***j* = **1***T***R**(*c*) is now decomposed into contributions to be determined by recursion. For convenience we introduced the one-vector **1** = (1, 1, . . .)*T* and inflate the vector **R** to a vector with each element being itself a vector **R***j*. Upon inserting the expansion in equation (16) one may regroup terms that obey the recursive equations and starts with the

From the construction of the recursion equation system it is evident that other schemes are possible. The specific choice made here allows us to solve the recursion initialisation using the procedure described in references Moreira et al. (2006a; 2009a), where a stationary **K** was assumed. For this reason the time dependence enters as a known source term from the first

*∂t***R** + *u*¯*∂x***R** − *∂z* (*K*¯ *z∂z***R**) + *Ky***ΛR** = *∂z* (*κz∂z***R**) (16)

*∂t***R**0 + *u*¯*∂x***R**0 − *∂z* (*K*¯ *z∂z***R**0) + *Ky***ΛR**0 = 0 (17)

*κz∂z***R***j*−1

. (18)

+ *Ky***ΛR***j* = *∂z*

structure of the equation that defines the recursive decomposition scheme.

 *K*¯ *z∂z***R***j*

The extension to the closed form recursion is then given by

*∂t***R***j* + *u*¯*∂x***R***j* − *∂z*

**B**0 = *Ly* 0 *dy***Y**[**Y**] = *Ly* 0 **Y***T* ∧ **Y** *dy* , **Z** = *Ly* 0 *dy***Y**[ˆ *∂y***Y**] = *Ly* 0 ˆ *∂y***Y***T* ∧ **Y** *dy* , **Ω**1 = *Ly* 0 *dy***Y**[(∇*T* 2 **K**)(∇2**R***T***Y**)] = *Ly* 0 (∇*T* 2 **K**)(∇2**R***T***Y**) *T* ∧ **Y** *dy* , (10) **Ω**2 = *Ly* 0 *dy***Y**[(**K**∇2)*T*(∇2**R***T***Y**)] = *Ly* 0 (**K**∇2)*T*(∇2**R***T***Y**) ∧ **Y** *dy* , **T**1 = *Ly* 0 *dy***Y**[((ˆ *∂T y* **K**)(ˆ *∂y***Y**)] = *Ly* 0 ((ˆ *∂T y* **K**)(ˆ *∂y***Y**) *T* ∧ **Y** *dy* , **T**2 = *Ly* 0 *dy***Y**[(**K**ˆ *∂y*)*T*(ˆ *∂y***Y**)] = *Ly* 0 (**K**ˆ *∂y*)*T*(ˆ *∂y***Y**) *T* ∧ **Y** *dy* .

Here, **B**0 = *Ly* 2 **I**, where **I** is the identity, the elements (**Z**)*mn* = 2 1−*n*2/*m*2 *δ*1,*j* with *δi*,*j* the Kronecker symbol and *j* = (*m* + *n*)mod2 is the remainder of an integer division (i.e. is one for *m* + *n* odd and zero else). Note, that the integrals **Ω***i* and **T***i* depend on the specific form of the eddy diffusivity **K**. The integrals (10) are general, but for practical purposes and for application to a case study we truncate the eigenfunction space and consider *M* components in **R** and **Y** only, though continue using the general nomenclature that remains valid. The obtained matrix equation determines now together with initial and boundary condition uniquely the components *Ri* for *i* = 1, . . . , *M* following the procedure introduced in reference Moreira et al. (2009a).

$$(\partial\_l \mathbf{R}^T) \mathbf{B} + \mathbf{\bar{U}} \left(\nabla\_2 \mathbf{R}^T \mathbf{B} + \mathbf{R}^T \mathbf{Z}\right) = \mathbf{\Omega}\_1(\mathbf{R}) + \mathbf{\Omega}\_2(\mathbf{R}) + \mathbf{R}^T(\mathbf{T}\_1 + \mathbf{T}\_2) \tag{11}$$

#### **3.2 A specific case for application**

In order to discuss a specific case we recall the convention already adopted above, that considered the average wind velocity aligned with the *x*-axis. Since we consider length scales *Lx*, *LY*, *h* typical for micro-meteorological scenarios we simplify **U**¯ = (*u*¯, 0, 0)*T*. By comparison of physically meaningful cases, one finds for the operator norm ||*∂xKx∂x*|| *<<* |*u*¯|, which can be understood intuitively because eddy diffusion is observable predominantly perpendicular to the mean wind propagation. As a consequence we neglect the terms with *Kx* and *∂xKx*.

The principal aspect of interest in pollution dispersion is the vertical concentration profile, that responds strongly to the atmospheric boundary layer stratification, so that the simplified eddy diffusivity **K** → **K**1 = diag(0, *Ky*, *Kz*) depends in leading order approximation only on the vertical coordinate **K**1 = **K**1(*z*). For this specific case the integrals **Ω***i* reduce to

$$\begin{aligned} \mathbf{0}\_1 &\rightarrow (\partial\_z K\_z)(\partial\_z \mathbf{R}^T)\mathbf{B} \ , \\ \mathbf{0}\_2 &\rightarrow K\_z(\partial\_z^2 \mathbf{R}^T)\mathbf{B} \ , \end{aligned} \tag{12}$$

$$\begin{aligned} \mathbf{T}\_1 &\rightarrow \mathbf{0}\_{\text{-}1} \\ \mathbf{T}\_2 &\rightarrow -K\_y \mathbf{A} \mathbf{B} \end{aligned} \tag{13}$$

where **Λ** = diag(*λ*2 1, *λ*2 2,...). The simplified equation system to be solved is then,

$$
\partial\_l \mathbf{R}^T \mathbf{B} + \bar{u} \partial\_\mathbf{x} \mathbf{R}^T \mathbf{B} = (\partial\_z K\_z) \partial\_{\bar{z}} \mathbf{R}^T \mathbf{B} + K\_{\bar{z}} \partial\_{\bar{z}}^2 \mathbf{R}^T \mathbf{B} - K\_y \mathbf{R}^T \mathbf{A} \mathbf{B} \tag{14}
$$

which is equivalent to the problem

6 Will-be-set-by-IN-TECH

**Y***T* ∧ **Y** *dy* ,

*∂y***Y***T* ∧ **Y** *dy* ,

 *Ly* 0

 *Ly* 0

2 **I**, where **I** is the identity, the elements (**Z**)*mn* = 2

In order to discuss a specific case we recall the convention already adopted above, that considered the average wind velocity aligned with the *x*-axis. Since we consider length scales *Lx*, *LY*, *h* typical for micro-meteorological scenarios we simplify **U**¯ = (*u*¯, 0, 0)*T*. By comparison of physically meaningful cases, one finds for the operator norm ||*∂xKx∂x*|| *<<* |*u*¯|, which can be understood intuitively because eddy diffusion is observable predominantly perpendicular to the mean wind propagation. As a consequence we neglect the terms with *Kx* and *∂xKx*. The principal aspect of interest in pollution dispersion is the vertical concentration profile, that responds strongly to the atmospheric boundary layer stratification, so that the simplified eddy diffusivity **K** → **K**1 = diag(0, *Ky*, *Kz*) depends in leading order approximation only on

the vertical coordinate **K**1 = **K**1(*z*). For this specific case the integrals **Ω***i* reduce to

**Ω**2 → *Kz*(*∂*2

**T**1 → **0** ,

**Ω**1 → (*∂zKz*)(*∂z***R***T*)**B** ,

Kronecker symbol and *j* = (*m* + *n*)mod2 is the remainder of an integer division (i.e. is one for *m* + *n* odd and zero else). Note, that the integrals **Ω***i* and **T***i* depend on the specific form of the eddy diffusivity **K**. The integrals (10) are general, but for practical purposes and for application to a case study we truncate the eigenfunction space and consider *M* components in **R** and **Y** only, though continue using the general nomenclature that remains valid. The obtained matrix equation determines now together with initial and boundary condition uniquely the components *Ri* for *i* = 1, . . . , *M* following the procedure introduced in

 *Ly* 0

> *Ly* 0

 ((ˆ *∂T y* **K**)(ˆ *∂y***Y**) *T*

 (**K**ˆ *∂y*)*T*(ˆ *∂y***Y**) *T*

 (∇*T*

2 **K**)(∇2**R***T***Y**)

(**K**∇2)*T*(∇2**R***T***Y**)

*T*

∧ **Y** *dy* ,

∧ **Y** *dy* .

= **Ω**1(**R**) + **Ω**2(**R**) + **R***T*(**T**1 + **T**2) (11)

*z***R***T*)**B** , (12)

**T**2 → −*Ky***ΛB** , (13)

∧ **Y** *dy* ,

∧ **Y** *dy* , (10)

1−*n*2/*m*2 *δ*1,*j* with *δi*,*j* the

**B**0 =

**Z** =

**Ω**1 =

**Ω**2 =

**T**1 =

**T**2 =

reference Moreira et al. (2009a).

**3.2 A specific case for application**

(*∂t***R***T*)**B** + **U**¯

∇2**R***T***B** + **R***T***Z**

Here, **B**0 = *Ly*

 *Ly* 0

 *Ly* 0

 *Ly* 0

 *Ly* 0

 *Ly* 0

 *Ly* 0

*dy***Y**[**Y**] =

*dy***Y**[(∇*T*

*dy***Y**[((ˆ *∂T y* **K**)(ˆ

*dy***Y**[(**K**ˆ

*∂y***Y**] =

*dy***Y**[ˆ

 *Ly* 0

> *Ly* 0 ˆ

*dy***Y**[(**K**∇2)*T*(∇2**R***T***Y**)] =

*∂y*)*T*(ˆ

2 **K**)(∇2**R***T***Y**)] =

*∂y***Y**)] =

*∂y***Y**)] =

$$
\partial\_l \mathbf{R} + \bar{u} \partial\_{\bar{x}} \mathbf{R} = (\partial\_{\bar{z}} K\_{\bar{z}}) \partial\_{\bar{z}} \mathbf{R} + K\_{\bar{z}} \partial\_{\bar{z}}^2 \mathbf{R} - K\_{\mathcal{Y}} \mathbf{A} \mathbf{R} \tag{15}
$$

by virtue of **B** being a diagonal matrix.

The specific form of the eddy diffusivity determines now whether the problem is a linear or non-linear one. In the linear case the **K** is assumed to be independent of *c*¯, whereas in more realistic cases, even if stationary, **K** may depend on the contaminant concentration and thus renders the problem non-linear. However, until now no specific law is known that links the eddy diffusivity to the concentration so that we hide this dependence using a phenomenologically motivated expression for **K** which leaves us with a partial differential equation system in linear form, although the original phenomenon is non-linear. In the example below we demonstrate the closed form procedure for a problem with explicit time dependence, which is novel in the literature.

The solution is generated making use of the decomposition method (Adomian (1984; 1988; 1994)) which was originally proposed to solve non-linear partial differential equations, followed by the Laplace transform that renders the problem a pseudo-stationary one. Further we rewrite the vertical diffusivity as a time average term *K*¯ *z*(*z*) plus a term representing the variations *κz*(*z*, *t*) around the average for the time interval of the measurement *Kz*(*x*, *z*, *t*) = *K*¯ *z*(*z*) + *κz*(*z*, *t*) and use the asymptotic form of *Ky*, which is then explored to set-up the structure of the equation that defines the recursive decomposition scheme.

$$
\partial\_t \mathbf{R} + \bar{u} \partial\_{\mathbf{x}} \mathbf{R} - \partial\_{\mathbf{z}} \left( \bar{K}\_2 \partial\_{\mathbf{z}} \mathbf{R} \right) + K\_{\mathcal{Y}} \mathbf{A} \mathbf{R} = \partial\_{\mathbf{z}} \left( \kappa\_2 \partial\_{\mathbf{z}} \mathbf{R} \right) \tag{16}
$$

The function **R** = ∑*j* **R***j* = **1***T***R**(*c*) is now decomposed into contributions to be determined by recursion. For convenience we introduced the one-vector **1** = (1, 1, . . .)*T* and inflate the vector **R** to a vector with each element being itself a vector **R***j*. Upon inserting the expansion in equation (16) one may regroup terms that obey the recursive equations and starts with the time averaged solution for *Kz*.

$$
\partial\_t \mathbf{R}\_0 + \mathbf{i}\partial\_x \mathbf{R}\_0 - \partial\_z \left( \bar{K}\_z \partial\_{\bar{z}} \mathbf{R}\_0 \right) + K\_y \mathbf{A} \mathbf{R}\_0 = 0 \tag{17}
$$

The extension to the closed form recursion is then given by

$$
\left(\partial\_l \mathbf{R}\_j + \vec{u} \partial\_{\mathbf{x}} \mathbf{R}\_j - \partial\_z \left(\bar{K}\_z \partial\_{\bar{z}} \mathbf{R}\_j\right) + K\_y \mathbf{A} \mathbf{R}\_j = \partial\_z \left(\kappa\_z \partial\_{\bar{z}} \mathbf{R}\_{j-1}\right) \,. \tag{18}
$$

From the construction of the recursion equation system it is evident that other schemes are possible. The specific choice made here allows us to solve the recursion initialisation using the procedure described in references Moreira et al. (2006a; 2009a), where a stationary **K** was assumed. For this reason the time dependence enters as a known source term from the first recursion step on.

where **P**˜(*s*,*r*) denotes the Laplace Transform of **P**(*x*,*r*). Here **X**(−1) is the (inverse) matrix of

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 47

(*s***I** + **D**)*ii* = *s* + *di*. After performing the Laplace transform inversion of equation (23), we

where **G**(*x*,*r*) is the diagonal matrix with components (**G**)*ii* = *e*−*dix*. Further the still

The analytical time dependence for the recursion initialisation (19) is obtained upon applying

 *γ*+*i*∞ *γ*−*i*∞

To overcome the drawback of evaluating the line integral appearing in the above solution, we perform the calculation of this integral by the Gaussian quadrature scheme, which is exact if

where **a** and **p** are respectively vectors with the weights and roots of the Gaussian quadrature

in the *k*-th row of **pR**0. Note, *k* is a component from contraction with **a**. Concerning the issue of possible Laplace inversion schemes, it is worth mentioning that this approach is exact, however we are aware of the existence of other methods in the literature to invert the Laplace transformed functions (Valkó & Abate (2004), Abate & Valkó (2004)), but we restrict our attention in the considered problem to the Gaussian quadrature scheme, which is sufficient in the present case. Finally, the solution of the remaining equations of the recursive system (18) are attained in an analogue fashion expanding the source term in series and solving the resulting linear first order differential matrix equation with known and integrable source by

To establish the solution of the stationary, three and two-dimensional, ADE under Fickian flow regime, with and without longitudinal diffusion, one only need to take the limit *t* → ∞ of the previously derived solutions, which is equivalent to take *r* → 0 (Buske et al. (2007a)). Further models that make use of simplifications are not discussed in this chapter because they can be determined with the present developments and applying them to the models of references Buske et al. (2007b); Moreira et al. (2005; 2006b; 2009b; 2010); Tirabassi et al. (2008; 2009; 2011);

The solution procedure discussed in the previous section was coded in a computer program and will be available in future as a program library add-on. In order to illustrate the suitability of the discussed formulation to simulate contaminant dispersion in the atmospheric boundary

*t* **a***T*

unknown arbitrary constant matrix is given by **Ω** = **X**−1**P**(0,*r*).

**R**0(*x*, *z*, *t*) = 1

the integrand is a polynomial of degree 2*M* − 1 in the 1

scheme (Stroud & Secrest (1966)), and the argument (*x*, *z*, **p**

the Laplace transform technique as shown above.

**4. Solution validation against experimental data**

**3.4 Reduction to solutions of simpler models**

Wortmann et al. (2005) and others.

2*πi*

**R**0(*x*, *z*, *t*) = 1

1 **B**2 with diagonal eigenvalue matrix **D** and the entries of matrix

**P**(*x*,*r*)**C**(*z*)*e*

**pR**0(*x*, *z*,

*r* variable.

**p** *t* )

**P**(*x*,*r*) = **XG**(*x*,*r*)**X**−1**Ω** , (24)

*rt dr* . (25)

, (26)

*t* ) signifies the *k*-th component of **p**

the eigenvectors of matrix **B**−1

the inverse Laplace transform definition.

come out with

#### **3.3 Recursion initialisation**

The boundary conditions are now used to uniquely determine the solution. In our scheme the initialisation solution that contains **R**0 satisfies the boundary conditions (equations (4)-(6)) while the remaining equations satisfy homogeneous boundary conditions. Once the set of problems (18) is solved by the GILTT method, the solution of problem (3) is well determined. It is important to consider that we may control the accuracy of the results by a proper choice of the number of terms in the solution series.

In references Moreira et al. (2009a;b), a two dimensional problem with advection in the *x* direction in stationary regime was solved which has the same formal structure than (18) except for the time dependence. Upon rendering the recursion scheme in a pseudo-stationary form problem and thus matching the recursive structure of Moreira et al. (2009a;b), we apply the Laplace Transform in the *t* variable, (*t* → *r*) obtaining the following pseudo-steady-state problem.

$$
\sigma \tilde{\mathbf{R}}\_0 + \overline{\boldsymbol{n}} \partial\_{\mathbf{x}} \tilde{\mathbf{R}}\_0 = \partial\_{\mathbf{z}} \left( \mathcal{K}\_{\mathbf{z}} \partial\_{\mathbf{z}} \tilde{\mathbf{R}}\_0 \right) - \Lambda \, \mathcal{K}\_{\mathcal{Y}} \tilde{\mathbf{R}}\_0 \tag{19}
$$

The *x* and *z* dependence may be separated using the same reasoning as already introduced in (7). To this end we pose the solution of problem (19) in the form:

$$
\vec{\mathbf{R}}\_0 = \mathbf{P} \mathbf{C} \tag{20}
$$

where **C** = (*ζ*1(*z*), *ζ*2(*z*),...)*T* are a set of orthogonal eigenfunctions, given by *ζi*(*z*) = cos(*γlz*), and *γi* = *iπ*/*h* (for *i* = 0, 1, 2, . . .) are the set of eigenvalues.

Replacing equation (20) in equation (19) and using the afore introduced projector (9) now for the *z* dependent degrees of freedom *h* 0 *dz***C**[**F**] = *h* 0 **F***T* ∧ **C** *dz* yields a first order differential equation system.

$$\mathbf{^0 \partial\_{\lambda} \mathbf{P} + \mathbf{B}\_1^{-1} \mathbf{B}\_2 \mathbf{P} = \mathbf{0}} \, \, \, \, \, \, \, \, \tag{21}$$

where **P**0 = **P**0(*x*,*r*) and **B**−1 1 **B**2 come from the diagonalisation of the problem. The entries of matrices **B**1 and **B**2 are

$$\begin{split} (\mathbf{B}\_{1})\_{i,j} &= -\int\_{0}^{h} \overline{u} \zeta\_{i}(z) \zeta\_{j}(z) \, dz \\ (\mathbf{B}\_{2})\_{i,j} &= \int\_{0}^{h} \partial\_{z} K\_{z} \partial\_{z} \zeta\_{i}(z) \zeta\_{j}(z) \, dz - \gamma\_{i}^{2} \int\_{0}^{h} K\_{z} \zeta\_{i}(z) \zeta\_{j}(z) \, dz \\ &\quad - r \int\_{0}^{h} \zeta\_{i}(z) \zeta\_{j}(z) \, dz - \lambda\_{i}^{2} K\_{y} \int\_{0}^{h} \zeta\_{i}(z) \zeta\_{j}(z) \, dz \, . \end{split}$$

A similar procedure leads to the source condition for (21).

$$\mathbf{P}(0,r) = Q\mathbf{B}\_1^{-1} \int dz \mathbf{C}[\delta(z - H\_S)] \int dy \mathbf{Y}[\delta(y - y\_0)] = Q\mathbf{B}\_1^{-1} \left(\mathbf{C}(H\_S) \wedge \mathbf{1}\right) \left(\mathbf{1} \wedge \mathbf{Y}(y\_0)\right) \tag{22}$$

Following the reasoning of Moreira et al. (2009b) we solve (21) applying Laplace transform and diagonalisation of the matrix **B**−1 1 **B**2 = **XDX**−1 which results in

$$\mathbf{P}(s,r) = \mathbf{X}(s\mathbf{I}+\mathbf{D})^{-1}\mathbf{X}^{-1}\mathbf{P}(0,r) \tag{23}$$

8 Will-be-set-by-IN-TECH

The boundary conditions are now used to uniquely determine the solution. In our scheme the initialisation solution that contains **R**0 satisfies the boundary conditions (equations (4)-(6)) while the remaining equations satisfy homogeneous boundary conditions. Once the set of problems (18) is solved by the GILTT method, the solution of problem (3) is well determined. It is important to consider that we may control the accuracy of the results by a proper choice

In references Moreira et al. (2009a;b), a two dimensional problem with advection in the *x* direction in stationary regime was solved which has the same formal structure than (18) except for the time dependence. Upon rendering the recursion scheme in a pseudo-stationary form problem and thus matching the recursive structure of Moreira et al. (2009a;b), we apply the Laplace Transform in the *t* variable, (*t* → *r*) obtaining the following pseudo-steady-state

> *Kz∂z***R**˜ 0

The *x* and *z* dependence may be separated using the same reasoning as already introduced in

where **C** = (*ζ*1(*z*), *ζ*2(*z*),...)*T* are a set of orthogonal eigenfunctions, given by *ζi*(*z*) =

Replacing equation (20) in equation (19) and using the afore introduced projector (9) now for

0 *dz***C**[**F**] = *h*

*∂x***P** + **B**−1

*uζi*(*z*)*ζj*(*z*) *dz*

*∂zKz∂zζi*(*z*)*ζj*(*z*) *dz* − *γ*2

*ζi*(*z*)*ζj*(*z*) *dz* − *λ*2

− Λ *Ky***R**˜ 0 (19)

**R**˜ 0 = **PC** (20)

0 **F***T* ∧ **C** *dz* yields a first order differential

1 **B**2**P** = **0** , (21)

*Kzζi*(*z*)*ζj*(*z*) *dz*

1 (**C**(*HS*) ∧ **1**) (**1** ∧ **Y**(*y*0)) (22)

*ζi*(*z*)*ζj*(*z*) *dz* .

**P**˜(*s*,*r*) = **X**(*s***I** + **D**)−1**X**−1**P**(0,*r*) (23)

1 **B**2 come from the diagonalisation of the problem. The entries of

*i h* 0

*i Ky h* 0

*dy***Y**[*δ*(*y* − *y*0)] = *Q***B**−1

1 **B**2 = **XDX**−1 which results in

Following the reasoning of Moreira et al. (2009b) we solve (21) applying Laplace transform

*r***R**˜ 0 + *u∂x***R**˜ 0 = *∂z*

(7). To this end we pose the solution of problem (19) in the form:

cos(*γlz*), and *γi* = *iπ*/*h* (for *i* = 0, 1, 2, . . .) are the set of eigenvalues.

 *h* 0

 *h* 0

−*r h* 0

A similar procedure leads to the source condition for (21).

*dz***C**[*δ*(*z* − *HS*)]

**3.3 Recursion initialisation**

problem.

equation system.

of the number of terms in the solution series.

the *z* dependent degrees of freedom *h*

(**B**1)*i*,*j* = −

(**B**2)*i*,*j* =

where **P**0 = **P**0(*x*,*r*) and **B**−1

matrices **B**1 and **B**2 are

**P**(0,*r*) = *Q***B**−1

1

and diagonalisation of the matrix **B**−1

where **P**˜(*s*,*r*) denotes the Laplace Transform of **P**(*x*,*r*). Here **X**(−1) is the (inverse) matrix of the eigenvectors of matrix **B**−1 1 **B**2 with diagonal eigenvalue matrix **D** and the entries of matrix (*s***I** + **D**)*ii* = *s* + *di*. After performing the Laplace transform inversion of equation (23), we come out with

$$\mathbf{P}(\mathbf{x},r) = \mathbf{X}\mathbf{G}(\mathbf{x},r)\mathbf{X}^{-1}\mathbf{D} \;/\tag{24}$$

where **G**(*x*,*r*) is the diagonal matrix with components (**G**)*ii* = *e*−*dix*. Further the still unknown arbitrary constant matrix is given by **Ω** = **X**−1**P**(0,*r*).

The analytical time dependence for the recursion initialisation (19) is obtained upon applying the inverse Laplace transform definition.

$$\mathbf{R}\_0(\mathbf{x}, z, t) = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathbf{P}(\mathbf{x}, r) \mathbf{C}(z) e^{rt} \, dr \,. \tag{25}$$

To overcome the drawback of evaluating the line integral appearing in the above solution, we perform the calculation of this integral by the Gaussian quadrature scheme, which is exact if the integrand is a polynomial of degree 2*M* − 1 in the 1 *r* variable.

$$\mathbf{R}\_0(\mathbf{x}, z, t) = \frac{1}{t} \mathbf{a}^T \left( \mathbf{p} \mathbf{R}\_0(\mathbf{x}, z, \frac{\mathbf{p}}{t}) \right) \; , \tag{26}$$

where **a** and **p** are respectively vectors with the weights and roots of the Gaussian quadrature scheme (Stroud & Secrest (1966)), and the argument (*x*, *z*, **p** *t* ) signifies the *k*-th component of **p** in the *k*-th row of **pR**0. Note, *k* is a component from contraction with **a**. Concerning the issue of possible Laplace inversion schemes, it is worth mentioning that this approach is exact, however we are aware of the existence of other methods in the literature to invert the Laplace transformed functions (Valkó & Abate (2004), Abate & Valkó (2004)), but we restrict our attention in the considered problem to the Gaussian quadrature scheme, which is sufficient in the present case. Finally, the solution of the remaining equations of the recursive system (18) are attained in an analogue fashion expanding the source term in series and solving the resulting linear first order differential matrix equation with known and integrable source by the Laplace transform technique as shown above.

#### **3.4 Reduction to solutions of simpler models**

To establish the solution of the stationary, three and two-dimensional, ADE under Fickian flow regime, with and without longitudinal diffusion, one only need to take the limit *t* → ∞ of the previously derived solutions, which is equivalent to take *r* → 0 (Buske et al. (2007a)). Further models that make use of simplifications are not discussed in this chapter because they can be determined with the present developments and applying them to the models of references Buske et al. (2007b); Moreira et al. (2005; 2006b; 2009b; 2010); Tirabassi et al. (2008; 2009; 2011); Wortmann et al. (2005) and others.

#### **4. Solution validation against experimental data**

The solution procedure discussed in the previous section was coded in a computer program and will be available in future as a program library add-on. In order to illustrate the suitability of the discussed formulation to simulate contaminant dispersion in the atmospheric boundary

*u* (115 *m*) *u* (10 *m*) *u*∗ *L w*∗ *h* Run (*ms*−1) (*ms*−1) (*ms*−1) (*m*) (*ms*−1) (*m*) 1 3.4 2.1 0.36 -37 1.8 1980 2 10.6 4.9 0.73 -292 1.8 1920 3 5.0 2.4 0.38 -71 1.3 1120 4 4.6 2.5 0.38 -133 0.7 390 5 6.7 3.1 0.45 -444 0.7 820 6 13.2 7.2 1.05 -432 2.0 1300 7 7.6 4.1 0.64 -104 2.2 1850 8 9.4 4.2 0.69 -56 2.2 810 9 10.5 5.1 0.75 -289 1.9 2090

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 49

Table 1. Meteorological conditions of the Copenhagen experiment Gryning & Lyck (1984).

for large diffusion travel times with *cv* = 0.36 and (*f* ∗

where *L* is the Obukhov length in the surface layer.

**4.2 Numerical results**

scale and *u*∗ is the friction velocity.

*ψ*1

3 = [(1 − *z*

the discussed solutions against the experimental centreline concentrations.

following the findings of reference Panosfsky & Dutton (1988).

Moreover, for the lateral eddy diffusivity we used the asymptotic behaviour of equation (27)

The measurements of the contaminant dispersion in the atmospheric boundary layer consist typically from a sequence of samples over a time period. The experiment used to validate the previously introduced solution was carried out in the northern part of Copenhagen and is described in detail by Gryning & Lyck (1984). Several runs of the experiment with changing meteorological conditions were considered as reference in order to simulate time dependent contaminant dispersion in the boundary layer and to evaluate the performance of

The essential data of the experiment are reported in the following. This experiment consisted of a tracer released without buoyancy from a tower at a height of 115*m*, and a collection of tracer sampling units were located at the ground-level positions up to the maximum of three crosswind arcs. The sampling unit distances varied between two to six kilometers from the point of release. The site was mainly residential with a roughness length of the 0.6*m*. Table 1 summarizes the meteorological conditions of the Copenhagen experiment where *L* is the Obukhov length, *h* is the height of the convective boundary layer, *w*∗ is the convective velocity

The wind speed profile used in the simulations is described by a power law expressed

Here *uz* and *u*1 are the horizontal mean wind speed at heights *z* and *z*1 and *n* is an exponent that is related to the intensity of turbulence (Irwin (1989)). As is possible to see in Irwin (1989)

*n*

*uz u*1 = *z z*1

dissipation function *ψ* according to H*φ*jstrup (1982) and Degrazia et al. (1998) has the form

*h* )2(− *z L* )− 2

*m*)*v* = 0.66 *z*

3 + 0.75]

1 2 *h* (Degrazia et al. (1997)). The

(28)

layer, we evaluate the performance of the new solution against experimental ground-level concentrations.

#### **4.1 Time dependent vertical eddy diffusivity**

As already mentioned in the derivation of the solution, one has to make the specific choice for the turbulence parametrisation. From the physical point of view the turbulence parametrisation is an approximation to nature in the sense that we use a mathematical model as an approximated (or phenomenological) relation that can be used as a surrogate for the natural true unknown term, which might enter into the equation as a non-linear contribution. The reliability of each model strongly depends on the way turbulent parameters are calculated and related to the current understanding of the ABL (Degrazia (2005)).

The present parametrisation is based on the Taylor statistical diffusion theory and a kinetic energy spectral model. This methodology, derived for convective and moderately unstable conditions, provides eddy diffusivities described in terms of the characteristic velocity and length scales of energy-containing eddies. The time dependent eddy diffusivity has been derived by Degrazia et al. (2002) and can be expressed as the following formula.

$$\frac{K\_{\rm a}}{w\_{\ast}h} = \frac{0.583c\_{i}\psi^{\frac{2}{3}}\left(\frac{z}{h}\right)^{\frac{4}{3}}X^{\*}\left[0.55\left(\frac{z}{h}\right)^{\frac{2}{3}} + 1.03c\_{i}^{\frac{1}{2}}\psi^{\frac{1}{3}}\left(f\_{m}^{\*}\right)\_{i}^{\frac{2}{3}}X^{\*}\right]}{\left[0.55\left(\frac{z}{h}\right)^{\frac{2}{3}}\left(f\_{m}^{\*}\right)\_{i}^{\frac{1}{2}} + 2.06c\_{i}^{\frac{1}{2}}\psi^{\frac{1}{3}}\left(f\_{m}^{\*}\right)\_{i}X^{\*}\right]^{2}}\tag{27}$$

Here *α* = *x*, *y*, *z* are the Cartesian components, *i* = *u*, *v*, *w* are the longitudinal and transverse wind directions,

$$c\_i = \alpha\_i (0.5 \pm 0.05)(2\pi\kappa)^{-\frac{2}{5}}\text{ .}$$

*αi* = 1 and 4/3 for *u*, *v* and *w* components, respectively (Champagne et al. (1977)), *κ* = 0.4 is the von Karman constant, (*f* ∗ *m*)*i* is the normalized frequency of the spectral peak, *h* is the top of the convective boundary layer height, *w*∗ is the convective velocity scale, *ψ* is the non-dimensional molecular dissipation rate function and

$$X^\* = \frac{\mathfrak{x}w\_\*}{\overline{\mathfrak{u}}h}.$$

is the non-dimensional travel time, where *u* is the horizontal mean wind speed.

For horizontal homogeneity the convective boundary layer evolution is driven mainly by the vertical transport of heat. As a consequence of this buoyancy driven ABL, the vertical dispersion process of scalars is dominant when compared to the horizontal ones. Therefore, the present analysis will focus on the vertical time dependent eddy diffusivity only. This vertical eddy diffusivity can be derived from equation (27) by assuming that *cw* = 0.36 and

$$(f\_m^\*)\_{\overline{w}} = 0.55 \frac{\frac{z}{\overline{h}}}{1 - \exp(-\frac{4z}{\overline{h}}) - 0.0003 \exp(\frac{8z}{\overline{h}})}$$

for the vertical component (Degrazia et al. (2000)). Furthermore, considering the horizontal plane, we can idealize the turbulent structure as a homogeneous one with the length scale of energy-containing eddies being proportional to the convective boundary layer height *h*.

Table 1. Meteorological conditions of the Copenhagen experiment Gryning & Lyck (1984).

Moreover, for the lateral eddy diffusivity we used the asymptotic behaviour of equation (27) for large diffusion travel times with *cv* = 0.36 and (*f* ∗ *m*)*v* = 0.66 *z h* (Degrazia et al. (1997)). The dissipation function *ψ* according to H*φ*jstrup (1982) and Degrazia et al. (1998) has the form

$$\psi^{\frac{1}{3}} = \left[ (1 - \frac{z}{h})^2 (-\frac{z}{L})^{-\frac{2}{3}} + 0.75 \right]^{\frac{1}{2}}$$

where *L* is the Obukhov length in the surface layer.

#### **4.2 Numerical results**

10 Will-be-set-by-IN-TECH

layer, we evaluate the performance of the new solution against experimental ground-level

As already mentioned in the derivation of the solution, one has to make the specific choice for the turbulence parametrisation. From the physical point of view the turbulence parametrisation is an approximation to nature in the sense that we use a mathematical model as an approximated (or phenomenological) relation that can be used as a surrogate for the natural true unknown term, which might enter into the equation as a non-linear contribution. The reliability of each model strongly depends on the way turbulent parameters are calculated

The present parametrisation is based on the Taylor statistical diffusion theory and a kinetic energy spectral model. This methodology, derived for convective and moderately unstable conditions, provides eddy diffusivities described in terms of the characteristic velocity and length scales of energy-containing eddies. The time dependent eddy diffusivity has been

Here *α* = *x*, *y*, *z* are the Cartesian components, *i* = *u*, *v*, *w* are the longitudinal and transverse

*ci* = *αi*(0.5±0.05)(2*πκ*)− 2

*αi* = 1 and 4/3 for *u*, *v* and *w* components, respectively (Champagne et al. (1977)), *κ* = 0.4 is

top of the convective boundary layer height, *w*∗ is the convective velocity scale, *ψ* is the

*X*∗ = *xw*∗ *uh*

For horizontal homogeneity the convective boundary layer evolution is driven mainly by the vertical transport of heat. As a consequence of this buoyancy driven ABL, the vertical dispersion process of scalars is dominant when compared to the horizontal ones. Therefore, the present analysis will focus on the vertical time dependent eddy diffusivity only. This vertical eddy diffusivity can be derived from equation (27) by assuming that *cw* = 0.36 and

1 − exp(−4*z*

for the vertical component (Degrazia et al. (2000)). Furthermore, considering the horizontal plane, we can idealize the turbulent structure as a homogeneous one with the length scale of energy-containing eddies being proportional to the convective boundary layer height *h*.

*z h*

*h* ) − 0.0003 exp( 8*z*

*h* )

is the non-dimensional travel time, where *u* is the horizontal mean wind speed.

3 + 1.03*c*

3 ,

*m*)*i* is the normalized frequency of the spectral peak, *h* is the

1 2 *i ψ*1 3 (*f* ∗ *m*)*iX*∗

1 2 *i ψ*1 3 (*f* ∗ *m*) 2 3 *i X*∗

2 (27)

and related to the current understanding of the ABL (Degrazia (2005)).

0.583*ciψ*2 3 *z h* 4 3 *X*∗ 0.55 *z h* 2

non-dimensional molecular dissipation rate function and

(*f* ∗

*m*)*w* = 0.55

 0.55 *z h* 2 3 (*f* ∗ *m*) 1 3 *i* + 2.06*c*

derived by Degrazia et al. (2002) and can be expressed as the following formula.

concentrations.

**4.1 Time dependent vertical eddy diffusivity**

*Kα w*∗*h* =

the von Karman constant, (*f* ∗

wind directions,

The measurements of the contaminant dispersion in the atmospheric boundary layer consist typically from a sequence of samples over a time period. The experiment used to validate the previously introduced solution was carried out in the northern part of Copenhagen and is described in detail by Gryning & Lyck (1984). Several runs of the experiment with changing meteorological conditions were considered as reference in order to simulate time dependent contaminant dispersion in the boundary layer and to evaluate the performance of the discussed solutions against the experimental centreline concentrations.

The essential data of the experiment are reported in the following. This experiment consisted of a tracer released without buoyancy from a tower at a height of 115*m*, and a collection of tracer sampling units were located at the ground-level positions up to the maximum of three crosswind arcs. The sampling unit distances varied between two to six kilometers from the point of release. The site was mainly residential with a roughness length of the 0.6*m*. Table 1 summarizes the meteorological conditions of the Copenhagen experiment where *L* is the Obukhov length, *h* is the height of the convective boundary layer, *w*∗ is the convective velocity scale and *u*∗ is the friction velocity.

The wind speed profile used in the simulations is described by a power law expressed following the findings of reference Panosfsky & Dutton (1988).

$$\frac{\overline{u}\_z}{\overline{u}\_1} = \left(\frac{z}{z\_1}\right)^n \tag{28}$$

Here *uz* and *u*1 are the horizontal mean wind speed at heights *z* and *z*1 and *n* is an exponent that is related to the intensity of turbulence (Irwin (1989)). As is possible to see in Irwin (1989)

Run Recursion *c*(*x*, *y*, *z*, *t*) Run Recursion *c*(*x*, *y*, *z*, *t*)

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 51

1 2 10.21 3.65 6 2 7.78 2.20 1.46

2 2 7.09 2.28 7 2 6.61 3.32 1.65

3 2 9.84 7.88 4.50 8 2 5.99 4.87 1.91

4 2 15.23 9 2 5.65 3.04 2.19

0 7.14 4.11 2.47 1 9.11 4.77 3.78

3 6.65 4.47 3.73 4 6.65 4.48 3.73 5 6.68 4.48 3.73

Table 2. Pollutant concentrations for nine runs at various positions of the Copenhagen experiment and model prediction by the 3D-GILTT approach with time dependent eddy

To perform statistical comparisons between GILTT results against Copenhagen experimental data we consider the set of statistical indices described by Hanna (1989) and defined in by

*σoσp* ,

0,5(*σ*0+*σp* ),

*Co Cp* ,

*Co* ≤ 2,

5 2 6.96 4.46 3.73

• NMSE (normalized mean square error) = (*Co*−*Cp* )2

• FA2 = fraction of data (%, normalized to 1) for 0, 5 ≤ *Cp*

0,5(*Co*+*Cp* )

,

• COR (correlation coefficient) = (*Co*−*Co* )(*Cp*−*Cp* )

• FS (fractional standard deviations) = *σ*0−*σp*

• FB (fractional bias) = *Co*−*Cp*

diffusivity.

0 8.68 4.49 0 3.65 2.10 1.31 1 10.45 3.66 1 7.28 2.23 1.47

3 10.21 3.65 3 7.78 2.20 1.46 4 10.21 3.65 4 7.78 2.20 1.46 5 10.21 3.65 5 7.78 2.20 1.46 0 8.62 2.01 0 5.77 3.07 2.27 1 6.60 2.30 1 6.73 3.33 1.65

3 7.09 2.28 3 6.60 3.32 1.65 4 7.09 2.28 4 6.60 3.32 1.65 5 7.09 2.28 5 6.60 3.32 1.65 0 5.58 6.38 3.72 0 6.38 3.94 2.76 1 10.76 7.97 4.51 1 6.40 4.91 1.91

3 9.79 7.88 4.50 3 5.97 4.87 1.91 4 9.79 7.88 4.50 4 5.97 4.87 1.91 5 9.80 7.88 4.50 5 5.97 4.87 1.91 0 10.73 0 5.01 2.83 2.08 1 15.32 1 5.73 3.05 2.20

3 15.24 3 5.64 3.04 2.19 4 15.24 4 5.64 3.04 2.19 5 15.24 5 5.64 3.04 2.19

*m*3 )

depth (10−7 *g*

Fig. 1. The vertical eddy diffusivity *Kz*(*z*, *t*) for three different travel times (*t* = 500*s*, 1500*s*, 2500*s*) using run 8 of the Copenhagen experiment.

*n* = 0.1 is valid for a power low wind profile in unstable condition. Moreover, U.S.EPA suggests for rural terrain (as default values used in regulatory models) to use *n* = 0.15 for neutral condition (class D) and *n* = 0.1 for stability class C (moderately unstable condition). In Figure 1 we present a plot of *Kz*(*z*, *t*) for three different travel times (*t* = 500*s*, 1500*s*, 2500*s*) using run 8 of the Copenhagen experiment.

In order to exclude differences due to numerical uncertainties we define the numerical accuracy 10−4 of our simulations determining the suitable number of terms of the solution series. As an eye-guide we report in table 2 on the numerical convergence of the results, considering successively one, two, three and four terms in the solution series. One observes that the desired accuracy, for the solved problem solved is attained including only four terms in the truncated series, which is valid for all distances considered. Once the number of terms in the series solution is determined numerical comparisons of the 3D-GILTT results against experimental data may be performed and are presented in table 3.

Figure 2 shows the scatter plot of the centreline ground-level observed concentrations versus the simulated the 3D-GILTT model predictions, normalized by the emission rate and using two points in the time Gaussian Quadrature inversion (Moreira et al. (2006a), Stroud & Secrest (1966)). In the scatter diagram analysis, the closer the data are from the bisector line, the better are the results. The lateral lines indicate a factor of two (FA2), meaning if all the obtained data are between these lines FA2 equals to 1 (the maximum value for stochastically distributed data). From the scatter diagram (figure 2 one observes the fairly good simulation of dispersion data by the 3D-GILTT model, even for the two-point Gaussian quadrature scheme considered here.

0 50 100 150 200 250

 Kz (500,z)

Kz

Kz

(1500,z)

(2500,z)

Kz instável

*n* = 0.1 is valid for a power low wind profile in unstable condition. Moreover, U.S.EPA suggests for rural terrain (as default values used in regulatory models) to use *n* = 0.15 for neutral condition (class D) and *n* = 0.1 for stability class C (moderately unstable condition). In Figure 1 we present a plot of *Kz*(*z*, *t*) for three different travel times (*t* = 500*s*, 1500*s*, 2500*s*)

In order to exclude differences due to numerical uncertainties we define the numerical accuracy 10−4 of our simulations determining the suitable number of terms of the solution series. As an eye-guide we report in table 2 on the numerical convergence of the results, considering successively one, two, three and four terms in the solution series. One observes that the desired accuracy, for the solved problem solved is attained including only four terms in the truncated series, which is valid for all distances considered. Once the number of terms in the series solution is determined numerical comparisons of the 3D-GILTT results against

Figure 2 shows the scatter plot of the centreline ground-level observed concentrations versus the simulated the 3D-GILTT model predictions, normalized by the emission rate and using two points in the time Gaussian Quadrature inversion (Moreira et al. (2006a), Stroud & Secrest (1966)). In the scatter diagram analysis, the closer the data are from the bisector line, the better are the results. The lateral lines indicate a factor of two (FA2), meaning if all the obtained data are between these lines FA2 equals to 1 (the maximum value for stochastically distributed data). From the scatter diagram (figure 2 one observes the fairly good simulation of dispersion data by the 3D-GILTT model, even for the two-point Gaussian quadrature scheme considered

Fig. 1. The vertical eddy diffusivity *Kz*(*z*, *t*) for three different travel times (*t* = 500*s*, 1500*s*, 2500*s*) using run 8 of the Copenhagen experiment.

experimental data may be performed and are presented in table 3.

0,0

using run 8 of the Copenhagen experiment.

here.

0,2

0,4

z / h

0,6

0,8

1,0

Table 2. Pollutant concentrations for nine runs at various positions of the Copenhagen experiment and model prediction by the 3D-GILTT approach with time dependent eddy diffusivity.

To perform statistical comparisons between GILTT results against Copenhagen experimental data we consider the set of statistical indices described by Hanna (1989) and defined in by

0 2 4 6 8 10 12 14 16 18 20

sm-3)

Cp (10 -7

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 53

Fig. 2. Observed (*Co*) and predicted (*Cp*) scatter plot of centreline concentration using the Copenhagen dataset. Data between dotted lines correspond to ratio *Cp*/*Co* ∈ [0.5, 2].

In the present contribution we focused on an analytical description of pollution related phenomena in a micro-scale that allows to simulate dispersion in an computationally efficient procedure. The reason why adopting an analytical procedure instead of using the nowadays available computing power resides in the fact that once an analytical solution to a mathematical model is found one can claim that the problem has been solved. We provided a closed form solution that may be tailored for numerical applications such as to reproduce the solution within a prescribed precision. As a consequence the error analysis reduces to model validation only, in comparison to numerical approaches where in general it is not straight

Our starting point, i.e. the mathematical model, is the advection-diffusion equation, which we solved in 3 ⊕ 1 space-time dimensions for a general eddy diffusivity. The closed form solution is obtained using a hybrid approach by spectral theory together with integral transforms (in the present case the Laplace transform) and the decomposition method. The general formalism was simplified in order to attend the meteorological situation of the Copenhagen experiment. By comparison the present approach was found to yield an acceptable solution for the time dependent three dimensional advection-diffusion equation and moreover predicted tracer concentrations closer to observed values compared to other approaches from the literature. Although K-closure is known to have its limitations comparison of measurements and theoretical predictions corresponded on a satisfactory level and thus supported the usage of such an approach for micro-scale dispersion phenomena. Note, that the dispersion of the experimental data in comparison to the predictions might suggest a considerable discrepancy between theory and experiment, but it is worth mentioning that

0

forward to disentangle model errors from numerical ones.

2

4

6

8

Co (10 -7sm-3)

**5. Conclusion**

10

12

14

16

18

Table 3. Numerical convergence of the 3D-GILTT model with time dependent eddy diffusivity for the 9 runs of the Copenhagen experiment.

Table 4. Statistical comparison between 3D-GILTT model results and the Copenhagen data set, changing the number of terms in equation (18).

where the subscripts *o* and *p* refer to observed and predicted quantities, respectively, and the bar indicates an averaged value. The best results are expected to have values near zero for the indices NMSE, FB and FS, and near 1 in the indices COR and FA2. Table 4 shows the findings of the statistical indices that show a fairly good agreement between the 3D-GILTT predictions and the experimental data. Moreover, the splitting proposed for the eddy diffusivity coefficient as a sum of the averaged eddy diffusivity coefficient plus time variation, appears to be a valid assumption, since we got compact convergence of the solution, in the sense that we attained results with accuracy of 10−4 with only a few terms in the solution series for all the distances considered.

Fig. 2. Observed (*Co*) and predicted (*Cp*) scatter plot of centreline concentration using the Copenhagen dataset. Data between dotted lines correspond to ratio *Cp*/*Co* ∈ [0.5, 2].

#### **5. Conclusion**

14 Will-be-set-by-IN-TECH

Run Distance Observed (Co) Predictions (Cp) (m) (10−7*sm*−3) (10−7*sm*−3)

1 1900 10.5 10.21 3700 2.14 3.65 2 2100 9.85 7.09 4200 2.83 2.28 3 1900 16.33 9.80 3700 7.95 7.88 5400 3.76 4.50 4 4000 15.71 15.24 5 2100 12.11 6.68 4200 7.24 4.48 6100 4.75 3.73 6 2000 7.44 7.78 4200 3.47 2.20 5900 1.74 1.46 7 2000 9.48 6.60 4100 2.62 3.32 5300 1.15 1.65 8 1900 9.76 5.97 3600 2.64 4.87 5300 0.98 1.91 9 2100 8.52 5.64 4200 2.66 3.04 6000 1.98 2.19

Table 3. Numerical convergence of the 3D-GILTT model with time dependent eddy

Recursion NMSE COR FA2 FB FS

Table 4. Statistical comparison between 3D-GILTT model results and the Copenhagen data

where the subscripts *o* and *p* refer to observed and predicted quantities, respectively, and the bar indicates an averaged value. The best results are expected to have values near zero for the indices NMSE, FB and FS, and near 1 in the indices COR and FA2. Table 4 shows the findings of the statistical indices that show a fairly good agreement between the 3D-GILTT predictions and the experimental data. Moreover, the splitting proposed for the eddy diffusivity coefficient as a sum of the averaged eddy diffusivity coefficient plus time variation, appears to be a valid assumption, since we got compact convergence of the solution, in the sense that we attained results with accuracy of 10−4 with only a few terms in the solution

0 0.38 0.83 0.83 0.32 0.59 1 0.16 0.90 1.00 0.11 -0.13 2 0.14 0.91 1.00 0.15 -0.07 3 0.14 0.91 1.00 0.15 -0.07 4 0.14 0.91 1.00 0.15 -0.07

diffusivity for the 9 runs of the Copenhagen experiment.

set, changing the number of terms in equation (18).

series for all the distances considered.

depth

In the present contribution we focused on an analytical description of pollution related phenomena in a micro-scale that allows to simulate dispersion in an computationally efficient procedure. The reason why adopting an analytical procedure instead of using the nowadays available computing power resides in the fact that once an analytical solution to a mathematical model is found one can claim that the problem has been solved. We provided a closed form solution that may be tailored for numerical applications such as to reproduce the solution within a prescribed precision. As a consequence the error analysis reduces to model validation only, in comparison to numerical approaches where in general it is not straight forward to disentangle model errors from numerical ones.

Our starting point, i.e. the mathematical model, is the advection-diffusion equation, which we solved in 3 ⊕ 1 space-time dimensions for a general eddy diffusivity. The closed form solution is obtained using a hybrid approach by spectral theory together with integral transforms (in the present case the Laplace transform) and the decomposition method. The general formalism was simplified in order to attend the meteorological situation of the Copenhagen experiment. By comparison the present approach was found to yield an acceptable solution for the time dependent three dimensional advection-diffusion equation and moreover predicted tracer concentrations closer to observed values compared to other approaches from the literature. Although K-closure is known to have its limitations comparison of measurements and theoretical predictions corresponded on a satisfactory level and thus supported the usage of such an approach for micro-scale dispersion phenomena. Note, that the dispersion of the experimental data in comparison to the predictions might suggest a considerable discrepancy between theory and experiment, but it is worth mentioning that

Concluding, analytical solutions of equations are of fundamental importance in understanding and describing physical phenomena, since they might take into account all the parameters of a problem, and investigate their influence. Moreover, when using models, while they are rather sophisticated instruments that ultimately reflect the current state of knowledge on turbulent transport in the atmosphere, the results they provide are subject to a considerable margin of error. This is due to various factors, including in particular the uncertainty of the intrinsic variability of the atmosphere. Models, in fact, provide values expressed as an average, i.e., a mean value obtained by the repeated performance of many experiments, while the measured concentrations are a single value of the sample to which the ensemble average provided by models refer. This is a general characteristic of the theory of atmospheric turbulence and is a consequence of the statistical approach used in attempting to parametrise the chaotic character of the measured data. An analytical solution can be useful in evaluating the performances of numerical models (that solve numerically the advection-diffusion equation) that could compare their results, not only against experimental data but, in an easier way, with the solution itself in order to check numerical errors without the uncertainties presented above. Finally, the program of providing analytical solutions for realistic physical problems, leads us to future problems with different closure hypothesis considering full space-time dependence in the resulting dynamical equation, which we will

Analytical Model for Air Pollution in the Atmospheric Boundary Layer 55

The authors thank to CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico)

Abate, J. & Valkó, P.P. (2004). Multi-precision Laplace transform inversion. *Int. J. for Num.*

Adomian, G. (1984). A New Approach to Nonlinear Partial Differential Equations. *J. Math.*

Adomian, G. (1988). A Review of the Decomposition Method in Applied Mathematics. *J. Math.*

Adomian, G. (1994). *Solving Frontier Problems of Physics: The Decomposition Method*, Kluwer,

Blackadar, A.K. (1997). *Turbulence and diffusion in the atmosphere: lectures in Environmental*

Bodmann, B.; Vilhena, M.T.; Ferreira, L.S. & Bardaji, J.B. (2010). An analytical solver for

Boichenko, V.A.; Leonov, G.A. & Reitmann, V. (2005). *NDimension theory for ordinary equations*,

Buske, D.; Vilhena, M.T.; Moreira, D.M. & Tirabassi, T. (2007a). An analytical solution of the

techniques. *Il Nuovo Cimento C*, Vol. 33, page numbers (199-206).

the multi-group two dimensional neutron-diffusion equation by integral transform

advection-diffusion equation considering non-local turbulence closure. *Environ. Fluid*

*Methods in Engineering*, Vol. 60, page numbers (979-993).

*Anal. Appl.*, Vol. 102, page numbers (420-434).

*Anal. Appl.*, Vol. 135, page numbers (501-544) .

*Mechanics*, Vol. 7, page numbers (43-54).

also approach by the proposed methodology.

for the partial financial support of this work.

**6. Acknowledgements**

Boston, MA. .

*Sciences*, Springer-Verlag.

Teubner, Stuttgart.

**7. References**

the measurements are a unique sample of a distribution around an average value, whereas the prediction of an average value is evaluated from a deterministic equation, where the stochastic character is hidden in the turbulence closure hypothesis, so that a spread of data along the bisector is to be expected.

A data dispersion due to numerical uncertainties may be excluded using convergence criteria to control the numerical precision. The quality of the solution is controlled by a genuine mathematical convergence criterion. Note, that for the *t* and *x* coordinate the Laplace inversion considers only bi-Lipschitz functions, which defines then a unique relation between the original function and its Laplace-transform. This makes the transform procedure manifest exact and the only numerical error comes from truncation, which is determined from the Sturm-Liouville problem. In order to determine the truncation index of the solution series we introduced a carbon-copy of the Cardinal theorem of interpolation theory.

Recalling, that the structure of the pollutant concentration is essentially determined by the mean wind velocity **U**¯ and the eddy diffusivity **K**, means that the quotient of norms *�* = ||**K**|| ||**U**¯ || defines a length scale for which the pollutant concentration is almost homogeneous. Thus one may conclude that with decreasing length ( *� m* and *m* an increasing integer number) variations in the solution become spurious. Upon interpreting *�*−1 as a sampling density, one may now employ the Cardinal Theorem of Interpolation Theory (Torres (1991)) in order to find the truncation that leaves the analytical solution almost exact, i.e. introduces only functions that vary significantly in length scales beyond the mentioned limit.

The square integrable function *χ* = *r c dt dx d* ¯ *η* ∈ *L*2 (*η* = *y* or *z*) with spectrum {*λi*} which is bounded by *m�*−1 has an exact solution for a finite expansion. This statement expresses the *Cardinal Theorem of Interpolation Theory* for our problem. Since the cut-off defines some sort of sampling density, its introduction is an approximation and is related to convergence of the approach and Parseval's theorem may be used to estimate the error. In order to keep the solution error within a prescribed error, the expansion in the region of interest has to contain *n* + 1 terms, with *n* = int *mLy*,*z* 2*π�* + 1 2 . For the bounded spectrum and according to the theorem the solution is then exact. In our approximation, if *m* is properly chosen such that the cut-off part of the spectrum is negligible, then the found solution is almost exact.

Further, the Cauchy-Kowalewski theorem (Courant & Hilbert (1989)) guarantees that the proposed solution is a valid solution of the discussed problem, since this problem is a special case of the afore mentioned theorem, so that existence and uniqueness are guaranteed. It remains to justify convergence of the decomposition method. In general convergence by the decomposition method is not guaranteed, so that the solution shall be tested by an appropriate criterion. Since standard convergence criteria do not apply in a straight forward manner for the present case, we resort to a method which is based on the reasoning of Lyapunov (Boichenko et al. (2005)). While Lyapunov introduced this conception in order to test the influence of variations of the initial condition on the solution, we use a similar procedure to test the stability of convergence while starting from an approximate (initial) solution **R**0 (the seed of the recursive scheme). Let |*δZn*| = � ∑∞ *i*=*n*+1 **R***i*� be the maximum deviation of the correct from the approximate solution Γ*n* = ∑*n i*=0 **R***i*, where �·� signifies the maximum norm. Then strong convergence occurs if there exists an *n*0 such that the sign of *λ* is negative for all *n* ≥ *n*0. Here, *λ* = 1 �Γ*n*� log |*δZn*| |*δZ*0| .

Concluding, analytical solutions of equations are of fundamental importance in understanding and describing physical phenomena, since they might take into account all the parameters of a problem, and investigate their influence. Moreover, when using models, while they are rather sophisticated instruments that ultimately reflect the current state of knowledge on turbulent transport in the atmosphere, the results they provide are subject to a considerable margin of error. This is due to various factors, including in particular the uncertainty of the intrinsic variability of the atmosphere. Models, in fact, provide values expressed as an average, i.e., a mean value obtained by the repeated performance of many experiments, while the measured concentrations are a single value of the sample to which the ensemble average provided by models refer. This is a general characteristic of the theory of atmospheric turbulence and is a consequence of the statistical approach used in attempting to parametrise the chaotic character of the measured data. An analytical solution can be useful in evaluating the performances of numerical models (that solve numerically the advection-diffusion equation) that could compare their results, not only against experimental data but, in an easier way, with the solution itself in order to check numerical errors without the uncertainties presented above. Finally, the program of providing analytical solutions for realistic physical problems, leads us to future problems with different closure hypothesis considering full space-time dependence in the resulting dynamical equation, which we will also approach by the proposed methodology.

#### **6. Acknowledgements**

The authors thank to CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) for the partial financial support of this work.

#### **7. References**

16 Will-be-set-by-IN-TECH

the measurements are a unique sample of a distribution around an average value, whereas the prediction of an average value is evaluated from a deterministic equation, where the stochastic character is hidden in the turbulence closure hypothesis, so that a spread of data

A data dispersion due to numerical uncertainties may be excluded using convergence criteria to control the numerical precision. The quality of the solution is controlled by a genuine mathematical convergence criterion. Note, that for the *t* and *x* coordinate the Laplace inversion considers only bi-Lipschitz functions, which defines then a unique relation between the original function and its Laplace-transform. This makes the transform procedure manifest exact and the only numerical error comes from truncation, which is determined from the Sturm-Liouville problem. In order to determine the truncation index of the solution series

Recalling, that the structure of the pollutant concentration is essentially determined by the mean wind velocity **U**¯ and the eddy diffusivity **K**, means that the quotient of norms *�* = ||**K**||

defines a length scale for which the pollutant concentration is almost homogeneous. Thus one

in the solution become spurious. Upon interpreting *�*−1 as a sampling density, one may now employ the Cardinal Theorem of Interpolation Theory (Torres (1991)) in order to find the truncation that leaves the analytical solution almost exact, i.e. introduces only functions that

is bounded by *m�*−1 has an exact solution for a finite expansion. This statement expresses the *Cardinal Theorem of Interpolation Theory* for our problem. Since the cut-off defines some sort of sampling density, its introduction is an approximation and is related to convergence of the approach and Parseval's theorem may be used to estimate the error. In order to keep the solution error within a prescribed error, the expansion in the region of interest has to

the theorem the solution is then exact. In our approximation, if *m* is properly chosen such that