1. Introduction

The air quality of a region is essential for the welfare of the population and the environment, and is directly influenced by the levels of air pollution. Industrial and technological developments generate excessive emission of pollutants, decreasing air quality. Therefore, there is a necessity to develop mathematical models and computer simulations in order to understand and predict impact of dispersion of pollutants in the environment and in case of incidents or accidents evaluate its risks on habitats.

Turbulent diffusion in the atmosphere is usually modelled by the advection-diffusion equation, even though does not explain all the observed phenomena. This equation is considered to be deterministic and its solution describes mean values of substance concentrations whereas the atmospheric dispersion is stochastic because of natural fluctuations that evidently cannot be reproduced by a purely deterministic model.

Thus, the main objective of this work is to investigate the behaviour of the solution to the deterministic advection-diffusion equation whilst introducing effects that shall mimic stochastic properties. For this purpose we modify the originally purely geometrical boundary conditions, i.e. the ground level and the boundary layer height, respectively. More specifically turbulent mixture is believed to take place in various scales, where the largest scale is limited by the boundary layer height, but also smaller scales shall be present.

One could think of the boundary layer as a superposition of various boundary layers, however with different ground and upper layer heights. Such a construction could model the escape of pollutant substances across the boundary layer horizon on the one side and the surface boundary on the other side and are modelled by probabilities to quantify the fraction of pollutant that returns into the boundary layer from above and the process of adsorption or deposition on the ground layer. These effects are represented by reflective and distributed boundary conditions that together with advection-diffusion dispersion define the model in consideration. The consequences of the reflections are analysed using the meteorological conditions and data from the Hanford and Copenhagen experiments.

2. A Locally Gaussian Model

Starting from the continuity equation one can obtain the advection-diffusion equation through the use of the Reynolds decomposition to separate the mean components for the concentration and the velocity fields. Upon taking averages and substitution of the average fluctuations by Fick’s closure, the desired equation for mean concentrations and an a priori known wind field and with all turbulent characteristics parametrised in a time dependent eddy diffusivity matrix coefficient K is attained. For the current study the eddy diffusion is simplified as locally constant coefficients, which may be justified by the fact that the coefficients vary softly only with changing coordinates and is typical for homogeneous turbulence. For details of the derivation see for instance the textbook by Arya [1].

(1)

where c represents the mean concentration of a contaminant (g/m³), u = (u, v, w) are the mean wind speeds (in m/s) in the longitudinal, vertical and cross wind directions, the Nabla symbol ∇ signifies the usual vector differential operator, the eddy diffusivity coefficient is represented by a diagonal matrix K = diag(K_x, K_y, K_z) and S is a source term.

Considering a point source at height H_s that releases instantaneously a pollutant at a time t = 0, in a fixed emission rate Q, then the source term can be cast in an initial condition and equation (1) is written as the initial value problem, disregarding further the slowly varying terms ∇K∇c.

(2)

(3)

This approximate problem can be solved analytically by separation of variables [16] and Fourier transform [19] and have the following solution.

(4)

However, most dispersion problems are due to continuous emissions, which can be idealized by the superposition of instantaneous emissions. Considering a small time interval dτ with an instantaneous emission, then the continuous emission is

(5)

where c is the concentration for the instantaneous and C for the continuous emission. Making use of equation (4) in equation (5)

(6)

Since the solutions (4) and (6) were obtained by Fourier transform, they are valid for the infinite ranges x ∈ (−∞, ∞), y ∈ (−∞, ∞), z ∈ (−∞, ∞), although the dispersion of contaminants is limited at the vertical domain by the ground (z = 0) and the top of the atmospheric boundary layer (z = z_i) thus the infinite range has to be mapped into a finite range.

3. Reflective Boundary Conditions

