InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Earth and Planetary Sciences » "Atmospheric Hazards - Case Studies in Modeling, Communication, and Societal Impacts", book edited by Jill S. M. Coleman, ISBN 978-953-51-2516-7, Print ISBN 978-953-51-2515-0, Published: August 17, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 6

Computational Tools for the Simulation of Atmospheric Pollution Transport During a Severe Wind Event in Argentina

By César Augusto Aguirre and Armando Benito Brizuela
DOI: 10.5772/63552

Article top


Scheme of Gong’s experiments [13, 34].
Figure 1. Scheme of Gong’s experiments [13, 34].
Profiles of normalized concentrations upwind hill. Left-hand side: x – xs = 150 mm (foot of the hill). Right-hand side: x – xs = 250 mm (halfway on the hill).
Figure 2. Profiles of normalized concentrations upwind hill. Left-hand side: x – xs = 150 mm (foot of the hill). Right-hand side: x – xs = 250 mm (halfway on the hill).
Profiles of mean concentration at the crest of the hill at x – xs = 350 mm.
Figure 3. Profiles of mean concentration at the crest of the hill at x – xs = 350 mm.
Profiles of normalized concentration downwind hill. Left-hand side: x – xs = 450 mm (halfway on the hill). Right-hand side: x – xs = 550 mm (foot of the hill).
Figure 4. Profiles of normalized concentration downwind hill. Left-hand side: x – xs = 450 mm (halfway on the hill). Right-hand side: x – xs = 550 mm (foot of the hill).
Normalized mean concentration at ground level (left) and standard deviation from the height of the gas plume centre (right).
Figure 5. Normalized mean concentration at ground level (left) and standard deviation from the height of the gas plume centre (right).
Mean gas concentration values simulated with LES-STO. Left-hand side: on the axial plane containing the source. Right-hand side: at ground level.
Figure 6. Mean gas concentration values simulated with LES-STO. Left-hand side: on the axial plane containing the source. Right-hand side: at ground level.
Location of study area and sample measurements.
Figure 7. Location of study area and sample measurements.
Horizontal (a) and vertical (b) configuration of the grid for LES computation and the location of two chimneys of smelter.
Figure 8. Horizontal (a) and vertical (b) configuration of the grid for LES computation and the location of two chimneys of smelter.
Calculation scheme of the initial vertical velocity of solid particles.
Figure 9. Calculation scheme of the initial vertical velocity of solid particles.
Copper concentration at ground level. Solid colour: simulated with 100 × 100 m grid mesh (part/m2). Lines: measurements (mgCu/kgsoil) with 100 × 100 m by (Fernández-Turiel [14]).
Figure 10. Copper concentration at ground level. Solid colour: simulated with 100 × 100 m grid mesh (part/m2). Lines: measurements (mgCu/kgsoil) with 100 × 100 m by (Fernández-Turiel [14]).
Copper concentration at ground level. Solid colour: simulated with 10 × 10 m grid mesh (part/m2). Lines: measurements (mgCu/kgsoil) with 100 × 100 m by (Fernández-Turiel [14]).
Figure 11. Copper concentration at ground level. Solid colour: simulated with 10 × 10 m grid mesh (part/m2). Lines: measurements (mgCu/kgsoil) with 100 × 100 m by (Fernández-Turiel [14]).
Copper concentration at ground level simulated with 10 × 10 m grid mesh on a satellite image.
Figure 12. Copper concentration at ground level simulated with 10 × 10 m grid mesh on a satellite image.

Computational Tools for the Simulation of Atmospheric Pollution Transport During a Severe Wind Event in Argentina

César Augusto Aguirre1, 2 and Armando Benito Brizuela2
Show details


This chapter details the theoretical aspects of numerical methods for the simulation of atmospheric phenomena, such as severe thunderstorms and turbulent transport of the dangerous gases and solid particles into the atmospheric boundary layer. Numerical methods are included in computational algorithms to solve large turbulent scales using large eddy simulation (LES) techniques to obtain acceptable results of turbulent flows. However, microphysics processes involving evaporation, condensation and precipitation water using LES techniques are parameterized. These atmospheric processes are simulated using the advanced regional prediction systems (ARPS) code. On the contrary, atmospheric transport of pollutants is simulated using ARPS code coupled with a Lagrangian stochastic one-particle method. The theoretical details of this coupling are presented. Later, we show some laboratory experiments of plume dispersion emitted from gaseous sources, and the results of the computational simulation tool are compared after obtaining good agreement of the gas concentrations on the stream-wise vertical plane and over the ground. Finally, we present a simulation of a pollution event of copper solid particles at San Miguel de Tucumán city, Argentina. The geographical distributions of copper particle concentrations are in good agreement with the measurements carried out experimentally.

Keywords: large eddy simulation, solid particle dispersion, smelter, severe wind events, Argentina

1. Introduction

In the last few years, diverse governmental and civic organizations have expressed the need to evaluate, reduce and legislate atmospheric emanations from factories, industries and motor vehicles so as to mitigate respiratory illnesses in residents living in large cities close to industrial zones. One of the tools in use to evaluate the adverse effects of industrial installations at the local level is the simulation of air pollution episodes of low, average and extreme operating conditions. Advances in hardware and software computing have led to implement codes of pollution simulation to aim at solving the balance equations of fluid mechanics in conjunction with chemical reaction models. These computation codes aim to estimate the concentration values and spatial distribution of chemical species dispersed in the atmosphere under different weather conditions. However, there are still limitations of computer power due to the large number of calculations needed to solve turbulent flows of a high Reynolds number (ℜe) as those normally are present in the atmospheric boundary layer. By using a large eddy simulation (LES), the instantaneous evolution of large turbulent structures can be computed with the balance equations of fluid mechanics. This technique of discretization in space and time of turbulent transport phenomena has been developed from the early work of Deardorff [1, 2] and Schumann [3].

