Overview of air dispersion modelling

Overview of air dispersion modelling

10.2.1. The transport equation

In a selected V₁ volume of the fluid, mass conservation of the component described with c concentration can be expressed as:

,

(10.1)

Where is the wind vector, S_c is the source term and D_c is the diffusion coefficient. Equation (10.1) represents the change of the total mass of the material within volume V₁ as the sum of the advective flux through the borders of the volume, source terms inside the volume, and the diffusive flux. In dispersion models, wind field and other meteorological data is obtained from measurement or a numerical weather prognostic (NWP) model, thus the only unknown term in equation (10.1) is the c concentration field. We can transform equation (10.1) into a differential form using Gauss’ formula and generalizing the integrates to any V₁ volumes:

.

(10.2)

This is the dispersion equation that describes advection, source and molecular diffusion processes. Dry and wet deposition, chemical or radioactive decay is part of the S_c term, while gravitational settling can be added as an extra advection component. Turbulent diffusion, however, is not represented in equation (2).
Turbulence is usually taken into account with Reynold’s theory that splits the wind and concentration field into time-averaged and turbulent perturbation values:

	(10.3)
	(10.4)

From equations (10.2) and (10.3, 10.4) we can construct a dispersion equation that represents both time-averaged and perturbation components:

(10.5)

Time averaging equation (10.5) will eliminate all the components that contain a single perturbation term, as the Reynolds model is based on the assumption that turbulent perturbations’ time average is zero. However, not all perturbations will disappear as the covariance term’s time average is not necessarily 0:

(10.6)

We can here conclude that dispersion equation (10.5) for turbulent flows can be written in the same form as equation (10.2) with the addition of three eddy covariance terms. Writing the turbulent components explicitly and assuming an isotropic molecular diffusion, equation (10.6) can be rewritten in a form that is widely used in atmospheric dispersion modelling:

(10.7)

The right side of equation (10.7) describes advection, source terms, molecular diffusion and horizontal and vertical turbulent fluxes. In the atmosphere, turbulent mixing is magnitudes more efficient than molecular diffusion thus the third component of equation (10.6) can usually be neglected. However, in the laminar layer within 0.1 − 3 cm of the ground, turbulence is very weak, therefore molecular diffusion gets a large importance in the investigation of soil-atmosphere fluxes and deposition processes, often treated as a resistance term.
Equation (10.7) involves four new variables in the equation. There are two ways for the closure of the turbulent dispersion equation: either we construct new transport equations for the turbulent fluxes or we use parameterization to express the turbulent fluxes with time-averaged concentration and wind values. The former approach leads to the Reynolds Stress Models (RSM), while the latter is the widely used gradient transport theory, or K-theory. These are presented in details in Stull (1988), here we provide a brief outline of their results.
As an analogy of Fick’s law for molecular diffusion, gradient transport theory is based on the assumption that the x directional turbulent flux is proportional to the first component of the gradient of the concentration field:

,

(10.8)

where K_x is the x directional turbulent diffusion coefficient or eddy diffusivity. Using this approach equation (10.7) results in a form where turbulent fluxes are expressed as an additional diffusion term:

,

(10.9)

where K is a diagonal matrix of the K_x, K_y, K_z eddy diffusivities. Due to the different atmospheric turbulent processes in horizontal and vertical direction, K cannot be assumed to be isotropic. Furthermore, while D_c is a property of the chemical species, K is a property of the flow, thus it varies in both space and time. Assuming an incompressible fluid and isotropic horizontal turbulence and neglecting the molecular diffusion, the dispersion equation can be written in a form:

,

(10.10)

where is the horizontal divergence operator, K_h is the horizontal and K_z is the vertical eddy diffusivity. These two parameters need to be estimated at each grid point and timestep through various parameterizations.
Reynolds Stress Models (RSM) construct new transport equations for the turbulent fluxes based on the turbulent kinetic energy and dissipation values (Launder et al., 1975). This approach has larger computational cost than gradient transport models, however, its results proved to be more accurate in several microscale cases (Rossi and Iaccarino, 2009, Chen, 1996). While RSMs have been successfully applied in CFD-based microscale atmospheric models (Riddle et al., 2004), on meso- and macroscale, computationally more efficient eddy diffusivity approach is used (Draxler and Hess, 1998, Lagzi et al., 2009).
The dispersion equation (10.10) can be solved numerically with spatial discretization of variables on a grid, which is often referred to as the Eulerian approach. Under some assumptions, (10.10) can also be solved analytically and provides a Gaussian distribution that is widely used in Gaussian dispersion models. A stochastic solution also exists where instead of solving the partial differential equation (PDE), equation (10.10), the concentration field is given as a superposition of a large number of drifting particles (Lagrangian approach).