To justify the mapping of the infinite range z ∈ (−∞, ∞) to the finite z ∈ [0, z_i] we first consider a cut of the distribution at z = 0 and z = z_i, respectively. Usually the boundary conditions for the pollutant dispersion problem are zero flux or zero concentration at the boundaries, however Fick’s hypothesis suggests there should be a flux across the boundaries. In that sense we copy from observation that the layer until the hight where temperature inversion occurs may be considered at least partially decoupled from the wind flux system beyond. Hence, in an ideally decoupled system the lost contributions should be recovered adopting reflecting boundaries, which intuitively agrees with a simple particle ensemble picture where the pollutant that reaches the ground or the top of the atmospheric boundary layer bounces completely back into the domain. For the distributions that means that even after reflections the Gaussian tails exceeding the allowed domain are mirrored back into the finite range z ∈ [0, zi].

Formally, the reflection on the ground and in the atmospheric boundary layer may be viewed as contributions due to a virtual source in some effective heights to both sides below ground and above the boundary layer [2], those heights are the center of the gaussians formed at the ground and at the top of the atmospheric boundary layer, as exemplified in Figures 1, 2 and 3. The sequences that represent the mirror maxima are

(7)

Figure 1. Scheme for the dispersion as if there were no limits in the vertical domain

Figure 2. Scheme for the reflection starting on the ground

Figure 3. Scheme for the reflection starting on the top of the ABL

Substituting those two sequences (7) in the solution for the continuous emission (6), the solution for continuous emission with complete reflection is obtained

(8)

and is now valid for x ∈ (−∞, ∞), y ∈ (−∞, ∞), z ∈ [0, zi].

So far the model still does not represent any property that might be associated to an effect from a stochastic feature. As already argued before, instead of a boundary layer with rigid limits one could mimic a sample of a distribution with different boundary layer heights upon changing the position of the mirror images that compose the total distributions. To this end we introduce the reduction factor ω_b and ω_g in the sequences (7). Note, that the system still maintains its deterministic character, but finite sample of boundary layer configurations with different heights and centre could be interpreted as a manifestation of stochasticity that are used to study the behaviour of the new solution.

(9)

4. Turbulent Diffusivity Parametrisation

To validate the proposed model, more specifically to analyse the impact of reflections on the results, turbulent diffusivity was parametrised to represent meteorological conditions of the Hanford [10] and Copenhagen experiment [11].

The Hanford experiment is a low source experiment (the height of the source H_s was 2 m) with stable to quasi-neutral conditions. A non-depositing tracer was released with an average rate of Q = 0.3 g/s and release time interval of 30 minutes, except for experiment run 05, where the release time was 22 minutes. The measurements were performed at distances 100 m, 200 m, 800 m, 1600 m and 3200 m from the source. The necessary micro-meteorological data for the parametrisation were provided by the experiment and are presented in Table 1. The height of the stable boundary layer (z_i,s) was calculated using the relation z_i,s=0.4(u_∗L/f_c)^1/2, where f_c=1.46x10^-4 is the Coriolis parameter.

Table 1. Micro-meteorological data for the Hanford experiment [10]

The Copenhagen campaign is a high source experiment (the height of the source H_s was 115 m) with convective conditions. The tracer sulphurhexafluoride (SF₆) was released without buoyancy at a constant flow rate ranging from 2.4 to 4.7 g/s and release time interval of 60 minutes. The measurements were performed at distances ranging from 1900 m to 6100 m from the source, depending on the experiment. The necessary micro-meteorological data for the parametrisation were provided by the experiment and are presented in Table 2. The convective velocity scale (w_∗) was calculated using the relation w_∗ = u_∗(−z_i/kL)^1/3, where k is the Von Kármán constant.

Table 2. Micro-meteorological data for the Copenhagen experiment [11]

4.1. Stable Conditions

The eddy diffusion coefficient for stable conditions proposed by Degrazia and Moraes [7] is based on the diffusion theory of Taylor [20] and the turbulent kinetic energy spectrum [17] and can be computed using the micrometeorological data set from table 1.