The large-scale resolution by LES allows a three-dimensional description of the wind field with a spatial resolution limited by the size of the cell that has subdivided the computational domain. While no major complications arise due to this limitation, certain meso- and macro-scale applications require a small cell grid size for computing molecular diffusion and chemical reactions with a high ℜe number, such as near pollution emission sources. One option to avoid too much reduce the size of the cell grid is considered fluid particles leading concentrations of chemical species (e.g. concentration of CO2, CO, N2O, NO, O3). This approach requires designing a model to get the trajectories of these particles by computing at each time step of the simulation, the position and velocity of the particles. These models are known as Lagrangian models (LM). When fluid dynamics are described by LES, the particles are driven following the movement of large scales of turbulence. Some authors added to the movement induced by LES one-random motion component that simulates the behaviour of smaller scales according to Brownian motion [412]. Models based on this technique are called Lagrangian stochastic models (STOs).

This chapter describes the assumptions and theoretical considerations of the coupling between these STO models with LES. Next, simulation results are compared with CO2 measurements emitted from an upwind source made by Gong [13]. Finally, model results are compared with measurements of copper concentration made by Fernández-Turiel [14] in Lastenia town, province of Tucuman, Argentina, caused by pollution from a smelter during a severe wind event.

2. Numerical methods

2.1. The large eddy simulation (LES)

In a fluid turbulent flow, many eddies of different sizes are induced. A technique for obtaining the time-space evolution of the most important eddies based on the energy they carry is known as large eddy simulation of turbulence (LES). The LES technique is an important tool for the simulation of wind turbulence into the atmosphere because the technique allows one to obtain a three-dimensional description of the wind field and its temporal evolution. The code used for LES is the ARPS (version 5.2.12), a mesoscale model of the non-hydrostatic and fully compressible type, developed by Center of Analysis and Prediction of Storms (CAPS) at the University of Oklahoma (USA). The model numerically integrates the time-dependent equations of mass balance, quantity of movement and energy of the largest turbulent scales. This model not only simulates the wind field but also has sub-models of heat and vapour flow, cloud formation and rainfall. For this, the orography and land cover are considered as well as the initial conditions of the ground and the atmospheric boundary layer.

The continuity, momentum equations and energy are resolved using the scheme of finite differences centred on an Arakawa C-grid cell type. A fully three-dimensional curvilinear coordinate system is used where the size of horizontal grids is constant, but the vertical coordinate follows the terrain elevation, and a stretching in size is applied to obtain more accuracy near the ground. The atmospheric model takes into account the compressibility of the flow. The numerical scheme used to obtain the solution of the differential equations is of 4th order centred of the explicit type, while that used to integrate the equations of pressure and vertical component of the air speed is implicit of Crank-Nicholson type. The prognostic variables of the model are Cartesian (0xyz) wind components u, v, w and scalars potential temperature θ, pressure p, air density ρ, mixing ratio of water vapour qv, cloud water qc, ice qi, rainwater qr, snow qs and hail qh. Initially, the states of these variables are included according to Reynolds decomposition (1) in a base-state ā and a perturbation a’.

{u(x,y,z,t)=u¯(z)+u'(x,y,z,t)v(x,y,z,t)=v¯(z)+u'(x,y,z,t)w(x,y,z,t)= w'(x,y,z,t)θ(x,y,z,t)=θ¯(z)+θ'(x,y,z,t)p(x,y,z,t)=p¯(z)+p'(x,y,z,t)ρ(x,y,z,t)=ρ¯(z)+ρ'(x,y,z,t)q(x,y,z,t)=q¯(z)+q'(x,y,z,t)

In this decomposition, the base state is assumed to be horizontally homogeneous. For this reason, the vertical component of the base state of the wind velocity is zero. Further, the base state is time invariant and hydrostatically balanced, so that the perturbation is integrated numerically in every time step for the filtered continuity, filtered momentum equation of wind velocity and filtered momentum equation of the scalars ψ. The filtered operation is carried out to obtain large scales of turbulent flow. This involves the application of a convolution spatial filter Gxi) where ∫Gxi)dxi = 1, a low pass filter, and Δxi is the size of grid elements of computational spatial domain. This filtered operation is applied to velocity champs of wind given the large scale for this variable:


where D is the entire spatial domain. The non-resolved scales are designed as ui=uiui.

Filtered continuity (3), filtered momentum of fluid velocity (4) and filtered momentum of scalars (5) are described as


The Leonard identity [15] is applied to Eqs. 3–5 as follows:

(u˜iuj)=(u˜iuj)+(u˜iuj+u˜iuj)+(u˜iuj)           =(u˜iuj)+L˜ij+C˜ij+R˜ij           =(u˜iuj)+τ˜ij
(u˜iψ)=(u˜iψ)+(u˜iψ+u˜iψ)+(u˜iψ)           =(u˜iψ)+τ˜iψ

In Eq. (4), B is a buoyancy force, Sija is the anisotropic deformation tensor and ν is the molecular viscosity. In Eq. (5), Φψ is the sink and source of the scalar variable ψ. The variables with tilde indicate that they have been weighted by the density state base u˜i=ρ¯(z)ui. The pressure equation is obtained by taking the material derivative of the state equation for moist air and replacing the time derivative of density by velocity divergence using the continuity equation. The correlation terms containing unsolved scales τ˜ij and τ˜iψ are modelled using Smagorinsky formulation [16, 17] with a dynamic scheme [18].