10.2.2. Turbulence parameterization

Turbulence is a key process in dispersion simulations: while downwind dispersion is usually dominated by advection, crosswind or even upwind turbulent mixing is magnitudes more efficient than molecular diffusion. Turbulence is treated on different scales: macroscale turbulence, with a scale larger than the numerical weather prediction (NWP) model’s grid resolution, is computed explicitly within the NWP thus it is taken into account in the advection term. Subgrid-scale turbulence causes a velocity and concentration fluctuation, often referred to as turbulent diffusion. We note that subgrid-scale velocity fluctuation also has an effect on the large scale flow through turbulent viscosity. This effect is treated with various turbulence parameterizations in the NWP models.
Atmospheric dispersion can be regarded as a sum of two main effects: the mechanical turbulence caused by wind shear, and the thermal turbulence caused by buoyancy. Their characteristics and dependence on measurable variables is very different (Table 10.1).
Mechanical turbulence estimates rely on 3D wind field measurement data to obtain wind shear values. Surface roughness also has an important role in generating mechanical turbulence through friction in the ground layer. Surface roughness is usually measured with the z₀ roughness length, a characteristic length of surface obstacles. Typical roughness lengths of different surfaces are shown in table 10.2. The models make the difference between complex terrain and surface roughness upon the scale of the problem: while the flow around large scale geometry covered with multiple grid points is computed explicitly in the model, sub-grid scale geometry is treated as roughness and parameterized through its effect on turbulence.
Table 10.1: Parameters that affect turbulence patterns in the planetary boundary layer (PBL)

Mechanical turbulence (wind shear)	Thermal turbulence (buoyancy)
3D PBL wind field	3D PBL (potential) temperature field
Surface roughness	Sun elevation
Complex terrain	Cloud cover
	Albedo
	Surface land use (sensible heat)
	Surface evapotranspiration (latent heat)

Table 10.2: Typical roughness length and albedo values of different surfaces

Surface	Roughness length	Albedo
Open water	1 mm	0.1
Fresh snow	5 mm	0.9
Flat terrain with low grass	5 cm	0.25
High crops	25 cm	0.2
Forest	1 m	0.15
Buildings	several meters	0.05

Thermal turbulence depends largely on the atmospheric stability and the surface’s radiation budget. Under stable conditions, thermal turbulence is low, and mixing is determined by the wind shear strength. On the other hand, in a convective boundary layer (CBL), thermal turbulence has a significant role in dispersion processes. Surface parameters, like albedo (Table 10.2) or potential evapotranspiration have a large importance in the estimation of sensible and latent heat fluxes that determine the thermal turbulence intensity.
The strength of the turbulence can be described using various measures. One of them is eddy diffusivities (m²s^-1) that represent a diffusion coefficient in the dispersion equation (10.10). Another approach is to estimate the u’, v’, w’ wind fluctuation components (equation 10.3), and calculate the turbulent kinetic energy (TKE, Jkg^-1) that describe the time-averaged energy of subgrid-scale turbulent eddies (Stull, 1988):

.

(10.11)

While eddy diffusivities and turbulent kinetic energy explicitly describe turbulence intensity, there are dispersion-oriented characteristics that try to estimate the mixing efficiency instead of describing the turbulence itself. Mixing efficiency is often treated as the deviation of a Gaussian plume, which is widely used in Gaussian and puff models. Lagrangian models use a stochastic random-walk simulation for mixing.
As both wind and temperature changes in space and time, atmospheric turbulence is a non-homogenous, non-stationary field. Furthermore, thermal turbulence is highly anisotropic, which means a difficult challenge for turbulence models especially in a convective boundary layer and leads to the separated treatment of horizontal and vertical turbulence in equation (R10.10). While horizontal dispersion is usually dominated by advection, vertical mixing of the planetary boundary layer (PBL) is caused by turbulence because of the large vertical wind shear and temperature gradients, together with the fairly low vertical wind speeds. This means that vertical turbulence is a key process in atmospheric dispersion simulations, and requires sophisticated methods to estimate its strength.
The vertical profile of vertical eddy diffusivity is presented in Figure 10.1 after Kumar and Sharan (2012). It can be seen that far from the ground, turbulence intensity decreases with height thus the h planetary boundary layer height can be defined as the elevation where turbulence becomes neglectable. On the other hand, near-ground turbulence intensity fast increases with height, and reaches its maximum value at the elevation approximately 30–40 % of the PBL height. Near-zero eddy diffusivity at the top of the PBL means that there is no vertical mixing upward, thus pollutants released from the surface will stay in the boundary layer. The PBL height and the advection speed together determine the volume in which the pollutant can dilute within a specified time, thus they have a very significant effect on concentration estimates.