(10)

Here u_∗ is the friction velocity, z is the observation height, z_i,s is the height of the stable boundary layer, X is the dimensionless distance, n′ is the dimensionless frequency of the turbulent kinetic energy spectrum and (f_m)_w = (f_m)_n,w (1 + 3.7 z/Λ) is the frequency of the spectral peak in the vertical eddy spectrum, where Λ = L(1 − z/z_i,s)^1.5α1−α2 is the local Obukhov length, α₁ = 1.5, α₂ = 1.0 [14] and (f_m)_n,w = 0.35 is the frequency of the spectral peak in neutral stratification [15].

4.2. Convective Conditions

The eddy diffusion coefficient for convective conditions proposed by Degrazia and Moreira [8] is also based on the diffusion theory of Taylor and the turbulent kinetic energy spectrum and can be computed using the micro-meteorological data set from Table 2. Considering α = x, y, z and i = u, v, w, the eddy diffusion coefficient for convective conditions is

(11)

Here w_∗ is the convective velocity scale, z is the observation height, z_i is the inversion height, X is the dimensionless distance, n′ is the dimensionless frequency of the turbulent kinetic energy spectrum, c_i = α_i(0.5 ± 0.05)(2πk)^−2/3 is a constant [6], (fm_∗)_i is the normalized frequency of the spectral peak regardless of stratification and ψ is the dissipation function and has the form [13; 9]

(12)

where L is the Obukhov length in the surface layer. The Copenhagen campaign is a two-dimensional experiment thus the three-dimensional solution (9) is integrated laterally and the eddy coefficient K_y vanishes, requiring only to calculate K_x and K_z. The values for the normalized frequency of the spectral peak are (f_m^∗)_u = 0.67 [15] and (f_m ^∗)_w = z/(λ_m)_w with (λ_m)_w = 1.8z_i [1 − exp(−4z/z_i) − 0.0003 exp(8z/z_i)] [5].

4.3. Wind Speed Profile

In the further, we introduce a simplification, without imposing restrictions on our numerical findings. We assume, that our coordinate system has its x-axis aligned with the average wind speed, which is to a good approximation horizontal with respect to the Earth’s surface. In order to determine the velocity field u=U(z)x with x a unit vector, we need to fix the vertical wind speed profile. The latter has been parametrised following Monin-Obukhov’s similarity theory manifest in the so-called OML-model [3], where close to the surface and because of its roughness, there is a raising profile, whereas sufficiently far from the surface the wind speed remains approximately constant. If z_b = min(|L|, 0.1z_i), then

(13)

where z₀ is the roughness length and Ψ_m is the stability function. For stable conditions the stability function is Ψ_m = (1 + 4.7 z /L) [4] and for convective condition is [18]

(14)

with A = [1 − (16z/L)]^1/4.

5. Validation of the Model

To simulate the results, it was used the complete data set of the Hanford experimental data, except those for the distances x = 100 m and x = 200 m and the complete data set of the Copenhagen experimental data. The comparison of observed (C_o) against predicted (C_p) concentrations for a variety of reflection parameter ω_b,g and number of reflections is shown in figures 4, 5, 6 and 7 for the Hanford experiment and in figures 8, 9, 10 and 11 for Copenhagen.

The corresponding statistical indices [12], i.e. the normalized mean square error (NMSE), the correlation coefficient (COR), the fractional bias (FB) and the fractional standard deviation (FS) are shown in the tables 3, 4, 5 and 6 for the Hanford experiment 7, 8, 9 and 10 for Copenhagen.

A general comment is in order here, although the afore mentioned statistical evaluations are similar to those from parametric inference procedures, their interpretations are different in the present context. In parametric inference the best estimates for parameters were attained for NMSE → 0, but in the present consideration a deterministic model is compared to a relatively small data set from a stochastic phenomenon, so that one does not expect vanishing values for this error. Further, the correlation coefficient does not converge to unity, recalling that the observed data are one sample out of a distribution for a specific situation, that are parametrised using their specific micro-meteorological data. The fractional bias may be interpreted in terms of model fidelity, significant deviations from zero indicate, that the model lacks some relevant physical features.