2.2. The dynamic Smagorinsky model

The unsolved scales can be viewed as the sub-grid scale viscosity due to the effect of small vortices whose size is less than Δ. The effect of small scales is to transfer the kinetic energy from large scales for dissipation following the theatrical energy cascade of Kolmogorov; then this behaviour can be modelled as a sub-grid stress (SGS) tensor τ˜ij, τ˜iψ. Furthermore, while large-scale motions are strongly dependent on the external flow conditions, small-scale motions are expected to behave more universally. Hence, the intention is that numerical modelling can be feasible and/or require few adjustments when applied to various flows [19].

Eddy viscosity models parameterized the SGS stress tensor in a way to relate it with the resolved scales through the deformation tensor Sij:


where τija=τij13δijτkk is the anisotropic part of SGS stress tensor, and νT denotes the SGS viscosity due to small scales.

The Smagorinsky model [16] proposes the equilibrium state where the small scales dissipate the kinetic energy entirely and instantaneously they receive from the large scales [20]. Following this assumption

with |S|=(2SijSij)12 and Cs=C known as the Smagorinsky coefficient. The Cs coefficient takes values between 0.18 and 0.23 depending on the kinetic energy, Reynolds number, solid boundary proximity and time-space decay energy.

More sophisticated models allow to progressively estimate the Smagorinsky coefficient in time space. The dynamic Smagorinsky model introduces a test filter operation where the filter size is taken as ΔT = αΔ. Typically α = 2 [17, 18]. The model coefficient is calculated to apply the Germano identity:


where LijT denotes the Leonard test filter stress, and Tij denotes the SGS stress to the test filter scale. The last term can be modelled following the eddy viscosity model (Eq. (8)):


Based on Eqs. (8) and (11), the relation in Eq. (10) can be estimated as

LijT=2CTΔ2[α2|S|TSijT(|S|Sij)T]Mij     =2CTΔ2Mij

where the scale invariance of the coefficients (CT=C) has been assumed, and CT is found by optimizing through least squares minimization of error function εrr=(LijT2CTΔ2Mij)2 when (εrr)CT=0 [17]:


The brackets  . H in Eq. (13) denote the average over homogeneous directions and need to be introduced in order to guarantee numerical stability of the procedure because when CT<0, the numerical stability is lost. This operation creates an important limitation, restricting the simulations to flat terrains.

In order to avoid this constraint, a modified dynamic Smagorinsky model is proposed using another test filter size ΔT=3Δ (α = 3) and computing the test filter coefficient [8] as