Figure 10.1: Profiles of the normalized vertical eddy diffusivity for different normalized planetary boundary layer heights under (a) unstable conditions and (b) stable conditions from Ref. (Kumar and Sharan, 2012), where h is the planetary boundary layer height and L is the Monin–Obukhov-length, respectively.

Table 10.3: Typical PBL height values

	Night	Day
Spring	150 m	1700 m
Summer	150 m	1900 m
Autumn	150 m	1200 m
Winter	100 m	500 m

PBL height has a strong diurnal and annual variability: it extends over 2000 m on convective summer days, however, it can shrink to a few 10 meters on clear nights (Table 10.3). A key process is the collapse of the mixed layer approximately one hour before sunset, when heat fluxes from the surface are stopped and thermal turbulence is ceased (Figure 10.2). The night time stable boundary layer keeps the surface-based pollutants close to the ground, while the residual layer is detached from the surface until approximately one hour after sunrise.

Figure 10.2: A typical daily cycle of the planetary boundary layer in fair weather (after Stull, 1988).

10.2.2.1. Stability classes

Atmospheric turbulence is determined by wind shear and thermal stratification. In the early years of dispersion modelling, it was an obvious solution to define categories based on easily measurable wind and radiation characteristics like albedo, cloud cover and sun elevation. The most popular stability classification is Pasquill’s method that defines six categories: from the very unstable A to the very stable F class (Table 10.4). Pasquill took into account wind speed, sun elevation and cloud cover data to determine the stability class that provided pre-defined mixing efficiency values. Despite the fact that this classification can handle only six discrete values of turbulence intensity, Pasquill’s method was used for scientific and regulatory purposes for many decades (Turner, 1997, Sriram et al., 2006).
Table 10.4: Pasquill’s stability classes from the very unstable (A) to the very stable (F) atmosphere. Note: neutral (D) class has to be used for overcast conditions and within one hour after sunrise / before sunset. [Source: http://ready.arl.noaa.gov/READYpgclass.php]

	Daytime insolation			Night time
Surface wind speed	Strong	Moderate	Slight	Thin overcast or max. 3 octas low cloud	Min. 4 octas low cloud
< 2 m/s	A	A – B	B	E	F
2 – 3 m/s	A – B	B	C	E	F
3 – 5 m/s	B	B – C	C	D	E
5 – 6 m/s	C	C – D	D	D	D
> 6 m/s	C	D	D	D	D

10.2.2.2. Richardson number

The difference between thermal and mechanical turbulence can be regarded as a separation of the subgrid-scale turbulent energy to potential (buoyant) and kinetic energy. Their ratio is described with the gradient Richardson number:

,

(10.12)

where is density and g is the gravitational acceleration. The numerator of equation (10.12) describes buoyancy, while the denominator is the wind shear term. It can be observed that the buoyancy can take both positive and negative values, representing both the unstable stratification that generates turbulence, and the stable stratification that destroys turbulence. The wind shear term, however, can take only positive values as it is always a generator of turbulence. The buoyancy term is the square of the Brunt-Vaisala-frequency, and is often written with more practical potential temperature gradients instead of density (Stull, 1988).
Richardson number gives a first-guess estimate on atmospheric stability and turbulence intensity. In the Ri < 0 case, unstable stratification develops thermal turbulence, which describes a convective boundary layer (CBL). Zero Richardson number means a neutral stratification, thus only mechanical turbulence is present. If the Richardson number is positive, but less than a critical Ri_c value (0 < Ri < Ri_c), mechanical turbulence overrides stable stratification that leads to a mechanic turbulence dominated boundary layer. For Richardson numbers larger than Ri_c, atmospheric stability destroys mechanical turbulence and results in a stable boundary layer (SBL). Estimates for the critical Richardson number are in a range 0.15-0.5, but it was also showed that taking qualitative difference between turbulence patterns based on a pre-defined critical value can lead to unrealistic results. However, in many cases, when surface characteristics and fluxes are not available to perform more sophisticated simulations, the Richardson number is used to estimate the turbulence intensity and/or PBL height.

10.2.2.3. Monin–Obukhov-theory

Monin and Obukhov developed a theory in 1954 for atmospheric turbulence that fully described vertical profiles of turbulence intensities. Their basic idea was to combine the vertical turbulent momentum and heat fluxes into a single number with dimension of length, the so called Monin–Obukhov-length

,

(10.13)

where T’, u’, w’ are the temperature, horizontal and vertical wind turbulent fluctuations, T is the temperature and is the von Karman constant, usually taken equal to 0.4 (Stull, 1988). Equation (10.12) is valid for dry air, however, for moist air, virtual temperature or virtual potential temperature is used (Stull, 1988). Hardly measurable turbulent fluxes can be estimated with the definition of two parameters: the u^* friction velocity and the q turbulent heat flux. Using these new parameters, the Monin–Obukhov-length can be expressed as a function of measurable variables (Foken, 2006)

,

(10.14)

where c_p is the specific heat and is the density of the air. There are different methods for the determination of u^* and q using net radiation components and (potential) temperature gradient, which largely differ in the convective and in the stable boundary layer.
The Monin–Obukhov-length in (10.13) is a reference scale for boundary layer turbulence, thus any z height can be given with the dimensionless z/L parameter. The z/L is not only a new vertical coordinate, but also a stability parameter: at any given z height, turbulence intensity can be described with the z/L value (Figure 10.1).
Monin and Obukhov defined universal functions for neutral, stable and unstable cases that describe vertical profiles of wind and temperature as a function of the z/L height (Foken, 2006). The u^* friction velocity is often also computed as a function of z/L using an iterative solution (Cimorelli et al., 2005). The turbulent momentum and heat fluxes are given explicitly from the u^* and q values (as in equations 10.13-10.14), however, dispersion models require eddy diffusivity profiles as well. These are obtained through model-dependent parameterizations using the u^* value, the vertical wind and temperature profile provided by the Monin–Obukhov-theory and the universal functions of z/L that describe stability. The Monin–Obukhov-theory was validated against numerous measurement and advanced computational fluid dynamics datasets, and although some weaknesses were shown, it still means the most reliable basis for all fields of atmospheric modelling (Foken, 2006).

10.2.3. Chemical reactions and radioactive decay

Chemical transformation of air pollutants can be diversified and complex especially in the troposphere or stratosphere, where thermal as well as photochemical reactions can occur. Photochemical reactions can produce radicals, which have high affinity to react with other chemical compounds in the atmosphere. This driving force and rich variety of atmospheric components can really provide a complex network of chemical reactions. Formation of air pollutants can be described by reaction mechanisms. A reaction mechanism is a set of elementary reactions by which overall chemical change occurs. Variation of concentrations of air pollutants in the atmosphere can be constructed from the rate equations and can be calculated by a set of ordinary differential equations (ODEs)

,

(10.15)

where is the chemical rate constant for n^th reaction and is the order of reaction in respect for a given chemical species. Usually, for a thermal reaction the rate constant is a “real” constant, however, for a photochemical reaction the rate constant depends on light intensity, and this dependence can be expressed through the annual and diurnal variation of the solar zenith angle ()

,

(10.16)

where a and b are parameters depending on the given photolytic reaction.
Radioactive decay in environmental models can be interpreted as a first order chemical reaction

.

(10.17)

Solution of this equation (10.17) is an exponential function in time. Radioactive decay constant (k) is an inherent property of the radionuclides. Models can take into account consecutive decays with different decay constants, which is mathematically not a complicated task. Equations (10.16, 10.17) are added to the source term (S_c) of atmospheric transport equations (10.10). Models using concentration field (e.g., Eulerian model) can easily handle this mean field description of the concentrations described by ODEs.
However, if we consider dispersion of air pollutants as a dispersion of individual particles, using deposition, chemical and decay rate constants could be meaningless, because these quantities represent “bulk” properties of the system. Here we deal with individual particles, where macroscopic parameters can hardly be used. We can overcome this problem using probabilistic nature for first order reaction, radioactive decay and deposition. The probability that individual particles during a time step will transform, decay or deposit can be estimated by the following relation

,

(10.18)

where k is either the “macroscopic” radioactive decay constant (first order rate constant) or wet or dry deposition constant. Thus we can connect microscopic event to macroscopic properties and constants.