5.1. Hanford Experiment

Figure 4. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 1 reflection with the parameters ωg = 1.0 and ωb = 1.0 for the Hanford experiment

Table 3. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 1.0 and ωb = 1.0 for the Hanford experiment

Figure 5. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 27 reflections with the parameters ωg = 0.01 and ωb = 0.01 for the Hanford experiment

Table 4. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 0.01 and ωb = 0.01 for the Hanford experiment

Figure 6. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 85 reflections with the parameters ωg = 0.003 and ωb = 0.003 for the Hanford experiment

Table 5. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 0.003 and ωb = 0.003 for the Hanford experiment

Figure 7. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 93 reflections with the parameters ωg = 0.01 and ωb = 0.003 for the Hanford experiment

Table 6. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 0.01 and ωb = 0.003 for the Hanford experiment

5.2. Copenhagen Experiment

Figure 8. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 1 reflection with the parameters ωg = 1.0 and ωb = 1.0 for the Copenhagen experiment

Table 7. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 1.0 and ωb = 1.0 for the Copenhagen experiment

Figure 9. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 4 reflections with the parameters ωg = 0.2 and ωb = 0.2 for the Copenhagen experiment

Table 8. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 0.2 and ωb = 0.2 for the Copenhagen experiment

Figure 10. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 5 reflections with the parameters ωg = 0.15 and ωb = 0.2 for the Copenhagen experiment

Table 9. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 0.15 and ωb = 0.2 for the Copenhagen experiment

Figure 11. Scatter plot for observed (Co) and predicted (Cp) concentrations for none and 7 reflections with the parameters ωg = 0.1 and ωb = 0.3 for the Copenhagen experiment.

Table 10. Statistical evaluation for observed (Co) and predicted (Cp) concentrations with the parameters ωg = 0.1 and ωb = 0.3 for the Copenhagen experiment

One observes in the presented cases, that after inclusion of reflections the correlation of experimental and predicted data improve, which may be interpreted as an indication that with our reasoning we have at least made some point, even though we have not solved the important question of how the best values for ω_b and ω_s shall be obtained. The reflections are repeated until there is no more contribution to the solution, although in all cases the statistical indices seem to remain the same after only a few reflections, this happens as a result of rounding.

6. Conclusions

The present work can be considered an attempt to incorporate stochastic aspects in a deterministic model, i.e. the dispersion of pollutants in the atmospheric boundary layer modelled by the deterministic advection-diffusion equation. The effective boundary layer height may vary according to the turbulent flow dynamics it incorporates hence the boundary layer should have stochastic aspects. A more realistic flow may be thought of as a superposition of various boundary layer problems but with different effective boundary layer heights. So far the intent was to evaluate the contribution of virtual sources superimposing them such as to mimic a finite size sample of a distribution from different boundary layer height realisations.

A variety of trials have shown that partial reflection on the boundary layer horizon and the ground obtain the most significant correlations between model and data suggesting that effects on boundary are essential to model dispersion processes in the atmospheric boundary layer, even though the values for the reflection parameters were established ad hoc.

This improvement in the solution can be associated with the fact that the deterministic equation predicts only mean values of an unknown distribution and is not capable to predict all stochastic properties, which in our case were modelled by the effects of the considered reflections. Moreover, the model does not consider deposition and adsorption on the ground, but the fact that concentration and vertical concentration fluxes are different from zero on the ground level one may reason, that the parameter ω_g is somehow incorporating these properties.

ACKNOWLEDGEMENTS

This work was supported by Linhares Geração SA, Coordination for the Improvement of the Higher Level Personnel (CAPES - Brazil) and National Council for Scientific and Technological Development (CNPq - Brazil).