{12Δ2LijTMijMijMij>0    then    CT=  12Δ2LijTMij(MijMij)12Δ2LijTMijMijMij0    then    CT=0

Therefore, another problem of this model is the inaccurate results near the solid wall (e.g. near of terrain) because in this region the viscous layer is present, and the sub-grid stress is a significant fraction of the total stress. These models are incapable of computing the small scales present in this region [21] and lead to the underestimation of coefficient CT, resulting in an overestimation of the velocity flow near the solid wall [22]. This has been accomplished by various types of ad hoc corrections such as Van Driest damping [23] and intermittence functions [24]. The Van Driest damping functions applied to the first row of equation (14) yield

CT=12Δ2LijTMijMijMij[1ez+/25] 2Van Driest damping      when     CT>0

where z+=zu*ν denotes the vertical coordinate in wall units and u* the friction velocity of flow.

The SGS stress tensor to scalar variables is computed in some way:


where Pr is the Prandtl number, a variable computed dynamically. The Germano identity applied to scalar ψ:


where Iiψdenotes the SGS stress tensor according to the test filter. This tensor can be modelled as


The relationship in Eq. (17) can be estimated by Eqs. (16) and (18):

QiT=CTΔ2Pr[α2|S|TψTxi(|S|ψxi)T]Ni      =CTΔ2PrNi

The least squares minimization of error function εrr=(QiTCTΔ2PrNi)2 yield


The computation of the SGS stress tensor νT in the eddy viscosity models (Eq. (8)) and the Prandtl number Pr in the turbulent diffusivity model (Eq. (16)) has been introduced into the ARPS code by Aguirre [8].

2.3. The parametrization of microphysical processes

In addition, the ARPS code has been designed specifically to describe thunderstorms. So, sub-models of heat flow and water steam, cloud formation and precipitation are included. Ground relief, vegetation types, soil types and initial conditions of the atmospheric state are taken into account for the simulation of these phenomena. The parameterization model is called warm-rain microphysics and is based on the descriptions of Klemp [25] and Soong [26]. This model considers three water categories such as water vapour (qv), cloud water (qc) and rainwater (qr), where the latter two are characterized by its size. At the beginning of the process, cloud water droplets are formed when the air becomes saturated and condensation occurs. Then, if the water mixing ratio exceeds a critical threshold in the cloud interior, raindrops form and the collision-coalescence process begins. If after crossing the base of the cloud they meet air below the saturation point, then the process of evaporation occurs, gradually diminishing the size of the raindrops. The time rate of evaporation and condensation will depend on certain parameters and on the water mixing ratio in the air. As all these processes involve transference of energy in heat and latent vapour, it is necessary to fit the potential temperature of the air. The theoretical formulation will not be described here because it is not the aim of this chapter. For more details on the microphysical processes in the ARPS code, see Xue [27].

2.4. Stochastic Lagrangian one-fluid particle model (STO)

The SGS models simulate the effect of unsolved scales as an energy sink of the solved LES scales. However, atmospheric dispersion phenomena are well simulated when small unsolved scales trajectory are computed. To obtain a more realistic description of the trajectories of fluid particles that transport the chemical species in a turbulent regime, the small-scale simulation unsolved by LES must be carried out. A coupled model between the LES and stochastic Lagrangian models (STO) for the small scales is proposed. The fluid particle Lagrangian velocity Ui is computed by solving the LES-STO coupled model by using a Langevin equation:


where hij is the dynamic deterministic coefficient, and qij is the dynamic random coefficient (in analogy with the Brownian movement), and it is linked to the statistical properties of the turbulence. ηj denotes a random variable whose mean value is null and covariance:


This property suggests that ηj is correlated neither in space nor in time.

To obtain the Langevin equation terms, the following hypotheses are proposed:

  1. The small scales of the turbulence are statistically isotropic; that is, they lose the memory of eddy geometry that originated them.

  2. They are located far from the Kolmogorov inertial range energy spectrum, then away from the production area of kinetic energy.

  3. There is transfer of energy from large scales towards small ones, which is dissipated by molecular viscosity.

The first is not completely sustainable in closeness to solid walls or to the ground. For this reason, an anisotropic model will be used. The last hypothesis suggests that a relationship exists between the results of the LES and the coefficients of the Langevin model of the equation (21). In this way, these terms are calculated dynamically in each cell and for each time step from the results of the LES. Consequently, the fluid particle moves inside the cell of the calculation grid by following the evolution of the large scales resolved by LES to which is added a fluctuation that simulates the behaviour of the small scales of movement produced in its interior. The fluid particle velocity decomposition can be expressed as


where ui is a Lagrangian fluctuation velocity of fluid particle due to sub-grid scale (SGS) turbulence. To obtain the coefficients of Eq. (21), the transport equation of probability density function (PDF) applied to velocities field, called the Fokker-Plank equation, can be used [4]:


The decomposition Eq. (23) is used into the PDF transport Eq. (24):

PLt=vi[hij(uj,t)PL]vi[hij(uj,t)PL]           +122vivj[qik(ui,t)qjk(uj,t)PL]+122vivj[qik(ui,t)qjk(uj,t)PL]

Also the same decomposition Eq. (23) is used in the Langevin Eq. (21):


In this analysis, some simplifications in Eq. (26) with reference to the hypothesis mentioned above can be performed. In particular, the hypothesis (b) proposes that small scales are far from the inertial range energy transfer, then it can be assumed that the random term qijis completely defined by uj small scales (sub-grid scales); in other words qij(uj,t)=0. Therefore, separating large scales and sub-grid scales of Eq. (26) can be assigned the following equivalences:


The first row of Eq. (27) is the material derivative of large scales. It is computed according to Eq. (4):


Therefore, the first term of the second member of Eq. (25) is obtained from LES:


For the other terms of Eq. (25), a turbulence model is proposed because these are terms dependent on the unsolved scales. For deterministic term, which takes account of unsolved scales, Gicquel [7] suggested that this value is proportional at velocity of small scales and a tensor that note the turbulent kinetic energy of sub-grid scales. The idea of this proposal is that the greater the turbulent kinetic energy of sub-grid scales, the greater is the number of small vortices and therefore the greater importance of this term. The other terms express the random component sub-grid. They can be modelled taking into account the assumption (a); in other words, the small scales have an isotropic behaviour. Pope [6] proposed an expression for these terms which takes into account also the rate of dissipation in homogeneous isotropic sub-grid turbulence ε.

With these considerations, the transport equation for the probability density function of velocity field in Eq. (25) can be summarized as


where C0=2.1 denotes the Kolmogorov coefficient.

With the hypotheses and using the PDF transport Eq. (30), the dynamic deterministic coefficient and dynamic random coefficient of Eq. (21) can be expressed as

{hij(Uj,t)=dujdt+αijuijqij(Uj,t)=C0ε δij

Since the material derivative of large scales is obtained using LES in Eq. (28) and the rate of dissipation of turbulent kinetic energy ε is computed with a gradient model [28] by the ARPS code, the remainder proposes an expression for the deterministic tensor αij. These tensors are related to the statistical properties of sub-grid turbulence [29]. For the more complex case, these features are not steady, inhomogeneous and anisotropic flow, as developed in the atmospheric boundary layer over heterogeneous rough terrain near the soil surface:


where K=12(uiuj)δij denotes the sub-grid turbulent kinetic energy; it is obtained with transport Eq. (28). Rij=(uiuj)denotes the Reynolds tensor of sub-grid scales; it is obtained using the turbulent diffusivity models (Eq. 8). More details of this resolution can be found in [9, 10] and [29].

2.5. Solid particle simulation model

The solid particle simulation addressed in this work is on the order of tens of micron size. Therefore, the following assumptions may be raised:

  1. The forces due to gravity and viscosity are dominant regarding the other forces on the solid particles.

  2. The shape of the solid particles is approximately spherical.

  3. The concentration of solid particles is not sufficiently large to influence kinematic or thermodynamic properties of air transports.

The motion equation of a solid particle immersed in fluid flow can be written as follows [30]:

where Vi is the solid particle velocity, and Fi denotes the forces per unit mass.

The assumption (a) has been proven when the density of the solid particles is more than 1000 times of the air density which transports. Then, the second member (Eq. (33)) can be summarized as


in which the first term on the right-hand side denotes the drag force, while the second term is the force of gravity acting in the vertical direction indicated with the index 3, both per unit mass of solid particle. τs denotes a time scale of particle acceleration. It depends on the rate of densities between the solid particle and the air, solid particle diameter ds and a drag coefficient CD.


The hypothesis (b) suggests that the drag coefficient CD is proportional to the diameter of the solid particle ds and also depends on the fluid viscosity. This condition is usually assessed by calculating the Reynolds number of the solid particlees.


where ν denotes the dynamic viscosity of fluid.

In laminar flow (or Stokes flow), ℜes ≤ 1 and the drag coefficient can be calculated as

In turbulent flows, the authors propose many different expressions. Sommerfeld [31] uses the following:

CD={24Res(1+0.15Res0.687)     if         1<Res<10000.44                              if               Res1000

In Eqs. (34), (35) and (36), Ui is the velocity of fluid particle on the position of solid one. Aguirre [32, 33] denotes it as `the fluid particle velocity seen by the solid particle´, and it is estimated with stochastic equation similar to Eq. (21) but using a Reynolds decomposition (RANS). These authors corrected the trajectories of fluid particles seen by the solid particle using a scale of weighted characteristic time taking into account the density and diameter of the solid particles.

In the decomposition of large and small scales, Ui is computed as Eq. (23), and ui can be obtained with the simulation of fluid particles emitted from the same sources as the solid particles and dispersed in the atmospheric boundary layer following the equations of motion described in the previous section. However, it is very difficult that after a few time steps, the positions of the pairs of particles (fluid-solid) launched from the same source, are still in matching positions. An alternative to this problem is to `detect´ the fluid particle closer to the solid particle after each time step, to assign the speed Ui. Firstly, this requires that both sets of particles are computed in the simulation, and secondly, the use of a greater amount of fluid particles that of solid particles because of their different trajectories. Each new position in a solid particle must find a nearby fluid. In this case, it is desirable that the source of fluid particles is much larger than that of the solid. Vinkovic [11] followed this method. An alternative, which uses less computer resources, is to calculate Ui in the exact position of the solid particle starting from the sub-grid turbulent kinetic energy K, then using a stochastic Eq. (21). Thus, it is not necessary to compute the trajectories of fluid particles or use search algorithms proximity to the solid particle. This method is used in this chapter and detailed below.

2.6. Numerical method to compute the solid particle trajectory

The discretized equations for the position and velocity of the solid particles can be written as

{Vi (n+1)=Vi (n)+Δtτs(Ui (n)Vi (n))giΔtδi3Xi (n+1)=Vi (n+1)Vi (n)2Δt

where the subscript in parentheses denotes the number of time instant simulation, and Δt=t(n+1)t(n)is the time step.

The velocity of fluid particle seen by the solid particle is computed following Eq. (23):

Ui (n)=ui (n) +ui (n)

The speed of the large scales to the position of the solid particle in the four nodes of the grid closest to it is weighed. Then, the component representing the speed of sub-grid fluid particle at solid particle position is computed by using the second row of Eq. (27) in discretized form as:

ui (n)=ui (n1)+αij (n)ui (n1)Δt+C0εΔt χ(n)

where χ(n) denotes an independent random variable with zero mean and unit variance at time t(n), and αij (n) denotes the tensor computed with Eq. (32) at solid particle position and to estimate the sub-grid scale of fluid particle velocity in the previous time step at solid particle position, is proposed isotropic turbulence (Pope, 1994):

ui (n1)=23K(n) χ(n)

3. Comparison with experimental laboratory measurements

3.1. Description of Gong’s experiment

Gong [34] carried out measurements of mean velocity and air flow fluctuation in the neutral turbulent layer produced in the wind tunnel of the Department of Agriculture of the University of Reading (UK) using the methodology of generation of turbulent flow of Counihan [35]. The authors installed a rubber sheet on the floor of tunnel to simulate a rough floor and a hill with a slight slope. Details of the geometry of the tunnel and the simulated hill can be found in Gong [34]. Gong [13] measured concentrations of a passive gas (carbon dioxide) incorporating a nearer point source upwind of a bi-dimensional symmetrical hill placed transversally to the air flow direction as shown in Figure 1. The characteristic data of the boundary layer generated in the laboratory, the diameter and height of the gas emission source as well as its position are as follows:

  • - Thickness of boundary layer: D = 300.00 mm.

  • - Mean velocity of flow upper of the boundary layer: Ue = 8.00 m/s.

  • - Friction velocity: u* = 0.44 m/s.

  • - Parameter of ground roughness: z0 = 0.17 mm.

  • - Height of emission source of CO2: zs 13.00 mm.

  • - Diameter of the emission source: ds =1.35 mm.

  • - Distance between the emission source and hill crest: xs = 350.00 mm.

  • - Height of hill crest: hc = 31.00 mm.

  • - Distance between crest and foot of hill: Lc = 200.00 mm.

  • - Mean CO2 concentration at the exit of the emission source: C0 = 400.00 ppm.

The measurements of gas concentration were made on the plane of axial symmetry containing the emission source, obtaining mean value profiles in five positions:

xxs={ 150.00 mm: foot, upwind 250.00 mm: halfway uphill, upwind 350.00 mm: crest 450.00 mm: halfway uphill, downwind 550.00 mm: foot, downwind 


Figure 1.

Scheme of Gong’s experiments [13, 34].

Below is the comparison between the experimental measurements and the results of the simulation using the coupled model LES-STO considering the hypothesis of statistically inhomogeneous and anisotropic turbulence for the αijtensor (32). Besides, other numeric simulations were carried out without using the stochastic model for sub-grid turbulence; that is, in this case, the particles only followed the trajectories imposed by the LES.

3.2. Numerical simulation of Gong’s experiment

For statistical calculations, results were seen after 4 seconds of the start simulation to get a permanent state of the particles within the computational domain. The total physical time simulated was 100 seconds. The rate of injection of particles was fixed at 50,000 particles per second. The time step of the simulation has been proposed to ensure the numerical stability of the computation at 0.02 seconds. The temporal evolution of the quantity of fluid particles in the domain of the calculations shows that close to 130,000 are present in a permanent state. The total quantity of particles injected for the entire simulation time was 10,125,151. The author of the experiment presented the results of gas concentration normalized:

Cn=CD Ue2zsus2C0


Figure 2.

Profiles of normalized concentrations upwind hill. Left-hand side: x – xs = 150 mm (foot of the hill). Right-hand side: x – xs = 250 mm (halfway on the hill).


Figure 3.

Profiles of mean concentration at the crest of the hill at x – xs = 350 mm.

where C is the gas concentration, and us is the air velocity to the height of the emission source. Figures 24 show the profiles of the mean concentration using LES and LES-STO. Figure 5 shows concentration at ground level (left) and the standard deviation of the height of the centre of the gas plume on the axial plane that contains the emission source (right). Figure 6 shows the mean concentration levels of the gas plume on the axial plane containing the source (left) and at the same at ground level (right), both of which correspond to the case of non-stationary, inhomogeneous and anisotropic sub-grid turbulence (LES-STO) that presents better results in comparison with the experimental measurements.


Figure 4.

Profiles of normalized concentration downwind hill. Left-hand side: x – xs = 450 mm (halfway on the hill). Right-hand side: x – xs = 550 mm (foot of the hill).


Figure 5.

Normalized mean concentration at ground level (left) and standard deviation from the height of the gas plume centre (right).


Figure 6.

Mean gas concentration values simulated with LES-STO. Left-hand side: on the axial plane containing the source. Right-hand side: at ground level.

4. Atmospheric dispersion of solid particles during severe wind

4.1. Description of study case

Lastenia town, northwest of San Miguel de Tucuman city, has operated a foundry of metals for 24 years. It is located on a flat residential area (430 m) and was operated until its closure in the mid-1990s. This plant supplied machinery to large sugar refineries since this is the most important crop in the area. The smelter has two chimneys of a height of 45 and 3 m in diameter. The outlet temperature of the gases was Tg = 220°C. In addition to these combustion gases from the foundry, the smelter emitted particles of different metals that are highly toxic when inhaled directly or ingested indirectly through local crops (e.g. citrus, leafy vegetables). According to Fernández-Turiel [14], the Lastenia region had dangerous concentrations of metals as such as silver (Ag), cadmium (Cd), cooper (Cu), nickel (Ni), plumb (Pb), tin (Sn), zinc (Zn), among others. These measurements were performed in the laboratory with very specific equipment samples taken from in situ in square areas of 10,000 m2 in both soil and plants, for a sector near the smelter. Lastenia town location and area of study are shown in Figure 7.


Figure 7.

Location of study area and sample measurements.

The study area has a rectangular shape of 3.5 km in north-south and 2.4 km east-west. The prevailing winds in the area are in north and south-west direction. However, the most damaging wind direction to the residential area is the eastern sector.

4.2. Simulation details

Calculation grid used for simulation consists of regular prismatic cells of varying height. They have a horizontal dimension of 100 × 100 m (of dimensions equal to those used for in situ sampling by Fernández-Turiel [14]) and height ranging from 3 m to adjacent cells to the ground, up to 42 m for which are the top of the computational domain. The law of height size of the cells follows a hyperbolic tangent function. Thus, the grid consists of 37 cells in the east-west direction, 54 cells in north-south direction and 52 cells in the vertical direction. Figure 8 shows the vertical arrangement (a) and horizontal (b) a sector of the grid to LES computation. The location of the chimneys of the smelter is appreciated too.


Figure 8.

Horizontal (a) and vertical (b) configuration of the grid for LES computation and the location of two chimneys of smelter.

In this case, it has been simulated that the copper particle dispersion has a density ρ = 8900 Kg/m3 and diameter ds = 46.5 μm. A time step Δt = 0.05 second was used. This value has been imposed by the acceleration time of the copper particles (τs = 0.0503 second) calculated using the terminal freefall speed of it (Δt < τs) to ensure numerical stability.

The initial velocity at chimney outlet of solid particles is conditioned to gas temperature Tg. This velocity at outlet chimney w(0) can be estimated with


where U0 is the air velocity at height of chimney h, Ta denotes the air temperature, f is the friction coefficient of wall chimney and D is the diameter.

This formula is derived using the Bernoulli theorem vertically and horizontally to consider the thermal draft and depression generated by the wind to the outlet chimney as shown in Figure 9.

A case of severe eastern wind is simulated. The boundary conditions of LES are considered a probability density function type Weibull with two parameters that serve to force the wind velocity:


where k is a shape parameter and c denotes the scale parameter. These are computed considering the average wind velocity [36]:


In this case, the average wind velocity has been obtained from register of meteorological stations during a severe eastern wind. This value at 10 m over ground is u¯10 = 2.5 m/s. Weibull parameters obtained with Eq. (46) are k = 1.486 and c = 2.766.

The logarithmic law has been used to estimate the average velocity of wind under 100 m height. It is used with a friction velocity u* = 0.276 m/s, a roughness parameter z0 = 0.228 m and a von-Kármán coefficient kv = 0.4. Above 100 m height, the potential law of wind is used with power coefficient n = 4.6.


Figure 9.

Calculation scheme of the initial vertical velocity of solid particles.

4.3. Results of the comparison between the measurements and the numerical simulation

The copper concentrations obtained at 4000 seconds of the simulation with those presented by Fernández-Turiel [14] are qualitatively compared. In Figure 10, the copper concentrations at ground level (in number of particles per square meters) computed using grids of 100 × 100 m and the concentration (in mgCu/Kgsoil) published in Fernández-Turiel [14] are shown. Figure 11 compares copper concentrations at ground level (in number of particles per square meters), but these have been computed using grids of 10 × 10 m. Figure 12 shows the same concentrations as in Figure 11 on a satellite image with the purpose of observing the affected sites.


Figure 10.

Copper concentration at ground level. Solid colour: simulated with 100 × 100 m grid mesh (part/m2). Lines: measurements (mgCu/kgsoil) with 100 × 100 m by (Fernández-Turiel [14]).


Figure 11.

Copper concentration at ground level. Solid colour: simulated with 10 × 10 m grid mesh (part/m2). Lines: measurements (mgCu/kgsoil) with 100 × 100 m by (Fernández-Turiel [14]).


Figure 12.

Copper concentration at ground level simulated with 10 × 10 m grid mesh on a satellite image.

4.4. Discussion

In Figure 10, it is noted that concentrations simulated by the model are calculated by counting the number of copper particles falling in grids of 100 × 100 m similar to those samples that were used in Fernández-Turiel [14]. Concentration contours indicate the soil contamination that has been produced by the smelter during operation; simulating this only considers pollution with the eastern wind. Considering this situation, it can be seen that there is good agreement between the maximum value that indicates the model with closed contour 45 mgCu/Kgsoil concentration. Other major peaks of concentration are studied by Fernández-Turiel [14], but they are associated with other wind directions. However, when the concentrations simulated by the model are calculated using grids 10 × 10 m, as shown in Figure 11, they can be individualized with the plumes of each chimney.

The model shows that concentrations of copper particles have little lateral dispersion. This implies that if concentrations are calculated using smaller grid cells, they will be much higher in the area where they fall. This is the reason why both figures legends differ by a factor of 10. It is likely that if the particle size was smaller, lateral dispersion increases.

The concentrations shown by Fernández-Turiel [14] to the west from the location of the smelter still indicate the existence of copper particles beyond those shown in the simulation. If the copper particles are simulated with different diameters, it is likely that smaller particles travel on greater distance before reaching the ground.

5. Conclusions

We have emphasized the use of computer simulation tools for the atmospheric dispersion of gases and solid particles in this chapter.

The first part of the chapter presented the large eddy simulation (LES) approach in order to numerically solve the turbulent flows of the great Reynolds number as those presented in the atmospheric boundary layer. Next, the coupling between the large eddy simulation (LES) and Lagrangian Stochastic one particle model (STO) was presented to detail how to dynamically calculate the model coefficients based on the kinetic energy of turbulent fluid flow.

The validation of these tools has relied on experimental measurements in a wind tunnel and in situ measurements. The stochastic Lagrangian one-particle method enables simulations without much computational cost and with good results. The results of the numerical simulation using the ARPS code (LES-STO) of the concentrations of CO2 emitted by a chimney up wind uphill are in good agreement with experimental measurements performed in a wind tunnel [13].

Fernández-Turiel [14] has found high levels of trace of copper in soil and plants in the vicinity of smelter. As a consequence, it is important to determine the extent of the contaminated area and the concentration of these elements that might be potential hazards due to inhalation and ingestion. The model has been used to simulate the atmospheric dispersion of these particles emitted from the smelter chimneys during a severe eastern wind event. The results closely matched field measurements made by these authors.

Finally, this approach has been able to predict the dispersion and level concentration of pollutants (gases and/or solid particles) in the atmosphere with the aim of predicting the impact on the populations and preventing environments problems. This information can be used to determine the affected areas close to industrial factories and emission control protocols that need to be used during specific weather conditions. By understanding the local pollution concentrations and their movement, future industrial factories sites could be installed in locations that would minimize impacts to the surrounding population.


1 - Deardorff J. W. A numerical study of three-dimensional turbulent channel flow. Journal of Fluid Mechanics. 1970; 41:453–480.
2 - Deardorff J. W. The use of sub-grid transport equations in a three-dimensional model of atmospheric turbulence. Journal of Fluid Engineering. 1973; 4:429–438.
3 - Schumann U. Sub-grid scale model for finite difference simulations of turbulent flow in plane channels and annuli. Journal of Computational Physics. 1975;18:376–404.
4 - Pope S. B. PDF methods for turbulent reactive flows. Energy Combustion Science. 1985;11:119–192.
5 - Haworth D. C. and Pope S. B. A generalized Langevin model for turbulent flow. Physics of Fluids. 1986;29(2):378–405.
6 - Pope S. B. Lagrangian PDF methods for turbulent flows. Annual Review of Fluid Mechanics. 1994;26:23–36.
7 - Gicquel L. Y. M., Givi P., Jaberi F. A. and Pope S. B. Velocity filtered density function for large-eddy simulation of turbulent flow. Physics of Fluids. 2002;14(3):1196–1213.
8 - Aguirre C. A. Dispersión et Mélange Atmosphérique Euléro-lagrangien de Particules Fluides Réactives. Application à des cas simples et complexes [thesis]. Université Claude Bernard, Lyon, France: Ecole Doctorale MEGA; 2005. 323 p. Available from:
9 - Aguirre C. A., Brizuela A. B., Vinkovic I. and Simoëns S. A sub-grid Lagrangian stochastic model for turbulent passive and reactive scalar dispersion. International Journal of Heat and Fluid Flow. 2006;27(4):627–635.
10 - Aguirre C. A., Brizuela A. B., Vinkovic I. and Simoëns S. Eulero-Lagrangian coupled model for the Simulation of atmospheric dispersion of chemically reactive species into the boundary layer. Serie Mecánica Computacional. 2006, 25 (2): 185-205.
11 - Vinkovic I., Aguirre C. A., Ayrault M. and Simoëns S. Large-eddy simulation of the dispersion of solid particles in a turbulent boundary layers. Journal of Boundary Layers Meteorology. 2006;121:283–311. DOI: 10.1007/s10546-006-9072-6
12 - Vinkovic I., Aguirre C. A., Simoëns S. and Gorokhovski M. Large-eddy simulation of droplet dispersion for inhomogeneous turbulent wall flow. International Journal of Multiphase Flow. 2006;32(3):344–364.
13 - Gong W. A wind tunnel study of turbulent dispersion over two – and three- dimensional gentle hills from upwind point sources in neutral flow. Journal of Boundary Layer Meteorology. 1991;54:211–230.
14 - Fernández-Turiel J. L., Aceñolaza P., Medina M. D., Llorens J. F. and Sardi F. Assessment of a smelter impact area using surface soils and plants. Journal of Environmental Geochemistry and Health. 2001;23:65–78.
15 - Leonard A. Energy cascade in large eddy simulation of turbulent fluid flow. Advances in Geophysics. 1974;18(A):237–248.
16 - Smagorinsky J. General circulation experiments with the primitive equations. I. The basic experiments. Monthly Weather Review. 1963;91:99–164.
17 - Lilly D. K. A. Proposed modification of the Germano subgrid-scale closure method. Physics of Fluids. 1992;4(A):633–635.
18 - Germano M., Piomelli U., Moin P. and Cabot W. H. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids. 1991;A(3):1760–1765.
19 - Lévêque E., Toschi F., Shao L. and Bertoglio J-P. Shear-improved Smagorinsky model for large-eddy simulation of wall-bounded turbulent flows. Journal of Fluid Mechanics. 2007;570:491–502.
20 - Piomelli U. and Balaras E. Wall-layer models for large-eddy simulations. Annual Review of Fluid Mechanics. 2002;34:349–374.
21 - Porté-Agel F., Meneveau C. and Parlange M. A scale-dependent dynamic model for large-eddy simulation: application to a neutral atmospheric boundary layer. Journal of Fluid Mechanics. 2000;415:261–284.
22 - Cui G. X., Xu C. X., Fang L., Shao L. and Zhang Z. S. A new subgrid eddy-viscosity model for large-eddy simulation of anisotropic turbulence. Journal of Fluid Mechanics. 2007;582:377–397.
23 - Van Drieas E. R. On turbulent flow near wall. Journal of Aero Sciences. 1956;23:1007–1011.
24 - Piomelli U., Ferziger J., Moin P. and Kim J. New approximate boundary conditions for large-eddy simulation of wall-bounded flows. Physics of Fluids. 1989;A(1):1061–1068.
25 - Klemp J. B. and Wilhelmson R. B. The simulation of three-dimensional convective storm dynamics. Journal of Atmospheric Science. 1978;35:1070–1096.
26 - Soong S-T. and Oruga Y. A comparison between axi-symmetric and slab-symmetric cumulus cloud models. Journal of Atmospheric Science. 1973;30:879–893.
27 - Xue M., Drogemeier K. and Wong V. The Advanced Regional Prediction System (ARPS). A multi-scale non-hydrostatic atmospheric simulation and prediction model. Part I: model dynamics and verification. Meteorology Atmospheric Physics. 2000;75:161–193.
28 - Deardorff J. W. Stratocumulus-capped mixed layer derived from a three dimensional model. Boundary Layer Meteorology. 1980;18:495–527.
29 - Aguirre C. A. and Brizuela A. B. Numerical Simulation of Atmospheric Dispersion Gas Passive on a hill using a coupled model. Serie Mecánica Computacional. 2008; 27 (4): 217-237.
30 - Kosinski P., Kosinska A. and Hoffmann A. C. Simulation of solid particles behaviour in a riven cavity flow. Powder Technology. 2009;191:327–339.
31 - Sommerfeld M. Analysis of collision effects for turbulent gas–particle flow in a horizontal channel: Part I. Particle transport. International Journal of Multiphase flow. 2003;29:675–699.
32 - Aguirre C. A., Y. Guo and Ayrault M. Dispersion of solid particles in saltation movement in a turbulent flow. Journal Comptes Rendus Mécanique. 2004; 332: 627-632.
33 - Aguirre C. A., Simoëns S. and Ayrault M. Dispersed of solid heavy particles in a homogeneous turbulence. In: Brebbia C. A. and Martin Duque J. F., editors. Air Pollution. 10th ed. Southampton, UK: Wessex Institute of Technology Press; 2002. p. 6.
34 - Gong W. and Ibbetson A. A wind tunnel study of turbulent flow over models hill. Boundary Layer Meteorology. 1989;49:113–148.
35 - Counihan J. An improved method of simulating an atmospheric boundary layers in a wind tunnel. Journal of Atmospheric Environment. 1969;3:197–214.
36 - Justus C. G. Wind and system performance. Philadelphia, Penn, USA: The Franklin Institute Press; 1978. 120 p.